Crout matrix decomposition
 | This article relies largely or entirely on a single source. (September 2024) |
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. [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.
Implementations
[edit]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
[edit]- ^ Press, William H. (2007). Numerical Recipes 3rd Edition: The Art of Scientific Computing. Cambridge University Press. pp. 50–52. ISBN 9780521880688.
- Implementation using functions In Matlab