Threaded Mode | Linear Mode

Valentin Albillo · 08-26-2023, 10:39 PM

Hi, all,

It seems this thread has run its course, without attracting much attention to be honest but that's life. Thanks for your interest and comments, Fernando del Rey, ctrclckws, EdS2, Werner (who managed to shorten my inversion code by two steps and further posted a related routine to solve large systems,) J-F Garnier, Namir and last-but-not-least, newcomer Bert, who registered just to let me know of a one-character but significant typo, much appreciated.

Now it's time to provide some closure and so I'm posting this Epilogue of sorts which includes a short 42-step timing-helper program which can be used to obtain arbitrarily accurate timings for both my matrix inverse and determinant routines for large N, namely 9 ≤ N ≤ 16, depending on the particular device's available memory.

Here you'll find the commented Program Listing, the Usage Instructions, the Timings I've obtained using it on my HP-15C CE/192, as well as ancillary sections for Things To Do and Miscellanea.

The program runs on the CE/192 for 9 ≤ N ≤ 12 (the CE/96 lacks the necessary memory) as well as any other devices physical or virtual which have enough memory. The DM15 with firmware M1B might be able to also do N=13, but no current device has enough memory to do 14 ≤ N ≤ 16 though both this timing program and the routines themselves would work unchanged.

Program listing:

    001  LBL A      018  RCL MATRIX A   035  LBL 0
    002  CF 8       019  GSB 0          036  STO I
    003   1         020  RCL MATRIX B   037  MATRIX 1
    004  DIM (i)    021  GSB 0          038  LBL 1
    005  X<>Y       022  RCL MATRIX C   039   RAN#
    006  MATRIX 0   023  GSB 0          040u  STO (i)  (*)
    007   8         024  RCL MATRIX D   041  GTO 1
    008   -         025  GSB 0          042  RTN
    009   8         026   2             043  LBL E        043  LBL D
    010  DIM C      027   0             ...  ...     or   ...  ...
    011  X<>Y       028  R/S            071  RTN          057  RTN
    012  DIM B      029  STO I
    013  ENTER      030  LBL 2
    014  DIM D      031   GSB E  or  GSB D
    015   8         032   DSE I
    016  ENTER      033  GTO 2
    017  DIM A      034  RTN

             (*) Step 040u STO (i) must be entered in USER mode.

Notes

Steps 001-017 initialize the program (disable the complex stack, allocate all pool registers for matrices, reset all matrices to 0x0) and dimension all 4 blocks A, B, C, D to their respective sizes: 8x8, 8x(N-8), (N-8)x8 and (N-8)x(N-8).
Steps 018-028 fill up all 4 blocks with random values, then stop with the default number of loops (20) in the display, ready for the user to start the loops and the timing.
Steps 029-034 perform the specified number of loops, calling the partitioned matrix inversion routine (LBL E) or the determinant routine (LBL D) that many times, then execution ends.
Steps 035-042 implement a subroutine that fills up with random values the matrix whose descriptor is passed to it in the X stack register.
Steps 043-071 or 043-057 implement my partitioned matrix inversion routine (LBL E) or my determinant routine (LBL D), respectively, whose listings can be found in the first post of this thread.

Usage Instructions:

Note: This step is not necessary, but if you want to check that the program has been correctly entered and works as intended, do this:

1, STO RAN#, 9, GSB A -> 20

and if now you recall the contents of the first element of matrices A and D you should get (in FIX 4):

A1,1 = 0.2018, D1,1 = 0.1746

This program works strictly for NxN matrices with 9 ≤ N ≤ 16 (subject to available memory; N ≤ 12 for the CE/192). Inputting N outside this range will generate an error. To time the inversion/determinant of an NxN matrix, do as follows:

First, run the program to dimension and quickly fill up all NxN elements with random values, then it stops with the default number of loops (20) in the display:

FIX 4, N, GSB A -> 20

The user can either accept this default (20 loops i.e. 20 NxN matrix inversions or determinants will be computed) or else key in another number of loops, say 50 or whatever. The more loops, the more time it takes to perform them all but the greater the timing accuracy. Assuming the user can get the final time accurate to the nearest second, the timing accuracy will be 1/#loops i.e. 1/20 = 0.05" for the default.
Now the user must simultaneously press R/S and start the timing, closely watching for the program to end and taking note of the final time.

