A Matrix Multi-Tool for the HP 35s Programmable Calculator

October 17, 2007

The HP 35s programmable calculator.

The HP 35s programmable calculator.

This is my second attempt at a large program for the new HP 35s programmable scientific calculator from Hewlett-Packard. My previous program addressed the curve fitting shortcomings of the 35s (compared to the HP-41C/CV/CX with Advantage Pac, or the HP-42S). This program is a start at doing the same for matrix functionality.

The HP-15C was the first calculator to introduce comprehensive support for matrix operations. HP quickly added equivalent functionality to the 41C series with the Advantage Pac, and the 42S came with these operations built in as well. The HP 35s does not have any general purpose matrix operations built in, although it can solve 2×2 and 3×3 linear systems.

The program I present here is a long way from being as powerful as the facilities provided by the 15C, 41C/Advantage, or 42S, but it does provide several useful matrix operations in one simple tool. It is loosely modelled after the Advantage Pac matrix program, which provides an easy-to-use subset of the functionality of the matrix library.

What it Does

Given an N×N matrix A, and an N-element vector b, the matrix multi-tool will do the following:

  • Compute A-1, the inverse of A.
  • Compute the determinant of A.
  • Solve the system of linear equations, Ax=b, giving column vector x.
  • Quickly solve additional Ax=b systems for different b vectors.
  • Perform matrix-vector multiplication of A and b.

The matrix multi-tool uses Gaussian elimination with partial pivoting to intially compute A-1 and simultaneously solve Ax=b. The determinant is also computed during this operation. The inverse is left in A, and x is left in b. When solving additional systems of equations for the same A and different b vectors, the program evaluates A-1b to yield x.

If you’re only interested in A-1 and/or A‘s determinant, you can omit entering values for b and ignore the solution of x.

You can also use the built-in matrix-vector multiplication routine directly if you just want to multiply a matrix by one or more vectors.

Thanks to the HP 35s’ support for complex numbers, the matrix multi-tool works with complex matrices too.

Using the Program

Before running the matrix multi-tool, make sure flag zero is clear by pressing yellow shift FLAGS CF 0 (the program doesn’t clear this flag itself because it is used to maintain information between invocations).

Start the program by pressing XEQ M ENTER. This will display the main menu, which looks like this:

 
1A2b 3SoL4Un5×

The menu choices are:

  1. Go to dimension/edit/view menu for matrix A.
  2. Go to edit/view menu for column vector b.
  3. Compute inverse and determinant of A, and solve Ax=b for x.
  4. Unsolve, restoring original A and most recent b.
  5. Multiply matrix A by vector b.

Pressing R/S without selecting a choice exits the matrix program (this is important, to ensure flag 10 is cleared so equations will work in other programs).

Entering Matrix A

Press 1 R/S from the main menu to display the matrix A menu:

 
1DIM 2ED 3VIEW

The menu choices are:

  1. Specify the dimension N for N×N matrix A and N-vector b.
  2. Edit the entries in matrix A.
  3. View the entries in A.

Pressing R/S without selecting a choice will return you to the main menu.

Specifying the Matrix Dimension

Press 1 R/S to specify the dimensions of N×N matrix A, and consequently N-vector b. The program will display:

N=            
3.0000

The currently specified dimension is displayed. Enter the desired dimension and press R/S (or just press R/S to keep the existing dimension). The lowest allowed dimension is 2 and the highest depends on the amount of free memory. If the matrix multi-tool program is the only one in memory, you can work with a matrix of dimension up to 19×19.

If you enter a dimension that is too small, the program will display N TO SMALL. Pressing R/S will return you to the matrix A menu. If you enter a dimension that is too large, the calculator will display INVALID (I) and the program will halt. You’ll have to restart it by pressing XEQ M ENTER.

After you’ve specified the dimension, the program will automatically proceed to the matrix editing mode described below.

Editing the Matrix A

To edit the entries of A, press 2 R/S from the matrix A menu. The program will begin in the upper left corner of the matrix and proceed left to right, top to bottom (the same way you read English text). For each entry, it will first briefly display the indices of that entry (the number above the indices is N):

2.0000        
[1.0000,1.0000]

After a brief pause, execution will automatically continue, and the program will display the existing value and prompt for the new value of the entry whose indices were just displayed:

A?            
1.2345

To keep the existing value, just press R/S. Otherwise, enter the new value before pressing R/S. If you’ve forgotten which entry you are about to change, you can press the xy key to retrieve the row and column indices to the display. If you’ve already performed calculations that have altered the stack, just press R a few times. The indices are probably still on the stack (as a bracketed vector).

At any time during editing, you can change which entry is being edited. When the program is stopped awaiting a new value for an entry, you can store new row and column indices into memory registers R and C respectively. When you then enter a new value and press R/S, it will be stored in the newly specified location and editing will resume from there.

After editing the last entry of A, the program will return to the matrix A menu.

Viewing the Matrix A or A-1

To examine A (or result matrix A-1), press 3 R/S from the matrix A menu. The program will proceed in the same way as during editing, first briefly displaying the indices of each entry. After a brief pause, the program will display the value of the entry whose indices were just displayed:

A=            
-3.1416

Press R/S to proceed to the next entry. If you store new row and column indices into memory registers R and C while the program is stopped, the next entry displayed will be the one whose indices you specified, and viewing will continue from there.

When you’ve examined the last entry in A, the program will return to the matrix A menu.

Entering the Vector b

Pressing 2 R/S from the main menu will display the vector b menu:

              
2ED 3VIEW

The menu choices are:

  1. Edit the entries in vector b.
  2. View the entries in b.

Note that there is no choice number 1. The choices are numbered 2 and 3 for consistency with the menu for matrix A.

Pressing R/S without selecting a choice will return you to the main menu.

Editing the Vector b

Press 2 R/S from the vector b menu to edit the entries of b. The program will begin at the top of b. Editing b proceeds the same way as editing A, except that there is only a single index (the row). During editing, you can change which entry is being edited by storing a new row number in memory register R.

After editing the last entry of b, the program will return to the vector b menu.

Viewing the Vector b or x

To examine b (or solution vector x after solving the linear system), press 3 R/S from the vector b menu. Viewing b (or x) proceeds the same way as viewing A, except that there is only a single index. During viewing, you can change which entry is to be viewed next by storing a new row number in memory register R.

