Matrix Operations for the HP-41

This program is supplied without representation or warranty of any kind.
Jean-Marc Baillard and The Museum of HP Calculators therefore assume no
responsibility and shall have no liability, consequential or otherwise, of
any kind arising from the use of this program material or any part thereof.

-All these programs deal with real matrices only.
-Some of them use synthetic registers.
-Don't interrupt a routine which uses registers P and/or Q (
no risk of "crash" but their contents could be altered )

-The elements aij are to be stored in column order into successive
data registers
and each matrix is identified by a control number of the form
bbb.eeerr

where Rbb is the first register
Ree is the last register
and rr
is the number of rows of the matrix.

-Key in:
1.01203 ENTER^
15.02204 ENTER^
26 XEQ "M*" >>>>
26.03103 = the control number of the matrix
C and the result is:

47 39 R26 R29
C = 71 51
= R27 R30
52 32 R28 R31

-The number of columns of the first matrix must equal the number of
rows of the second matrix.

b) Program#2 ( product
tA.B
)

-It is sometimes useful to multiply the transpose of a matrix A by another
matrix B directly.
-This routine is used to solve underdetermined and overdetermined systems
( cf "Linear and non-Linear Systems for the HP-41" )

-Key in:
1.01204 ENTER^
15.02204 ENTER^
26 XEQ "TM*" >>>>
26.03103 = the control number of the matrix
C and the result is:

47 39 R26 R29
C = 71 51
= R27 R30
52 32 R28 R31

3°) Matrix Inversion

"INV" can invert up to a 16x16 matrix, or even a 17x17 matrix
if it's executed directly from extended memory.
Gauss-Jordan elimination ( also called the "exchange method"
) is used.

Here, the first element of the matrix must be stored into R06.
( R01 thru R05 are used for control numbers of different rows and columns.)
You put the order of the matrix into X-register and XEQ "INV"
the determinant is in X-register and in R00 and the inverse matrix
has replaced the original one ( in registers R06 ... etc ... )

If flag F01 is clear, Gaussian elimination with partial pivoting
is used.
If flag F01 is set, the pivots are the successive elements of
the diagonal.

The XTOA AROT ATOX and REGSWAP functions of the X-Functions
module are used.
If you don't have an X-Functions module, you can delete lines
148 to 177 and lines 25 to 69 but your matrix will have to be "completely
regular"
I mean that every diagonal pivot must be non-zero. However, if
the first element is 0 , you can swap row n°1 and row n°2
for instance ,
and after the whole calculation , you will have to swap column
n°1 and column n°2

In fact, the HP-41 has to remember the different exchanges of
rows in order to do the same exchanges of columns in the reverse order
after the calculation is achieved, and XTOA is used
for that purpose.

-If you try to invert a Pascal's matrix of order 16 for instance with
flag F01 cleared, you'll obtain great roundoff errors ( even with an HP48
)
because these matrices are very ill-conditioned. But if you set
F01, all the coefficients of the inverse will be computed exactly!
( Execution time = 42mn )

-This program also works if p = 0.
-"MPOW2" is actually slower than "MPOW" for n = 3 ( the order of the
matrix ) and p = 7
-But if n = 3 and p = 100 , execution time is only 2mn31s ( instead
of 21mn29s with "MPOW"! )

-The elements of exp(A) will be accurately computed if all the elements
of A are non-negative.
-Otherwise - especially if some elements are large negative numbers
- the results may even be meaningless because of cancellation of leading
digits!

-If it's possible, try to diagonalize A ( i-e find a regular matrix
B so that B-1A.B is diagonal ) and apply the formula exp
( B-1AB ) = B-1exp(A).B
-The exponential of a diagonal matrix A = [ aij ] is a diagonal
matrix too where the diagonal elements are exp(aii)

11°) Logarithm of a Square Matrix

-If || A - I || < 1 , the logarithm of a nxn matrix A is defined
by the series expansion:

-You enter the different elements in column order.-When the first column is stored, simply press R/S without any
digit entry,
the next column indexes will be automatically incremented.
-When all the elements are stored, press R/S and you'll get the control
number of the matrix!