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).^[1]

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.

being L a lower triangular matrix, D a diagonal matrix and U a normalised upper triangular matrix

C implementation:

void crout (double **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

References

• Implementation using functions In Matlab