Systems of Linear Equations



Doolittle Decomposition

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 k = 0, … , n-1;
Ukj = Akj - (Lk0 U0j + · · · + Lk,k-1 Uk-1,j),   for j = k, … , n-1.
and
Lik = ( Aik - (Li0 U0k + · · · + Li,k-1 Uk-1,k)) / Ukk,   for i = k+1, … , n-1,
Note that Lik = 0 for k > i, Uik = 0 for k < i, and Lkk = 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 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.  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( ).