Householder

Householder Transformation for the HP-41

Overview

1°) Householder Algorithm
2°) A Modified Algorithm

-The 1st program uses the Householder method to triangularize one - or several - linear system(s)
-Then "TRS" is called to compute the solutions.

-In the 2nd program, the method is modified to diagonalize the principal matrix - if it is possible - so that the diagonal only contains "1"s
-No subroutine is needed in this case.

-In both cases, the elements are to be stored in column order, the right-hand side first.

B = A.X

1°) Householder Algorithm

-The Householder method uses orthogonal symmetries to gradually triangularize a matrix.

-If ( a₁ , ...... , a_n ) is the first column of the matrix A, the other columns ( x₁ , ...... , x_n ) are replaced by ( y₁ , ...... , y_n ) defined as

( y₁ , ...... , y_n ) = ( x₁ , ...... , x_n ) - C [ a₁ + sgn(a₁) ( a₁² + ....... + a_n² )^1/2 , a₂ , ...... , a_n ]

where C = [ SUM_i=1,...,n a_i x_i + x₁ sgn(a₁) ( SUM_i=1,...,n a_i² )^1/2 ] / [ SUM_i=1,...,n a_i² + a₁ sgn(a₁) ( SUM_i=1,...,n a_i² )^1/2 ]

-The first column becomes ( A₁ , 0 , .......... , 0 ) with A₁ = -sgn(a₁) ( SUM_i=1,...,n a_i² )^1/2
-Then,the same process is repeated with the 2nd column, without taking into account the 1st column and the 1st raw, and so on...

-Actually, the program listed below starts with the last column and the last element, and then it continues with the next-to-last column and so on...

Data Registers: R00 = Det A ( Registers R01 thru Rn.m are to be initialized before executing "HLS" )

• R01 = b₁ • Rn+1 = a_1,1 ..................... • Rnm-n+1 = a_1,n
• R02 = b₂ • Rn+2 = a_2,1 ..................... • Rnm-n = a_2,n

................ ...............................................................................

• Rnn = b_n • R2n = a_n,1 ..................... • Rnm = a_n,n

>>> At the end x₁ , .......... , x_n are in registers R01 to Rnn and the square matrix has become a lower triangular matrix.

• Of course, this program may be used to solve several linear systems with the same "principal" matrix,
especially if you want to invert the square-matrix itself.

Flag: F00

CF 00 = we solve the linear system(s) except if det A = 0
SF 00 = we stop when the array has only zeros above the "diagonal"

F00 is automatically set if m <= n ( overdetermined systems )

Subroutine: "TRS" ( cf "Triangular Systems for the HP-41" )

-Lines 51 & 143 are three-byte GTOs
-Line 72 is a TEXT0 or any other NOP instruction ( STO X .... )

01 LBL "HLS"
02 X<=Y?
03 SF 00
04 X<>Y
05 .1
06 %
07 +
08 STO N
09 ST* Y
10 FRC
11 -
12 STO O
13 1
14 STO 00
15 LASTX
16 %
17 RCL Z
18 INT
19 +
20 LBL 01
21 STO M
22 INT
23 RCL O
24 FRC
25 +
26 E-6
27 RCL O
28 INT
29 RCL N
30 INT
31 -
32 STO P

33 CLX
34 LBL 02
35 RCL IND Z
36 X^2
37 +
38 DSE Z
39 GTO 02
40 SQRT
41 RCL IND M
42 SIGN
43 *
44 STO Q
45 ABS
46 X<=Y?
47 CLX
48 X=0?
49 STO 00
50 X=0?
51 GTO 07
52 LBL 03
53 RCL P
54 X=0?
55 GTO 06
56 RCL O
57 INT
58 ST- P
59 RCL M
60 INT
61 ST+ P
62 RCL O
63 FRC
64 ST+ Y

