LU - decomposition

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).

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).

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

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
  if Max = 0 then
    call ERROR "LUDCMP: Error - singular matrix"
  Vv.I = 1 / Max
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
    A.I.J = Sum
  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
    A.I.J = Sum
    W = Vv.I * ABS(Sum)
    if W >= Max then do; Max = W; Imax = I; 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
    D = -D; Vv.Imax = Vv.J
  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
return D
ERROR: say ARG(1); exit


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

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

last modified 8th August 2001
Copyright 2000-2001 Vladimir Zabrodsky
Czech Republic.