Complex Cholesky Function

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("&#37;.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.

From a programming perspective, what might help is to write a small but helpful library dealing with complex numbers. For example, you could write a struct representing a complex number with the real and imaginary part, and then a collection of functions to perform basic math operations on such a struct. From there, in your actual program, you would simply need to call these math functions. This separates the complex numbers code from the logic of your program, and might make it much easier to think about. If you write the complex numbers code in a separate source file complete with a header, you're able to reuse your code at a later date in other programs, which is a good thing.

Complex numbers are built into C as of C99, so if you have a reasonable compiler you shouldn't have to actually do anything.

... except use a valid algorithm .....

The simplest thing would be to read up on the Cholesky decomposition for complex matrices (or, more precisely complex hermitian positive definite matrices, which is roughly as far as its applicability goes). A number of the operations that are multiplications of real numbers when working with real positive definite matrices turn out to be involve multiplying one number by the complex conjugate of the other when working with complex matrices.

As to your NaNs when computing sqrt of complex values, or when dividing complex values ..... the usual reason is either that you are dividing by zero or doing an intermediate step of computing sqrt(real*real + imag*imag) for some complex value. The computation real*real + imag*imag will overflow for moderately large values of real or imag. The simple solution is to reorder the computation to avoid overflow. I'll leave finding such a solution as an exercise.

Incidentally, complex values (other than zero) have two distinct square roots. In many numerical algorithms that work with complex values or matrices, it is important to get one of those square roots and not the other.