The Museum of HP Calculators

Eigenvalues & Eigenvectors for the HP-41

This program is supplied without representation or warranty of any kind. Jean-Marc Baillard and The Museum of HP Calculators therefore assume no responsibility and shall have no liability, consequential or otherwise, of any kind arising from the use of this program material or any part thereof.

Overview

1°) The Characteristic Polynomial
   a) Program#1   ( improved )
   b) Program#2   ( new )
2°) The Power Method: Dominant Eigenvalue
   a) Program#1   ( improved )
   b) Program#2   ( new )
3°) Deflation   ( improved )
4°) The Jacobi's Method
   a) Symmetric ( and some non-symmetric ) Matrices   ( improved )
   b) Symmetric Matrices only ( new )
5°) Complex Matrices ( new )

-The following programs - except the last one - only deal with real n x n matrices A

1°) "CP" calculates the characteristic polynomial of A ( n < 11 )
        "CP2" ------------------------------------------ ( n < 17 )
2°) "EGV" & "EGV2" give the dominant real eigenvalue of Aand a corresponding eigenvector by the power method ( n < 17 )
3°) "DFL" uses a deflation method to produce the real eigenvalues k_i and the eigenvectors, provided that | k₁ | > | k₂ | > ....... > | k_n | > 0   ( n < 17 )
4°) "JCB" uses the Jacobi's method to produce all the n real eigenvalues and eigenvectors of a symmetric matrix of order n ( n < 13 )
        -However, it also works if the matrix is not symmetric provided all the eigenvalues are real and distinct. ( n < 13 )
        "JCB2" runs faster and uses less data registers but it works for symmetric matrices only ( n < 15 )
5°) "JCBZ" uses a generalization of the Jacobi's method to compute all the eigenvalues and eigenvectors of a Hermitian matrix. ( n < 9 )
        -It also works for a general complex matrix provided all the eigenvalues are distinct. ( n < 9 )

-All these programs use the REGMOVE function of the X-functions module.

1°) The Characteristic Polynomial

a) Program#1

-If there exist a number k and a non-zero vector V such that A.V = k.V k is called an eigenvalue and V an eigenvector of the matrix A.
-The characteristic polynomial p(x) = xⁿ + c₁.x^n-1 + c₂.x^n-2 + ......... + c_n-2.x² + c_n-1.x + c_n is defined by p(x) = det ( x.I - A ) where I is the n x n unit matrix.
thus, its roots are actually the n eigenvalues of A.

-"CP" uses the following formulae: we define a sequence of matrices M_k by: M₀ = I and M_k+1 = M_k.A - trace(M_k.A) I /k ( the trace of a matrix is the
then we have: c_k = - trace(M_k-1.A)/k sum of its diagonal elements )

-The elements of A are to be stored into contiguous registers from R01 to Rn² , n being stored in R00
-When the program stops, R01 = c₁ ; R02 = c₂ ; ............... ; Rnn = c_n ( c₁ = -trace(A) and c_n = (-1)ⁿ det(A) )

Data Registers:

                            •   R00 = n
                            •   R01 = a₁₁    ......................    • Rn²-n+1 = a_1n
                              ...................................................................                          ( the n² elements of A )
   INPUTS
                            •   Rnn = a_n1    ......................    • Rn² = a_nn

---------

                          R00 = n
                   R01 = c₁       Rn+1 = a₁₁          Rn²+1 = a_1n
                         .............           ..........................................
OUTPUTS
                   Rnn = c_n         R2n = a_n1           Rn²+n = a_nn

Flag: F10
Subroutine: "M*" ( cf "Matrix Operations for the HP-41" )

01 LBL "CP"
02 RCL 00
03 ENTER^
04 X^2
05   E3
06 /
07 ISG X
08 +
09   E3
10 /
11 1
12 +
13 REGMOVE
14 RCL 00
15 .1
16 %
17 STO IND 00
18 ISG X
19 *
20 +
21 REGMOVE
22 RCL 00
23 .1
24 %
25 +
26 RCL 00
27 X^2
28 .1
29 %
30 ST+ X
31 +
32 +
33 ISG X
34 RCL 00
35   E5
36 /
37 +
38 CF 10
39 LBL 01
40 ENTER^
41 ISG IND 00
42   E-5
43 +
44 STO M
45 X<>Y
46 0
47 LBL 02
48 RCL IND Z
49 -
50 ISG Z
51 GTO 02
52 RCL IND 00
53 INT
54 ST/ Y
55 RDN
56 STO IND L
57 ISG L
58 FS? 30
59 GTO 04
60 LBL 03
61 ST+ IND M
62 ISG M
63 GTO 03
64 X<>Y
65 STO Y
66 RCL 00
67 X^2
68 .1
69 %
70 +
71 FS? 10
72 ST+ X
73 -
74 ENTER^
75 INT
76 RCL 00
77 X^2
78 FC? 10
79 ST+ X
80 +
81 XEQ "M*"
82 FC?C 10
83 SF 10
84 GTO 01
85 LBL 04
86 CF 10
87 CLA
88 END

( 138 bytes / SIZE 3n²+n+1 )

STACK INPUT OUTPUT

X / c_n

                                                                                                   1 2 4
Example:   Find the characteristic polynomial of the matrix A = 4 3 5
                                                                                                   7 4 7
                                                                                                                       R01 R04 R07         1 2 4
-First we store these 9 elements in registers R01 thru R09                                R02 R05 R08   =    4 3 5    respectively and n = 3 STO 00
                                                                                                                       R03 R06 R09         7 4 7

-Then XEQ "CP" >>>> c₃ = 5 = R03 ( in 31 seconds )
and R01 = -11 , R02 = -25 , R03 = 5 Thus p(x) = x³ - 11.x² - 25.x + 5

Notes:

-A is saved in registers Rn+1 thru Rn²+n and n is saved in R00
-The elements of A may be stored either in column order or in row order since A and ^tA ( the transpose of A ) have the same characteristic polynomial.

-Execution time = 1mn25s for n = 4
--------------- = 3mn18s for n = 5

-Any root-finding program yields the 3 eigenvalues k₁ = 12.90692994 , k₂ = -2.092097589 , k₃ = 0.185167648
-Then, subtracting ( for instance ) k₁ to the diagonal elements of A leads to a singular system with an infinite set of solutions ( the eigenvectors )
-One solution is ( 0.451320606 ; 0.686921424 ; 1 ) which is an eigenvector corresponding to the first eigenvalue.

-This method could be used to obtain all the eigenvalues and eigenvectors of a n x n matrix A.

b) Program#2 ( HP-41CX )

-"CP" uses 3 blocks of n² registers to store the matrices A , M_k and the products M_k.A so that n cannot exceed 10.
-"CP2" creates a temporary data file to store the matrix A and uses another method to perform the multiplicationof 2 square-matrices:
-Thus, it can produce the characteristic polynomial of a 16x16 matrix!

Data Registers:

                            •   R00 = n
                            •   R01 = a₁₁    ......................    • Rn²-n+1 = a_1n
                              ..................................................................                ( the n² elements of A )
   INPUTS
                            •   Rnn = a_n1    ......................    • Rn² = a_nn

---------

R00 = n
OUTPUTS R01 = c₁ ...................... Rnn = c_n

Flags: /
Subroutine: "TRN2" Transpose of a square matrix ( cf "Matrix Operations for the HP-41" )

