Cholesky decomposition

[Bogy007]

I am looking for a way to write a code implementing the Cholesky decomposition with only one loop (on K), utilizing outer product. Not sure how to go about this. MATLAB can do it, but i have to use c++. I have looked at parallelism but that is over my head. Can someone help point my in the right direction.

for(int k = 0; k < N-1; k++)
    {
        a[k][k] = pow(a[k][k], 0.5);
        for(int i = k+1; i < N; i++)
        {

            a[i][k] = a[i][k] / a[k][k];
        }
        for(int j = k + 1; j < N; j++)
        {
            for(int i = j; i < N; i++)
            {
                a[i][j] = a[i][j] - (a[i][k] * a[j][k]);
            }
        }
    }
    a[N-1][N-1] = pow(a[N-1][N-1], 0.5);

[mike_2000_17]

That makes no sense to me. First, utilizing the outer product doesn't make the inner loops disappear, it simply masks them inside a function (or product operation), but they still occur just the same. Second, it would be difficult to implement this outer product mechanism without losing some performance of the overall algorithm. And for such a simple algorithm as Cholesky (which fits in 10 lines of code or so), it doesn't make much sense to do it in layers like that (i.e., one outer loop that calls some outer product function, etc.). So, why exactly do you want to use an outer product inside the loop?

About parallelism, it doesn't really make sense to use parallelism unless the dimensions are very large, because of the overhead of parallelizing it. And if the dimensions were that big, you wouldn't be using a vanilla implementation of Cholesky and a vanilla storage strategy for the matrix anyways.

Btw, about your code, you should not use pow(a, 0.5); where you could instead use sqrt(a); because the latter is clearer and more efficient. Also, you need to detect zero elements on the diagonal, otherwise you will have divide-by-zero problems. Other than that, it looks good.

[Bogy007]

That is what I thought. This was given to me as an assignment for a numerical methods course. I could not get passed the idea that you were to only have one loop on k, since you still needed the function to do the same work.