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

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

mail