01 LBL "CP2"
02 RCL 00
03 ENTER^
04 X^2
05 "T"
06 CRFLD
07   E3
08 /
09 ISG X
10 SAVERX
11 +
12   E3
13 /
14 1
15 +
16 REGMOVE
17 RCL 00
18 .1
19 %
20 STO IND 00
21 1.001
22 +
23 *
24 ISG X
25 RCL 00
26   E5
27 /
28 +
29 STO M
30 XEQ "TRN2"
31 LBL 01
32 ISG IND 00
33 RCL M
34   E-5
35 +
36 ENTER^
37 ENTER^
38 CLX
39 LBL 02
40 RCL IND Y
41 -
42 ISG Y
43 GTO 02
44 RCL IND 00
45 INT
46 ST/ Y
47 X<>Y
48 STO IND L
49 ISG L
50 FS? 30
51 GTO 07
52 LBL 03
53 ST+ IND T
54 ISG T
55 GTO 03
56 R^
57 INT
58 STO Y
59 RCL 00
60 ST- Z
61 SIGN
62 -
63   E3
64 /
65 +
66 STO N
67 RCL 00
68 .1
69 %
70 ST+ X
71 +
72 ISG X
73 STO O
74 LBL 04
75 CLX
76 SEEKPT
77 LBL 05
78 RCL O
79 0
80 LBL 06
81 RCL IND Y
82 GETX
83 *
84 +
85 ISG Y
86 GTO 06
87 STO IND N
88 ISG N
89 GTO 05
90 RCL O
91 INT
92 RCL 00
93 ST- N
94 .1
95 %
96 +
97 ST+ O
98 X<> L
99 +
100   E3
101 /
102 RCL N
103 INT
104 +
105 REGMOVE
106 RCL N
107 RCL O
108 X#Y?
109 GTO 04
110 GTO 01                 ( a three-byte GTO )
111 LBL 07
112 "T"                          do not key in lines 112-113 if you want to save the matrix A in this data file.
113 PURFL
114 CLA
115 END

( 192 bytes / SIZE n²+2n+1 )

STACK INPUT OUTPUT

X / c_n

                                                                                                     1 2 4
Example1:   Find the characteristic polynomial of the matrix A = 4 3 5
                                                                                                     7 4 7

                                                                                                                R01 R04 R07         1 2 4
-We store these 9 elements in registers R01 thru R09                                R02 R05 R08   =    4 3 5    respectively and n = 3 STO 00
                                                                                                                R03 R06 R09         7 4 7

-Then XEQ "CP2" >>>> c₃ = 5 = R03 ( in 37 seconds )
and R01 = -11 , R02 = -25 , R03 = 5 whence p(x) = x³ - 11.x² - 25.x + 5

Example2: The 10x10 matrix A = [ a_ij ] is defined by a_ij = i^j mod 13

-Store these 100 elements into registers R01 thru R100 ( you can use for instance the "MATRIX" program listed in "Creating a Matrix for the HP-41"
and the subroutine LBL "U" X<>Y Y^X 13 MOD RTN then "U" ASTO 00 1.10010 XEQ "MATRIX" )

-Then 10 STO 00 XEQ "CP2" and you should find ( in about 54mn45s ):

p(x) = x¹⁰ - 43 x⁹ - 968 x⁸ - 2462 x⁷ + 40796 x⁶ - 488852 x⁵ - 10916340 x⁴ + 15630136 x³ + 441980832 x² - 1282786560 x + 155105280

-There is no roundoff error: all the coefficients are exact.
-Obviously, this is not a very fast routine for large n-values!
-For n = 16 , the execution time is probably greater than 5 hours, so a good emulator like V41 in "turbo mode" is quite useful...

-However, if you only need the real dominant eigenvalue , the following program may be interesting.

2°) The Power method: Dominant Eigenvalue

a) Program#1

-Here, we assume that the matrix A has n eigenvalues   k₁ ; k₂ ; ....... ; k_n corresponding to n independent eigenvectors V₁ ; V₂ ; ....... ; V_n
and that k₁ is the dominant eigenvalue | k₁ | > | k_i |    i = 2 ; 3 ; ....... ; n
-Any vector V can be written    V = a₁ V₁ + a₂ V₂ + ....... + a_n V_n whence AV = a₁ k₁ V₁ + a₂ k₂ V₂ + ....... + a_n k_n V_n
and A^pV = k₁^p ( a₁ V₁ + a₂ ( k₂/k₁ )^p V₂ + ....... + a_n ( k_n/k₁ )^p V_n ) which tends to k₁^p a₁ V₁ as p tends to infinity.

-Furthermore, if we define W_p+1 = ( A.W_p ) / || A.W_p || with W₀ = V ( almost ) arbitrary , then W_p tends to a unit eigenvector U
and ^tW_p AW_p tends to the corresponding eigenvalue k₁

Data Registers:

                            •   R00 = n
                            •   R01 = v₁           •   Rn+1 = a₁₁        • Rn²+1 = a_1n
                              .................                  ( the n² elements of A )                       where V ( v₁ ; ....... ; v_n )   is an initial guess of the eigenvector.
   INPUTS
                            •   Rnn = v_n          •   R2n = a_n1           • Rn²+n = a_nn

---------

                       R00 = n
                R01 = u₁       Rn+1 = a₁₁          Rn²+1 = a_1n
                       .............           .............................................      and k in X-register.   U ( u₁; ........ ; u_n) = a unit eigenvector
OUTPUTS
                Rnn = u_n         R2n = a_n1           Rn²+n = a_nn

Flags: /
Subroutine: "M*" ( cf "Matrix Operations for the HP-41" )

01 LBL "EGV"
02 LBL 01
03 RCL 00
04 X^2
05   E3
06 /
07 RCL 00
08 ST+ Y
09 1
10 %
11 +
12   E3
13 /
14 1
15 +
16 ST+ Y
17 LASTX
18 RCL 00
19 ST+ Y
20 X^2
21 +
22 XEQ "M*"
23 LASTX
24 RCL 00
25   E3
26 /
27 ISG X
28 LASTX
29 /
30 +
31 REGSWAP
32 INT
33 RCL 00
34 STO N
35 +
36 DSE X
37 STO O
38 0
39 LBL 02
40 RCL IND Y
41 RCL IND N
42 X^2
43 ST+ M
44 X<> L
45 *
46 +
47 DSE Y
48 DSE N
49 GTO 02
50 VIEW X
51 STO N
52 RCL 00
53 RCL 00
54 RCL M
55 SQRT
56 R^
57 SIGN
58 *
59 LBL 03
60 ST/ IND Y
61 DSE Y
62 GTO 03
63 RDN
64 LBL 04
65 RCL IND Y
66 RCL IND O
67 -
68 ABS
69 +
70 DSE O
71 DSE Y
72 GTO 04
73   E-8                    ( or another "small" number like E-7 or 3 E-9   RCL 00 *   in order to take the size of the matrix into account )
74 X<Y?
75 GTO 01
76 RCL N
77 CLA
78 END

( 122 bytes / SIZE n²+2n+1 )

STACK INPUTS OUTPUTS

X / k

                                                                                                                                                              1 2 4
Example:   Find the dominant eigenvalue k and the corresponding unit eigenvector U of the matrix A = 4 3 5
                                                                                                                                                              7 4 7
-Here, we have n = 3 therefore     3 STO 00
and if we choose ( 1 ; 1 ; 1 ) as initial guess, we store the 12 numbers

   1 1 2 4                              R01 R04 R07 R10
   1 4 3 5     into registers      R02 R05 R08 R11     respectively.
   1 7 4 7                              R03 R06 R09 R12

-Then XEQ "EGV" the HP-41 displays the successive k-approximations
and 104 seconds later, it yields k = 12.90692995 in X-register and U is in registers R01 R02 R03 : U ( 0.348663346 ; 0.530674468 ; 0.772540278 )

Notes:

-The greater | k₁/ k₂ | , the faster the convergence.
-If you calculate A^-1 ( provided A is non singular ) and apply "EGV" again , you'll find the reciprocal of the smallest eigenvalue and the corresponding eigenvector.
in the above example, 1/k₃ = 5.400511417 and the eigenvector ( -0.094824736 ; 0.897989404 ; -0.429678136 ) which is also an eigenvector of A.
-Subtracting a number N to the diagonal elements of A can make the other extreme eigenvalue dominant:
for instance choosing N = 13 and applying "EGV" leads to k₂ - 13 = -15.092... ( add 13 to the result ) and the last eigenvector is correctly found.
-Unfortunately in our example, the convergence is very slow since the new ratio | k'₁/ k'₂ | = 1.17... is very close to 1. More than 100 iterations are needed
to achieve a good accuracy! All the digits of k are correct , but the components of U are only obtained to 7 decimals.

-Thus, in this case, the convergence of k is faster ( I mean less slow ) than the convergence of U:
That's why the test ( line 74 ) compares 2 successive approximations of U instead of 2 successive approximations of k.
Otherwise, the program could have been both shorter and faster: see below!

b) Program#2

-The following program uses the same "power method" but the iteration stops when 2 successive k-values are (approximately ) equal.
-The eigenvector V( x₁ , ....... , x_n ) is choosen so that its first coordinate = the eigenvalue k.
-So, it will not work if x₁ = 0

01 LBL "EGV2"
02 LBL 01
03 RCL 00
04 X^2
05   E3
06 /
07 RCL 00
08 ST+ Y
09 1
10 %
11 +
12   E3
13 /
14 1
15 +
16 ST+ Y
17 LASTX
18 RCL 00
19 ST+ Y
20 X^2
21 +
22 XEQ "M*"
23 LASTX
24 RCL 00
25   E3
26 /
27 ISG X
28 LASTX
29 /
30 +
31 REGSWAP
32 RCL 00
33 RCL IND Y
34 LBL 02
35 ST/ IND Y
36 DSE Y
37 GTO 02
38 RCL 01
39 VIEW X
40 ST- Y
41 /
42 ABS
43 2 E-9
44 X<Y?
45 GTO 01
46 RCL 01
47 END

( 77 bytes / SIZE n²+2n+1 )

STACK INPUTS OUTPUTS

X / k

-With the same example as "EGV" it yields k = 12.90692994 in X-register in 90 seconds

and V is in registers R01R02R03 : V ( 12.90692994 ; 19.64467513 ; 28.59814010 )

-With this less reliable method, always check that the eigenvector V is also in registers Rn²+n+1 ...... Rn²+2n ( in this example, R13 R14 R15 )
-Of course, the last decimals may differ slightly.

3°) Deflation

-If k₁ is an eigenvalue , U₁ the corresponding unit eigenvector of A and U₁' the corresponding unit eigenvector of ^tA
then A and the new matrix A' = A - k₁ ( U₁ ^tU₁' ) / ( U₁.U₁' ) have the same eigenvalues and eigenvectors , except that k₁ has been replaced by 0.
-Thus, the 2^nd dominant eigenvalue of A is now the 1^st dominant eigenvalue of A'

-"DFL" uses this method to calculate all the eigenvalues and eigenvectors , provided they are real and verify | k₁ | > | k₂ | > ....... > | k_n | > 0

1°) "DFL" calculates k₁ and U₁ ( with flag F02 clear ) and the program stops ( line 87 )

>>>> k₁ is in X-register
>>>> U₁ is in registers R01 thru Rnn

2°) If you press R/S , the HP-41 sets flag F02 and computes k₁ once again & U₁'
the matrix A is replaced by A' in registers Rn+1 thru Rn²+n and when the program stops again:

>>>> k₂ is in X-register
>>>> U₂ is in registers R01 thru Rnn

... etc ...

Data Registers:

---------

R00 = n
R01 = u₁ Rn+1 = a'₁₁ Rn²+1 = a'_1n

............. A will be replaced with k in X-register. U ( u₁; ........ ; u_n ) = a unit eigenvector
OUTPUTS A' = A - k ( U.^tU' ) / ( U.U' ) when F02 is set.

Rnn = u_n R2n = a'_n1 Rn²+n = a'_nn

Flag: F02
Subroutine: "M*" ( cf "Matrix Operations for the HP-41" )

01 LBL "DFL"
02 LBL 12
03 CF 02
04 LBL 01
05 RCL 00
06 X^2
07 LASTX
08 ST+ Y
09   E3
10 ST/ Z
11 /
12 1
13 ST+ Z
14 %
15 ISG Y
16 ST+ Z
17 FS? 02
18 CLX
19 +
20 ISG Y
21 FS? 02
22 X<>Y
23 RCL 00
24 ENTER^
25 X^2
26 +
27 1
28 +
29 XEQ "M*"
30 LASTX
31 RCL 00
32   E3
33 /
34 ISG X
35 LASTX
36 /
37 +
38 REGSWAP
39 INT
40 RCL 00
41 STO N
42 +
43 DSE X
44 STO O
45 0
46 LBL 02
47 RCL IND Y
48 RCL IND N
49 X^2
50 ST+ M
51 X<> L
52 *
53 +
54 DSE Y
55 DSE N
56 GTO 02
57 VIEW X
58 STO N
59 RCL 00
60 RCL 00
61 RCL M
62 SQRT
63 R^
64 SIGN
65 *
66 LBL 03
67 ST/ IND Y
68 DSE Y
69 GTO 03
70 RDN
71 LBL 04
72 RCL IND Y
73 RCL IND O
74 -
75 ABS
76 +
77 DSE O
78 DSE Y
79 GTO 04
80   E-8                    ( or another "small" number like E-7 or something like 3 E-9   RCL 00 *   in order to take the size of the matrix into account )
81 X<Y?
82 GTO 01              ( a three-byte GTO )
83 RCL N
84 FS? 02
85 GTO 05
86 CLA
87 STOP                 ( or RTN )
88 SF 02
89 RCL 00
90 1
91 +
92 X^2
93   E3
94 /
95 ISG X
96 RCL 00
97 LASTX
98 X^2
99 /
100 +
101 REGMOVE
102 GTO 01              ( a three-byte GTO )
103 LBL 05
104 3
105 RCL 00
106 ST+ Y
107 ST* Y
108 STO M
109 CLX
110 LBL 06
111 RCL IND Y
112 RCL IND M
113 *
114 +
115 DSE Y
116 DSE M
117 GTO 06
118 ST/ N
119 CLX
120 .1
121 %
122 X<>Y
123 ST+ Y
124 RCL 00
125 ST- Y
126 ST+ Z
127 STO M
128 LBL 07
129 CLX
130 RCL IND Z
131 RCL IND M
132 *
133 RCL N
134 *
135 ST- IND Y
136 DSE Y
137 DSE Z
138 GTO 07
139 CLX
140 RCL 00
141 ST+ Z
142 DSE M
143 GTO 07
144 GTO 12              ( a three-byte GTO )
145 END

( 229 bytes / SIZE n²+3n+1 )

STACK INPUTS OUTPUTS

X / k₁ , k₂, ...

                                                                                                                                                        1 2 4
Example:   Find all the eigenvalues and the corresponding unit eigenvectors of the same matrix A = 4 3 5
                                                                                                                                                        7 4 7
1°) Once again,    n = 3 therefore     3 STO 00

and if we choose ( 1 ; 1 ; 1 ) as initial guess, we store the 12 numbers

   1 1 2 4                              R01 R04 R07 R10
   1 4 3 5     into registers      R02 R05 R08 R11     respectively.
   1 7 4 7                              R03 R06 R09 R12

- XEQ "DFL" the HP-41 displays the successive k₁-approximations ( with F02 clear )

- Eventually, k₁ = 12.90692995 is in X-register and U₁ is in registers R01 R02 R03 : U₁ ( 0.348663346 ; 0.530674468 ; 0.772540278 )