65 CLX
66 LBL 04
67 RCL IND Y
68 RCL IND P
69 *
70 +
71 DSE P
72 ""
73 DSE Y
74 GTO 04
75 RCL Q
76 /
77 RCL N
78 RCL P
79 +
80 RCL M
81 +
82 INT
83 RCL O
84 INT
85 -
86 STO P
87 X<>Y
88 RCL IND Y
89 +
90 RCL IND M
91 STO Z
92 RCL Q
93 ST+ IND M
94 +
95 /
96 RCL M

97 INT
98 RCL O
99 FRC
100 +
101 X<>Y
102 SIGN
103 LBL 05
104 CLX
105 RCL IND Y
106 LASTX
107 *
108 ST- IND P
109 DSE P
110 CLX
111 DSE Y
112 GTO 05
113 R^
114 STO IND M
115 GTO 03
116 LBL 06
117 RCL O
118 FRC
119 RCL M
120 INT
121 PI
122 INT
123 10^X
124 ST* Z
125 /
126 +
127 ISG X
128 CLRGX

129 RCL Q
130 ST* 00
131 CHS
132 STO IND M
133 LBL 07
134 RCL N
135 ST- O
136 RCL O
137 RCL M
138 1
139 -
140 X<=Y?
141 CLX
142 DSE X
143 GTO 01
144 RCL M
145 RCL N
146 INT
147 /
148 INT
149 RCL N
150 INT
151 RCL 00
152 CLA
153 FC? 00
154 X=0?
155 RTN
156 X<> Z
157 XEQ "TRS"
158 RCL 00
159 END

( 244 bytes / SIZE 1+n.m )

STACK	INPUTS	OUTPUTS
Y	n	/
X	m	Det A

where n = number of raws
and m = number of columns

Example1: The instructions are similar to those of "LS3" ( cf "Linear & Non-Linear Systems for the HP-41" )

   -3 = 2 x + 3 y - 4 z              R01   R04   R07   R10
   21 = 4 x - 5 y + 7 z      =     R02   R05   R08   R11
   38 = 4 x + 2 y + 6 z            R03   R06   R09   R12

         CF 00
    3   ENTER^
    4   XEQ "HLS"   >>>>   Det = -188                       ---Execution time = 26s---

   R01 = x = 2              R04   R07   R10          -3.6778         0            0
   R02 = y = 3    and    R05   R08   R11    =      1.8181    5.0864        0                   rounded to 4D
   R03 = z = 4              R06   R09   R12          -4.3782    3.4826   -10.0499

-Before calling "TRS", we had:

R01 = -7.3556 R02 = 18.8953 R03 = -38.5079

Example2: To invert the same matrix, store

   R01   R04   R07   R10   R13   R16         1    0    0    2    3   -4
   R02   R05   R08   R11   R14   R17   =    0    1    0    4   -5   7
   R03   R06   R09   R12   R15   R18         0    0    1    4    2    6

         CF 00
    3   ENTER^
    6   XEQ "HLS"   >>>>   Det = -188                      ---Execution time = 47s---

   R01   R04   R07         11/47    13/94    -1/188
   R02   R05   R08   =    -1/47    -7/47      15/94          with small roundoff-errors as usual
   R03   R06   R09         -7/47    -2/47      11/94

Note:

-With n = 7 & m = 8, Execution time = 2mn39s

2°) A Modified Algorithm

-A similar method may be used to diagonalize a matrix:

-After using the formulae mentionned above, the last raw of the matrix is divided by A_m So, the last column is now ( 0 , ........... , 0 , 1 )
-And we continue in the same way with the other complete columns.

Data Registers: R00 = Det A ( Registers R01 thru Rn.m are to be initialized before executing "HLS" )

• R01 = b₁ • Rn+1 = a_1,1 ..................... • Rnm-n+1 = a_1,n
• R02 = b₂ • Rn+2 = a_2,1 ..................... • Rnm-n = a_2,n

................ ...............................................................................

• Rnn = b_n • R2n = a_n,1 ..................... • Rnm = a_n,n

>>> At the end x₁ , .......... , x_n are in registers R01 to Rnn and the square matrix - if it's not singular - has become a Unit matrix.

• This program may also be used to solve several linear systems with the same "principal" matrix,
especially if you want to invert the square-matrix itself.

