c syntax for choleskey Decomposition

Question

double *cholesky(double *A, int n) {
    double *L = (double*)calloc(n * n, sizeof(double));
    if (L == NULL)
        exit(EXIT_FAILURE);
    for (int i = 0; i < n; i++)
        for (int j = 0; j < (i+1); j++) {
            double s = 0;
            for (int k = 0; k < j; k++)
                s += L[i * n + k] * L[j * n + k];
            L[i * n + j] = (i == j) ?
                sqrt(A[i * n + i] - s) :
                (1.0 / L[j * n + j] * (A[i * n + j] - s));
        }
        return L;
    }

Now my question with this code is a follows, I am trying to see the step by step but I am getting a bit confused.

when I write

for(condition)
    for(condition){
        For(k=0;k<j;k++)
            s += L[i * n + k] * L[j * n + k];
        L[i * n + j] = (i == j) ?
            sqrt(A[i * n + i] - s) :
            (1.0 / L[j * n + j] * (A[i * n + j] - s));

    }

This is what I see happening:

First i = 0 and j = 0; then I step further into the code and we get this: for(k) loop. Now my first question is this since j=0 in first instance this for loop does not get evaluated, since k is less than j.

for (int k = 0; k < j; k++)

But what of the kode below this for loop does get evaluated. Since L is an array of zeros, then s+=l[0]*l[0] should equal 0, but I dont get how the loop did run at all. Next does everything else below the for(k) get evaluated?

if so i==j is true so L[i * n + j] = sqrt(A[i * n + i] - s)(in this case equals 0).

now going back to the top my nest questions becomes since

for(condition i )
   for(condition j)

don't have brackets does for(i=1 ) get evaluated following that for(j) gets evaluated twice below before j=0 and j=i?

Really would appreciate all the help I could get.

Thank you,

[Casting calloc is bad](https://stackoverflow.com/questions/605845/do-i-cast-the-result-of-malloc) — Ed Heal, Jun 23 '17 at 07:54
K loop won't run on the first iteration as j == 0. L[0][0] is set. On the second iteration j is 0, 1, and the k loop runs once to modify s. — Malcolm McLean, Jun 23 '17 at 08:44
Compile with all warnings & debug info (e.g. `gcc -Wall -Wextra -g` with [GCC](http://gcc.gnu.org/)....) then **use the debugger** (e.g. `gdb`) to run your program step by step and examine its state. — Basile Starynkevitch, Jun 23 '17 at 09:13

score 1 · Accepted Answer · answered Jun 23 '17 at 19:38

This code may be easier to understand if you add all the curly brackets, and convert the ternary operator into an if-else statement:

double *cholesky(double *A, int n) {
    double *L = (double*)calloc(n * n, sizeof(double));
    if (L == NULL) {
        exit(EXIT_FAILURE);
    }
    for (int i = 0; i < n; i++) { // loop A
        for (int j = 0; j < (i+1); j++) { // loop B
            double s = 0;
            for (int k = 0; k < j; k++) { // loop C
                s += L[i * n + k] * L[j * n + k];
            }
            // still in loop B
            if (i == j) {
                L[i * n + j] = sqrt(A[i * n + i] - s);
            } else {
                L[i * n + j] = (1.0 / L[j * n + j] * (A[i * n + j] - s));
            }
        }
    }
    return L;
}

The code execution goes like this:

loop A starts and sets i = 0, and assuming i is less than n it runs
loop B starts and sets j = 0, j is less than i+1 (1) so it runs
s = 0
loop C starts and sets k = 0, but since k < j (0 < 0) is false, the code in loop C does not run

The rest of loop B then executes:

if (i == j) {
    L[i * n + j] = sqrt(A[i * n + i] - s);
} else {
    L[i * n + j] = (1.0 / L[j * n + j] * (A[i * n + j] - s));
}

Since i = 0 and j = 0, L[i * n + j] is set to sqrt(A[i * n + i] - s)

loop B increments j to 1, but since j < (i+1) (1 < 1) is false, it terminates
loop A increments i to 1, and assuming i is less than n it runs
loop B starts and sets j = 0, j is less than i+1 (2) so it runs
s is set to 0
loop C starts and sets k = 0, but since k < j (0 < 0) is false, the code in loop C does not run
The rest of loop B runs. i == j is false, so L[i * n + j] is set to (1.0 / L[j * n + j] * (A[i * n + j] - s))
loop B increments j to 1, and since j < (i+1) (1 < 2) is true, it runs
s is set to 0
loop C starts and sets k = 0, and since k < j (0 < 1) is true, it runs
s is incremented by L[i * n + k] * L[j * n + k]
loop C increments k to 1, but since k < j (1 < 1) is false, it terminates
The rest of loop B runs. i == j is true, so L[i * n + j] is set to sqrt(A[i * n + i] - s)
loop B increments j to 2, but since j < (i+1) (2 < 2) is false, it terminates
loop A increments i to 2, and assuming i is less than n it runs
...

c syntax for choleskey Decomposition

1 Answers1