Row reduced echelon form on HP-41

This program transforms a matrix into row reduced echelon form. This means the program will calculate determinants and inverses and will solve systems of equations.\(\) Bits and pieces of this code originates PPC ROM Module [1,2,3,4,5,6,7,8] released by the international group People Programming Computers on December 1981. There the routine is known as RRM.

Method

A single routine calculates the row reduced echelon form. It handles the three matrix problems individually or simultaneously. It uses a technique known as partial pivoting with helps reduce round-off error.

The only limitation on the size of the matrices is the number of available data registers. The size can be increased by running MA in extended memory [9].

Instructions

The subroutines of the program map to the top row of keys on the USER keyboard.

Keyboard overlay for MA

Key functions:

keystrokes
operation key description
NEW ∑+ Enter a new Matrix
VIEW 1/x View a matrix
(y,x) √x Recall (y,x)
SOLVE LOG Solve the matrix
Det LN Show the determinant of the square coefficient matrix

Start by GTO “MA” and enabling the USER keyboard. Enter a new matrix using NEW. It will prompt with (r,c)? for each of the coefficients at row r and column c. After entering the coefficients, you press

  • VIEW to review the matrix. The number of decimal places shown is controlled by R05.
  • row ENTER column (y,x) to inspect any particular coefficient. The program will show the register number, followed by the coefficient. Knowing the register number allows you to change the coefficient directly.
  • SOLVE, to solve the matrix. After solving the matrix, press VIEW to view the result, Det to view the determinant of the square matrix. If zero, the system of equations may have no solutions or an infinite number of solutions.

Examples

Three examples illustrate the use of MA. You may have to claim some memory using at least “SIZE 031”.

Solve a system of equations

We start with problem 1 from the PPC ROM Users Manual. Solve the system of equations as shown below.

$$ \left. \begin{array}{rcrcrcr} -5x&+&10y&+&15z&=&5\\ 2x&+&y&+&z&=&6\\ x&+&3y&-&2z&=&13 \end{array} \right\}\Rightarrow\\[30pt] \left[ \begin{array}{rrr|r} -5 & 10 & 15 &5\\ 2&1&1&6\\ 1&3&-2&13 \end{array} \right] \times \left( \begin{array}{c} x \\ y \\ z \end{array} \right) $$

Use the calculator to determine the row reduced echelon form.

keystrokes
Keystrokes Display
GTO “MA”
USER
NEW REG.?
10 R/S R^C
3 ENTER  4 R/S (1,1)=?
5 CHS R/S (1,2)=?
10 R/S :
: (3,4)=?
13 R/S 13
SOLVE beep
Det D=150
VIEW (1,1)=1.000
(1,2)=0.000
:
(3,4)=-1.000

In other words the solution is

$$ \left[ \begin{array}{rrr|r} 1 & 0 & 0 & 2.000 \\ 0 & 1 & 0 & 3.000 \\ 0 & 0 & 1 & -1.000 \end{array} \right] \times \left( \begin{array}{c} x \\ y \\ z \end{array} \right) \Rightarrow \\[30pt] \left\{ \begin{array}{c r} x&=&-2.000\\ y&=&3.000\\ z&=&-1.000 \end{array} \right. $$

Inverse of matrix

Moving on to problem 2: find the inverse matrix by augmenting the original matrix by an identity matrix.

$$ A=\left( \begin{array}{rrr} 2&-3&1\\ 3&2&-1\\ 5&-2&1 \end{array} \right) \Rightarrow\\[30pt] \left[ \begin{array}{rrr|ccc} 2&-3&1&1&0&0\\ 3&2&-1&0&1&0\\ 5&-2&1&0&0&1 \end{array} \right] $$

Enter the matrix similar to our first example. After solving it, the identify matrix will be on the left, and the inverse matrix in the right square.

$$ \left[ \begin{array}{ccc|rrr} 1 & 0 & 0 & 0.000 & 0.125 & 0.125\\ 0 & 1 & 0 & -1.000 & -0.375 & 0.625 \\ 0 & 0 & 1 & -2.000 & -1.375 & 1.625 \end{array} \right] $$

In other words the inverse matrix is as shown below.

$$ A^{-1}=\left( \begin{array}{rrr} 0.000 & 0.125 & 0.125\\ -1.000 & -0.375 & 0.625 \\ -2.000 & -1.375 & 1.625 \end{array} \right) $$

Solve system of equations, find the inverse of the coefficient matrix

Akin to problem 3 in the PPC ROM Manual.

$$ \left. \begin{array}{rcrcrcr} 14x&+&2y&-&6z&=&9\\ -4x&+&y&+&9z&=&3\\ 6x&-&4y&+&3z&=&-4 \end{array} \right\} $$

