I was wondering if it would be possible to get some help on writing a cholesky function that works for complex numbers. Thus far, I have one that works for real numbers. I figured that adding a few lines of code to handle complex numbers would be easy, but it has been very difficult.

Here is the code I have so far...

Code:

#include<stdio.h>
#include<math.h>
int main()
{
double A_real[4][4] ={{1.0000, -0.7275, 0.6726 , -0.3499},{-0.7275, 1.0000, -0.6288, 0.6726},
{0.6726, -0.6288, 1.0000, -0.7275},{-0.3499, 0.6726, -0.7275, 1.0000}};
double A_imag[4][4] ={{ 0, 0.4042, - 0.3450, 0.5229},{ - 0.4042, 0, - 0.0208, - 0.3450},
{0.3450, 0.0208, 0, 0.4042}, {-0.5229, 0.3450, 0.3450, 0}};
int i, j, k, n=4;
for (k = 0; k < n; k++)
{
if (A_real[k][k] < 0.0)
return;
if (A_real[k][k] < 1e-20)
A_real[k][k] = 1e-20;
A_real[k][k] = sqrt(A_real[k][k]);
for (i = k + 1; i < n; i++)
{
A_real[i][k] /= A_real[k][k];
for (j = k + 1; j <= i; j++)
A_real[i][j] -= A_real[i][k] * A_real[j][k];
}
}
for(i = 0; i < n; i++)
for(j = i+1; j < n; j++)
A_real[i][j] = 0.0;
for(i=0;i<4;i++)
{
for(j=0;j<4;j++)
{
printf("%.4lf\t", A_real[i][j]);
}printf("\n");
}
return;
}

I ended up messing around with that code, making additions to attempt to handle complex numbers, but was not successful at all. I understand that taking the square root of a complex number and dividing two complex numbers by each other are different from the method by which real numbers are calculated. I wrote up code for complex division and taking the complex square root, but keep getting Nan as answers everywhere except for the simplest part of the Cholesky decomposed matrix...the 0's in a triangle.

I am no math genius or good programmer by any means and was hoping someone more intelligent can help me out with this. I'm not looking for someone to do it for me, but I certainly need all the assistance I can get.

Thanks in advance.