LU decomposition: Difference between revisions

→‎{{header|REXX}}: simplified the REXX program.
m (→‎Library gonum/matrix: library churn)
(→‎{{header|REXX}}: simplified the REXX program.)
Line 2,809:
#=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 makeMat call bldAMat; call showMat 'A' /*makebuild theand display A matrix from the numbers.*/
call showMat 'A', N call bldPmat; call showMat 'P' /*display the " " A matrix." P " */
call manPmat call multMat; call showMat 'PA' /* " " /*manufacture P " PA (permutation). " */
call showMat 'P', N do y=1 for N; call bldUmat; call bldLmat /*display the build P matrix.U and L " */
end /*y*/
call multMat /*multiply the A and P matrices. */
call showMat 'PA', N call showMat 'L'; call showMat 'U' /*display the PA matrix.L and U " */
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 all done. */
/*──────────────────────────────────────────────────────────────────────────────────────*/
erbldAMat: say?=words(x); do N=1 say '***error***';for ? until N**2>=? say; say arg(1); /*find matrix say;size. exit 13*/
end /*N*/
/*──────────────────────────────────────────────────────────────────────────────────────*/
makeMat: ?=words(x); do N=1 for ?; if N**2\==? then leavedo; say end /'*N*/*error*** wrong # of elements entered:' ?;exit 9;end
if N**2 \==? then call er 'not correct number of elements entered: ' ?
 
do r=1 for N /*build the "A" matrix from the 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
/*──────────────────────────────────────────────────────────────────────────────────────*/
manLmatbldLmat: parse arg ? do r=1 for N /*manufacture L ( /*build lower) matrix.*/
do c=1 for N; if r==c then do; L.r.c=1; foriterate; N end
doif c\=1=y for N; if| r==c | then do; L.c>r.c=1; then iterate; end
if c\_==? | PA.r==.c | c>r then iterate
do k=1 for c-1; _=PA_ - U.rk.c * L.r.k
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
/*──────────────────────────────────────────────────────────────────────────────────────*/
manPmatbldPmat: c=N; do r=N by -1 for N; P.r.c=1; c=c+1 /*manufacture P build (permutation)perm. matrix*/
P.r.c=1; c=c+1; if c>N then c=N % 2; if c==N then c=1
end /*r*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
manUmatbldUmat: parse arg ? do r=1 for N; if r\==y then iterate /*manufacture U build (upper) matrix.*/
do do rc=1 for N; if c<r\==? then iterate
do c_=1 for N; if PA.r.c<r then iterate
do k=1 for r-1; _=PA_ - U.rk.c * L.r.k
do k=1end 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 & 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 end /*j*/
end /*i*/; return
/*──────────────────────────────────────────────────────────────────────────────────────*/
showMat: parse arg mat,rows,cols; w=0 say; rows=word(rows N,1); cols=word(cols rows,1); say
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</lang>
{{out|output|text=&nbsp; when using the input of: &nbsp; &nbsp; <tt> 1 3 5 &nbsp; 2 4 7 &nbsp; 1 1 0 </tt>}}
<pre>