2°) To obtain the second eigenvalue simply press R/S ( you may re-initialize R01 to Rnn but in most cases, it's not necessary )

the HP-41 sets flag F02 and displays the successive k₁-approximations converging to 12.90692995 again
then flag F02 is cleared and the HP-41 displays the successive k₂-approximations and

k₂ = -2.092097594 in X-register and U₂ in R01R02R03 : U₂ ( 0.800454175 ; -0.041651079 ; -0.597945065 )

3°) Once again R/S ( similar displays ) and finally,

k₃ = 0.185167649 in X-register and U₃ in R01R02R03 : U₃ ( -0.094824730 ; 0.897989404 ; -0.429678136 )

Warning:   Watch the successive k-numbers displayed by the HP-41 ( with flag F02 set and clear ): the 2 limits must be ( approximately ) the same.
                   Otherwise change the initial guess and start over again.
        ( It may happen that an initial vector V leads to a dominant eigenvalue of A but a non-dominant eigenvalue of ^tA.
          In such a case, A' will be miscalculated and the next eigenvalues will be wrong! )

-If k₁ and k₂ are real but k₃ ... are complex, the first 2 eigenvalues and eigenvectors will be correctly computed.

4°) The Jacobi's method

a) Symmetric ( and some non-symmetric ) Matrices

*** A non-symmetric real matrix may have complex eigenvalues, but all the eigenvalues of a real symmetric matrix A are necessarily real !

-In the Jacobi's method, the eigenvalues are the elements of a diagonal matrix P equal to the infinite product ...O_k^-1.O_k-1^-1....O₂^-1.O₁^-1.A.O₁.O₂....O_k-1.O_k....
where the O_k are rotation matrices. The eigenvectors are the columns of O₁.O₂....O_k-1.O_k....

-"JCB" finds the greatest off-diagonal element a_{i j} ( lines 22 to 46 )
-Then, O₁ is determined so that it makes a pair of off-diagonal elements zero in O₁^-1.A.O₁
-It happens if the rotation angle µ is choosen so that Tan 2µ = 2.a_{i j} / ( a_{i i} - a_{j j} )
-The process is repeated until the greatest off-diagonal element is smaller than E-8 in absolute value ( line 48 )

-The successive greatest | a_{i j} | ( i > j ) are displayed ( line 47 ) when the routine is running.

*** Though it's not really orthodox, we can try to apply the same program to non-symmetric matrices: of course, it cannot work if some eigenvalues are complex.

But if all the eigenvalues are real, the infinite product T = [ t_{i j} ] = ...O_k^-1.O_k-1^-1....O₂^-1.O₁^-1.A.O₁.O₂....O_k-1.O_k.... is an upper triangular matrix,
the diagonal elements of T ( i-e t_{i i} = k_i ) are the n eigenvalues and the first column of U = O₁.O₂....O_k-1.O_k.... is an eigenvector corresponding to k₁

-If all the eigenvalues are distinct, we can create an upper triangular matrix W = [ w_{i j} ] defined by

   w_{i j} = 0    if i > j
   w_{i i} = 1
   w_{i j} = Sum_{k = i+1, i+2 , .... , j}      t_{i k} w_{k j} / ( t_{j j} - t_{i i} )      if i < j

-Then, the n columns of U.W are the eigenvectors corresponding to the n eigenvalues t₁₁ = k₁ , ............... , t_nn = k_n

*** When the following program stops: the n eigenvalues are stored in registers R01 thru Rnn
                                                             the first eigenvector is stored ---------- Rn+1 thru R2n
                                                             ..................................................................................
                                                             the n^th eigenvector is stored ---------- Rn²+1 thru Rn²+n as shown below.

Data Registers:
                            •   R00 = n
                            •   R01 = a₁₁ .................      • Rn²-n+1 = a_1n
                              ..................................................................              ( the n² elements of A in column order )
   INPUTS
                            •   Rnn = a_n1 .................      • Rn² = a_nn

---------

   DURING            R01 thru Rn² = the "infinite" product   ...O_k^-1.O_k-1^-1....O₂^-1.O₁^-1.A.O₁.O₂....O_k-1.O_k....
       THE                Rn²+1 thru R2n² = O₁.O₂........O_k.....
ITERATIONS      R2n²+1 thru R2n²+n: temp ( for non-symmetric matrices only )

---------

                             R00 = n
                             R01 = k₁       Rn+1 = V₁⁽¹⁾       R2n+1 = V₁⁽²⁾ ....................    Rn²+1 = V₁⁽ⁿ⁾
                             .............
OUTPUTS
                             Rnn = k_n        R2n = V_n⁽¹⁾         R3n = V_n⁽²⁾    .....................   Rn²+n = V_n⁽ⁿ⁾

( n eigenvalues ; 1st eigenvector ; 2nd eigenvector ; .................. ; n^th eigenvector )

Flag: F02 Set flag F02 for a symmetric matrix
Clear flag F02 for a non-symmetric matrix
Subroutines: /

01 LBL "JCB"
02 RCL 00
03 X^2
04 .1
05 %
06 ST+ X
07 +
08 ISG X
09 CLRGX                  if you don't have an HP-41CX, replace line 09 by   ENTER^   SIGN CLX LBL 12 STO IND L ISG L GTO 12 CLX
10 RCL 00
11 1
12 +
13   E5
14 /
15 +
16 SIGN
17 LBL 00
18 STO IND L
19 ISG L
20 GTO 00
21 LBL 01
22 RCL 00
23   E3
24 /
25 STO M
26 2
27 +
28 ENTER^
29 ENTER^
30 CLX
31 LBL 02
32 RCL IND Z
33 ABS
34 X>Y?
35 STO Z
36 X>Y?
37 +
38 RDN
39 ISG Z
40 GTO 02
41 ISG M
42 RCL M
43 ST+ T
44 RDN
45 ISG Z
46 GTO 02
47 VIEW X              the HP-41 displays the successive greatest sub-diagonal elements | a_ij |   ( i > j )
48   E-8                    ( or another "small" number )
49 X>Y?
50 GTO 05              ( a three-byte GTO )
51 RDN
52 SIGN
53 -
54 INT
55 ENTER^
56 STO M
57 RCL 00
58 ST/ Z
59 MOD
60 ENTER^
61 STO N
62 LASTX
63 *
64 +
65 X<>Y
66 INT
67 RCL X
68 RCL 00
69 *
70 STO O
71 +
72 1
73 ST+ M
74 ST+ N
75 ST+ Z
76 +
77 RCL IND M
78 ST+ X
79 RCL IND Y
80 RCL IND T
81 -
82 R-P
83 CLX
84 2
85 /
86 1
87 P-R
88 X<> O
89 X<>Y
90 X<> N
91 RCL 00
92 .1
93 %
94 +
95 *
96 1
97 ST+ Z
98 +
99 RCL 00
100 ST- Y
101 X^2
102 ST+ Z
103 .1
104 %
105 +
106 +
107 XEQ 03
108 ST- Y
109 ST- Z
110 X^2
111 ST- Z
112 .1
113 %
114 +
115 -
116 XEQ 03
117 DSE Z
118 ST/ Z
119 /
120 INT
121 RCL 00
122 1
123 %
124 X<>Y
125 X^2
126 +
127   E3
128 /
129 ST+ Z
130 +
131 XEQ 03
132 GTO 01              ( a three-byte GTO )
133 LBL 03
134 STO M
135 LBL 04
136 CLX
137 RCL IND M
138 RCL O
139 *
140 RCL IND Y
141 RCL N
142 *
143 -
144 X<> IND M
145 RCL N
146 *
147 X<> IND Y
148 RCL O
149 *
150 ST+ IND Y
151 ISG Y
152 CLX
153 ISG M
154 GTO 04
155 X<> M
156 RCL 00
157 RTN
158 LBL 05
159 CLD
160 RCL 00
161 X^2
162 LASTX
163 1
164 +
165   E5
166 /
167 +
168 STO O
169 RCL 00
170 LBL 06
171 RCL IND Y
172 STO IND Y
173 DSE Y
174 RDN
175 DSE Y
176 GTO 06
177 FS? 02                   Lines 177 to 269 are unuseful if the matrix is symmetric
178 GTO 11                 ( a three-byte GTO )
179 RCL 00
180 DSE X
181   E3
182 /
183 ISG X
184 STO M
185 LBL 07
186 RCL 00
187 ENTER^
188 X^2
189 ST+ X
190 +
191 STO Y
192 RCL M
193 INT
194 -
195   E3
196 /
197 +
198 STO N
199 RCL O
200 RCL M
201 INT
202 -
203   E-5
204 -
205 1
206 STO IND Z
207 CLX
208 LBL 08
209 RCL IND Y
210 RCL IND N
211 *
212 +
213 DSE Y
214 DSE N
215 GTO 08
216 RCL IND O
217 RCL IND Z
218 -
219 /
220 STO IND N
221 ISG M
222 GTO 07
223 RCL 00
224 DSE M
225 LBL 09
226 RCL 00
227 X^2
228 LASTX
229 RCL M
230 INT
231 *
232 +
233   E3
234 /
235 RCL 00
236 X^2
237 +
238 +
239 RCL 00
240   E5
241 /
242 +
243 STO P                  ( synthetic )
244 CLX
245 RCL N
246 STO Q
247 CLX
248 LBL 10
249 RCL IND P
250 RCL IND Q
251 *
252 +
253 ISG Q
254 CLX
255 ISG P
256 GTO 10
257 ST+ IND P
258 X<>Y
259 DSE X
260 GTO 09
261 RCL M
262 FRC
263   E-3
264 -
265 STO M
266 DSE O
267 ISG M
268 GTO 07                ( a three-byte GTO )
269 LBL 11
270 RCL 00
271 X^2
272 LASTX
273 1
274 ST+ Z
275 +
276   E3
277 /
278 +
279 RCL 00
280 X^2
281   E6
282 /
283 +
284 REGMOVE
285 RCL 01
286 CLA
287 END

