Editing Crout matrix decomposition

{{Short description|Type of matrix factorization}}
{{One source|date=September 2024}}
In [[linear algebra]], the '''Crout matrix decomposition''' is an [[LU decomposition]] which decomposes a [[matrix (mathematics)|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]]. <ref name=Press2007>{{cite book|last=Press|first=William H.|title=Numerical Recipes 3rd Edition: The Art of Scientific Computing|year=2007|publisher=Cambridge University Press|isbn=9780521880688|pages=50–52}}</ref>

The Crout matrix decomposition [[algorithm]] differs slightly from the [[Doolittle decomposition|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:

<syntaxhighlight lang="c">
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];
		}
	}
}
</syntaxhighlight>

Octave/Matlab implementation:
<syntaxhighlight lang="matlab">
   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
</syntaxhighlight>

==References==
{{reflist}}
*[https://web.archive.org/web/20131003135132/http://tutorialfrenzy.com/matlab-crout-lu-factorization-linear-equations-systems.php  Implementation using functions ] In Matlab

{{DEFAULTSORT:Crout Matrix Decomposition}}
[[Category:Matrix decompositions]]