Solving the System Ax=b

If you’re at the matrix A or vector b menus, press R/S to return to the main menu. Then press 3 R/S to solve the system of equations. First, the program will display a time estimate (in seconds):

T=            
4.9680

Press R/S to continue. After a while (a long while for a really large system), the matrix multi-tool will display the computed determinant:

D=            
1278.3225

Press R/S to return to the main menu. The results can be found as follows:

  • A-1 will have replaced A. You can examine A-1 by pressing 1 R/S at the main menu, and then 3 R/S at the matrix A menu.
  • Solution vector x will have replaced b. You can examine x by pressing 2 R/S at the main menu, and then 3 R/S at the vector b menu.
  • The determinant is stored in memory register D, which you can examine by pressing RCL D.

Sometimes there is no solution. If this turns out to be the case, the program will display the message SINGULAR. Pressing R/S will then return you to the main menu.

Solving for Another Vector b

You can solve additional systems for the same matrix A and different vectors b by entering new values for the entries of b, and then running the solver again. So long as A-1 has not been overwritten by editing A, the matrix multi-tool will solve these additional systems by just computing x = A-1b, which is much faster (O(N2)) than the full solution process (O(N3)).

Restoring A and b

If you have not overwritten A-1 or x since the last system you’ve solved, you can undo the solution by pressing 4 R/S from the main menu. This will restore matrix A and the most recently entered vector b. If you try to perform the undo operation when A-1 or x have been modified, the program displays the message CANNOT UNDO and returns to the main menu.

Performing Matrix-Vector Multiplication

The built-in matrix-vector multiplication routine that the matrix multi-tool uses to solve additional systems when matrix A hasn’t changed is also available directly for your use. After entering A and b as described earlier, pressing 5 R/S from the main menu will compute the matrix-vector product of A and b.

The result vector will replace b, and can be examined the same way b can, by pressing 2 R/S at the main menu, and then 3 R/S at the vector b menu. The matrix A is unaltered. Another multiplication of A by a different b vector can be carried out by entering the new vector and invoking the multiplication.

Note that the undo operation (4 R/S from the main menu) does not apply to matrix-vector multiplication.

Example

Solve the following system of equations:

3.8x1 + 7.2x2 = 16.5

1.3x1 – 0.9x2 = -22.1

Keystrokes Display Description
 XEQ M ENTER   1A2b 3SoL4Un5×  Start the matrix multi-tool
 1 R/S   1DIM 2ED 3VIEW  Select matrix A
 1 R/S   N?  Select the DIMension command
 2 R/S   [1,1]
 A? 
Dimension is 2×2
 3.8 R/S   [1,2]
 A? 
A11=3.8
 7.2 R/S   [2,1]
 A? 
A12=7.2
 1.3 R/S   [2,2]
 A? 
A21=1.3
 0.9 +/- R/S   1DIM 2ED 3VIEW  A22=-0.9
Return to A menu
 R/S   1A2b 3SoL4Un5×  Return to main menu
 2 R/S   2ED 3VIEW  Select vector b
 2 R/S   [1]
 B? 
EDit command
 16.5 R/S   [2]
 B? 
b1=16.5
 22.1 +/- R/S   2ED 3VIEW  b2=-22.1
Return to b menu
 R/S   1A2b 3SoL4Un5×  Return to main menu
 3 R/S   T= 1.4720  Start solver and display time estimate
 R/S   D= -12.7800  Continue solver and display determinant 
 R/S   1A2b 3SoL4Un5×  Solution completed
 2 R/S   2ED 3VIEW  Select vector b (now x)
 3 R/S   [1]
 B= -11.2887 
Select the VIEW command
x1=-11.2887
 R/S   [2]
 B= 8.2496 
x2=8.2496
 R/S   2ED 3VIEW  Return to b menu
 R/S   1A2b 3SoL4Un5×  Return to main menu

Performance

Compared to the built-in matrix operations of an HP-15C or HP-42S, the matrix multi-tool is quite slow, but still quite usable for real problems. Furthermore, due to the larger built-in memory of the HP 35s, it can solve larger problems than the 15C or 42S can. The following table gives the approximate times for solving different sizes of systems using the matrix multi-tool:

Size (N) Solve Ax = b Compute x = A-1b
or Multiply Ab
2 4 sec 1 sec
3 8 sec 2 sec
4 14 sec 3 sec
5 26 sec 4 sec
6 43 sec 5 sec
7 1 min 7 sec
8 1½ min 8 sec
9 2 min 10 sec
10 3 min 13 sec
11 4 min 15 sec
12 5½ min 18 sec
13 7 min 21 sec
14 8½ min 25 sec
15 10½ min 28 sec
16 13 min 32 sec
17 15½ min 36 sec
18 18 min 41 sec

Program Listing