( 429 bytes / SIZE 2n²+1 if A is symmetric / SIZE 2n²+n+1 if A is not symmetric )

STACK INPUTS OUTPUTS

X / k₁

1 2 4
Example1: Let's compute all the eigenvalues and the eigenvectors of the matrix A = 2 7 3
4 3 9

n = 3 therefore 3 STO 00 and we store the 9 numbers

     1 2 4                              R01 R04 R07
     2 7 3     into registers      R02 R05 R08     respectively.
     4 3 9                              R03 R06 R09

SF 02 ( A is symmetric )

   XEQ "JCB" yields ( in 1m33s )                               k₁ = 12.81993499 in R01
                                                                                  k₂ = 4.910741214 in R02
                                                                                  k₃ = -0.730676199 in R03

and the 3 corresponding eigenvectors: V₁ ( 0.351369026 ; 0.521535689 ; 0.777521917 )   in   ( R04 ; R05 ; R06 )
                                                           V₂ ( -0.101146468 ; 0.846760701 ; -0.522269766 ) in   ( R07 ; R08 ; R09 )
                                                           V₃ ( -0.930757326 ; 0.104865823 ; 0.350276976 )   in   ( R10 ; R11 ; R12 )

                                  1 2 4
Example2:        A = 4 3 5
                                  7 4 7

n = 3 STO 00 and we store the 9 numbers

     1 2 4                              R01 R04 R07
     4 3 5     into registers      R02 R05 R08     respectively.
     7 4 7                              R03 R06 R09

CF 02 ( A is not symmetric )

XEQ "JCB" >>>> k₁ = 12.90692994 ( in 2mn44s )

and we have          k₁ = 12.90692994 in R01
                               k₂ = 0.185167649 in R02
                               k₃ = -2.092097593 in R03

and the 3 corresponding eigenvectors: V₁ ( 0.348663347 ; 0.530674467 ; 0.772540277 )   in   ( R04 ; R05 ; R06 )
                                                           V₂ ( -0.095420104 ; 0.903627529 ; -0.432375917 ) in   ( R07 ; R08 ; R09 )
                                                           V₃ ( -0.830063031 ; 0.043191757 ; 0.620063097 )   in   ( R10 ; R11 ; R12 )

-The eigenvectors are not unit vectors except the first one.
-Divide V₂ & V₃ by || V₂ || & || V₃ || respectively if you want unit eigenvectors.

Notes:

-If A is non-symmetric, this program only works if all the eigenvalues are real and distinct.
-If all the eigenvalues are real but not distinct, they will be correctly computed but not the eigenvectors ( DATA ERROR line 219 )

-The former version of this program used the subroutine "M*" to multiply matrices by the successive "rotation" matrices, but it was prohibitive!
-This new "JCB" uses a much simpler ( and faster ) way to perform these products: in fact, at each step, only 2 columns and 2 rows are modified in O_i^-1.M.O_i
and only 2 columns in M.O_i .

-Execution time = 13 seconds/iteration if n = 3
- -------------- = 16 ---------------- if n = 4
- -------------- = 27 ---------------- if n = 7

-Unfortunately, the number of required iterations may increase very much with n, especially if A is not symmetric.

Improvement?

-Actually, the Jacobi's rotations do not zero the element in position ( i , j ) in O_k^-1.A.O_k if the matrix is non-symmetric!
-We could use the formula of paragraph 5 below. Unfortunately, it involves complex arithmetic if the radicand is negative.
-But we can replace this radicand by 0 in this case!

-If you want to use this hybrid method, replace lines 77 thru 87 by

    STO P         RCL M                ST+ X           RCL IND P       X^2           SQRT        SIGN
    X<>Y          -                          STO Z           RCL IND Q      +                +                P-R
    STO Q        RCL IND X         *                    -                        X<0?         R-P
    +                 RCL IND M        ST+ X           STO Z               CLX          CLX                            ( 302 lines / 451 bytes )

-Numerical results are mixed: for the matrix of example 2 above, the problem is solved in 2mn14s instead of 2mn44s
but in other examples, the number of iterations may be greater and execution time longer!
-So I let you decide whether it's worthwhile or not!

                                                                                                                     1 2 4 7
Exercise:      Find all the eigenvalues and eigenvectors of the 4x4 matrix        2 3 7 1
                                                                                                                     4 7 2 4
                                                                                                                     7 1 4 9

Answer: k₁ = 16.97583168 , k₂ = 6.365547530 , k₃ = -3.301311094 , k₄ = -5.040068160

       V₁ ( 0.455772321 ; 0.346041152 ; 0.464075961 ; 0.676136537 )
       V₂ ( -0.142731960 ; 0.681492880 ; 0.448494335 ; -0.560399745 )
       V₃ ( 0.842568185 ; -0.247658954 ; 0.050153515 ; -0.475634862 )
       V₄ ( -0.248953877 ; -0.595388965 ; 0.762214511 ; -0.050625961 )

-These results are obtained in 4mn22s.

b) Symmetric Matrices only