Enter the coefficient matrix augmented by both the identity matrix and the final column of constants.

$$ \left[ \begin{array}{rrr|ccc|r} 14 & 2 & -6 & 1 & 0 & 0 & 9\\ -4 & 1 & 9 & 0 & 1 & 0 & 3\\ 6 & -4& 3 & 0 & 0 & 1 &-4 \end{array} \right] $$

The calculator will show the answer as shown below.

$$ \left[ \begin{array}{ccc|rrr|ccc|r} 1&0&0&0.063 & 0.029 & 0.039 & 0.500 \\ 0&1&0&0.107 & 0.126 & -0.165 & 2.000\\ 0&0&1&0.016 & 0.110 & 0.036 & 0.333 \end{array} \right] $$

The inverse of the original matrix is in the middle.

$$ A^{-1}=\left( \begin{array}{rrr} 0.063 & 0.029 & 0.039\\ 0.107 & 0.126 & -0.165\\ 0.016 & 0.110 & 0.036 \end{array} \right) $$

The last column contains the solutions of the system of equations.

$$ \left\{ \begin{array}{rr} x&=0.500\\ y&=2.000\\ z&=0.333 \end{array} \right. $$

Source code

  • Requires: X-Functions module on the HP-41cv
  • Available as source code, raw for the V41 emulator or bar code for the HP Wand (HP82153A)
  • The program takes 442 bytes, 4 magnetic tracks.
  • For line by line comments refer to [1,2,3,4,5,6,7,8].
  • Minimum SIZE 31.

Listing

Available through the repository

01      LBL "MA"
02      LBL D
03      CLA
04      AVIEW

05      .               ; store a 1 in R01 for the determinant
06      STO 03
07      STO 04
08      SIGN
09      STO 01
10      SF 10           ; for the BX routine
11      LBL 05          ; make R03 and R04 point to the next pivot position
12      ISG 03
13      LBL 06
14      ISG 04
15      ""
16      RCL 08          ; determine when the program ends by checking if
17      RCL 04          ;  either a row or column boundary has been crossed.
18      X>Y?
19      GTO E
20      RCL 09
21      RCL 03
22      X>Y?
23      GTO E
24      RCL 04          ; set up the block control word for BX
25      XEQ 17          ; matrix (i,j) to register address
26      X<> Z
27      XEQ 17          ; matrix (i,j) to register address
28       E3
29      /
30      +
31      RCL 08
32       E5
33      /
34      +
35      XEQ 18          ; find the pivot number and check if all the remaining
36      RCL IND M       ;  column entries are zero in which case the determinant
37      ST* 01          ;  must be zero and only the next column is incremented
38      X=0?            ;  by branching to LBL 06
39      GTO 06
40      1/X             ; make a 1 in the row containing the pivot number
41      RCL M
42      INT
43      XEQ 15          ; register address to (i,j)
44      RDN
45      STO 02
46      XEQ 19
47      RCL 02          ; check if the pivot number is already in the pivot
48      ST- 02          ;  position.
49      RCL 03
50      X=Y?
51      GTO 07
52      XEQ 12          ; row interchange to move the pivot to the true pivot
53      RCL 01          ;  position and adjust the sign of the determinant
54      CHS             ;  accordingly,
55      STO 01
56      LBL 07          ; make 0's in the current pivot column in all rows
57      ISG 02          ;  except the pivot row.
58      ""
59      RCL 09
60      RCL 02
61      X>Y?
62      GTO 05
63      RCL 03
64      X=Y?
65      GTO 07
66      RCL 02
67      RCL 04
68      XEQ 17          ; matrix (i,j) to register address
69      RDN
70      RCL IND T
71      CHS
72      XEQ 21
73      GTO 07
74      LBL A
75      "REG. ?"
76      PROMPT
77      STO 07
78      "R^C"
79      PROMPT
80      STO 08
81      X<>Y
82      STO 09
83      SF 09
84      GTO 01
85      LBL B
86      CF 09
87      LBL 01
88      CF 29
89      RCL 07
90      STO 04
91      RCL 08
92      RCL 09
93      *
94      STO 03
95      LBL 02
96      RCL 04
97      XEQ 15          ; register address to (i,j)
98      FIX 0
99      " ("
100     ARCL Y
101     >","
102     ARCL X
103     >")="
104     FC? 09
105     GTO 03
106     >"?"
107     TONE 89
108     PROMPT
109     STO IND 04
110     GTO 04
111     LBL 03
112     FIX IND 05
113     ARCL IND 04
114     AVIEW
115     LBL 04
116     ISG 04
117     ""
118     DSE 03
119     GTO 02
120     RTN
121     LBL C
122     XEQ 17          ; matrix (i,j) to register address
123     " R"
124     ARCL X
125     AVIEW
126     STO 04
127      E
128     STO 03
129     CF 09
130     GTO 02
131     LBL 16
132     X<>Y
133     STO O
134     X<>Y
135     MOD
136     ST- O
137     LASTX
138     ST/ O
139     CLX
140     X<> O
141     X<>Y
142     RTN

; MULTIPLY ROW BY CONSTANT (M2)
;
; on entry
;   the row number (i) in X
;   the constant (k) in Y

143     LBL 19
144     XEQ 13           ; compute the block control word (bbb.eee) for row i
145     X<>Y             ;
146     LBL 20           ; the constant k is then placed in X and M2 runs
147     ST* IND Y        ;  through a short loop to multiply each element of
148     ISG Y            ;  row i by the constant k
149     GTO 20
150     RTN

; ADD MULTIPLE OF ANOTHER ROW (M3)
;
; Add a constant multiple of one row in a matrix to another row.  The row
; that is multiplied does not change.

151     LBL 21
152     STO M           ; M is used for tmp storage
153     RDN
154     XEQ 13          ; compute the block control word (bbb.eee) for row i
155     X<>Y
156     XEQ 13          ; compute the block control word (bbb.eee) for row i
157     RCL M
158     SIGN            ; store k in LASTX
159     LBL 22          ; the main loop.
160     RDN
161     RCL IND Y
162     LASTX
163     *
164     ST+ IND Y
165     ISG Y
166     ""
167     ISG Z
168     GTO 22
169     RTN

; INTERCHANGE ROWS (M1)
;
; The matrix is assumed to be stored with each row occupying a consecutive
; block of registers.  The entire matrix is assumed to be stored row by row
; as one string of consecutive data registers.
; Like the other 4 matrix routines, M1 requires two stored values: R07 holds
; the starting register of the matrix, and R08 holds the number of columns in
; the matrix.  Both the row and column numbers start counting from 1.
; To interchange any two rows in the matrix enter the two row numbers in Y
; and X (the order is unimportant) and call M1.  M1 performs a block
; exchange of the two rows involved by dropping into the routine BE. 

170     LBL 12      ; Feeds into the block exchange routine BE after
171     XEQ 13      ;  setting up the two block control words for the
172     X<>Y        ;  two rows by calling LBL 13 twice.
173     XEQ 13

; BLOCK EXCHANGE (BE)
;
; On entry the two block control words are assumed to be in X and Y on the
; stack.

174     LBL 14
175     RCL IND Y
176     X<> IND Y   ; Perform the exchange on an element by element basis
177     STO IND Z   ;  as part of the loop.
178     RDN         ; Put the stack back in the correct config for the next
180     ""           ;  pass through the loop.
181     ISG X       ; Both X and Y are incremented, but note that only Y
182     GTO 14      ;  is tested.  (line 180 is a NOP.)
183     RTN

; COMPUTE THE BLOCK CONTROL WORD (BBB.EEE) FOR ROW I
;
;   bbb = s + c*(i-1)
;   eee = s + c*i - i

184     LBL 13        ; on entry
185     RCL 08        ;   the starting register (s) in R07
186     *             ;   the number of columns in the matrix (c) in R08
187     RCL 07        ;   the row number of the ith row (i) in X
188     +
189     RCL X
190     RCL 08
191     ST- Z
192     SIGN
193     -
194      E3
195     /
196     +
197     RTN

; M4 - register address to (i,j)
;
; Determine the row number i nd the column number j from the register
; number r by the following formulas using the starting register (s) of the
; matrix and the number of columns (c) in the matrix:
;   i = INT((r-s)/c) + 1
;   j = (r-s) MOD c + 1

198     LBL 15
199     RCL 07
200     -
201     RCL 08
202     XEQ 16          ; calculate i-1 and j-1
203     ISG Y           ; increment these values
204     ""                ; NOP
205     ISG X
206     ""                ; NOP
207     RTN

; MATRIX (I,J) TO REGISTER ADDRESS (M5)
;
; Determine the register number of the (i,j) element in a matrix (row i,
; column j).  On entry column number in X, and row number in Y.  Returns
; the register number in X

208     LBL 17          ; r = s + c*(i-1) + (j-1)
209     X<> 08
210     ST- 08
211     *
212     ST+ 08
213     X<> L
214     X<> 08
215      E
216     -
217     RCL 07
218     +
219     RTN

; BLOCK EXTREMA (BX)
;
; Finds the largest or smallest element of a block of registers.  By setting
; F10, only absolute values of the elements will be considered.  The maximum
; and minimum values as well as the their register numbers are returned:
;   Y: maximum value
;   X: minimum value
;   M: register of the maximum value (INT part)
;   N: register of the minimum value (INT part)

220     LBL 18
221     STO M
222     STO N
223     STO O
224     RCL IND X
225     FS? 10
226     ABS
227     ENTER^
228     ENTER^
229     RDN
230     LBL 08
231     CLX
232     RCL IND Z
233     FS? 10
234     ABS
235     X>Y?
236     GTO 10
237     R^
238     X>Y?
239     GTO 11
240     RDN
241     LBL 09
242     ISG Z
243     GTO 08
244     X<>Y
245     R^
246     RTN
247     LBL 10
248     X<>Y
249     CLX
250     RCL Z
251     STO M
252     GTO 09
253     LBL 11
254     CLX
255     RCL T
256     STO N
257     X<>Y
258     RDN
259     GTO 09
260     RTN
261     LBL E
262     "D="
263     ARCL 01
265     TONE 89
266     TONE  7
267     TONE 89
268     END

References

[1] M1 – Matrix, interchange rows
John Kennedy, December 1981
PPC ROM User’s Manual, page 260-265
[2] M2 – Matrix, multiply row by constant
John Kennedy, December 1981
PPC ROM User’s Manual, page 266-267
[3] M3 – Matrix, add multiple of another row
John Kennedy, December 1981
PPC ROM User’s Manual, page 268-269
[4] M4 – Matrix, add multiple of another row
John Kennedy, December 1981
PPC ROM User’s Manual, page 270-275
[5] M5 – Matrix, (i,j) to register address
John Kennedy, December 1981
PPC ROM User’s Manual, page 274-275
[6] QR – Quotent Remainder
John Kennedy & Roger Hill, December 1981
PPC-User Manual, pg. 372
[7] BE – Block Extrema
Richard Schwartz, December 1981
PPC ROM User’s Manual, pg. 54-55
[8] BX – Block
Richard Schwartz, December 1981
PPC ROM User’s Manual, pg. 68-69
[9] Running prgms in X-memory
Philip Karras, April 1982
PPC Journal, V9N3, pg.26

Leave a Reply

Your email address will not be published. Required fields are marked *

 

This site uses Akismet to reduce spam. Learn how your comment data is processed.