Line Instruction Comments
M001♦   LBL M  Display main menu
M002  CLx   
M003  SF 10   
M004♦   EQN 1A2b 3SoL4Un5×  Text of main menu (see note 1)
M005  CF 10   
M006  x=0?   
M007  RTN   
M008  1   
M009  x=y?   
M010  GTO M036  Go to dimension/enter/view A menu
M011  R↓   
M012  2   
M013  x=y?   
M014  GTO M057  Go to enter/view b menu
M015  R↓   
M016  3   
M017  x≠y?   
M018  GTO M021  Not the solve command
M019  XEQ M208  Solve Ax = b
M020  GTO M001  Return to main menu
M021♦   R↓   
M022  4   
M023  x≠y?   
M024  GTO M027  Not the undo command
M025  XEQ M183  Restore A and b
M026  GTO M001  Return to main menu
M027♦   R↓   
M028  5   
M029  x≠y?   
M030  GTO M034  Not the multiply command
M031  XEQ M396  Multiply A and b
M032  CF 2  Vector b has been overwritten
M033  GTO M001  Return to main menu
M034♦   XEQ M179  Invalid command
M035  GTO M001  Return to main menu
M036♦   CLx   
M037  SF 10   
M038  EQN 1DIM 2ED 3VIEW   
M039  CF 10   
M040  x=0?   
M041  GTO M001  Return to main menu
M042  SF 1   
M043  1   
M044  x=y?   
M045  GTO M074  Dimension A
M046  R↓   
M047  2   
M048  x=y?   
M049  GTO M095  Edit A
M050  R↓   
M051  CF 1   
M052  3   
M053  x=y?   
M054  GTO M095  View A
M055  XEQ M179  Invalid command
M056  GTO M036  Return to A menu
M057♦   CLx   
M058  SF 10   
M059  EQN 2ED 3VIEW   
M060  CF 10   
M061  x=0?   
M062  GTO M001  Return to main menu
M063  SF 1   
M064  2   
M065  x=y?   
M066  GTO M141  Edit b
M067  R↓   
M068  CF 1   
M069  3   
M070  x=y?   
M071  GTO M141  View b
M072  XEQ M179  Invalid command
M073  GTO M057  Return to b menu
M074♦   INPUT N  Prompt for new size (default is old size)
M075  IP   
M076  2   
M077  x>y?  Is dimension at least 2?
M078  GTO M081  Dimension is too small
M079  XEQ M085  Set up memory
M080  GTO M095  Proceed to matrix A editing mode
M081♦   SF 10  Dimension is too small
M082  EQN N TOO SMALL   
M083  CF 10   
M084  GTO M036  Return to A menu
M085♦   RCL N  Compute fence location 2N^2+N = N(2N+1)
M086  ENTER   
M087  ENTER   
M088  +  2N
M089  1   
M090  +  2N+1
M091  x  N(2N+1)
M092  STO I   
M093  STO (I)  Set up matrix end-of-memory fence
M094  RTN   
M095♦   1  Initialize row and column indices
M096  STO R   
M097  STO C   
M098♦   XEQ M133  Compute index of A[R,C]
M099  STO I   
M100  RCL R   
M101  RCL N   
M102  x<y?  Is it past the bottom of A?
M103  GTO M036  If so, return to A menu
M104  EQN [R,C]  Flash [row,col] on display
M105  PSE   
M106  FS? 1  Are we editing?
M107  GTO M119  Skip pre-increment of R,C
M108♦   1  Increment row/column before view
M109  STO+ C   
M110  RCL C  Is C still within [1..N]?
M111  RCL− N   
M112  x≤0?   
M113  GTO M116  If yes, don’t increment row
M114  STO C  Reset column to 1
M115  STO+ R  Compute next row index
M116♦   FS? 1  Were we called by edit’s post-increment?
M117  RTN  If so, return
M118  REGZ  Bring [R,C] back into x register (see note 2)
M119♦   RCL (I)  Fetch matrix element
M120  STO A  Make it the default
M121  FS? 1   
M122  GTO M125  Edit entry
M123  VIEW A   
M124  GTO M098  Skip editing code
M125♦   INPUT A  Prompt for A[R,C]
M126  CF 0  Previous solve is no longer undoable
M127  XEQ M133  Recompute index in case user changed it
M128  STO I   
M129  R↓  Retrieve value to store in A[R,C]
M130  STO (I)   
M131  XEQ M108  Increment row/column after edit
M132  GTO M098  Process next entry of A
M133♦   RCL C  Compute index of A[R,C]
M134  1   
M135  −   
M136  RCL× N  N(C-1)
M137  RCL+ R   
M138  1   
M139  −  N(C-1) + R-1
M140  RTN   
M141♦   1  Initialize row index
M142  STO R   
M143♦   XEQ M171  Compute index of b[R]
M144  STO I   
M145  RCL R   
M146  RCL N   
M147  x<y?  Is R past the end of b?
M148  GTO M057  If so, return to b menu
M149  EQN [R]  Flash [row] on display
M150  PSE   
M151  FS? 1  Are we editing?
M152  GTO M156  Skip pre-increment of R,C
M153  1  Increment row before view
M154  STO+ R   
M155  R↓  Bring [R] back into x register
M156♦   RCL (I)  Fetch vector element
M157  STO B  Make it the default
M158  FS? 1   
M159  GTO M162  Edit entry
M160  VIEW B   
M161  GTO M143  Skip editing code
M162♦   INPUT B  Prompt for b[R]
M163  CF 2  Previous solve no longer undoable
M164  XEQ M171  Recompute index in case user changed it
M165  STO I   
M166  R↓  Retrieve value to store in b[R]
M167  STO (I)   
M168  1   
M169  STO+ R  Increment row after edit
M170  GTO M143  Process next entry of b
M171♦   RCL N  Compute index of b[R]
M172  x2   
M173  ENTER   
M174  +  2N^2 – beginning of b
M175  RCL+ R   
M176  1   
M177  −  2N^2 + R-1
M178  RTN   
M179♦   SF 10   
M180  EQN INVALID CMD   
M181  CF 10   
M182  RTN   
M183♦   FS? 0  Flag 0 has to be set to be able to undo
M184  GTO M186  If it is, go check flag 2
M185  GTO M188  Can’t undo if flag 0 not set
M186♦   FS? 2  Flag 2 also has to be set for undo
M187  GTO M192  If it is, use Gaussian elimination to undo
M188♦   SF 10  Display error message
M189  EQN CANNOT UNDO   
M190  CF 10   
M191  RTN  Return to main menu
M192♦   XEQ M215  Call Gaussian elimination routine
M193  RCL D  Restore determinant to pre-undo value
M194  1/x   
M195  STO D   
M196  CF 0  Clear flag 0 since we’ve unsolved back to A
M197  CF 2  Clear flag 2 since we’re back to b
M198  RTN  Return to main menu
M199♦   RCL N  Routine to set up I, J, and E for solve or multiply
M200  x2  N^2
M201  STO I  Set I to N^2 for counting
M202  ENTER   
M203  +   
M204  STO J  Set J to 2N^2 for indexing
M205  RCL+ N  2N^2+N
M206  STO E  Set end of memory pointer to 2N^2 + N
M207  RTN   
M208♦   FS? 0  Do we still have A’ intact?
M209  GTO M392  Use matrix-vector multiplication to solve
M210  XEQ M215  Call Gaussian elimination
M211  SF 0  Set flag 0 indicating we have A’
M212  SF 2  Set flag 2 indicating we have x
M213  VIEW D  Display determinant
M214  RTN   
M215♦   RCL N  Estimate time for Gaussian elimination (see note 3)
M216  3   
M217  yx   
M218  0.184   
M219  x   
M220  STO T   
M221  VIEW T  Display estimate and make program stoppable
M222♦   XEQ M199  Set up I, J, and E for solver
M223  CLx  Fill temporary matrix with 0
M224♦   DSE J  Decrement J (will never skip)
M225  STO (J)   
M226  DSE I  Decrement count
M227  GTO M224  Fill in next 0 entry
M228♦   1  Fill diagonal of temporary matrix with 1
M229  STO (J)   
M230  RCL+ N  N+1
M231  STO+ J  Next diagonal entry of temporary matrix
M232  RCL J   
M233  RCL E   
M234  x>y?   
M235  GTO M228  Fill in next 1 entry
M236  1   
M237  STO D  Initialize determinant to 1
M238  CLx   
M239  STO K  Index of first diagonal entry is 0
M240  RCL N  N >= 2
M241  STO R  Number of rows remaining to process
M242♦   RCL R  Start of Gaussian elimination
M243  1   
M244  x≥y?   
M245  GTO M325  There’s nothing to do for the last row
M246  −   
M247  STO R   
M248  STO P  Find row with largest value in current column
M249  RCL K  Start at current diagonal entry
M250  STO I   
M251  STO J   
M252  RCL (I)  Save current entry as largest seen so far
M253  ABS   
M254  STO T   
M255♦   1  Point to next row
M256  STO+ I   
M257  RCL T   
M258  RCL (I)   
M259  ABS   
M260  x≤y?  Is new entry larger than largest so far?
M261  GTO M265  No, don’t save new entry
M262  STO T  Save new largest seen so far
M263  RCL I  Save index of largest seen so far
M264  STO J   
M265♦   DSE P   
M266  GTO M255  Check next row
M267  RCL J  Swap rows if necessary
M268  RCL K   
M269  x=y?  Is the largest in a row other than the current row?
M270  GTO M293  No, so don’t swap rows
M271  STO I  Point to diagonal entry of current row
M272  RCL T   
M273  x≠0?   
M274  GTO M279  Found a non-zero row; okay to swap
M275  SF 10   
M276  EQN SINGULAR   
M277  CF 10   
M278  RTN   
M279♦   RCL D  Swapping rows flips sign of determinant
M280  +/−   
M281  STO D   
M282♦   RCL E   
M283  RCL J   
M284  x≥y?  Is column index past end of memory?
M285  GTO M293  Yes, done swapping rows
M286  x↔ (I)  Swap entries between (I) and (J)
M287  x↔ (J)   
M288  x↔ (I)   
M289  RCL N  Increment pointers to next column
M290  STO+ I   
M291  STO+ J   
M292  GTO M282  Swap next column of rows
M293♦   RCL R  Subtract multiple of current row from remaining rows
M294  STO P   
M295♦   RCL K   
M296  STO I  Point to diagonal entry of current row
M297  RCL+ P   
M298  STO J  Point to same column in Pth row
M299  RCL (J)  Compute multiplier so entry becomes zero
M300  RCL÷ (I)   
M301  STO T   
M302  CLx  Clear entry in starting column
M303  STO (J)   
M304♦   RCL N  Increment to next column in both rows
M305  STO+ I   
M306  STO+ J   
M307  RCL E   
M308  RCL J   
M309  x≥y?  Is column index past end of memory?
M310  GTO M315  Yes, done subtraction
M311  RCL (I)  Compute multiple of current row
M312  RCL× T   
M313  STO− (J)  Subtract from same column of Pth row
M314  GTO M304  Proceed to next column
M315♦   DSE P   
M316  GTO M295  Proceed to next row
M317  RCL K  Update determinant for current row
M318  STO I   
M319  RCL (I)   
M320  STO× D   
M321  1  Point to next diagonal entry
M322  RCL+ N   
M323  STO+ K   
M324  GTO M242  Process next row
M325♦   RCL K  Update determinant for last row
M326  STO I   
M327  RCL (I)   
M328  STO× D   
M329  RCL N  Start of back substitution
M330  STO R  Begin with last (Nth) row
M331♦   RCL K  Divide row R by its diagonal entry
M332  STO J   
M333  RCL (J)   
M334  STO T   
M335  1  Diagonal entry becomes 1
M336  STO (J)   
M337  STO P  Initialize counter to be used later
M338♦   RCL N  Divide rest of row R by diagonal entry in T
M339  STO+ J  Increment pointer to next column
M340  RCL E   
M341  RCL J   
M342  x≥y?  Are we past the end of memory?
M343  GTO M347  Yes, done dividing this row
M344  RCL T  Divide entry by T
M345  STO÷ (J)   
M346  GTO M338  Proceed to next column of row R
M347♦   DSE R  Compute address of previous row
M348  GTO M363  Process previous row
M349  RCL N  Done back substitution; swap result into A
M350  x2   
M351  STO K  Have to move N^2 entries
M352  STO J  Point to beginning of temporary matrix
M353  CLx   
M354  STO I  Point to beginning of A
M355♦   RCL (J)  Copy entry from temporary to A
M356  STO (I)   
M357  1  Point to next entry in each
M358  STO+ I   
M359  STO+ J   
M360  DSE K   
M361  GTO M355  Copy next entry
M362  RTN   
M363♦   RCL K  Get index of diagonal entry of row R
M364  STO I   
M365  RCL− P   
M366  STO J  Index of corresponding entry in row R-P
M367  RCL (J)  Fetch and save to use as a multiplier
M368  STO T   
M369  CLx  Element in this column will become 0
M370  STO (J)   
M371♦   RCL N  Update indices to next column in both rows
M372  STO+ I   
M373  STO+ J   
M374  RCL E   
M375  RCL J   
M376  x≥y?  Are we past the end of memory?
M377  GTO M382  Done subtracting from row R-P
M378  RCL (I)  Subtract multiple of row R from row R-P
M379  RCL× T   
M380  STO− (J)   
M381  GTO M371  Proceed to next column
M382♦   1  Point to next row to subtract from
M383  STO+ P   
M384  RCL P   
M385  RCL R   
M386  x≥y?   
M387  GTO M363  Process next row to subtract from
M388  1  Point to previous row’s diagonal entry
M389  RCL+ N   
M390  STO− K   
M391  GTO M331  Process previous row
M392♦   XEQ M396  Call matrix multiplication routine
M393  SF 0  Set flag 0 indicating we have A’
M394  SF 2  Set flag 2 indicating we have x
M395  RTN   
M396♦   RCL N  Quick solve by multplication A’b = x
M397  x2  Estimate how long this will take (see note 3)
M398  0.125   
M399  x   
M400  STO T   
M401  VIEW T  Display estimate and make program stoppable
M402♦   XEQ M199  Set up I, J, and E for multiplication
M403  RCL N   
M404  STO R  Start one past the last row
M405♦   CLx   
M406  STO T  Accumulate each dot product in T
M407  RCL R   
M408  1   
M409  −   
M410  STO I  Point to row R of A
M411  RCL E   
M412  RCL− N   
M413  STO J  Point to first row of b
M414♦   RCL (I)  Inner multiplication loop
M415  RCL× (J)   
M416  STO+ T   
M417  RCL N   
M418  STO+ I  Next column of A
M419  1   
M420  STO+ J  Next row of b
M421  RCL E   
M422  RCL J   
M423  x<y?  Are we past the end of memory (and thus b)?
M424  GTO M414  Process next column of A and row of b
M425  RCL T   
M426  STO (I)  Store result in empty space between A and b
M427  DSE R   
M428  GTO M405  Process next row of A
M429  RCL N  Move result into b
M430  STO− J   
M431♦   RCL (I)   
M432  STO (J)   
M433  1  Increment row counters
M434  STO+ I   
M435  STO+ J   
M436  RCL E   
M437  RCL J   
M438  x<y?  Are we past the end of memory (and thus b)?
M439  GTO M431  If not, move next row
M440  SF 0  Can once again undo the calculation
M441  SF 2   
M442  RTN  Return to main menu

