hernlund
11-17-2004, 11:47 PM
I've modified ludcmp to make it a complex routine, and used it to calculate the determinant of a complex matrix. The magnitude of the real part of the determinant is good to 3 significant figures, but the imaginary part is way off, and the sign of the real part seems random and is incorrect half of the time. Looking at the LU decomposed elements, the random sign errors are completely in the imaginary parts, which explains the errors in the imaginary parts of the determinant. I've gone through this routine several times and can't quite figure out what is going wrong...any ideas?
hernlund
11-18-2004, 01:47 AM
Forgot to add: This all indicates that there must be a branch cut problem happening some where in Crout's algorithm, though I've not been able to identify where it occurs...
agezerlis
05-29-2007, 07:17 PM
...but I am posting it because it may be useful to those who search for "complex determinant" on Google (which is how I got to this page).
I just modified the ludcmp routine and checked it for a number of different matrices (agreement with Mathematica up to at least 5 decimal places).
It's actually pretty simple. One needs to:
1) change the declarations of "a(np,np)" and "sum" from double precision to COMPLEX(kind=8).
2) declare a new COMPLEX(kind=8) variable called (say) "dum2", and then just add a 2 at the end of the last four occurences of the "dum" variable.
3) assign the variable "d" to a new, complex, variable "newd" thus: newd = dcmplx(d,0.D0), after one is done calling ludcmp (but before the loop described in p. 41 of Numerical Recipes).
Everything else can be left as is. (The abs() function is no longer the absolute value but the complex modulus, but this is a question of interpretation, not code modification).
Hope this helps!
Alexandros Gezerlis