-If the matrix is symmetric, we can save data registers: it is unuseful to store a_{i j} and a_{j i} ( i # j ) in 2 registers since they are equal!
-So the program is approximately 20% faster than "JCB" and it can be used with a matrix of order 14

Data Registers: • R00 = n ( All these registers are to be initialized before executing "JCB2" )

                          • R01 = a₁₁    • Rn+1 = a_1,2   • R2n    = a_2,3   .................................. • R(n²+n)/2 = a_n-1,n
                          • R02 = a₂₂    • Rn+2 = a_1,3   • R2n+1 = a_2,4     ...................
   INPUTS          ...............         .................      ...................
                           ...............         .................    • R3n-2 = a_2,n
                           ................      • R2n-1 = a_1,n
                          • Rnn = a_nn

-----------------------------------------------------------------------------------------------------------------------------------------------------------------

When the program stops:

( n eigenvalues ; 1st eigenvector ; 2nd eigenvector ; .................. ; n^th eigenvector )

Flags: /
Subroutines: /

01 LBL "JCB2"
02 RCL 00
03 X^2
04 STO Y
05 3
06 *
07 RCL 00
08 ST+ Z
09 +
10 2
11 ST/ Z
12   E3
13 *
14 /
15 +
16 ISG X
17 CLRGX                  if you don't have an HP-41CX, replace line 17 by   ENTER^   SIGN CLX LBL 12 STO IND L ISG L GTO 12 CLX
18 RCL 00
19 1
20 +
21   E5
22 /
23 +
24 SIGN
25 LBL 08
26 STO IND L
27 ISG L
28 GTO 08
29 LBL 01
30 RCL 00
31 ENTER^
32 X^2
33 +
34 2
35 /
36 RCL 00
37   E3
38 /
39 +
40 ENTER^
41 ENTER^
42 CLX
43 LBL 02
44 RCL IND Z
45 ABS
46 X>Y?
47 STO Z
48 X>Y?
49 +
50 RDN
51 DSE Z
52 GTO 02
53 VIEW X
54   E-8                        or another "small" number
55 X>Y?
56 GTO 10                 ( a three-byte GTO )
57 CLX
58 ENTER^
59 R^
60 INT
61 STO M
62 RCL 00
63 -
64 X<>Y
65 LBL 03
66 ISG Z
67 CLX
68 RCL 00
69 +
70 R^
71 -
72 X<Y?
73 GTO 03
74 -
75 RCL 00
76 +
77 RCL IND Y
78 RCL IND Y
79 -
80 RCL IND M
81 ST+ X
82 X<>Y
83 R-P
84 CLX
85 2
86 /
87 1
88 P-R
89 STO N
90 RDN
91 STO O
92 X^2
93 RCL IND Y
94 *
95 STO P                    ( synthetic )
96 CLX
97 RCL IND Z
98 RCL N
99 X^2
100 *
101 ST+ P
102 CLX
103 X<> IND M
104 ST+ X
105 RCL N
106 *
107 RCL O
108 *
109 ENTER^
110 X<> P
111 +
112 X<> IND Z
113 RCL O
114 X^2
115 *
116 RCL P
117 -
118 X<> IND Y
119 RCL N
120 X^2
121 *
122 ST+ IND Y
123 CLX
124 RCL 00
125 ENTER^
126 X^2
127 +
128 2
129 /
130 X<>Y
131 STO Q
132 R^
133 STO P                    ( synthetic )
134 ENTER^
135 SIGN
136 ST- Z
137 -
138 RCL 00
139 ST* Z
140 *
141 ST+ T
142 RDN
143 +
144 RCL X
145 RCL 00
146 +
147 PI
148 INT
149 10^X
150 /
151 +
152 ISG X
153 ISG Y
154 LBL 09
155 XEQ 07
156 ISG Y
157 CLX
158 ISG X
159 GTO 09
160 RCL 00
161 DSE X
162 X<> P
163 RCL Q
164 LBL 04
165 RCL P
166 ST+ Z
167 +
168 RCL M
169 X=Y?
170 GTO 00
171 RDN
172 XEQ 07
173 DSE P
174 GTO 04
175 LBL 00
176 RDN
177 LBL 05
178 RCL P
179 ENTER^
180 SIGN
181 ST+ T
182 -
183 +
184 X<>Y
185 RCL M
186 X=Y?
187 GTO 00
188 X<> Z
189 XEQ 07
190 DSE P
191 GTO 05
192 LBL 00
193 X<> Z
194 LBL 06
195 DSE P
196 X=0?
197 GTO 01                 ( a three-byte GTO )
198 ISG X
199 CLX
200 ISG Y
201 CLX
202 XEQ 07
203 GTO 06
204 LBL 07
205 RCL IND Y
206 RCL O
207 *
208 SIGN
209 CLX
210 RCL IND Y
211 RCL N
212 ST* Y
213 X<> L
214 -
215 X<> IND Y
216 RCL O
217 *
218 X<> IND Z
219 RCL N
220 *
221 ST+ IND Z
222 RDN
223 RTN
224 LBL 10
225 RCL 00
226 ENTER^
227 X^2
228 +
229 2
230 /
231 RCL 00
232 1
233 ST+ Z
234 +
235   E3
236 /
237 +
238 RCL 00
239 X^2
240   E6
241 /
242 +
243 REGMOVE
244 RCL 01
245 CLA
246 CLD
247 END

( 358 bytes / SIZE (3n²+n+2)/2 )

STACK INPUTS OUTPUTS

X / k₁

                                1 2 4
Example:        A = 2 7 3
                                4 3 9

n = 3 therefore 3 STO 00 and we store the 6 numbers

     1                                         R01
     7 2        into registers          R02 R04                 respectively.
     9 4 3                                 R03 R05 R06

             XEQ "JCB2" yields ( in 1m12s )                   k₁ = 12.81993499 in R01
                                                                                  k₂ = 4.910741213 in R02
                                                                                  k₃ = -0.730676198 in R03

and the 3 corresponding eigenvectors: V₁ ( 0.351369026 ; 0.521535689 ;   0.777521917 ) in   ( R04 ; R05 ; R06 )
                                                           V₂ ( -0.101146468 ; 0.846760701 ; -0.522269766 ) in   ( R07 ; R08 ; R09 )
                                                           V₃ ( -0.930757326 ; 0.104865823 ; 0.350276976 )   in   ( R10 ; R11 ; R12 )

Notes:

-Here, the coefficients are stored differently.
-Store in successive registers starting with R01:

   -the diagonal elements a₁₁ a₂₂ ........ a_nn
   -then a₁₂   a₁₃   a₁₄ a₁₅ ......... a_1n
   -then a₂₃ a₂₄ a₂₅ .......... a_2n
      ............ and so on ......
   -until a_n-1,n

-Execution time = 10 seconds/iteration if n = 3
- -------------- = 12 ---------------- if n = 4
- -------------- = 22 ---------------- if n = 7

                                     1 2 4 7
-For the matrix   A =    2 3 7 1      execution time = 3mn25s    ( instead of 4mn22s with "JCB" )
                                     4 7 2 4
                                     7 1 4 9

-This program is very similar to "QUAD" ( cf "Quadratic Hypersurfaces for the HP-41" )

5°) Complex Matrices

-This program uses a variant of the Jacobi's algorithm:

*** The eigenvalues of A are the diagonal-elements of an upper triangular matrix T equal to the infinite product ...U_k^-1.U_k-1^-1....U₂^-1.U₁^-1.A.U₁.U₂....U_k-1.U_k....
where the U_k are unitary matrices. The eigenvectors are the columns of U = U₁.U₂....U_k-1.U_k.... if A is Hermitian ( i-e if A equals its transconjugate )
Actually, T is diagonal if A is Hermitian.

-"JCB" finds the greatest element a_{i j} below the main diagonal ( lines 23 to 60 )
-Then, U₁ is determined so that it places a zero in position ( i , j ) in U₁^-1.A.U₁

U₁ has the same elements as the Identity matrix, except that u_{i i} = u_{j j} = x and u_{i j} = y + i.z , u_{j i} = -y + i.z

with y + i.z = C.x , x = ( 1 + | C |² )^-1/2 and C = 2.a_{i j} / [ a_{i i} - a_{i j} + ( ( a_{i i} - a_{i j} )² + 4 a_{i j} a_{j i} )^1/2 ] ( line 156 R-P avoids a test if the denominator = 0 )

-The process is repeated until the greatest sub-diagonal element is smaller than E-8 in magnitude ( line 62 )

-The successive greatest | a_{i j} | ( i > j ) are displayed ( line 61 ) when the routine is running.