Length: 1455, Checksum: 9A4A

Note 1: The keystrokes to enter line M004 are: EQN  1  RCL A  2  yellow shift  L.R.  5  blue shift  SPACE  3  RCL S  blue shift  BASE  7  RCL L  4  RCL U  blue shift  SUMS  1  5  ×  ENTER

Note 2: To enter the REGZ instruction, press any function key that displays a menu of choices (such as yellow shift FLAGS), press the R key, and select z.

Note 3: There is a bug in the HP 35s (and HP 33s) firmware that if a program last stopped to display an EQN message, it will be uninterruptible until it stops for some other reason (a VIEW, INPUT, or STOP instruction). By computing a time estimate, we have an excuse to use a VIEW instruction (to display the time estimate) before entering into a potentially long-running computation.

Registers and Flags

Register Use
A Input/display of current entry of matrix A
B Input/display of current entry of vector b
C Column index (only during editing of A)
D Determinant
E Address of register beyond end of matrix memory
I, J General purpose indexing
K Loop control and diagonal traversal
N Matrix dimensions (NxN)
P Secondary row index
R Current row index
T Temporary accumulator and time estimate display
0…N2-1 Storage for matrix A and inverse matrix A-1
N2…2N2-1 Temporary matrix used for inversion
2N2…2N2+N-1 Column vector b and solution vector x
2N2+N End of memory fence (to ensure memory stays allocated if b[N] = 0)
Flag Meaning
0 Indicates memory contains A-1 instead of A. In conjunction with flag 2, indicates that undo is possible.
1 Set when editing, cleared when viewing
2 Indicates memory contains x instead of b. In conjunction with flag 0, indicates that undo is possible.