R/S -> E 1 8 (say, depends on N) for the inversion or some numeric value for the determinant.

and dividing the final time by the number of loops will give the time per NxN operation accurate to 1/#loops.

Timings:

These are the times I obtain on my CE/192 when using the above program with new, fresh lithium batteries (yours may vary slightly.) Note that I use the default 20 loops to time the inversion routine but change it to 50 loops when timing the determinant routine, as it is about 3x faster:

                                       Inverse        Determinant
         Dim   Block A  Block D   20 loops  time    50 loops  time
       -------------------------------------------------------------
         9x9     8x8      1x1       18"     0.90"      17"    0.34"
        10x10    8x8      2x2       24"     1.20"      22"    0.44"
        11x11    8x8      3x3       31"     1.55"      29"    0.58"
        12x12    8x8      4x4       40"     2.00"      37"    0.74"

Thus, when using my matrix inversion routine, inverting a 9x9 matrix takes my HP-15C CE less than one second, and inverting a sizable 12x12 matrix takes about two seconds flat.

As for the determinant routine, it can compute the determinant of a 9x9 matrix in 1/3 of a second and that of a 12x12 matrix in less than 3/4 of a second. Thus, it can compute three 9x9 determinants per second. That's fast !

Things To Do:

For those interested, there's a number of further things to do if you feel like it, e.g.:

In my OP I provided routines for matrix inversion and determinant computations but, as Fernando del Rey pointed out, I (purposefully) didn't provide one for system solving, which would complete the "matrix trilogy" and might be quite useful for large systems.

I may eventually post mine here but in the meantime Werner has recently posted a 212-step, 238-byte system solving program which doesn't use or rely on my code, and just a few days ago posted another version but this time in the spirit of my routines, i.e. dealing with suitably partitioned matrices, which is thus 5x shorter and 10x faster than his previous standard attempt. This further proves that when solving large systems in the memory-expanded HP-15C partitioning is the way, which was the whole point of this SRC#15 !

Be as it may, there are other ways to try and solve systems and as I said above, I may eventually post one or two different attempts, not doing it right now because documenting and formatting those routines per my standards does take a lot of time. Meanwhile, other people's routines are most welcome.
My routine for matrix inversion takes a 4-block partitioned NxN matrix as input and does the inversion in place, so the matrix inverse is also in the same form and can be manually output by the user or utilized by some other program as-is. However, it might be useful to write two supporting, complementary routines:
1. to convert a non-partioned NxN matrix to a 4-block partitioned matrix, as needed by the matrix inverse and determinant routines.
2. to convert the partitioned NxN matrix inverse to a single non-partitioned NxN matrix. This might simplify further operations with the matrix inverse and would also free 3 matrix descriptors for other purposes.
Call my inversion routine from one program of yours to perform, say, Nth-degree Polynomial Fit or Nth-degree Least-Squares Polynomial Regression or Multivariate Linear Regression, or any other tasks which would benefit from using an NxN matrix inverse or determinant for 9 ≤ N ≤ 12 (subject to available memory.)

Miscellanea:

Once R/S is pressed to start the timing of the matrix inversion routine, my CE immediately displays one or two running messages in very quick succession, then the display remains blank for the whole duration of the loops. This seems very similar to what happens with the dreaded LE's PSE (Pause) bug (except there's no PSE instruction whatsoever in my program/routines,) so perhaps the cause is somehow related.

On the other hand, if timing the determinant routine, upon pressing R/S a single running message is immediately displayed once and it stands still, never blinking or blanking till the timing process ends. All in all, I find the display of the running message to be quite haphazard, sort of yet "unfinished business".
While checking this timing-helper program, it accidentally tried to invert a 9x9 matrix using the buil-in 1/x function (which can't reliably work for matrices larger than 8x8,) and it immediately halted with an Error message and a blinking display. All my bad, the code I wrote to compute the blocks' dimensions was faulty. Fortunately nothing got corrupted AFAICT.

Regards.
V.