Flags: /
Subroutines: /

-Lines 49 & 143 are three-byte GTOs

01 LBL "HLS"
02 X<>Y
03 .1
04 %
05 +
06 STO N
07 ST* Y
08 FRC
09 -
10 STO O
11 1
12 STO 00
13 LASTX
14 %
15 RCL Z
16 INT
17 +
18 LBL 01
19 STO M
20 INT
21 RCL O
22 FRC
23 +
24 E-6
25 RCL O
26 INT
27 RCL N
28 INT
29 -
30 STO P

31 CLX
32 LBL 02
33 RCL IND Z
34 X^2
35 +
36 DSE Z
37 GTO 02
38 SQRT
39 RCL IND M
40 SIGN
41 *
42 STO Q
43 ABS
44 X<=Y?
45 CLX
46 X=0?
47 STO 00
48 X=0?
49 GTO 08
50 LBL 03
51 RCL P
52 X=0?
53 GTO 06
54 RCL O
55 INT
56 ST- P
57 RCL M
58 INT
59 ST+ P
60 RCL O

61 FRC
62 ST+ Y
63 CLX
64 LBL 04
65 RCL IND Y
66 RCL IND P
67 *
68 +
69 DSE P
70 ""
71 DSE Y
72 GTO 04
73 RCL Q
74 /
75 RCL N
76 INT
77 ST+ P
78 CLX
79 RCL P
80 RCL M
81 +
82 INT
83 RCL O
84 INT
85 -
86 X<>Y
87 RCL IND Y
88 +
89 RCL IND M
90 STO Z

91 RCL Q
92 ST+ IND M
93 +
94 /
95 RCL O
96 X<>Y
97 SIGN
98 LBL 05
99 CLX
100 RCL IND Y
101 LASTX
102 *
103 ST- IND P
104 DSE P
105 CLX
106 DSE Y
107 GTO 05
108 R^
109 STO IND M
110 GTO 03
111 LBL 06
112 RCL O
113 INT
114 LASTX
115 FRC
116 PI
117 INT
118 10^X
119 ST/ Z
120 *

121 +
122 ISG X
123 CLRGX
124 RCL M
125 RCL Q
126 ST* 00
127 CHS
128 STO IND Y
129 LBL 07
130 ST/ IND Y
131 DSE Y
132 GTO 07
133 LBL 08
134 RCL N
135 ST- O
136 RCL O
137 RCL M
138 1
139 -
140 X<=Y?
141 CLX
142 DSE X
143 GTO 01
144 RCL 00
145 CLA
146 END

( 223 bytes / SIZE 1+n.m )

STACK	INPUTS	OUTPUTS
T	/	n.nnn
Z	/	/
Y	n	/
X	m	Det A

where n = number of raws
and m = number of columns

Example1:

   -3 = 2 x + 3 y - 4 z              R01   R04   R07   R10
   21 = 4 x - 5 y + 7 z      =     R02   R05   R08   R11
   38 = 4 x + 2 y + 6 z            R03   R06   R09   R12

3 ENTER^
4 XEQ "HLS" >>>> Det = -188 ---Execution time = 25s---

   R01 = x = 2              R04   R07   R10           1    0    0
   R02 = y = 3    and    R05   R08   R11    =     0    1    0
   R03 = z = 4              R06   R09   R12           0    0    1

Example2: To invert the same matrix, store

3 ENTER^
6 XEQ "HLS" >>>> Det = -188 ---Execution time = 42s---

   R01   R04   R07         11/47    13/94    -1/188
   R02   R05   R08   =    -1/47    -7/47      15/94          with small roundoff-errors...
   R03   R06   R09         -7/47    -2/47      11/94

Notes:

-With n = 7 & m = 8, Execution time = 2mn50s
-The Householder method is very sound and no pivoting is needed.
-On the other hand, it's slower and roundoff-errors may be slightly larger than when using Gaussian elimination.
-"LS3" inverses the matrix above in 25 seconds instead of 42 !

hp41programs

Householder Transformation for the HP-41