Revision History

2007-Nov-07 — Initial release.

2007-Nov-18 — Made the matrix dimensioning function and solver callable so they can be used by other programs I might write.

2007-Nov-20 — Compute and display time estimate before invoking solvers. This works around the HP 35s bug that programs are not interruptible if the last thing displayed was an EQN message, which can be disconcerting during an 18 minute matrix inversion.

2008-Apr-03 — Made the matrix multiplication routine accessible from the main menu (5 R/S), and also callable as a subroutine from other programs.

Calling Matrix Multi-Tool Subroutines from Other Programs

Several subroutines in the matrix multi-tool can be readily used from other programs. Note that the line numbers of these routines might change as I post future revisions to this program.

M085: Set Up Matrix and Vector Storage

Sets up memory to hold an N×N matrix A, an N×N scratch matrix, and an N-element vector b. The dimension N is expected to be in register N.

On return, register I will contain 2N2+N, the indirect address of the first register after the allocated memory, and that register will also contain its own address (to ensure that all the required memory remains allocated even if the last entries in b are zero).

M133: Compute Index of Arc

Given row and column indices in registers R and C respectively, compute the index of Arc in the memory allocated above. Requires that register N still contains the dimension, N, of matrix A and vector b.

M171: Compute Index of br

Given a row index in register R, compute the index of br in the memory allocated above. Requires that register N still contains the dimension, N, of matrix A and vector b.

M222: Solve Ax=b or Compute A-1 by Gaussian Elimination

Solves the system Ax=b for matrix A and vector b. The dimension must be in register N. Matrix A is expected to be stored column-wise in indirect registers 0 through N2-1. Vector b is expected to be stored in indirect registers 2N2 through 2N2+N-1.

On return, A will have been replaced by its inverse (A-1) and b will have been replaced by x. The determinant of A will be in register D. This subroutine also uses registers E, I, J, K, P, R, and T, and indirect registers N2 through 2N2-1.

An alternate entry point for this subroutine is M215, which will first display a time estimate (T=…). The Gaussian elimination will commence only after the user presses R/S. This works around an HP 35s firmware bug where a program is uninterruptible if the last thing displayed was an EQN message.

M402: Compute Matrix-Vector Product Ab

Computes the matrix-vector product Ab for matrix A and vector b. The dimension must be in register N. Matrix A is expected to be stored column-wise in indirect registers 0 through N2-1. Vector b is expected to be stored in indirect registers 2N2 through 2N2+N-1.

On return, A will be unchanged and b will have been replaced by the product Ab. This subroutine also uses registers E, I, J, R, and T, and indirect registers N2 through N2+N-1.

An alternate entry point for this subroutine is M396, which will first display a time estimate (T=…). The matrix-vector multiplication will commence only after the user presses R/S. Again, this is useful for working around the uninterruptible program bug.

Why I Wrote this Program

