## Doolittle LU Decomposition

Doolittle's method decomposes a nonsingular*n×n*matrix

**A**into the product of an

*n×n*unit lower triangular matrix

**L**and an

*n × n*upper triangular matrix

**U**. A unit triangular matrix is a triangular matrix with 1's along the diagonal. Doolittle's algorithm proceeds as follows:

Evaluate the following pair of expressions for

and

Note that *k*= 0*, … , n-*1;*U*

_{kj}=

*A*

_{kj}- (

*L*

_{k0}

*U*

_{0j}+ · · · +

*L*

_{k,k-1}

*U*

_{k-1,j}), for

*j = k*, … ,

*n-*1.

*L*

_{ik}= (

*A*

_{ik}- (

*L*

_{i0}

*U*

_{0k}+ · · · +

*L*

_{i,k-1}

*U*

_{k-1,k})) /

*U*

_{kk}, for

*i = k*+1

*, … , n-*1,

*L*

_{ik}= 0 for

*k > i*,

*U*

_{ik}= 0 for

*k < i*, and

*L*

_{kk}= 1 for

*k =*0

*, … , n-*1. The matrix

**U**forms an upper triangular matrix, and the matrix

**L**forms a unit 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**.

Doolittle'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 Doolittle'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 Doolittle_LU_Decomposition( double* A, int n )

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

- int Doolittle_LU_Solve( double *LU, double* B, double* x, int n )

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

- int Doolittle_LU_Decomposition_with_Pivoting( double *A, int pivot[ ], int n )

Given the*n×n*matrix*A*, Doolittle_LU_Decomposition_with_Pivoting uses Doolittle's algorithm with pivoting to decompose a row interchanged version of*A*into the product of a unit lower triangular matrix and an upper triangular matrix. The non-diagonal lower triangular part of the unit lower triangular matrix is returned in the lower triangular part of*A*and the upper triangular part of the 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. Doolittle_LU_Decomposition returns 0 if the decomposition was successful and returns -1 if the matrix is singular.

- int Doolittle_LU_with_Pivoting_Solve( double *LU, double B[ ], int pivot[ ], double x[ ], int n )

Doolittle_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 Doolittle_LU_Decomposition_with_Pivoting and*pivot*is the array of rows determined by Doolittle_LU_Decompositon_with_Pivoting. Doolittle_LU_with_Pivoting_Solve returns 0 if the solution was found and returns -1 if the matrix is singular.

*C* Source Code

- The file, doolittle.c, contains the routines Doolittle_LU_Decomposition( ) and Doolittle_LU_Solve( ).

- The file, doolittle_pivot.c, contains the routines

Doolittle_LU_Decomposition_with_Pivoting( ) and

Doolittle_LU_with_Pivoting_Solve( ).