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