[VA] SRC #015b - HP-15C & clones: COMPLEX Matrix Inverse up to 8x8

[Valentin Albillo]

  
Hi all,

This is my 1,001^th post here so I'm celebrating it with a new, significant improvement to the inversion routine I included in my original post, but first a few comments are in order, as after a slow start this thread has been gathering momentum by the day. Let's see:

Werner Wrote:You make it increasingly harder for me to improve upon your code ;-) I did manage to squeeze off a puny byte [...] The one thing I miss in your extensive posts, however, is some explanation of the algorithm used. [...] Shamelessly stealing borrowing your idea, we can solve systems of equations the same way:

Thanks for your interest, Werner, I hope you're enjoying your new HP-15C CE and do not regret your purchase at all. It certainly seems to have boosted your inspiration. As for the algorithms, as I said I leave that math-heavy discussion for the future article out of consideration for people who just want to use the code and don't care for its innards.

As for squeezing off a puny byte, well, that's your specialty. See if you can oblige in my new code below. And as for stealing/borrowing my idea, I appreciate your interest in my contributions and I don't resent you doing it at all, quite the contrary, thanks for it.

Werner Again Wrote:Not quite what you asked for, Valentin, but here's a short routine to multiply two split complex matrices.

Very nice, and yes, I was expecting a simple piece of code that would use conventional nested loops to perform the multiplication by taking advantage of the 15C's native complex arithmetic, that's why I left it as an easy exercise for the reader. As it wouldn't involve matrix inversion, it should run fast and accurately.

Werner Once More Wrote:[Update: shaved off 2 bytes ;-) ] Apologies for monopolizing this thread ;-)] [...] to invert A+iB, do (GSB) C and (GSB) E [...] Routine C (Factor) - 21 bytes [...] Routine E (Invert) - 15 bytes

No need to apologize for monopolizing, I'm used to it by now. And congratulations for the 2-byte savings, it's no mean feat and I know you enjoy doing it.

Per your figures, it seems you need to call two routines to invert a complex matrix, totalling 36 bytes in all, which is nice taking into account that routine C is also used to solve a complex system if C is called first so C serves double duty. However, just inversion alone can be done by calling a single 21-byte routine, as I show below.

Gerson Wrote:The matrix capabilities in the HP-15C were nice, but the methods in the manual for solving complex linear systems were somewhat complicated and not easy to remember [...] Doing it programmatically is a great idea, I wonder why no one thought of this back then [...]

Thanks a lot for your great 4x4 EE examples for the original HP-15C, Gerson. The methods in the manual (both the OH and the AFH) are wholly inadequate: extremely cumbersome, very difficult to remember when you needed them, outrageously complicated and requiring several ad-hoc pre- and post-transformations, but the worst of all was their enormous inefficiency, turning an NxN complex matrix into a 2Nx2N real one, doubling the memory required and more than doubling the computing time. In short, a very poor performance by HP in the OH, and the AFH is no better in that regard.

And as for "I wonder why no one thought of this back then", I did, but when the HP-15C was introduced (1982) I had alreay left PPC (and PPC TN) so I had nowhere to submit any further articles, plus I didn't have the money for a 15C or a 71B or any other HP model with complex/matrix capabilities to implement my ideas, nor was I able to until 1985 and afterwards it took me several years to join the MoHPC forum and it's only now that the 15C CE has appeared (with its tremendous 180x speed and 3x the memory,) that I considered worthwhile and interesting to dust off my old notes, implement them and post the results here.

So here you are, the announced improvement ...

The improvement:  Complex Matrix Inversion up to 8x8 in 21 bytes

The complex matrix inversion routine I included in my original post is nice and fast and all that jazz, and even quite short at 31 steps (32 bytes, 5 regs), but it's a bit conservative and thus somewhat longish and I can do better by selecting a different variant, as I'll explain right now.

