Systems of Linear Equations



Crout Decomposition

Crout LU Decomposition

Crout's method decomposes a nonsingular n × n matrix A into the product of an n×n lower triangular matrix L and an n×n unit upper triangular matrix U. A unit triangular matrix is a triangular matrix with 1's along the diagonal. Crout's algorithm proceeds as follows:
Evaluate the following pair of expressions for k = 0, . . . , n-1;
Lik = ( Aik - (Li0 U0k + · · · + Li,k-1 Uk-1,k)),   for i = k, . . . , n-1,
and
Ukj = Akj - (Lk0 U0j + · · · + Lk,k-1 Uk-1,j) / Lkk,   for j = k+1, . . . , n-1.
Note that Lik = 0 for k > i, Uik = 0 for k < i, and Ukk = 1 for k = 0, … , n-1. The matrix U forms a unit upper triangular matrix, and the matrix L forms a lower triangular matrix. The matrix
A = LU.
After the LU decomposition of A is performed, the solution to the system of linear equations A x = L U x = B is solved by solving the system of linear equations L y = B by forward substitution for y, and then solving the system of linear equations U x = y by backward substitution for x.

Crout's LU decomposition with pivoting is similar to the above algorithm except that for each k a pivot row is determined and interchanged with row k, the algorithm then proceeds as before. Source code is provided for the two different versions of Crout's LU decomposition, one version performs pivoting and the other version does not. If the matrix A is positive definite symmetric or if the matrix is diagonally dominant, then pivoting is not necessary; otherwise the version using pivoting should be used.

Function List

  • int Crout_LU_Decomposition( double* A, int n )

    Given the n×n matrix A, Crout_LU_Decomposition uses Crout's algorithm to decompose A into the product of a lower triangular matrix and a unit upper triangular matrix. The lower triangular part of the lower triangular matrix is returned in the lower triangular part of A and the non-diagonal upper triangular part of the unit upper triangular matrix is returned in the upper triangular part of A.  Crout_LU_Decomposition returns 0 if the decomposition was successful and returns -1 if the matrix is singular.

  • int Crout_LU_Solve( double *LU, double* B, double* x, int n )

    Crout_LU_Solve solves the system of linear equations LU x = B for x, where LU is the LU decomposition of A returned from Crout_LU_Decomposition.  Crout_LU_Solve returns 0 if the solution was found and returns -1 if the matrix is singular.

  • int Crout_LU_Decomposition_with_Pivoting( double *A, int pivot[ ], int n )

    Given the n×n matrix A, Crout_LU_Decomposition_with_Pivoting uses Crout's algorithm with pivoting to decompose a row interchanged version of A into the product of a lower triangular matrix and a unit upper triangular matrix. The lower triangular part of the lower triangular matrix is returned in the lower triangular part of A and the non-diagonal upper triangular part of the unit upper triangular matrix is returned in the upper triangular part of A. Upon completion the ith element of the array pivot contains the row interchanged with row i when k = i, k being the k in the description of the algorithm above. The array pivot should be dimensioned at least n in the calling routine.  Crout_LU_Decomposition returns 0 if the decomposition was successful and returns -1 if the matrix is singular.

  • int Crout_LU_with_Pivoting_Solve( double *LU, double B[ ], int pivot[ ], double x[ ], int n )

    Crout_LU_with_Pivoting_Solve solves the system of linear equations LU x = B for x, where LU is the LU decomposition of A returned from Crout_LU_Decomposition_with_Pivoting and pivot is the array of rows determined by Crout_LU_Decompositon_with_Pivoting.  Crout_LU_with_Pivoting_Solve returns 0 if the solution was found and returns -1 if the matrix is singular.

C Source Code

  • The file, crout.c, contains the routines Crout_LU_Decomposition( ) and Crout_LU_Solve( ).

  • The file, crout_pivot.c, contains the routines
    Crout_LU_Decomposition_with_Pivoting( ) and Crout_LU_with_Pivoting_Solve( ).