(42S) Cholesky decomposition and solve - Printable Version +- HP Forums (https://www.hpmuseum.org/forum) +-- Forum: HP Software Libraries (/forum-10.html) +--- Forum: General Software Library (/forum-13.html) +--- Thread: (42S) Cholesky decomposition and solve (/thread-6704.html) (42S) Cholesky decomposition and solve - Werner - 08-19-2016 09:02 AM 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?