One of my first programmable calculators was a Texas Instruments SR-52 that I purchased at Eaton’s for $149 back in 1978 or so. At that time, the TI-59 had already been introduced and the store was clearing out the SR-52. This calculator had 224 steps of program memory, 20 registers, and a magnetic card reader. It was only slightly less powerful than the HP-67.

When I was in the 10th or 11th grade, we learned how to solve systems of linear equations using matrices, Gaussian elimination, and back substitution. After learning the method, it became clear to me that it would be easy to do on a computer, and I decided I’d try to program the SR-52 to do it. My teacher (who also taught computer science), was all in favour of the idea, and even said I could use the program during a test. He figured if I was smart enough to program it, especially in only 224 program steps, I obviously knew the algorithm inside out and backwards. Well, I managed to fit it onto the SR-52, although it was limited to 2×2 and 3×3 systems due to the limited number of storage registers. Of course with only 224 steps, there was no room for a user interface; I had to manually store the matrix elements into the appropriate registers.

Now, about 28 years later, an expanded version of this program came to mind as an ideal second project on the HP 35s, especially since my earlier HP calculators (41CX with Advantage Pac, and 42S) had this capability built in. This program will also prove useful as a subroutine in a polynomial curve fitting program I have planned for the future.

Related Articles

If you've found this article useful, you may also be interested in:

