Cholesky decomposition: Difference between revisions

← Older editNewer edit →
Content deleted Content added
m →‎{{header|REXX}}: added/changed whitespace and comments, made other cosmetic improvements.
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
<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> aLinez=aLinez right( format(@.row.col,, decPlacesdPlaces) /1, width) </code>
<lang rexx>/*REXX program to performperforms the Cholesky decomposition on a square matrix. */
exitniner = '25 15 -5' , /*stickdefine a fork in3x3 it, we'rematrix. done.*/▼
 
niner= '25 '15 -518 0' ,
'15 18'-5 0 11' ,
'-5 0 11' call Cholesky niner
hexer = 18 22 54 42, call Cholesky niner /*define a 4x4 matrix. */
hexer= 18 22 5470 86 4262,
22 7054 86 174 62134,
54 8642 174 62 134, 106
42 62 134 106 call Cholesky hexer
exit call Cholesky hexer /*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──────────────*/
else do row=1 r=1 for orderord▼
▲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 if r=1c for orderthen !.r.r=sqrt(!.r.r-$)
do col=1 for row; s else !.r.c=01/!.c.c*(@.r.c-$)
end do i=1 for col-1/*c*/
end s=s+$.row.i/*r*$.col.i/
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
/*─────────────────────────────────ERR subroutine─────────────────────────────*/
/*─────────────────────────────────────TELL subroutine───&find the order*/
err: say; say; say '***error***!'; say; say arg(1); say; say; exit 13</lang>▼
tell: parse arg hdr,x,y,sep; #=0; if sep=='' then sep='═'
/*─────────────────────────────────TELL subroutine────────────────────────────*/
decPlaces = 5 /*number of decimal places past the decimal point. */▼
widthtell: parse arg hdr,x,y,sep; #=0; 10 /*width of field to be usedif tosep=='' display thethen elements*/sep='═'
▲decPlaces dPlaces= 5 /*number of# decimal places past the decimal point. */
 
width =10 /*width of field used to display elements*/
if y=='' then $!.=0
▲ else do row=1 for order
else do row=1 for ord; do col=1 for ord; x=x do !.row.col=1; forend; orderend
x=x $.row.col▼
end /*row*/▼
end /*col*/
w=words(x)
do orderord=1 until orderord**2>=w; end /*a fast way to find the MATmatrix's order. */▼
 
say
▲ do order=1 until order**2>=w /*fast way to find the MAT order.*/
if orderord**2\==w then call err "matrix elements don't matchform itsa ordersquare matrix."▼
end /*order*/
say; say center(hdr, ((width+1)*w)%orderord, sep); say▼
 
say
▲if order**2\==w then call err "matrix elements don't match its order"
▲ do row=1 for ord; xz=x $.row.col
▲say; say center(hdr, ((width+1)*w)%order, sep); say
▲ do col=1 for ord; end /*row*/#=#+1
 
do row=1 for order; aLine @.row.col=word(x,#)
doif col<=1row for order;then #!.row.col=#+1@.row.col
z=z right( format(@.row.col,, dPlaces) / 1, @.row.col=word(x,#width)
end if col<=row then $.row.col=@.row./*col*/
say aLinez▼
aLine=aLine right( format(@.row.col,, decPlaces) /1, width)
end /*col*/
▲ say aLine
end /*row*/
return
/*─────────────────────────────────SQRT subroutine────────────────────────────*/
/*─────────────────────────────────────SQRT subroutine──────────────────*/
sqrt: procedure; parse arg x; if x=0 then return 0; d=digits(); p=d+d%4+2; m.=11
numeric form; numeric digits 11m.; g=.sqrtGuessparse value format(x,2,1,,0);'E0' dowith j=0g while'E' p>9;_ m.j=p; p=p%2+1; end
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
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}}
<pre>
Line 2,076 ⟶ 2,060:
3 3 0
-1 1 3
 
 
 
════════════════input array═════════════════
Retrieved from "https://rosettacode.org/wiki/Cholesky_decomposition"