0
$\begingroup$

Consider the following code, found on wikipedia, that implements the Crout decomposition algorithm:

    void crout(double const **A, double **L, double **U, int n) {
    int i, j, k;
    double sum = 0;

    for (i = 0; i < n; i++) {
        U[i][i] = 1;
    }

    for (j = 0; j < n; j++) {
        for (i = j; i < n; i++) {
            sum = 0;
            for (k = 0; k < j; k++) {
                sum = sum + L[i][k] * U[k][j];  
            }
            L[i][j] = A[i][j] - sum;
        }

        for (i = j; i < n; i++) {
            sum = 0;
            for(k = 0; k < j; k++) {
                sum = sum + L[j][k] * U[k][i];
            }
            if (L[j][j] == 0) {
                printf("det(L) close to 0!\n Can't divide by 0...\n");
                exit(EXIT_FAILURE);
            }
            U[j][i] = (A[j][i] - sum) / L[j][j];
        }
    }
}

Which part is responsible for pivoting?

The thing is, I've read that partial pivoting is still a rather non-trivial operation (on the order $O(n^2)$), but the above algorithm seems similar in terms of complexity to the "naive" LU decomposition without any pivoting.

To check, I've run the algorithm on my 28x28 test matrix

[ 0 0.291601633 0 -0.062262937 0 0 0 0 0 -0.22949092 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0.291601633 0 -0.062262937 0 0 0 0 0 -0.22949092 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0 -0.062262937 0 0.633860174 0 -0.203470344 0 -0.185648901 0 -0.182968732 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ -0.062262937 0 0.633860174 0 -0.203470344 0 -0.185648901 0 -0.182968732 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0 0 0 -0.203470344 0 0.386208493 0 -0.18368887 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0 0 -0.203470344 0 0.386208493 0 -0.18368887 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0 0 0 -0.185648901 0 -0.18368887 0 1.163032501 0 -0.044175498 0 0 0 -0.20912 0 0 0 -0.55618 0 0 0 0 0 0 0 0 0 0 ]
[ 0 0 -0.185648901 0 -0.18368887 0 1.163032501 0 -0.044175498 0 0 0 -0.20912 0 0 0 -0.55618 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0 -0.22949092 0 -0.182968732 0 0 0 -0.044175498 0 0.702145321 0 -0.25202 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ -0.22949092 0 -0.182968732 0 0 0 -0.044175498 0 0.702145321 0 -0.25202 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 0 0 0 -0.25202 0 0.883882188 0 0 0 0 0 0 0 0 0 -0.220414179 0 -0.28380561 0 -0.14610303 0 0 ]
[ 0 0 0 0 0 0 0 0 -0.25202 0 0.883882188 0 0 0 0 0 0 0 0 0 -0.220414179 0 -0.28380561 0 -0.14610303 0 0 0 ]
[ 0 0 0 0 0 0 0 -0.20912 0 0 0 0 0 0.49528 0 -0.17615 0 -0.11001 0 0 0 0 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 -0.20912 0 0 0 0 0 0.49528 0 -0.17615 0 -0.11001 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.17615 0 0.17615 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 0 0 0 0 0 0 -0.17615 0 0.17615 0 0 0 0 0 0 0 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 0 -0.55618 0 0 0 0 0 -0.11001 0 0 0 1.033363204 0 -0.090289125 0 0 0 0 0 0 0 -0.298767964 ]
[ 0 0 0 0 0 0 -0.55618 0 0 0 0 0 -0.11001 0 0 0 1.033363204 0 -0.090289125 0 0 0 0 0 0 0 -0.298767964 0 ]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.090289125 0 0.299090395 0 -0.208861407 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.090289125 0 0.299090395 0 -0.208861407 0 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 0 0 0 0 0 -0.220414179 0 0 0 0 0 0 0 -0.208861407 0 0.429181968 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 0 0 0 0 -0.220414179 0 0 0 0 0 0 0 -0.208861407 0 0.429181968 0 0 0 0 0 0 0 ]
[ 0 0 0 0 0 0 0 0 0 0 0 -0.28380561 0 0 0 0 0 0 0 0 0 0 0 0.570852385 0 -0.29792224 0 0 ]
[ 0 0 0 0 0 0 0 0 0 0 -0.28380561 0 0 0 0 0 0 0 0 0 0 0 0.570852385 0 -0.29792224 0 0 0 ]
[ 0 0 0 0 0 0 0 0 0 0 0 -0.14610303 0 0 0 0 0 0 0 0 0 0 0 -0.29792224 0 0.818338896 0 -0.387730558 ]
[ 0 0 0 0 0 0 0 0 0 0 -0.14610303 0 0 0 0 0 0 0 0 0 0 0 -0.29792224 0 0.818338896 0 -0.387730558 0 ]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.298767964 0 0 0 0 0 0 0 -0.387730558 0 0.68647389 ]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -0.298767964 0 0 0 0 0 0 0 -0.387730558 0 0.68647389 0 ]

which is invertible (you can check here), when it returned a number NaNs and Infinities.

$\endgroup$
1
$\begingroup$

That algorithm doesn't do any pivoting. It requires non-zero elements on the diagonal of A to work, otherwise you divide by zero. Note the check:

if (L[j][j] == 0) {
    printf("det(L) close to 0!\n Can't divide by 0...\n");
    exit(EXIT_FAILURE);
}
$\endgroup$

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.