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

October 17, 2007

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 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 x↔y 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:

2. Edit the entries in vector b.
3. 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(N²)) than the full solution process (O(N³)).

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.8x₁ + 7.2x₂ = 16.5

1.3x₁ – 0.9x₂ = -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    L.R.  5    SPACE  3  RCL S    BASE  7  RCL L  4  RCL U    SUMS  1  5  ×  ENTER

Note 2: To enter the REGZ instruction, press any function key that displays a menu of choices (such as 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 2N²+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 A_rc

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

M171: Compute Index of b_r

Given a row index in register R, compute the index of b_r 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 N²-1. Vector b is expected to be stored in indirect registers 2N² through 2N²+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 N² through 2N²-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 N²-1. Vector b is expected to be stored in indirect registers 2N² through 2N²+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 N² through N²+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:

• Op-Amp Gain and Offset Design with the HP-41C Programmable Calculator
• Op-Amp Oscillator Design with the HP-41C Programmable Calculator
• Low-Sensitivity Sallen-Key Filter Design with the HP-41C Programmable Calculator
• Low-Sensitivity Sallen-Key Filter Design with the HP-67 Programmable Calculator
• Op-Amp Gain and Offset Design with the HP-67 Programmable Calculator
• Op-Amp Oscillator Design with the HP-67 Programmable Calculator
• Resistor Network Solver for the HP-67 Programmable Calculator
• Curve Fitting for the HP 35s Programmable Calculator