Cholesky decomposition: Difference between revisions

Content added Content deleted
m (→‎{{header|REXX}}: added/changed whitespace and comments, made other cosmetic improvements.)
Line 1,996: Line 1,996:
If trailing zeroes are wanted after the decimal point for values of zero (0), &nbsp; the &nbsp; &nbsp; <big>'''/1'''</big> &nbsp; &nbsp; (division by one performs
If trailing zeroes are wanted after the decimal point for values of zero (0), &nbsp; the &nbsp; &nbsp; <big>'''/1'''</big> &nbsp; &nbsp; (division by one performs
<br>REXX number normalization) &nbsp; can be removed from the line &nbsp; (number 73) &nbsp; which contains the source statement:
<br>REXX number normalization) &nbsp; can be removed from the line &nbsp; (number 73) &nbsp; which contains the source statement:
<br> &nbsp; &nbsp; &nbsp; <code> aLine=aLine right( format(@.row.col,, decPlaces) /1, width) </code>
<br> &nbsp; &nbsp; &nbsp; <code> z=z right( format(@.row.col,, dPlaces) /1, width) </code>
<lang rexx>/*REXX program to perform the Cholesky decomposition on square matrix.*/
<lang rexx>/*REXX program performs the Cholesky decomposition on a square matrix. */
⚫
niner = '25 15 -5' , /*define a 3x3 matrix. */

niner= '25 15 -5' ,
'15 18 0' ,
'15 18 0' ,
'-5 0 11'
'-5 0 11'
call Cholesky niner
call Cholesky niner
hexer = 18 22 54 42, /*define a 4x4 matrix. */
hexer= 18 22 54 42,
22 70 86 62,
22 70 86 62,
54 86 174 134,
54 86 174 134,
42 62 134 106
42 62 134 106
call Cholesky hexer
call Cholesky hexer
exit /*stick a fork in it, we're all done. */
/*──────────────────────────────────Cholesky subroutine───────────────────────*/
⚫
exit /*stick a fork in it, we're done.*/
⚫
Cholesky: procedure; parse arg title; say; say; call tell 'input array',title
/*─────────────────────────────────────Cholesky subroutine──────────────*/
⚫
do r=1 for ord
⚫
Cholesky: procedure; arg !; call tell 'input array',!
⚫
do c=1 for r; $=0; do i=1 for c-1; $=$+!.r.i*!.c.i; end /*i*/

do row=1 for order
if r=c then !.r.r=sqrt(!.r.r-$)
do col=1 for row; s=0
else !.r.c=1/!.c.c*(@.r.c-$)
do i=1 for col-1
end /*c*/
s=s+$.row.i*$.col.i
end /*r*/
⚫
call tell 'Cholesky factor',,!.,'─'
⚫
end /*i*/
if row=col then $.row.row=sqrt($.row.row-s)
else $.row.col=1/$.col.col*(@.row.col-s)
end /*col*/
end /*row*/

⚫
call tell 'Cholesky factor',,$.,'─'
return
return
/*─────────────────────────────────ERR subroutine─────────────────────────────*/
/*─────────────────────────────────────TELL subroutine───&find the order*/
⚫
err: say; say; say '***error***!'; say; say arg(1); say; say; exit 13
tell: parse arg hdr,x,y,sep; #=0; if sep=='' then sep='═'
/*─────────────────────────────────TELL subroutine────────────────────────────*/
⚫
decPlaces = 5 /*number of decimal places past the decimal point. */
width = 10 /*width of field to be used to display the elements*/
tell: parse arg hdr,x,y,sep; #=0; if sep=='' then sep='═'
⚫
dPlaces= 5 /*# decimal places past the decimal point*/

width =10 /*width of field used to display elements*/
if y=='' then $.=0
if y=='' then !.=0
⚫
else do row=1 for order
do col=1 for order
else do row=1 for ord; do col=1 for ord; x=x !.row.col; end; end
⚫
x=x $.row.col
⚫
end /*row*/
end /*col*/
w=words(x)
w=words(x)
⚫
do ord=1 until ord**2>=w; end /*a fast way to find the matrix's order. */

say
⚫
do order=1 until order**2>=w /*fast way to find the MAT order.*/
⚫
if ord**2\==w then call err "matrix elements don't form a square matrix."
end /*order*/
⚫
say center(hdr, ((width+1)*w)%ord, sep)

say
⚫
if order**2\==w then call err "matrix elements don't match its order"
⚫
do row=1 for ord; z=
⚫
say; say center(hdr, ((width+1)*w)%order, sep); say
⚫
do col=1 for ord; #=#+1

do row=1 for order; aLine=
@.row.col=word(x,#)
do col=1 for order; #=#+1
if col<=row then !.row.col=@.row.col
@.row.col=word(x,#)
z=z right( format(@.row.col,, dPlaces) / 1, width)
if col<=row then $.row.col=@.row.col
end /*col*/
⚫
say z
aLine=aLine right( format(@.row.col,, decPlaces) /1, width)
end /*col*/
⚫
say aLine
end /*row*/
end /*row*/
return
return
/*─────────────────────────────────SQRT subroutine────────────────────────────*/
/*─────────────────────────────────────SQRT subroutine──────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits()
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); p=d+d%4+2; m.=11
numeric digits 11; g=.sqrtGuess(); do j=0 while p>9; m.j=p; p=p%2+1; end
numeric form; numeric digits m.; parse value format(x,2,1,,0)'E0' with g 'E' _ .
g=g*.5'E'_%2; do j=0 while p>9; m.j=p; p=p%2+1; end
do k=j+5 to 0 by -1;if m.k>11 then numeric digits m.k;g=.5*(g+x/g);end
do k=j+5 to 0 by -1; if m.k>11 then numeric digits m.k; g=.5*(g+x/g); end
numeric digits d; return g/1
numeric digits d; return g/1</lang>

.sqrtGuess: if x<0 then call err 'SQRT of negative #'; numeric form
m.=11; p=d+d%4+2; parse value format(x,2,1,,0) 'E0' with g 'E' _ .
return g*.5'E'_%2
/*─────────────────────────────────────ERR subroutine───────────────────*/
⚫
err: say; say; say '***error***!'; say; say arg(1); say; say; exit 13</lang>
{{out}}
{{out}}
<pre>
<pre>
Line 2,076: Line 2,060:
3 3 0
3 3 0
-1 1 3
-1 1 3




════════════════input array═════════════════
════════════════input array═════════════════
Retrieved from "https://rosettacode.org/wiki/Cholesky_decomposition"