Anonymous user
LU decomposition: Difference between revisions
→{{header|REXX}}: added the REXX language. -- ~~~~
(→{{header|Go}}: added library based solution) |
(→{{header|REXX}}: added the REXX language. -- ~~~~) |
||
Line 1,401:
[0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 0.0, 1.0]]</pre>
=={{header|REXX}}==
<lang rexx>/*REXX pgm makes a matrix from input, performs/shows LU decomposition. */
#=0; P.=0; PA.=0; L.=0; U.=0 /*initialize some variables to 0.*/
parse arg x /*get the matrix elements from CL*/
call makeMat /*make the A matrix from numbers.*/
call showMat 'A', N /*display the A matrix. */
call manPmat /*manufacture P (permutation).*/
call showMat 'P', N /*display the P matrix. */
call multMat /*multiply the A and P matrices.*/
call showMat 'PA', N /*display the PA matrix. */
do y=1 for N; call manUmat y /*manufacture U matrix, parts*/
call manLmat y /*manufacture L matrix, parts*/
end
call showMat 'L', N /*display the L matrix. */
call showMat 'U', N /*display the U matrix. */
exit /*stick a fork in it, we're done.*/
/*──────────────────────────────────er subroutine───────────────────────*/
er: say; say '***error!***'; say; say arg(1); say; exit 13
/*──────────────────────────────────makeMat subroutine──────────────────*/
makeMat: ?=words(x); do N=1 for ?; if N**2==? then leave; end
if N**2\==? then call er 'not correct number of elements entered: ' ?
do r=1 for N /*build the "A" matrix from input*/
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
/*──────────────────────────────────manLmat subroutine──────────────────*/
manLmat: arg ? /*manufacture L (lower) matrix.*/
do r=1 for N
do c=1 for N; if r==c then do; L.r.c=1; iterate; end
if c\==? | 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
/*──────────────────────────────────manPmat subroutine──────────────────*/
manPmat: c=N; do r=N by -1 for N /*manufacture P (permutation). */
P.r.c=1; c=c+1; if c>N then c=N%2; if c==N then c=1
end /*r*/
return
/*──────────────────────────────────manUmat subroutine──────────────────*/
manUmat: arg ? /*manufacture U (upper) matrix.*/
do r=1 for N; if r\==? then iterate
do c=1 for N; if c<r then iterate
_=PA.r.c
do k=1 for r-1; _=_-U.k.c*L.r.k; end /*k*/
U.r.c=_
end /*c*/
end /*r*/
return
/*──────────────────────────────────multMat subroutine──────────────────*/
multMat: do i =1 for N /*multiply matrix P & 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
end /*j*/
end /*i*/
return
/*──────────────────────────────────showMat subroutine──────────────────*/
showMat: parse arg mat,rows,cols; w=0; cols=word(cols rows,1); say
do r =1 for rows
do c=1 for cols; w=max(w,length(value(mat'.'r'.'c))); end
end
say center(mat 'matrix',cols*(w+1)+7,"─")
do r =1 for rows; _=
do c=1 for cols; _=_ right(value(mat'.'r'.'c),w+1); end; say _
end
return</lang>
'''output''' when using the input of: <tt> 1 3 5 2 4 7 1 1 0 </tt>
<pre style="overflow:scroll">
──A matrix───
1 3 5
2 4 7
1 1 0
──P matrix───
0 1 0
1 0 0
0 0 1
──PA matrix──
2 4 7
1 3 5
1 1 0
─────L matrix──────
1 0 0
0.5 1 0
0.5 -1 1
───────U matrix───────
2 4 7
0 1.0 1.5
0 0 -2.0
</pre>
'''output''' when using the input of: <tt> 11 9 24 2 1 5 2 6 3 17 18 1 2 5 7 1 </tt>
<pre style="overflow:scroll">
─────A matrix──────
11 9 24 2
1 5 2 6
3 17 18 1
2 5 7 1
───P matrix────
1 0 0 0
0 0 1 0
0 1 0 0
0 0 0 1
─────PA matrix─────
11 9 24 2
3 17 18 1
1 5 2 6
2 5 7 1
───────────────────────────L matrix────────────────────────────
1 0 0 0
0.272727273 1 0 0
0.0909090909 0.287500001 1 0
0.181818182 0.23125 0.00359712804 1
───────────────────────U matrix────────────────────────
11 9 24 2
0 14.5454545 11.4545455 0.45454545
0 0 -3.47500002 5.68750000
0 0 0 0.510791339
</pre>
=={{header|Ruby}}==
| |||