LU - decomposition

PROBLEM
We wish to write a square N by N matrix A. as product of two matrices L. (lower triangular - has elements only on the diagonal and below) and U. (upper triangular - has elements only on the diagonal and above).

ALGORITHM
Crout's algorithm is the decomposition in place. We suppose L.J.J=1 for J=1,...,N; the other elements of L. will in A. below the diagonal. Elements of U. will in A. on the diagonal and above. The following algorithm doesn't actually decompose the matrix A. into LU form; it rather decompose a rowwise permutation of input matrix. This permutation is recorded in the Indx. output vector. The LU decomposition in LUDCMP requires about (1/3)*N**3 executions of the inner loops (each with one multiply and one add).

IMPLEMENTATION
Unit: internal function

Global variables: input N by N matrix A., Indx. is an N element vector

Parameter: a positive integer N

Returns: +1 or -1 depending on whether the number of row interchanges was even or odd (for calculation of a determinant)

Result: Indx. records the row permutation; a rowwise permutation of A. is decomposed

Note:
W. H. Press et al. say: "If the pivot element is zero the matrix is singular. For some applications on singular matrices, it is desirable to substitute Tiny for zero." They recommend the value Tiny=1E-20.

 LUDCMP: procedure expose A. Indx. parse arg N D = 1; Tiny = 1E-20 do I = 1 to N   Max = 0   do J = 1 to N     W = ABS(A.I.J)     if W > Max then Max = W   end   if Max = 0 then     call ERROR "LUDCMP: Error - singular matrix"   Vv.I = 1 / Max end do J = 1 to N   do I = 1 to J - 1     Sum = A.I.J     do K = 1 to I - 1       Sum = Sum - A.I.K * A.K.J     end     A.I.J = Sum   end   Max = 0   do I = J to N     Sum = A.I.J     do K = 1 to J - 1       Sum = Sum - A.I.K * A.K.J     end     A.I.J = Sum     W = Vv.I * ABS(Sum)     if W >= Max then do; Max = W; Imax = I; end   end   if J <> Imax then do     do K = 1 to N       W = A.Imax.K; A.Imax.K = A.J.K; A.J.K = W     end     D = -D; Vv.Imax = Vv.J   end   Indx.J = Imax   if A.J.J = 0 then A.J.J = Tiny   if J <> N then do     W = 1 / A.J.J     do I = J + 1 to N       A.I.J = A.I.J * W     end   end end return D   ERROR: say ARG(1); exit

CONNECTIONS
Linear Algebraic Equations
Determinant of a matrix
Inverse of a matrix
Solution of linear algebraic equations

Literature
Press W.H., Teukolsky S.A., Vetterling W.T., Flannery B.P. Numerical Recipes in C : the art of scientific computing
- 2nd ed. University Press, Cambridge, 1992

 cover contents index main page