Download Here

1
A
=
x
b
Using a Gaussian elimination technique, A can be factored (decomposed) into
the following two matrices L and U:
2
L
=
1
0
3
0:4
42
05
1
2
3
0
and
U
=
1
40
0
1
2
3
3
05 05 5
0 08
After the matrix A of coefficients is factored into the lower and upper triangular
matrices L and U, values for the vector x of variables can be determined easily:
Since
A
1
x
=
b
and
A
=
L
1
U
, therefore
x
= (L
1 )01 1
U
b
or
x
=
U
01 1 L01 1 b
In effect, the application of U 01 and L01 is performed by the processes of
forward elimination (for L01 ) and back substitution (for U 01 ). Consequently,
the computation is easily done in two steps, by:
•
Calculating an intermediate vector equal to
elimination
•
Applying
x.
U
L
01 1 b using forward
01 to this vector by back substitution to yield the solution vector
Once A has been factored into L and U, this two-step procedure can be used
repeatedly to solve the system of equations for different values of b.
B.2 Coding the Algorithm
A standard algorithm for LU decomposition, described in Numerical Recipes,1
transforms a square matrix ‘‘in place,’’ storing the values for all the elements
of L and U in the same space in memory where the original square matrix was
stored. This can be done by overlapping the two arrays so that the mandatory
zeros on the opposite sides of both L and U, and the ones on the diagonal of L,
are not explicitly represented in memory. The algorithm transforms the array
A in the previous example into the following array:
2
1
42
3
1
2
05
0:4
3
3
05 5
08
William H. Press [et al.], Numerical Recipes in FORTRAN : The Art of Scientific
Computing. 2nd ed. Cambridge University Press, 1992.
B–2 HPF Tutorial: LU Decomposition