60 Comments

  1. Palmer O. Hanson, Jr.
    January 17, 2008

    The UNDO function of this program does NOT return the machine to the exact input values. The UNDO function in machines such as the HP-28S does return the machine to the exact input values.

    Matrix A and the inverse matrix are stored in registers 0 through N^2 – 1 as indicated in the text. However, although the elements of the input are entered row by row the elements are stored column by column, e.g., for an nth order matrix A element A(2,1) will be stored in register 1, element A(n,1) will be stored in register n-1, element A(1,2) will be stored in register n, and so on.

  2. Stefan Vorkoetter
    January 17, 2008

    Palmer, you are correct on both counts. The UNDO function re-inverts the inverse matrix, restoring a close approximation of the original matrix. This may not be the way a 28S does it, but it is the way the Advantage Pac matrix program works.

    The matrix elements are stored in column order in memory because it made the index arithmetic simpler (and therefore faster). Although I’m a C programmer and am used to storing things in row order, there’s nothing conceptually different in storing them in column order (which, for example, FORTRAN does).

  3. Palmer O. Hanson, Jr.
    January 30, 2008

    I need a routine for the HP-35s which will multiply an nxn matrix by a vector. It occurred to me that there must be such a routine buried in your matrix program which is used to multiply the inverse of the a matrix by the B vector. I found the notation at M383 of your program which states "Quick solve by multiplication A’b = x". So I installed my matrix in A and my vector in B using your input routine, but rather than solving I did XEQ M383. I got the first element of the product but I have not been able to find out how to get the remaining elements of the product. Can you help me? Even better, can you add to the capability of your program to provide this function as one of the options?

  4. Stefan Vorkoetter
    January 30, 2008

    Palmer, the product replaces the original b vector, so you can use the method for viewing b to examine the result (i.e. select B from the main menu, and then VIEW from the B menu).

    That’s a good idea about adding it as a function, since all the code is already there anyway. I’ll do that as soon as I get a chance.

  5. Stefan Vorkoetter
    February 05, 2008

    Palmer, I think I found the reason that a direct call to M383 isn’t working for you. Because the multiplication routine was designed as part of the whole program, it’s expecting certain values in a few registers, specifically:

    N – Number of rows (or columns) in A and length of vector b.

    E – End of memory pointer = 2N2+N.

  6. Palmer O. Hanson, Jr.
    March 05, 2008

    Stefan:

    Your discussion of why you wrote this program states that your program for the SR-52 was limited to 2×2 and 3×3 systems by memory limitations. Back in the heyday of the SR-52 some of the "superprogrammers" managed to get 4×4 in place. The following paragraph from my history of the matrix manipulation "friendly" competition gives references where the old issues of 52 Notes can be seen on Viktor Toth’s site.

    "The so-called friendly competition between users of HP and TI programmable calculators was announced for subscribers of 52 Notes in April 1977 (V2N4P1). The first problem which was agreed on was to mechanize a 5×5 determinant and inverse program where the competition would be between the SR-52/PC-100 and the HP-97. Work on calculating determinants and inverses had appeared in 52 Notes as early as August 1976 with the publication of a 4×4 determinant program by Dix Fulton (V1N3P4). The September 1976 issue published a program by Alan Trimble which would solve 4 simultaneous equations which would fit on one card (V1N4P2). The January 1977 issue presented a program by Rick Wenger which would solve 4 simultaneous equations, which was faster than that by Alan Trimble and which used fewer steps (V2N1P4/5). That issue also published a 4×4 determinant and inverse program by Barbara Osofsky. The February 1977 issue presented a program by Barbara Osofsy which would calculate the determinant and inverse of a 5×5 matrix, and also calculate a 4×4 determinant and solve a system of four simultaneous equations (V2N2P3/4). The May 1977 issue presented the latest version of Barbara Osofsky’s 5×5 matrix program (V2N5P5)."

    Palmer

  7. Rick
    April 01, 2008

    Can the HP-33s calculator use this program?

  8. Stefan Vorkoetter
    April 03, 2008

    Rick, the HP 33s doesn’t have the 801 indirectly-addressable registers that the 35s has. This program uses those registers to store the matrix in, so no, it won’t work on the 33s.

  9. Carlos T.
    February 01, 2009

    Tried several times to insert program. Checked all lines twice and checksum doesn’t match the one indicated. Lines are OK. Matrix multiplication OK, but determinant when resolving Ax=b is always -1 or 1. Checked subrutine from line 222 and all are OK. Can you tell me if porobably there is one thing that I am doing incorrectly or there is an error in the list herein indicated. Many thanks for your help. Regards\\carlos

  10. Stefan Vorkoetter
    February 01, 2009

    Unfortunately, the checksum isn’t all that useful. There’s a bug in the HP 35s firmware which makes the checksum inconsistent.

    When you try the example, does it solve the system correctly? If not, then the problem with the determinant is probably the result of a larger problem. If it does solve correctly, then probably only the determinant code itself is wrong. Is line M319 correct in your calculator (i.e. "RCL (I)", not "RCL I"). Also double check lines M321 to M323.

    I’m quite sure there are no problems in the program listing itself, since several people have entered the program from this listing, and had no problems.

  11. Pepe Fdez
    June 13, 2009

    Documentation speak about flags 0, 1 & 2. Program uses 0, 2 & 10 ones. Program running very well for me but my checksum (DD37) don´t concorde. Its possible that fix checksum problem change flag 10 to 1?

    Thanks & Regards.

  12. Stefan Vorkoetter
    June 14, 2009

    Pepe, the program uses flags 0, 1, 2, and 10. For example, line M042 sets flag 1.

    Flag 10 isn’t listed in the flags section because it’s not used to control the flow of the program. It’s just used to tell the calculator to display EQN messages instead of evaluating them (see the manual for details).

    The checksum mismatch can occur even when the program has been entered correctly, due to a bug in the HP-35s firmware.

  13. Dennis A
    January 29, 2010

    From the first step, my calculator says nonexistent. So I wonder if the program is right or my calculator is faulty.

  14. Stefan Vorkoetter
    January 29, 2010

    It sounds as if you haven’t typed the program into the calculator yet.

  15. juan
    March 02, 2010

    how put instruction m075 IP in a program

  16. Stefan Vorkoetter
    March 02, 2010

    Left-shift, INTG, 6. It’s in the manual on page G-9.

  17. mattia Tosi
    May 12, 2010

    I tried with N=2 and it works; with N=3,4 after the "time estimate", I press R/S and the program continues "RUNNING" and never ends…

  18. Stefan Vorkoetter
    May 12, 2010

    It sounds like you’ve made a mistake entering the program. It works for me and many other people who have tried it.

  19. Pablo
    May 22, 2010

    Hi Stefan, thanks for the program. You have made a good work.

    You can visit my web page where there are some simple programs I have written for my new HP 35s.

    blog.sortingthecomputer.com

  20. Turbato Tomas
    June 29, 2010

    Thx for the program Stefanv!

    I haven’t understand only one thing: i must use RPN mode for running the prgm…

    i’ve tried in alg mode for fun but don’t run…(have i wrong the prgm or is right? i’ve read 1.000 times all

  21. Stefan Vorkoetter
    June 29, 2010

    Yes, you need to use RPN mode.

  22. John
    July 17, 2010

    After having Maple crunch the closed forms, I wrote a 2×2, 3×3, adjoint, determinant, inverse, and multiply programs for my 35s. It’s plug dimension N, coefs A,B,C… the formulas chug out answers. It works fine, but I can’t go any higher than N=3 as I ran out of labels for input/output, and they total over 500 lines. Also, it doesn’t directly solve systems. Then I stumbled upon this marvelous program of yours. You do it all in one program, using an index counter for Row,Col and indirect labels is so much better, and shorter. It gave me a base to write a different program, a big hat’s off and thanks.

  23. Stefan Vorkoetter
    July 17, 2010

    What a coincidence! I actually used Maple during the development of this program. I first wrote it in Maple, then started to change the Maple code to use single letter variables, then hand-optimized the Maple code, until it was in a simple enough form to translate to HP-35s code. Of course, I do nearly everything in Maple, since I work for Maplesoft.

  24. Umer Jamshaid
    January 23, 2011

    Hey Stefan, I was trying to code the program into my calculator when I came across line M093. I’m not sure how to enter that, and was hoping you could help me.

    Thanks for the help, and thanks for the program!

  25. Stefan Vorkoetter
    January 26, 2011

    Umer, press STO and then the (I) key (i.e. the zero key).

  26. abel gonzalez
    June 20, 2011

    estas chido pero echemen la mano con un programa para mi hp 35s para calcular una poligonal abierta calculando puntos sobre un mismo azimut

  27. Chris Copela
    February 15, 2012

    I’m a new owner of a 35s. I found your program and it looks great. I wanted to ask one question though. What effect would inserting CF 0 into line 2 have? I ask because there could be long stretches of time between uses of the program and I could forget to do it. (If it matters, I wouldn’t mind losing anything I keyed into the program before I turned the calculator off last.) Thanks very much.

  28. Francesco Baldrighi
    May 08, 2012

    I think I found a bug.
    If you want to multiply a matrix (not solve a system of equations)for a vector your program set the flags 0 and 2 (line M440 and M441) and so the program think A is the inverse matrix A’ but it’s not true.
    A is the original matrix, untouched.
    Now, if you want to solve a system of equations with the precedent matrix A the program executes the product thinking to use A’ because the flag 0 is on; obviously the solution is to clear the flag 0 manually before to solve the system or delete the two lines M440 and M441.
    I hope it helps.
    Thanks for your program.

  29. luca
    August 29, 2012

    How put instruction M286 M287 and M288? Please

  30. lucas
    September 05, 2012

    Is there a way to download this to my calculator so I don’t have to manually enter 442 lines of code?

    Thanks for your help.

  31. Bernardo
    September 09, 2012

    Lucas, you need to write all lines manually.

  32. lucas
    September 10, 2012

    Yikes.. that’s what I was afraid of!

    Thanks for getting back to me.

  33. Joss
    November 28, 2012

    How to insert line M339

  34. Joss
    November 28, 2012

    is the + sign really next to J and not to STO

  35. Stefan Vorkoetter
    November 28, 2012

    Joss, no, it was just a space in the wrong place (not that it matters, since you don’t enter the space anyway).

  36. Tom
    January 06, 2013

    Luca, regarding to M286, M287, M288 lines question:

    Press:
    “Right Shift” RCL 0
    “Right Shift” RCL .
    “Right Shirt” RCL 0

  37. JKiD
    January 16, 2013

    I have a problem whenever I’m at the step to enter the data in the b matrix. When I press 2 R/S in the b matrix menu I’m being redirected to the same menu, it won’t let me write anything on the matrix I’m just stuck at the menu, I tried to view so I’m pressing 3 R/S in the b matrix menu, I have the same reaction as before: stuck in the b matrix menu.

  38. Stefan Vorkoetter
    January 16, 2013

    Check that lines M061 through M073 of the program were entered correctly. It sounds like the code that’s supposed to send you to the editing routines isn’t doing so.

  39. Philalethes
    January 19, 2013

    I’d like to be able to row reduce NxM matrices into RE form using Gauss-Jordan, per this page: http://vale.thus.ch/software/rref.html

    I would be interested to know you think your program could easily be modified to do this. Thanks!

  40. Dan Lewis
    March 03, 2013

    Very nice program. After I had checked it over a few times, I found that I had made some mistakes. After fixing them, it works like a charm! Even after fixing y mistakes, the checksums still didn’t match, but I guess there will always be some bugs. One thing I didn’t like about the 35s was its inability to work with matrices, and this program fits the bill. I can’t imagine ever having to solve an equation with 19 variables, but it’s nice to know that I could if the need ever arose. The undo feature is also a nice touch.

    Thanks,
    Dan

  41. Stefan Vorkoetter
    March 03, 2013

    I’m glad you like the program Dan. Regarding the checksum, it’s commonly acknowledged (for example, in the forums at that the checksums can be inconsistent between calculators, so it doesn’t necessarily mean there’s something wrong with the program you entered (perhaps that’s what you meant by “bugs”).

  42. Dan Lewis
    April 12, 2013

    Furthermore,

    I’m sure you know what mesh and nodal analysis are. This involves multiplying the inverse of a matrix by a vector. So, just fake solve a system of linear equations (to get the inverse of the matrix), then go back and put in the vector and voila! I had to use these analyses on the Fundamentals of Engineering Exam (on which the only HP allowed is the 35s) and I would have had to solve them by hand had I not had your program.

    Thanks again,
    Dan

  43. al penero
    December 14, 2013

    Mr stefanv could you please give an example on how to enter the program for N=3.tnx for help

  44. Stefan Vorkoetter
    December 14, 2013

    Al, the program is the same for all N. The value of N is specified when you run the program. Please look at the “Using the Program” section.

  45. jason foose
    April 01, 2014

    Great programming Stefan. A great addition to the 35s. Can you build a rectangular / polar conversion key too 😉 Thanks for sharing.

  46. Jason Foose
    April 04, 2014

    Stefan,

    I am writing an article for land surveyors regarding surface modeling. I would like to reference your Matrix multi tool in the article. Do you have any objections providing the link?

    I also am working to modify MMT (matrix multi tool) to accomdate the six term surface modeling equation and matrix found at http://www.ahinson.com/algorithms_general/Sections/InterpolationRegression/FitPolynomialSurface.pdf

    I have played with this to include more data points which requires a rectangular “A” matrix of “n” rows x 6 columns. The two modifications I would like to make are 1.) input program for x,y,z coordiantes that arranges the data in the six term format (i.e. X^2 XY Y^2 ect) then stores the refined data in the appropriate matrix designated registers. I think I can pull this off with my skill level. ANy thoughts or comments??? Thanks again for the wonderful tool. I really opens the door for the 35s.

  47. Marco
    June 18, 2014

    You did a great job. Thanks for that nice little and handy program. 🙂

  48. Carolina
    August 15, 2014

    Thanks 😀

  49. Martin
    September 03, 2014

    Stefan,

    nice work, reviving my HP calculator enthusiasm. Talking about Maple – wouldn’t it be a nice hobby project for you to create a CodeGeneration function which outputs HP 35s or better RPL code for HP48/49/50 calculators … ? Would this be feasible?

    Martin

  50. Mandla Skhosana
    October 06, 2014

    What are the keystrokes for the following? (M038, M059,M082, M180, M189 and M104)

  51. Stefan Vorkoetter
    October 06, 2014

    See page 13-16 of the HP 35s manual, “Using Equations to Display Messages”.

  52. Mandla Skhosana
    October 07, 2014

    Thanks Stefan, it took me almost 20hrs to programme, test and recheck the programm. Enventually it worked, its magic.My worry is that during exams they sometimes reset all the prommable calculators and I will loose the programme at the time i need. Is there a way this cannot be deleted

  53. Stefan Vorkoetter
    October 08, 2014

    Unfortunately, there’s no way to prevent a program from being deleted.

  54. mandla Skhosana
    October 08, 2014

    is the a way I can load an engineering equation like
    CTerror = Ie/(I’+Ie)

  55. Gerard Guerra
    February 26, 2015

    Stefan
    Thanks for the hard the program,the code is new to me I am using the manual as a reference, but still have a couple of questions. In M075 IP, what is the appropriate input for this, also are the EQN titles space sensitive? Is there a space key on my calculator? Thanks for your help.

    GGUERR

  56. Pedro Daniel Leiva
    October 31, 2015

    Dear Gerard:
    To insert IP in step M075, just press yellow keys; left arrow INTG 6
    To get space during equation typing press (blue keys): right arrow 0

    I hope it helps
    Pedro

  57. Raymond
    June 20, 2017

    Hi, Im going by the instructions. Xeq m, 1 rs, 1 rs, I choose the dimension, enter my matrix, after that when I go to 3 rs to view the inverse, I can see my first entry, after that its shows “invalid (I)”, any suggestion?

  58. Raymond
    June 20, 2017

    When I go to programing, it takes me to M093 = STO(I)

  59. Stefan Vorkoetter
    June 20, 2017

    All I can suggest is that you double check that the program was entered correctly.

  60. Stefan Vorkoetter
    June 20, 2017

    Somehow you’ve ended up back in the matrix entry routine with the index set up incorrectly, which suggests that the code branched to the wrong place from somewhere. Again, check that the program was entered correctly.

Leave a Comment

Want to see your picture next to your comments on this site and others? Visit gravatar.com to register your own globally recognized avatar.

Buy Stefan a coffee! If you've found this article
useful, consider leaving a donation to help support
stefanv.com

Disclaimer: Although every effort has been made to ensure accuracy and reliability, the information on this web page is presented without warranty of any kind, and Stefan Vorkoetter assumes no liability for direct or consequential damages caused by its use. It is up to you, the reader, to determine the suitability of, and assume responsibility for, the use of this information. Links to Amazon.com merchandise are provided in association with Amazon.com. Links to eBay searches are provided in association with the eBay partner network.

Copyright: All materials on this web site, including the text, images, and mark-up, are Copyright © 2017 by Stefan Vorkoetter unless otherwise noted. All rights reserved. Unauthorized duplication prohibited. You may link to this site or pages within it, but you may not link directly to images on this site, and you may not copy any material from this site to another web site or other publication without express written permission. You may make copies for your own personal use.