The thing with complex inversion routines based on this specific split-matrix approach (M = A + iB) is that they require inverting some real matrices (say A or B or A+B or A-B or a random linear combination of A and B, etc.) and thus said matrices must be invertible (i.e. non-singular i.e. det#0).

The good point is that we get to choose which matrix must be invertible, and if the complex matrix includes that component or combination and it happens to be singular, we can choose another that isn't, which would of course require a variant of the inversion routine.

Needless to say, all of this only applies if complex M is invertible even if real A or B (say) are not, but even if both are singular it might be the case that M is not and thus can be inverted using the particular routine which requires A+B (for instance) to be invertible, which is the 31-step routine featured in my original post.

But that routine can be significantly shortened (by some 11 steps no less, 33% shorter) by assuming that A is invertible (which for real-life problems is usually the case and for random matrices is always the case,) i.e. det(A)#0, and here's the resulting improved complex inversion routine (which, as detailed below, can also be used when A is singular but B is not by means of a simple trick):

Program listing

This 20-step (21-byte, 3-reg) routine can be used to invert an NxN complex matrix M = A + iB when A is invertible, and also when A is singular but B is invertible using a simple trick, see details below.

    LBL C           001- 42,21,13
    RESULT A        002- 42,26,11
    RCL MATRIX E    003- 45,16,15
    RCL MATRIX B    004- 45,16,12
    RCL MATRIX A    005- 45,16,11
    1/x             006-       15
    RESULT E        007- 42,26,15
     x              008-       20
    RESULT A        009- 42,26,11
    CHS             010-       16
    RCL MATRIX B    011- 45,16,12
    LASTX           012-    43 36
    1/x             013-       15
    R▼              014-       33
    MATRIX 6        015- 42,16, 6
    1/x             016-       15
    RESULT B        017- 42,26,12
    X<>Y            018-       34
     x              019-       20
    RTN             020-    43 32

Changes in the documentation

I could have edited the docs in my original post to cater for this new routine but that's not my personal policy; unlike most people I only edit a post to correct typos or blatant errors, never to change, remove or add significant contents. Besides, the original documentation fully applies to my original 31-step routine which is the one to use if both A and B are singular but A+B is not, so its docs will remain as they currently are.

What I'll do here is to specify the changes in the docs pertaining to this new 20-step routine, including for context the abridged paragraphs where they appear, with the changed parts in bold red. You should still refer to the original docs for everything else (Notes, Requirements, Worked Example, etc.)

Changes in Notes:

• It runs sequentially from its first to its last step, executing each just once, which means it executes exactly 20 user-code instructions in all (not hundreds or thousands like other approaches,) so it runs very fast.
   
  Comment: As 11 fewer instructions are executed in the new routine (33% less,) it runs somewhat faster.

Changes in Requirements:

• The maximum size NxN complex matrix M you can invert depends on the memory available in your physical or virtual device, as per this table:
  
         M    A,B,E   Regs  +Prog    Comments
      -------------------------------------------------------------
        1x1    1x1      3      6        -
        2x2    2x2     12     15        -
        3x3    3x3     27     30        -
        4x4    4x4     48     51    Max. size w/ 15C/64 but see (*)
        5x5    5x5     75     78    ditto CE/96
        6x6    6x6    108    111        -
        7x7    7x7    147    150    ditto CE/192 but see (**)
        8x8    8x8    192    195    ditto DM15/M1B
  
  so e.g. if you want to invert an 8x8 complex matrix you need 195 regs available, which includes matrices A, B (and E,) plus 3 regs to hold the routine itself.
   
  Comment: The new routine is only 21 bytes long, so it occupies exactly 3 regs of program memory. As the initial routine occupied 5 regs, this means that now 2 extra regs are available for additional data or up to 14 additional program steps.
   
• The matrix A must be invertible, i.e. det(A)#0. For most real-life uses this will be the case but if this condition isn't met there are slight variations to this routine that would work Ok in those cases, namely for any of the following invertible matrices: B, A+B (dealt with in my original post in this thread,) A-B, etc.
  

Changes in (*) Observation re the original HP-15C and 4x4 complex matrices:

• The original HP-15C could invert a 4x4 complex matrix but it took all 64 regs available and it was a complicated, completely manual process as there wasn't any memory left for program code. On the other hand, running my routine is a fast, automated process that leaves out all the drudgery (transformations, etc.) and still leaves 13 regs (up to 91 extra program steps) free for additional code or data. [...]
   
• Call my complex inversion routine. When it ends, the complex inverse is stored in A, B and there's still 13 regs free. [...]
   
• Get rid of auxiliary matrix E (redimension it to 0x0.) This frees another 16 regs, so there's now 29 regs available for what follows.
   
  Comment: There's now 29 regs still available in the original HP-15C after inverting a 4x4 complex matrix, which can be used to store further data (e.g. the constant matrix and the solution matrix, when programmatically solving a 4x4 system of complex equations (ignore the abstruse manual methods featured in the OH and AFH) or up to 203 program steps for additional code (e.g. your own program which uses the inversion routine,) or a combination thereof.
   
• Dimension both the constant matrix (say C) and the solution matrix (say D) to be 4x2 and populate the constant matrix. This will leave 13 regs (i.e. as many as 91 program steps) still available for the matrix-multiplication code, which is left as a fairly easy exercise for the reader (just a loop which multiplies each row of the inverse by the constant matrix using the HP-15C's native complex arithmetic.)
  

Changes in (*) Observation re the HP-15C CE and 8x8 complex matrices:

• Though there's not enough room in the HP-15C CE in 192-regs Mode to invert an 8x8 complex matrix M by running the routine featured here (195 regs would be needed; the routine itself wouldn't fit) there's just enough room for the split matrices A, B and the auxiliary matrix E (192 regs in all,) so the 8x8 complex matrix can be inverted in a pinch if the user executes the 18 program instructions manually in sequential order. The procedure would be like this: [...]
   
  • Carefully execute manually the 18 instructions from 002 RESULT A  to  019  x  in sequential order. Assuming you're reasonably proficient with the HP-15C, this should take 2 min or less [...]
   
  Comment: As now there's only 18 instructions to manually execute instead of the 29 required by the original routine, this makes the process all the more expeditious, less tiring to perform and much more likely to be completed without error.
  
  Also, it's a pity that the index registers R0, R1 and RI are permanent and can't be allocated to the common pool. If it were possible, this inversion routine (which doesn't use the index registers,) could be stored in the 3 regs (21 bytes) they would provide, thus allowing for inverting 8x8 complex matrices programmatically. Alas, no such luck, HP didn't see fit to allow for the possibility.

Using this routine when A is singular but B is not

As I explained above, usually different variants of my complex inversion routine are needed depending on whether (1) A is invertible, or (2) A is singular but B is invertible, or (3) both A and B are singular but A+B (or A-B) is invertible.

My original 32-byte routine somewhat conservatively addressed the third case, where A and B can both be singular as long as A+B isn't (a slightly-different, same-length variant would handle the invertible A-B case.)

Now, my new 21-byte routine caters for the first case (A is invertible.) The second case (A is singular but B is invertible) can be dealt with using a near-symmetric same-length variant of the routine but actually there's no need to have two almost-identical routines in memory at once, as a simple trick will allow the present one to also handle the second case (A singular, B invertible), just proceed like this:

To compute using this routine the complex inverse of M = A + iB when A is singular but B is invertible:

1. Initialize and dimension matrices A, B as described in my OP's Worked example.
    
2. Swap the roles of A and B, i.e. store M's real parts in B and the imaginary parts in A.
    
3. Compute the complex inverse of this "swapped" matrix:  GSB C
    
4. Now swap back the contents of A and B while also changing the signs of their elements on the fly. From the keyboard, execute:
   
       RCL MATRIX A        { swap back A and B using auxiliary matrix E ... }
       CHS                 { ... while also changing their signs at the same time }
       STO MATRIX E        { done; matrix E now holds the negated contents of A }
   
       RCL MATRIX B        { now matrix B is recalled to the X stack register ... }
       CHS                 { ... its contents are negated ... }
       STO MATRIX A        { ... and stored back in A }
   
       RCL MATRIX E        { finally, the negated contents of A held in matrix E ... }
       STO MATRIX B        { ... are stored back in B and the swap is completed. }
   
   Now A contains the real parts of M's true complex inverse, M', and B contains the imaginary parts, as desired.
     
   Note: Of course, all those 8 instructions executed manually can be included as program lines to be executed programmatically after calling my routine, in case A in singular but B is not. Using the ad-hoc variant of my routine for such case wouldn't need those 8 steps but then you'd have to include said routine's 20 steps as well. Oh, and there's a trick to swap the contents of two NxN matrices without needing a third auxiliary one.
   

That's all. I've missed a few usual suspects posting here but that's life, no hard feelings .  As always, hope you enjoyed it.

V.
Edit: a few typos.