*** If the matrix is non-Hermitian, and if all the eigenvalues are distinct, we define an upper triangular matrix W = [ w_{i j} ] by

   w_{i j} = 0    if i > j
   w_{i i} = 1
   w_{i j} = Sum_{k = i+1, i+2 , .... , j}      t_{i k} w_{k j} / ( t_{j j} - t_{i i} )      if i < j

-Then, the n columns of U.W are the eigenvectors corresponding to the n eigenvalues t₁₁ = k₁ , ............... , t_nn = k_n

Data Registers:

                            •   R00 = n
                            •   R01,R02 = a₁₁     .................      • R2n²-2n+1,R2n²-2n+2 = a_1n
                              .................................................................................................                         ( the n² elements of A in column order )
   INPUTS
                            •   R2n-1,R2n = a_n1 .................      • R2n²-1,R2n² = a_nn

--------- odd registers = real part of the coefficients , even registers = imaginary part of the coefficients

   DURING            R01 thru R2n² = the "infinite" product   ...U_k^-1.U_k-1^-1....U₂^-1.U₁^-1.A.U₁.U₂....U_k-1.U_k....
       THE                R2n²+1 thru R4n² = U₁.U₂........U_k.....
ITERATIONS      R4n²+1 thru R4n²+2n: temp ( for non-Hermitian matrices only )

---------

                             R00 = n
                             R01,R02 = k₁          R2n+1,R2n+2 = V₁⁽¹⁾      ....................    R2n²+1,R2n²+2 = V₁⁽ⁿ⁾
                             .......................            ........................                ....................     ..................................
OUTPUTS
                             R2n-1,R2n = k_n        R4n-1,R4n =    V_n⁽¹⁾      .....................   R2n²+2n-1,R2n²+2n = V_n⁽ⁿ⁾

( n eigenvalues ; 1st eigenvector ; .................................. ; n^th eigenvector )

Flag: F02 Set flag F02 for a Hermitian matrix
Clear flag F02 for a non-Hermitian matrix
Subroutines: /

01 LBL "JCBZ"
02 RCL 00
03 X^2
04 ST+ X
05 .1
06 %
07 ST+ X
08 +
09 ISG X
10 CLRGX                  if you don't have an HP-41CX, replace line 10 by   ENTER^   SIGN CLX LBL 13 STO IND L ISG L GTO 13 CLX
11 RCL 00
12 1
13 +
14 ST+ X
15   E5
16 /
17 +
18 SIGN
19 LBL 00
20 STO IND L
21 ISG L
22 GTO 00
23 LBL 01
24 RCL 00
25 ST+ X
26   E3
27 /
28 ISG X
29 STO M
30 2
31 +
32 ENTER^
33 ENTER^
34 CLX
35 LBL 02
36 RCL IND Z
37 X^2
38 SIGN
39 ST+ T
40 CLX
41 RCL IND T
42 ST* X
43 ST+ L
44 X<> L
45 X>Y?
46 STO Z
47 X>Y?
48 +
49 RDN
50 ISG Z
51 GTO 02
52 2
53 ST+ M
54 CLX
55 RCL M
56 ST+ T
57 RDN
58 ISG Z
59 GTO 02
60 SQRT
61 VIEW X
62   E-8                      or another "small" number
63 X>Y?
64 GTO 06                a three-byte GTO
65 X<> Z
66 INT
67 STO M
68 2
69 ST- Y
70 /
71 STO Y
72 RCL 00
73 ST/ Z
74 MOD
75 ST+ X
76 RCL X
77 LASTX
78 *
79 +
80 X<>Y
81 INT
82 ST+ X
83 RCL X
84 RCL 00
85 *
86 +
87 2
88 ST+ Z
89 +
90 STO N
91 X<>Y
92 STO O
93 +
94 RCL M
95 -
96 RCL IND X
97 DSE Y
98 RCL IND Y
99 R-P
100 4
101 *
102 RCL IND M
103 RCL M
104 DSE X
105 RDN
106 RCL IND T
107 R-P
108 X<>Y
109 ST+ T
110 RDN
111 *
112 P-R
113 RCL IND N
114 RCL IND O
115 -
116 DSE N
117 DSE O
118 RCL IND N
119 STO N
120 CLX
121 RCL IND O
122 ST- N
123 RDN
124 STO O
125 RCL N
126 R-P
127 X^2
128 X<>Y
129 ST+ X
130 X<>Y
131 P-R
132 X<>Y
133 ST+ T
134 RDN
135 +
136 R-P
137 SQRT
138 X<>Y
139 2
140 /
141 X<>Y
142 P-R
143 RCL O
144 ST+ Z
145 X<> N
146 +
147 R-P
148 RCL IND M
149 DSE M
150 RCL IND M
151 R-P
152 ST+ X
153 R^
154 ST- Z
155 X<> T
156 R-P
157 CLX
158 SIGN
159 ST- M
160 P-R
161 STO O
162 RDN
163 P-R
164 STO N
165 X<>Y
166 X<> M
167 2
168 /
169 STO Y
170 RCL 00
171 ST/ Z
172 MOD
173 X<>Y
174 INT
175 RCL 00
176 ST+ X
177 ST* Y
178 ST* Z
179 +
180 .1
181 %
182 +
183 RCL 00
184 ST+ X
185 -
186 1
187 ST+ Z
188 +
189 RCL 00
190 X^2
191 ST+ X
192 ST+ Z
193 .1
194 %
195 +
196 +
197 XEQ 03
198 RCL 00
199 ST+ X
200 ST- Z
201 -
202 RCL 00
203 X^2
204 ST+ X
205 ST- Z
206 .1
207 %
208 +
209 -
210 XEQ 03
211 RCL 00
212 ST/ Z
213 /
214 INT
215 X<>Y
216 INT
217 RCL 00
218 X^2
219 ST+ X
220   E3
221 /
222 +
223 RCL 00
224 ST+ X
225 1
226 CHS
227 ST* M
228 ST* N
229 ST+ Z
230 ST+ T
231 +
232   E5
233 /
234 ST+ Z
235 +
236 XEQ 03
237 GTO 01                a three-byte GTO
238 LBL 03
239 STO P                  synthetic
240 LBL 04
241 CLX
242 RCL IND P
243 RCL O
244 *
245 RCL IND Y
246 RCL N
247 *
248 +
249 PI
250 SIGN
251 ST+ T
252 CLX
253 RCL IND T
254 RCL M
255 *
256 -
257 STO Q
258 CLX
259 RCL IND Y
260 RCL O
261 *
262 RCL IND P
263 RCL N
264 *
265 -
266 PI
267 SIGN
268 ST+ P
269 CLX
270 RCL IND P
271 RCL M
272 *
273 -
274 X<> IND Y
275 RCL M
276 *
277 RCL IND P
278 RCL O
279 *
280 +
281 PI
282 SIGN
283 ST+ Z
284 CLX
285 RCL IND Z
286 RCL N
287 *
288 +
289 X<> IND P
290 RCL N
291 *
292 RCL IND Y
293 RCL O
294 *
295 X<>Y
296 -
297 RCL P
298 SIGN
299 ST- L
300 X<> Q
301 X<> IND L
302 RCL M
303 *
304 +
305 STO IND Y
306 ISG Y
307 CLX
308 ISG P
309 GTO 04
310 X<> P
311 RTN
312 LBL 05
313 RCL P
314 SIGN
315 ST- L
316 RCL IND P
317 RCL IND L
318 R-P
319 RCL IND 04
320 DSE 04
321 RCL IND 04
322 R-P
323 X<>Y
324 ST+ T
325 RDN
326 *
327 P-R
328 RTN
329 LBL 06
330 CLD
331 FS? 02
332 GTO 11                a three-byte GTO          Lines 333 thru 471 ( and lines 312 thru 328 ) are unuseful if the matrix is Hermitian
333 RCL 00
334 X^2
335 ST+ X
336 LASTX
337 1
338 +
339 ST+ X
340   E5
341 /
342 +
343 STO O
344 RCL 00
345 DSE X
346   E3
347 /
348 ISG X
349 STO M
350 LBL 07
351 RCL 00
352 ST+ X
353 ENTER^
354 X^2
355 +
356 STO Y
357 RCL M
358 INT
359 ST+ X
360 STO 03                 we can use registers R03 & R04 for temporary data storage after the triangularization
361 -
362   E3
363 /
364 +
365 STO 04
366 RCL O
367 RCL 03
368 -
369 2 E-5
370 -
371 STO P                  synthetic
372 CLX
373 STO 03
374 STO N
375 STO IND Y
376 SIGN
377 ST- Y
378 STO IND Y
379 LBL 08
380 XEQ 05
381 ST+ 03
382 X<>Y
383 ST+ N
384 DSE P
385 DSE 04
386 GTO 08
387 RCL IND O
388 RCL IND P
389 -
390 PI
391 SIGN
392 ST- 04
393 ST- O
394 ST- P
395 CLX
396 RCL IND O
397 RCL IND P
398 -
399 R-P
400 RCL N
401 RCL 03
402 R-P
403 R^
404 ST- Z
405 X<> T
406 /                            there will be a DATA ERROR line 406 if all the eigenvalues are not distinct
407 P-R
408 STO IND 04
409 CLX
410 SIGN
411 ST+ 04
412 ST+ O
413 X<>Y
414 STO IND 04
415 ISG M
416 GTO 07
417 RCL 00
418 STO 03
419 DSE M
420 LBL 09
421 RCL 00
422 X^2
423 ENTER^
424 STO 04
425 ST+ 04
426 LASTX
427 ST+ 04
428 RCL M
429 INT
430 ST- 04
431 *
432 +
433 STO N
434   E3
435 /
436 +
437 RCL 03
438 ST+ N
439 +
440 RCL 00
441   E5
442 /
443 +
444 2
445 ST* 04
446 ST* N
447 *
448 STO P                      synthetic
449 LBL 10
450 XEQ 05
451 RCL N
452 DSE X
453 RDN
454 ST+ IND T
455 X<>Y
456 ST+ IND N
457 PI
458 INT
459 ST+ 04
460 ISG P
461 GTO 10
462 DSE 03
463 GTO 09
464 RCL M
465 FRC
466   E-3
467 -
468 STO M
469 DSE O
470 ISG M
471 GTO 07                             a three-byte GTO
472 LBL 11
473 RCL 00                              Lines 473 thru 511 are not essential:
474 X^2                                    if you delete these lines, the eigenvalues will be the diagonal elements of the matrix in registers R01 thru R2n²
475 ST+ X                                and the eigenvectors will be the columns of the matrix in registers R2n²+1 thru R4n²
476 STO M
477 LASTX
478 ST+ X
479 STO Z
480 1
481 +
482   E5
483 /
484 +
485 LBL 12
486 RCL IND X
487 STO IND Z
488 CLX
489 SIGN
490 ST- Z
491 -
492 RCL IND X
493 STO IND Z
494 DSE Y
495 RDN
496 DSE Y
497 GTO 12
498 CLX
499 RCL M
500 R^
501 1
502 ST+ Z
503 +
504   E3
505 /
506 +
507 RCL M
508   E6
509 /
510 +
511 REGMOVE
512 RCL 02
513 RCL 01
514 CLA
515 END

