How to resolve the algorithm LU decomposition step by step in the REXX programming language
How to resolve the algorithm LU decomposition step by step in the REXX programming language
Table of Contents
Problem Statement
Every square matrix
A
{\displaystyle A}
can be decomposed into a product of a lower triangular matrix
L
{\displaystyle L}
and a upper triangular matrix
U
{\displaystyle U}
, as described in LU decomposition. It is a modified form of Gaussian elimination. While the Cholesky decomposition only works for symmetric, positive definite matrices, the more general LU decomposition works for any square matrix. There are several algorithms for calculating L and U. To derive Crout's algorithm for a 3x3 example, we have to solve the following system: We now would have to solve 9 equations with 12 unknowns. To make the system uniquely solvable, usually the diagonal elements of
L
{\displaystyle L}
are set to 1 so we get a solvable system of 9 unknowns and 9 equations. Solving for the other
l
{\displaystyle l}
and
u
{\displaystyle u}
, we get the following equations: and for
l
{\displaystyle l}
: We see that there is a calculation pattern, which can be expressed as the following formulas, first for
U
{\displaystyle U}
and then for
L
{\displaystyle L}
We see in the second formula that to get the
l
i j
{\displaystyle l_{ij}}
below the diagonal, we have to divide by the diagonal element (pivot)
u
j j
{\displaystyle u_{jj}}
, so we get problems when
u
j j
{\displaystyle u_{jj}}
is either 0 or very small, which leads to numerical instability. The solution to this problem is pivoting
A
{\displaystyle A}
, which means rearranging the rows of
A
{\displaystyle A}
, prior to the
L U
{\displaystyle LU}
decomposition, in a way that the largest element of each column gets onto the diagonal of
A
{\displaystyle A}
. Rearranging the rows means to multiply
A
{\displaystyle A}
by a permutation matrix
P
{\displaystyle P}
: Example: The decomposition algorithm is then applied on the rearranged matrix so that
The task is to implement a routine which will take a square nxn matrix
A
{\displaystyle A}
and return a lower triangular matrix
L
{\displaystyle L}
, a upper triangular matrix
U
{\displaystyle U}
and a permutation matrix
P
{\displaystyle P}
, so that the above equation is fulfilled. You should then test it on the following two examples and include your output.
Let's start with the solution:
Step by Step solution about How to resolve the algorithm LU decomposition step by step in the REXX programming language
Source code in the rexx programming language
/*REXX program creates a matrix from console input, performs/shows LU decomposition.*/
#= 0; P.= 0; PA.= 0; L.= 0; U.= 0 /*initialize some variables to zero. */
parse arg x /*obtain matrix elements from the C.L. */
call bldAMat; call showMat 'A' /*build and display A matrix.*/
call bldPmat; call showMat 'P' /* " " " P " */
call multMat; call showMat 'PA' /* " " " PA " */
do y=1 for N; call bldUmat; call bldLmat /*build U and L " */
end /*y*/
call showMat 'L'; call showMat 'U' /*display L and U " */
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
bldAMat: ?= words(x); do N=1 for ? until N**2>=? /*find matrix size. */
end /*N*/
if N**2\==? then do; say '***error*** wrong # of elements entered:' ?; exit 9
end
do r=1 for N /*build A matrix.*/
do c=1 for N; #= # + 1; _= word(x, #); A.r.c= _
if \datatype(_, 'N') then call er "element isn't numeric: " _
end /*c*/
end /*r*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
bldLmat: do r=1 for N /*build lower matrix.*/
do c=1 for N; if r==c then do; L.r.c= 1; iterate; end
if c\==y | r==c | c>r then iterate
_= PA.r.c
do k=1 for c-1; _= _ - U.k.c * L.r.k
end /*k*/
L.r.c= _ / U.c.c
end /*c*/
end /*r*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
bldPmat: c= N; do r=N by -1 for N; P.r.c= 1; c= c + 1 /*build perm. matrix.*/
if c>N then c= N%2; if c==N then c= 1
end /*r*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
bldUmat: do r=1 for N; if r\==y then iterate /*build upper matrix.*/
do c=1 for N; if c
_= PA.r.c
do k=1 for r-1; _= _ - U.k.c * L.r.k
end /*k*/
U.r.c= _ / 1
end /*c*/
end /*r*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
multMat: do i=1 for N /*multiply matrix P and A ──► PA */
do j=1 for N
do k=1 for N; pa.i.j= (pa.i.j + p.i.k * a.k.j) / 1
end /*k*/
end /*j*/ /*÷ by one does normalization [↑]. */
end /*i*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
showMat: parse arg mat,rows,cols; say; rows= word(rows N,1); cols= word(cols rows,1)
w= 0; do r=1 for rows
do c=1 for cols; w= max(w, length( value( mat'.'r"."c ) ) )
end /*c*/
end /*r*/
say center(mat 'matrix', cols * (w + 1) + 7, "─") /*display the header.*/
do r=1 for rows; _=
do c=1 for cols; _= _ right( value(mat'.'r"."c), w + 1)
end /*c*/
say _
end /*r*/; return
You may also check:How to resolve the algorithm Variadic function step by step in the ActionScript programming language
You may also check:How to resolve the algorithm Mandelbrot set step by step in the Huginn programming language
You may also check:How to resolve the algorithm Horner's rule for polynomial evaluation step by step in the Java programming language
You may also check:How to resolve the algorithm Sum digits of an integer step by step in the Dart programming language
You may also check:How to resolve the algorithm 99 bottles of beer step by step in the Dart programming language