1
$\begingroup$

I am using my own Routine to calculate the Cholesky-Factorization of a complex, positive definite symmetric Matrix.

My Code Looks like this:

void CholeskyDecomposition(Complex **L, Complex **A, const int dim)
{
    Complex sum, t;
    Complex sum1;

    for (int k = 0; k < dim; k++)
    {
        // First, calculate the diagonal element of a column
        sum1 = A[k][k];
        for (int j = 0; j < k; j++)
        {
            Complex product;
            Cconjmul(product, L[j][k], L[j][k]);
            Csub(sum1, sum1, product);
        }
        Csqrt(&sum1, sum1);
        L[k][k] = sum1;

        // Now, calculate the elements right the diagonal element
        for (int i = k + 1; i < dim; i++)
        {
            sum = A[k][i];
            for (int j = 0; j < k; j++)
            {
                Cconjmul(t, L[j][k], L[j][i]);
                Csub(sum, sum, t);
            }
            Cdiv(L[k][i], sum, sum1);
        }
    }
}

I tested this Code against the MATLAB Cholesky-Factorization which results in

A =

  13.6393 + 0.0000i   1.8844 + 3.4319i  -4.8788 - 4.0195i
   1.8844 - 3.4319i   5.7265 + 0.0000i   1.1568 - 0.4937i
  -4.8788 + 4.0195i   1.1568 + 0.4937i   6.6512 + 0.0000i

K>> L = chol(A)

L =

   3.6931 + 0.0000i   0.5102 + 0.9293i  -1.3210 - 1.0884i
   0.0000 + 0.0000i   2.1454 + 0.0000i   1.3248 - 0.5435i
   0.0000 + 0.0000i   0.0000 + 0.0000i   1.2927 + 0.0000i

The result of my C-Implementation seems not to be totally wrong. The first row is identical to the MATLAB-solution. The second and third row contains Errors in the last column:

MATLAB         <-->       C
1.3248-0.5435i <--> 1.3248+0.08324i
1.2927+0.0000i <--> 1.3998+0.00000i

Can anyone help me to identify the error source?

$\endgroup$
1
$\begingroup$

Sorry for annoying you guys. The error was quite simple:

// Now, calculate the elements right the diagonal element
for (int i = k + 1; i < dim; i++)
{
    sum = A[k][i];
    for (int j = 0; j < k; j++)
    {
        Cconjmul(t, L[j][k], L[j][i]);
        Csub(sum, sum, t);
    }
    Cdiv(L[k][i], sum, sum1);
}

Just swap the L's of the complex conjugate multiplication. Has to be:

Cconjmul(t, L[j][i], L[j][k]);
$\endgroup$
  • 1
    $\begingroup$ Glad you got it sorted! To stop your question from getting bumped by the community bot, you can either accept your own answer (enough time has passed) or delete your question (if you think it won't be useful for other people). $\endgroup$ – Christian Clason Nov 12 '15 at 16:23

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.