( 780 bytes / SIZE 4n²+1 if A is Hermitian / SIZE 4n²+2n+1 if A is non-Hermitian )

      STACK        INPUTS      OUTPUTS

           Y             /            y₁

           X             /            x₁

where k₁ = x₁ + y₁ = the first eigenvalue.

                                  1      4 -7.i   3-4.i
Example1:     A = 4+7.i      6      1-5.i             A is equal to its transconjugate so A is Hermitian
                               3+4.i   1+5.i      7

                   1 0   4 -7 3 -4           R01 R02    R07 R08    R13 R14
-We store   4 7   6 0   1 -5   into   R03 R04    R09 R10    R15 R16    respectively     and n = 3 STO 00
                  3 4   1 5   7 0            R05 R06    R11 R12    R17 R18

SF 02 ( the matrix is Hermitian ) XEQ "JCBZ" yields in 5mn07s

           k₁ = 15.61385271 + 0.i = ( X , Y ) = ( R01 , R02 )
           k₂ = 5.230678474 + 0.i = ( R03 , R04 )
           k₃ = -6.844531162 + 0.i = ( R05 , R06 )   and the 3 corresponding eigenvectors

    V₁( 0.481585119 + 0.108201123 i , 0.393148652 + 0.527305194 i , -0.142958847 + 0.550739892 i )    = ( R07 R08 , R09 R10 , R11 R12 )
    V₂( -0.436763638 + 0.075843695 i , 0.210271354 - 0.463555612 i , -0.516799072 + 0.526598642 i )   = ( R13 R14 , R15 R16 , R17 R18 )
    V₃( -0.706095475 + 0.247553486 i , 0.326530385 + 0.449069516 i , 0.363126603 - 0. i )                      = ( R19 R20 , R21 R22 , R23 R24 )

-Due to roundoff errors, registers R02 R04 R06 contain very small numbers instead of 0
but actually, the eigenvalues of a Hermitian matrix are always real.

                               1+2.i   2+5.i   4+7.i
Example2:     A = 4+7.i   3+6.i   3+4.i        A is non-Hermitian
                               3+4.i   1+7.i   2+4.i

                   1 2   2 5   4 7           R01 R02    R07 R08    R13 R14
-We store   4 7   3 6   3 4   into   R03 R04    R09 R10    R15 R16    respectively     and n = 3 STO 00
                  3 4   1 7   2 4            R05 R06    R11 R12    R17 R18

CF 02 ( the matrix is non-Hermitian ) XEQ "JCBZ" produces in 6mn45s

           k₁ = 7.656606009 + 15.61073835 i = ( X , Y ) = ( R01 , R02 )
           k₂ = 1.661248138 - 1.507335315 i   = ( R03 , R04 )
           k₃ = -3.317854151 - 2.103403073 i = ( R05 , R06 )   and the 3 corresponding eigenvectors

    V₁( 0.523491429 + 0.008555748 i , 0.650242685 - 0.046353000 i , 0.546288477 + 0.049882590 i )      = ( R07 R08 , R09 R10 , R11 R12 )
    V₂( -0.4777850326 - 0.196368365 i , 0.660547554 - 0.265888145 i , -0.356842635 + 0.318267519 i ) = ( R13 R14 , R15 R16 , R17 R18 )
    V₃( -0.567221101 - 0.574734664 i , 0.141603481 + 0.549379028 i , 0.480533312 - 0.090337847 i )     = ( R19 R20 , R21 R22 , R23 R24 )

Notes:

-If n = 4 , a typical execution time is 19 minutes for a non-Hermitian matrix.
-If you don't want to get the eigenvectors, but you do want to compute the eigenvalues:

Delete lines 498 to 511 , line 476 , lines 331 to 472 , lines 312 to 328 , lines 189 to 209 and lines 02 to 22 ( 299 lines / 450 bytes / SIZE 2n²+1 )

-The eigenvalues of the 2 matrices above are now obtained in 3mn59s & 4mn39s respectively.
-Since the SIZE is now 2n²+1 , you might calculate the eigenvalues of a 12x12 complex matrix! Not very quickly however:
-For a complex matrix of order 7 , the eigenvalues are computed in about 2h30mn, but with a good emulator ...

References: Francis Scheid "Numerical Analysis" McGraw-Hill ISBN 07-055197-9
J-M Ferrard "Mastering your HP-48" D3I Diffusion ISBN 2-908791-04-8

Go back to the HP-41 software library
Go back to the general software library
Go back to the main exhibit hall