(42S) Cholesky decomposition and solve

[Werner]

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?