# HP Forums

Full Version: (42S) Cholesky decomposition and solve
You're currently viewing a stripped down version of our content. View the full version with proper formatting.
0.Abstract

The listed routines allow for decomposing a positive definite matrix A into:

A = U'*U

Where ' means real transpose or complex hermitian.
Subsequently, a system of linear equations A*X=B may be solved.

The routines use only the stack, and only the upper triangular part of A and U
is referenced.
A and U are NxN matrices, real or complex. In fact, U overwrites A.
For accuracy reasons, matrix products have been used when possible, as
they use higher intermediate precision on a real 42S.

1."u'u" - decomposition A=U'U, U upper triangular

usage:
X Y
In : A
Out: INFO U

A and U are NxN real or complex matrices. U overwrites A (so A is lost).
INFO=0 when ok, INFO=i when A was found not to be positive definite at
step i.
Stack levels Z and T are lost, but Y is preserved.

2."u/" and "u'/" - solve U*X=B and U'*X=B, respectively

usage:
X Y
In : U B
Out: U X

U is the result of the decomposition with "u'u".
B and X are NxM matrices. If U is complex, then so must B.

3."PD/" - all in one solve A*X=B with A positive definite

usage:
X Y
In : A B
Out: X

A is NxN, real or complex.
B and X are NxM. If A is complex and B is real, X will be complex.
When A is not positive definite, execution stops with the message "Not Pos Def"

4.Source code

Code:
```00 { 332-Byte Prgm } 01>LBL "u'u" 02 EDIT 03 Rv 04>LBL 10 05 I- 06 RCLIJ 07 I+ 08 CLX 09 RCLIJ 10 STO- ST Z 11 X<> ST Z 12 + 13 GETM 14 RCLIJ 15 XEQ 09 16 +/- 17 X<>Y 18 EDIT 19 Rv 20 DIM? 21 STO ST Y 22 STOIJ 23 X<> ST T 24 LASTX 25 R^ 26 x 27 PUTM 28 RCLEL 29 REAL? 30 GTO 00 31 COMPLEX 32 X<>Y 33 Rv 34 X!=0? 35 GTO 14 36 X<> ST T 37>LBL 00 38 X<=0? 39 GTO 14 40 SQRT 41 / 42 PUTM 43 X<> ST L 44 STOEL 45 RCLIJ 46 SIGN 47 X<>Y 48 STOIJ 49 R^ 50 J+ 51 FC? 76 52 GTO 10 53 0 54 GTO 00 55>LBL 14 56 RCL ST Z 57 RCLIJ 58 Rv 59>LBL 00 60 RCLEL 61 EXITALL 62 X<>Y 63 RTN 64>LBL "u/" 65 X<>Y 66 ENTER 67 DIM? 68 Rv 69>LBL 11 70 X<> ST Z 71 DIM? 72 X<> ST L 73 EDIT 74 R^ 75 STO ST Y 76 STOIJ 77 Rv 78 STO- ST Y 79 SIGN 80 STO+ ST Y 81 X<>Y 82 GETM 83 X<>Y 84 DIM? 85 Rv 86 CLX 87 RCLEL 88 EXITALL 89 X<>Y 90 EDIT 91 CLX 92 -1 93 EXITALL 94 X<> ST L 95 EDIT 96 Rv 97 X<>Y 98 1 99 STOIJ 100 Rv 101 X<> ST Z 102 LASTX 103 DIM? 104 X<>Y 105 X<> ST L 106 X<> ST Z 107 GETM 108 x 109 RCLIJ 110 XEQ 13 111 - 112 RCLEL 113 EXITALL 114 X<>Y 115 X>0? 116 GTO 11 117 X<> ST Z 118 RTN 119>LBL "u'/" 120 X<>Y 121 1 122>LBL 12 123 X<> ST Z 124 EDIT 125 1 126 R^ 127 STOIJ 128 XEQ 09 129 R^ 130 DIM? 131 Rv 132 X<> ST L 133 EDIT 134 Rv 135 DIM? 136 Rv 137 X<> ST L 138 Rv 139 Rv 140 GETM 141 STOx ST Y 142 DIM? 143 XEQ 13 144 + 145 I+ 146 RCLEL 147 EXITALL 148 X<>Y 149 FC? 76 150 GTO 12 151 X<> ST Z 152 RTN 153>LBL 13 154 CLX 155 RCLEL 156 EXITALL 157 R^ 158 EDIT 159 Rv 160 X<>Y 161 ENTER 162 STOIJ 163 CLX 164 RCLEL 165 +/- 166 STO/ ST T 167 +/- 168 EXITALL 169 X<> ST Z 170 EDIT 171 CLX 172 SIGN 173 STOIJ 174 R^ 175 PUTM 176 Rv 177 RTN 178>LBL 09 179 X<>Y 180 GETM 181 TRANS 182 RCLEL 183 REAL? 184 GTO 00 185 Rv 186 COMPLEX 187 +/- 188 COMPLEX 189 RCLEL 190>LBL 00 191 EXITALL 192 X<>Y 193 EDIT 194 <- 195 -1 196 EXITALL 197 END```

And PD/:

Code:
```00 { 51-Byte Prgm } 01>LBL 14 02 "Not Pos Def" 03 PROMPT 04>LBL "PD/" 05 EDIT 06 RCL- ST X 07 STO+ ST Y 08 X<> ST L 09 EXITALL 10 XEQ "u'u" 11 X>0? 12 GTO 14 13 Rv 14 XEQ "u'/" 15 XEQ "u/" 16 Rv 17 END```

5.Thoughts

In case anyone is wondering:
- the routines are slower than the built-in solving
- they use up more space (if you keep the original A)
- the results are less accurate, and that is actually a bit of a surprise.

So why did I do it? Because I could ;-) And they are a nice example of how to code matrix
routines on the 42S, and a prelude of things to come.

As to the triangular solving, since the 42S factors a general matrix A into a lower triangular matrix L
and an upper unit triangular matrix U, we can just call the builtin solve for U*X=B, and using a 'Mirror' permutation matrix M,
we can do the same for U'X=B. The results are no better though - in fact, ever so slightly worse.

6.Attachment .raw code, annotated excel

Cheers, Werner
How can I align text (the In/out parts) without using Code?
Reference URL's
• HP Forums: https://www.hpmuseum.org/forum/index.php
• :