Crout matrix decomposition

In linear algebra, the Crout matrix decomposition is an LU decomposition which decomposes a matrix into a lower triangular matrix (L), an upper triangular matrix (U) and, although not always needed, a permutation matrix (P). It was developed by Prescott Durand Crout.

The Crout matrix decomposition algorithm differs slightly from the Doolittle method. Doolittle's method returns a unit lower triangular matrix and an upper triangular matrix, while the Crout method returns a lower triangular matrix and a unit upper triangular matrix.

So, if a matrix decomposition of a matrix A is such that:

A = LDU

being L a unit lower triangular matrix, D a diagonal matrix and U a unit upper triangular matrix, then Doolittle's method produces

A = L(DU)

and Crout's method produces

A = (LD)U.

Implementations

C implementation:

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];		}	}}

Octave/Matlab implementation:

   function [L, U] = LUdecompCrout(A)                [R, C] = size(A);        for i = 1:R            L(i, 1) = A(i, 1);            U(i, i) = 1;        end        for j = 2:R            U(1, j) = A(1, j) / L(1, 1);        end        for i = 2:R            for j = 2:i                L(i, j) = A(i, j) - L(i, 1:j - 1) * U(1:j - 1, j);            end                        for j = i + 1:R                U(i, j) = (A(i, j) - L(i, 1:i - 1) * U(1:i - 1, j)) / L(i, i);            end        end   end

References

• Implementation using functions In Matlab