Matrix Division on the 28C/S (and later)

[Werner]

Solving A*X=B for square A is done factorizing A into P*L*U, with P a permutation matrix and L and U lower and upper triangular matrices, respectively, of which one has only 1's on the diagonal (it depends on the algorithm used which one has the 1's). L and U are kept in A, overwriting it.
Then, A*X=B is solved by solving, in order:
(1) P*Z=B, effectively applying the same permutations to B as had been done to A
(2) L*Y=Z
(3) U*X=Y

Solving A*X=B as X=inv(A)*B requires a number of additional steps:
- invert A from A=P*L*U
1.inverting L: R=inv(L)
2.solving U*C=R to obtain C=inv(U)*inv(L)
3.applying P from the right: C=C*P
- multiply C and B
these steps require a lot more calculations, and are somewhat less accurate as a result.

BTW SYS on the 71B solves the system as I set out above, but it performs an additional 'iterative refinement' step to obtain a more accurate result:
Solving A*x=b does not give you x but x', an approximation of the true x.
Let x = x' + dx
Let's try and determine dx:
A*(x' + dx) = b
A*dx = b - A*x'
If we can calculate c = b-A*x' with greater accuracy, we can determine dx, solving A*dx = c.
In the 71B, b-A*x' is determined with 15-digit arithmetic (it's just a number of DOT products).
Alas, the 42S does not do this, but I have written a routine that will perform an iterative refinement step. It is only useful in the real 42S and emu42, not in Free42 or the DM42, as these do not carry out dot products in higher precision (sadly). I already posted it somewhere here before, I think.

Code:

                L       X       Y       Z       T
00 { 58-Byte Prgm }
01 LBL "/IR"            A       b
02 X<>Y
03 ENTER
04 TRANS                bt      b       A
05 RCL ST Z             A       bt      b       A
06 STO/ ST Z            A       bt      x       A
07 TRANS
08 EDIT
09 INSR                 0       bt      x       A
10 Rv                   bt      x       A       .
11 PUTM
12 CLX
13 RCLEL
14 EXITALL
15 TRANS                bA      x       A
16 X<>Y
17 LBL 02               x       bA      A
18 STO ST L    x        x       bA      A
19 EDIT
20 INSR
21 -1
22 EXITALL     x        -1x     bA      A
23 X<> ST L  -1x        x       bA      A
24 RCL ST Y  -1x        bA      x       bA      A
25 RCLx ST L            Ax-b    x       bA      A
26 RCL/ ST T           -dx      x       bA      A
27 -                    x        bA        A
28 RTN
29 GTO 02
30 END

Provide b ENTER A, then run "/IR".
It will solve the system and perform a single refinement step, leaving x on the
stack. If you're careful, you can examine x with EDIT and the arrow functions,
leaving the stack undisturbed. Exit x, and press R/S to carry out another
refinement step.
It will only work with vectors b (matrices with a single column).
Usually, a single refinement step suffices.

Cheers, Werner