Complex Cholesky Function

This is a discussion on Complex Cholesky Function within the C Programming forums, part of the General Programming Boards category; I was wondering if it would be possible to get some help on writing a cholesky function that works for ...

  1. #1
    Registered User
    Join Date
    Jun 2008
    Posts
    45

    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.
    Last edited by magda3227; 07-18-2008 at 08:36 PM.

  2. #2
    Deathray Engineer MacGyver's Avatar
    Join Date
    Mar 2007
    Posts
    3,211
    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.

  3. #3
    and the Hat of Guessing tabstop's Avatar
    Join Date
    Nov 2007
    Posts
    14,185
    Complex numbers are built into C as of C99, so if you have a reasonable compiler you shouldn't have to actually do anything.

  4. #4
    Registered User
    Join Date
    Jun 2005
    Posts
    6,252
    ... 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.
    Last edited by grumpy; 07-18-2008 at 08:58 PM.

Popular pages Recent additions subscribe to a feed

Similar Threads

  1. Getting an error with OpenGL: collect2: ld returned 1 exit status
    By Lorgon Jortle in forum C++ Programming
    Replies: 6
    Last Post: 05-08-2009, 08:18 PM
  2. Troubleshooting Input Function
    By SiliconHobo in forum C Programming
    Replies: 14
    Last Post: 12-05-2007, 06:18 AM
  3. Screwy Linker Error - VC2005
    By Tonto in forum C++ Programming
    Replies: 5
    Last Post: 06-19-2007, 02:39 PM
  4. Replies: 28
    Last Post: 07-16-2006, 11:35 PM
  5. Why am I getting 'undelcared identifier' ???
    By Bill83 in forum C++ Programming
    Replies: 2
    Last Post: 02-15-2006, 12:00 PM

Tags for this Thread


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21