Threaded Mode | Linear Mode

Werner · (This post was last modified: 08-21-2023 01:34 PM by Werner.)

Solve a system of linear equations

The matrix routines to solve, invert or calculate the determinant do not work for matrices of order greater then 8; (eg the determinant of a 9x9 identity matrix is calculated as 0..)
Valentin has already shared a short and elegant routine to go beyond order 8 for determinant and inverse. My routine is neither short nor elegant, but it gets the job done ;-)

Use: exactly like the built-in 15C functionality, ie.
to solve A*C=B (so C=inv(A)*B) enter

RCL MATRIX B
RCL MATRIX A
RESULT C
GSB B

to subsequently solve for a different right-hand side, eg D; this time overwriting D with the result:

RCL MATRIX A
RCL MATRIX D
RESULT D
GSB C

With the bigMem version of the 15C CE , we can solve a 12x12 system. It takes about 10.5s, only slightly longer than a 42S running built-in code. Not bad for usercode ;-)
(it won't work for higher dimensions, even if we had room, as the pivotation scheme may fail).
Incidentally, the LU-factorisation produces the same results as the built-in routine, apart from reduced accuracy (the built-in routines follow the same algorithm, but the calculations are done in a different order to make maximum use of extended precision dot products).

Register usage:

I : matrix descriptor
0 : row index
1 : column index
2 : counter
3 : pivot array = SUM(p*(n-i)!), p is relative distance to pivot

So you need at least 3 DIM (i).

routine overview:
LBL .0 : LU-factorisation
LBL 1 : main loop
LBL 2 : search for largest element in column
LBL 3 : scale column
LBL 4 : rank-1 matrix modification

LBL .1 : solve LU.X=P.B
LBL 5 : Z = P.B
LBL 6 : solve L.Y=Z
LBL .6 : solve U.X=Y

subroutines
LBL 7 : rank-1 modification (2 matrices)
LBL 8 : swap two rows
LBL 9 : y := y - a*x subroutine

LBL B : user routine to factor and solve
LBL C : user routine to solve with an already factored matrix

212 lines, 238 bytes
001 - 42,21,.0  LBL .0
002 - 45,23,25  RCL DIM I
003 -       26  E
004 -        3  3
005 -       10  /
006 -    44  2  STO 2
007 -    43 35  CLX
008 -    44  3  STO 3
009 - 42,21, 1  LBL 1
010 - 45,23,25  RCL DIM I
011 -    45  2  RCL 2
012 -    43 44  INT
013 -       30  -
014 - 44,20, 3  STOx 3
015 - 42, 6, 2  ISG 2
016 -    45  2  RCL 2
017 -    44  0  STO 0
018 -    44  1  STO 1
019 -        0  0
020 -       36  ENTER
021 - 42,21, 2  LBL 2
022 -    45  0  RCL 0
023 -    43 33  R^
024 -    45 24  RCL (i)
025 -    43 16  ABS
026 -    43 33  R^
027 - 43,30, 8  TEST 8
028 -       33  Rv
029 - 42, 6, 0  ISG 0
030 -    22  2  GTO 2
031 -       33  Rv
032 -       33  Rv
033 -    44  0  STO 0
034 - 45,30, 2  RCL- 2
035 -    43 20  X=0?
036 -    22  0  GTO 0
037 - 44,40, 3  STO+ 3
038 - 45,23,25  RCL DIM I
039 -    44  1  STO 1
040 -    32  8  GSB 8
041 - 42,21, 0  LBL 0
042 -    45  2  RCL 2
043 -    44  0  STO 0
044 - 42, 6, 0  ISG 0
045 -    43 20  X=0?
046 -    43 32  RTN
047 -       36  ENTER
048 -    44  1  STO 1
049 - 45,43,24  RCL g (i)
050 -    43 20  X=0?
051 -    22  1  GTO 1
052 - 42,21, 3  LBL 3
053 - 44,10,24  STO/ (i)
054 - 42, 6, 0  ISG 0
055 -    22  3  GTO 3
056 -    45  2  RCL 2
057 -    44  0  STO 0
058 - 42, 6, 0  ISG 0
059 - 42,21, 4  LBL 4
060 -    45  2  RCL 2
061 -    44  1  STO 1
062 - 42, 6, 1  ISG 1
063 -    45  0  RCL 0
064 -    45  2  RCL 2
065 - 45,43,24  RCL g (i)
066 -    43 16  ABS
067 -    32  9  GSB 9
068 - 42, 6, 0  ISG 0
069 -    22  4  GTO 4
070 -    22  1  GTO 1
071 - 42,21,.1  LBL .1
072 -    45 26  RCL RESULT
073 - 42, 4,25  X<> I
074 -    44 26  STO RESULT
075 - 45,23,25  RCL DIM I
076 -    43 35  CLX
077 -       48  .
078 -        1  1
079 -    43 14  %
080 -    44  2  STO 2
081 -    45  3  RCL 3
082 - 42,21, 5  LBL 5
083 -       36  ENTER
084 - 42, 6, 2  ISG 2
085 - 45,23,25  RCL DIM I
086 -    44  1  STO 1
087 -    43 35  CLX
088 -    45  2  RCL 2
089 -    43 44  INT
090 -    44  0  STO 0
091 -       30  -
092 -    42  0  x!
093 -       10  /
094 -    43 36  LASTX
095 -       34  X<>Y
096 -    43 44  INT
097 - 44,40, 0  STO+ 0
098 -       20  x
099 -       30  -
100 -    32  8  GSB 8
101 - 43,30, 1  TEST 1
102 -    22  5  GTO 5
103 - 45,23,25  RCL DIM I
104 -    43,35  CLX
105 -        1  1
106 -       30  -
107 -       26  E
108 -        3  3
109 -       10  /
110 -    44  2  STO 2
111 -    22  0  GTO 0
112 - 42,21, 6  LBL 6
113 -    45  2  RCL 2
114 -        1  1
115 -       48  .
116 -        0  0
117 -        0  0
118 -        1  1
119 -       40  +
120 -    44  0  STO 0
121 -    32  7  GSB 7
122 - 42,21, 0  LBL 0
123 - 42, 6, 2  ISG 2
124 -    22  6  GTO 6
125 - 45,23,25  RCL DIM I
126 -       34  X<>Y
127 -    44  2  STO 2
128 - 42,21,.6  LBL .6
129 - 45,23,25  RCL DIM I
130 -    44  1  STO 1
131 -    45 26  RCL RESULT
132 - 42, 4,25  X<> I
133 -    45  2  RCL 2
134 -       36  ENTER
135 -    44  0  STO 0
136 - 45,43,24  RCL g (i)
137 -       34  X<>Y
138 -    44 25  STO I
139 -       34  X<>Y
140 - 42,21,.7  LBL .7
141 - 44,10,24  STO/ (i)
142 - 42, 5, 1  DSE 1
143 -    22 .7  GTO .7
144 -    45  2  RCL 2
145 -        1  1
146 - 43,30, 5  TEST 5
147 -    43 32  RTN
148 -       30  -
149 -       26  E
150 -        3  3
151 -       10  /
152 -    44  0  STO 0
153 - 42, 6, 0  ISG 0
154 -    32  7  GSB 7
155 - 42, 5, 2  DSE 2
156 -    22 .6  GTO .6
157 - 42,21, 7  LBL 7
158 - 45,23,25  RCL DIM I
159 -       26  E
160 -        3  3
161 -       10  /
162 -    44  1  STO 1
163 - 42, 6, 1  ISG 1
164 -    45 26  RCL RESULT
165 - 42, 4,25  X<> I
166 -    45  0  RCL 0
167 -    45  2  RCL 2
168 - 45,43,24  RCL g (i)
169 -    43 16  ABS
170 -       34  X<>Y
171 -    44 25  STO I
172 -    32  9  GSB 9
173 - 42, 6, 0  ISG 0
174 -    22  7  GTO 7
175 -    43 32  RTN
176 - 42,21, 8  LBL 8
177 -    45  2  RCL 2
178 -    45  1  RCL 1
179 - 45,43,24  RCL g (i)
180 - 42, 4,24  X<> (i)
181 -    45  2  RCL 2
182 -    45  1  RCL 1
183 - 44,43,24  STO g (i)
184 -       34  X<>Y
185 - 42, 5, 1  DSE 1
186 -    22  8  GTO 8
187 -    43 32  RTN
188 - 42,21, 9  LBL 9
189 -    45  2  RCL 2
190 -    45  1  RCL 1
191 - 45,43,24  RCL g (i)
192 -    43 36  LASTX
193 -       20  x
194 - 44,30,24  STO- (i)
195 - 42, 6, 1  ISG 1
196 -    22  9  GTO 9
197 -    43 32  RTN
198 - 42,21,12  LBL B
199 -    44 25  STO I
200 -    43 35  CLX
201 -       40  +
202 -    32 .0  GSB .0
203 -    22  0  GTO 0
204 - 42,21,13  LBL C
205 -    44 25  STO I
206 -    43 35  CLX
207 -       40  +
208 - 42,21, 0  LBL 0
209 -    32 .1  GSB .1
210 -    45 25  RCL I
211 -    44 26  STO RESULT
212 -    43 32  RTN

Comments welcome!
Cheers, Werner

Werner · 10-12-2023, 11:18 AM

New, shorter version, with comments:
(I also prefix 2-byte instructions with '_' to easily determine the byte count)

I : descriptor
0 : r row index
1 : c col index
2 : i counter
3 : P pivots
RESULT : matrix

207 bytes

  001 LBL B        Factor and Solve
  002 STO I             X Y RESULT
  003 CLX          In:  X Y R
  004 +            Out: R = inv(X).Y
  005 RCL DIM I    PA=LU factorisation
  006 E3
  008 /
  009 STO 2
  010 CLX
  011 STO 3
  012 LBL 1        Outer loop i=1..n
  013 RCL DIM I
  014 RCL 2
  015 INT
  016 -
  017_STOx 3       P := P*(n-i);
  018_ISG 2
  019 RCL 2
  020 STO 0
  021 STO 1
  022 0            max |Ari|, r=i..n
  023 ENTER        X   Y   Z   T
+>024 LBL 2        G   -   p   -
| 025 RCL 0
| 026 R^
| 027 RCL (i)
| 028 ABS
| 029 R^           G |Ari| p   r
| 030 TEST 8       X<Y?
| 031 Rv
| 032 ISG 0
+-033 GTO 2
  034 RCL 2
  035 R^
  036 STO 0
  037 -
  038_STO- 3       Save pivot
  039 RCL DIM I    n  n
  040 STO 1
  041 -            0
  042 TEST 6       X#Y?
  043 GSB 8        Aic :=: Apc; c=n..1
  044 RCL 2
  045 STO 0
  046 STO 1
  047 RCL (i)
  048 ENTER
  049 ISG 0
  050 TEST 6       X#Y? - AFF
  051 GTO 0        end
  052 X=0?
  053 GTO 1
+>054 LBL 3        Ari := Ari/Aii; r=i+1..n
| 055_STO/ (i)
| 056 ISG 0
+-057 GTO 3
  058 RCL 2
  059 STO 0
  060 ISG 0
+>061 LBL 4        Rank-1 update
| 062 RCL 2        for r=i+1 to n
| 063 STO 1          t := Ari;
| 064uRCL (i)        for c=i+1 to n
| 065 ABS              Arc := Arc - t*Aic;
| 066 GSB 9          next c;
| 067 ISG 0        next r;
+-068 GTO 4
  069 GTO 1

  070 LBL C        Solve only
  071 STO I
  072 CLX
  073 +            Set RESULT
  074 LBL 0        Solve LUX = PB
  075 RCL RESULT
  076 X<> I        Work on B
  077 STO RESULT
  078 RCL DIM I
  079 CLX
  080 E3
  082 /
  083 STO 2        i=0.00n
  084 RCL 3
+>085 LBL 5        Z=PB
| 086 ENTER        for i=1 to n
| 087_ISG 2          pi := P/(n-i)!;
| 088 RCL DIM I      pi := INT(pi);
| 089 STO 1          P := P - pi*(n-i)!;
| 090 CLX            r := i + pi;
| 091 RCL 2          Bic :=: Brc; c=m..1
| 092 INT          next i;
| 093 STO 0
| 094 -
| 095 x!
| 096 /
| 097 LASTX
| 098 X<>Y
| 099 INT          pi
| 100_STO+ 0
| 101 x
| 102 -            Pnew Pold
| 103 TEST 6       X#Y?
| 104 GSB 8
| 105 TEST 1       X>0?
+-106 GTO 5
  107 RCL 2
  108 FRAC
  109 1
  110 +
  111 STO 2
+>112 LBL 6        LY=Z
| 113 RCL 2        for i=1 to n-1
| 114 STO 0          for r=i+1 to n
| 115 ISG 0            t := Ari;
| 116 GSB 7            Brc := Brc - t*Bic; c=1..m
| 117_ISG 2          next r;
+-118 GTO 6        next i;
  119 RCL DIM I
  120 X<>Y
  121 STO 2
+>122_LBL .6       UX=Y
| 123 RCL 2        for i=n to 1
| 124 STO 0          Bic := Bic/Aii; c=1..m
| 125 GSB 0          for r=1 to i-1
| 126_LBL .7  <-+      t := Ari;
| 127_STO/ (i)  |      Brc := Brc - t*Bic; c=1..m
| 128 ISG 1     |    next r;
| 129_GTO .7  --+  next i;
| 130 RCL 2
| 131 1
| 132 -
| 133 E3
| 135 /
| 136 STO 0        r=1..i-1
| 137 ISG 0
| 138 GSB 7
| 139_DSE 2
+-140_GTO .6
  141 RCL I
  142 STO RESULT
  143 RTN
                   Subroutines
  144 LBL 7        Rank-1 modification
  145 GSB 0
  146 ABS
  147 GSB 9
  148 ISG 0
  149 GTO 7
  150 RTN

  151 LBL 0        recall Rri
  152 RCL DIM I
  153 E3
  155 /
  156 STO 1
  157 ISG 1        c=1.00m
  158 RCL RESULT
  159 X<> I
  160 RCL 0
  161 RCL 2
  162_RCL g (i)
  163 X<>Y
  164 STO I
  165 X<>Y
  166 RTN

  167 LBL 8        swap subroutine
  168 RCL 2        Aic :=: Arc; DSE c
  169 RCL 1             X Y Z T 0 1 2
  170_RCL g (i)    In:  X       r c i
  171 X<> (i)      Out: X       r   i
  172 RCL 2
  173 RCL 1
  174_STO g (i)
  175 X<>Y
  176 DSE 1
  177 GTO 8
  178 RTN

  179 LBL 9        axpy subroutine
  180 RCL 2        Arc := Arc - t*Aic; ISG c
  181 RCL 1            L X Y Z T 0 1 2
  182_RCL g (i)    In: t         r c i
  183 LASTX
  184 x
  185_STO- (i)
  186 ISG 1
  187 GTO 9
  188 RTN

Cheers, Werner