## 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

and

Note that *k*= 0*, . . . , n-*1;*L*

_{ik}= (

*A*

_{ik}- (

*L*

_{i0}

*U*

_{0k}+ · · · +

*L*

_{i,k-1}

*U*

_{k-1,k})), for

*i = k, . . . , n-*1,

*U*

_{kj}=

*A*

_{kj}- (

*L*

_{k0}

*U*

_{0j}+ · · · +

*L*

_{k,k-1}

*U*

_{k-1,j}) /

*L*

_{kk}, for

*j = k+*1

*, . . . , n-*1.

*L*

_{ik}= 0 for

*k > i*,

*U*

_{ik}= 0 for

*k < i*, and

*U*

_{kk}= 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*i*element of the array^{th}*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( ).