The Museum of HP Calculators


Eigenvalues & Eigenvectors for the HP-41

This program is Copyright © 2003-2006 by Jean-Marc Baillard and is used here by permission.

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 A and a corresponding eigenvector by the power method  ( n < 17 )
  3°) "DFL" uses a deflation method to produce the real eigenvalues ki and the eigenvectors, provided that | k1 | >  | k2 | > ....... >  | kn | > 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) = xn + c1.xn-1 + c2.xn-2 + ......... + cn-2.x2 + cn-1.x + cn 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 Mk by:      M0 = I   and  Mk+1 = Mk.A - trace(Mk.A) I /k      ( the trace of a matrix is the
                                                                                             then we have:       ck  = - trace(Mk-1.A)/k                                         sum of its diagonal elements )

-The elements of A are to be stored into contiguous registers from R01 to Rn2 , n being stored in R00
-When the program stops,  R01 = c1 ; R02 = c2 ; ............... ; Rnn = cn            ( c1 = -trace(A)  and cn = (-1)n det(A) )
 

Data Registers:

                               R00 =  n
                               R01 = a11    ......................      Rn2-n+1 = a1n
                              ...................................................................                          ( the n2 elements of A )
   INPUTS
                               Rnn = an1    ......................      Rn2 = ann

   ---------

                          R00 =  n
                   R01 = c1       Rn+1 =  a11          Rn2+1 = a1n
                         .............           ..........................................
  OUTPUTS
                   Rnn = cn         R2n = an1           Rn2+n  = ann
 

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 3n2+n+1 )
 
 
      STACK        INPUT      OUTPUT
           X             /            cn

                                                                                                   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"   >>>>  c3 = 5 = R03   ( in 31 seconds )
   and  R01 = -11  ,  R02 = -25  ,  R03 = 5     Thus p(x) = x3 - 11.x2 - 25.x + 5
 

Notes:

-A is saved in registers Rn+1 thru Rn2+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    k1 = 12.90692994  ,    k2 = -2.092097589  ,   k3 =  0.185167648
-Then, subtracting ( for instance ) k1 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 n2 registers to store the matrices A , Mk and the products Mk.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 multiplication of 2 square-matrices:
-Thus, it can produce the characteristic polynomial of a 16x16 matrix!
 

Data Registers:

                               R00 =  n
                               R01 = a11    ......................      Rn2-n+1 = a1n
                              ..................................................................                ( the n2 elements of A )
   INPUTS
                               Rnn = an1    ......................      Rn2 = ann

   ---------

                          R00 =  n
  OUTPUTS      R01 = c1  ......................  Rnn = cn

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 n2+2n+1 )
 
 
      STACK        INPUT      OUTPUT
           X             /            cn

                                                                                                     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"   >>>>  c3 = 5 = R03   ( in 37 seconds )
   and  R01 = -11  ,  R02 = -25  ,  R03 = 5     whence    p(x) = x3 - 11.x2 - 25.x + 5
 

Example2:      The 10x10 matrix A = [ aij ] is defined by  aij = 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) = x10 - 43 x9 - 968 x8 - 2462 x7 + 40796 x6 - 488852 x5 - 10916340 x4 + 15630136 x3 + 441980832 x2 - 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   k1 ; k2 ; ....... ; kn  corresponding to n independent eigenvectors  V1 ; V2 ; ....... ; Vn
  and that  k1 is the dominant eigenvalue | k1 | > | ki |    i = 2 ; 3 ; ....... ; n
-Any vector V can be written    V =  a1 V1 + a2 V2 + ....... + an Vn  whence AV  =  a1 k1 V1 + a2 k2 V2 + ....... + an kn Vn
 and  ApV = k1p ( a1 V1 + a2 ( k2/k1 )p V2 + ....... + an ( kn/k1 )p Vn )  which tends to  k1p  a1 V1 as p tends to infinity.

-Furthermore,  if we define  Wp+1 = ( A.Wp ) / || A.Wp ||   with  W0 = V ( almost ) arbitrary , then Wp tends to a unit eigenvector U
  and  tWp AWp tends to the corresponding eigenvalue  k1
 

Data Registers:

                               R00 =  n
                               R01 = v1              Rn+1 = a11          Rn2+1 = a1n
                              .................                  ( the n2 elements of A )                       where V (  v1 ; ....... ;  vn )   is an initial guess of the eigenvector.
   INPUTS
                               Rnn =  vn             R2n = an1             Rn2+n = ann

   ---------

                       R00 =  n
                R01 = u1       Rn+1 =  a11          Rn2+1 = a1n
                       .............           .............................................      and k in X-register.   U ( u1; ........ ; un) = a unit eigenvector
  OUTPUTS
                Rnn = un         R2n = an1           Rn2+n  = ann
 

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 n2+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 | k1/ k2 | , 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/k3 = 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  k2 - 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'1/ k'2 | = 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(  x1 , ....... ,  xn )  is choosen so that its first coordinate = the eigenvalue k.
-So, it will not work if x1 = 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 n2+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  Rn2+n+1 ...... Rn2+2n  ( in this example,  R13  R14  R15 )
-Of course, the last decimals may differ slightly.
 

3°) Deflation
 

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

-"DFL" uses this method to calculate all the eigenvalues and eigenvectors , provided they are real and verify   | k1 | >  | k2 | > ....... >  | kn | > 0

  1°) "DFL" calculates k1 and U1  ( with flag F02 clear )  and the program stops ( line 87 )

     >>>>     k1 is in X-register
     >>>>     U1  is in registers R01 thru Rnn

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

     >>>>     k2 is in X-register
     >>>>     U2  is in registers R01 thru Rnn

   ... etc ...
 

Data Registers:

                               R00 =  n
                               R01 = v1              Rn+1 = a11          Rn2+1 = a1n
                              .................                  ( the n2 elements of A )                       where V (  v1 ; ....... ;  vn )   is an initial guess of the eigenvector.
   INPUTS
                               Rnn =  vn             R2n = an1             Rn2+n = ann

   ---------

                      R00 =  n
               R01 = u1       Rn+1 =  a'11          Rn2+1 = a'1n

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

               Rnn = un         R2n = a'n1           Rn2+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 n2+3n+1 )
 
 
      STACK        INPUTS      OUTPUTS
           X             /         k1 , k2, ...

                                                                                                                                                        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 k1-approximations ( with F02 clear )

- Eventually,  k1 = 12.90692995 is in X-register and U1 is in registers R01 R02 R03 :    U1 ( 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 k1-approximations converging to 12.90692995 again
     then flag F02 is cleared and the HP-41 displays the successive k2-approximations and

          k2 = -2.092097594 in X-register and  U2 in  R01R02R03 :    U2 ( 0.800454175 ; -0.041651079 ; -0.597945065 )

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

          k3 =   0.185167649 in X-register and U3 in R01R02R03 :    U3 ( -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  k1 and  k2 are real but  k3 ...  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   ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
 where the Ok are rotation matrices. The eigenvectors are the columns of  O1.O2....Ok-1.Ok....

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

-The successive greatest | ai 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 = [ ti j ] =  ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....  is an upper triangular matrix,
 the diagonal elements of T ( i-e  ti i = ki ) are the n eigenvalues and the first column of  U = O1.O2....Ok-1.Ok....   is an eigenvector corresponding to k1

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

   wi j = 0    if  i > j
   wi i = 1
   wi j = Sumk = i+1, i+2 , .... , j      ti k wk j / ( tj j - ti i )      if  i < j

-Then, the n columns of  U.W  are the eigenvectors corresponding to the n eigenvalues  t11 = k1 , ............... , tnn = kn

*** 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 nth eigenvector is stored ---------- Rn2+1 thru Rn2+n  as shown below.
 
 

Data Registers:
                               R00 =  n
                               R01 = a11  .................        Rn2-n+1 = a1n
                              ..................................................................              ( the n2 elements of A in column order )
   INPUTS
                               Rnn = an1  .................        Rn2 = ann

   ---------
 

   DURING            R01 thru Rn2 = the "infinite" product   ...Ok-1.Ok-1-1....O2-1.O1-1.A.O1.O2....Ok-1.Ok....
       THE                Rn2+1 thru R2n2 =  O1.O2........Ok.....
 ITERATIONS      R2n2+1 thru R2n2+n: temp  ( for non-symmetric matrices only )

  ---------

                             R00 = n
                             R01 = k1       Rn+1 = V1(1)       R2n+1 = V1(2)  ....................    Rn2+1 = V1(n)
                             .............
  OUTPUTS
                             Rnn = kn        R2n =  Vn(1)         R3n =  Vn(2)    .....................   Rn2+n = Vn(n)
 

                        ( n eigenvalues  ;  1st eigenvector ;  2nd eigenvector ; .................. ;  nth 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 | aij |   ( 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 2n2+1 if A is symmetric / SIZE 2n2+n+1 if A is not symmetric )
 
 
      STACK        INPUTS      OUTPUTS
           X             /            k1

                                                                                                                                  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 )                               k1 =  12.81993499  in R01
                                                                                  k2 =  4.910741214  in R02
                                                                                  k3 = -0.730676199  in R03

and the 3 corresponding eigenvectors:  V1 (  0.351369026 ; 0.521535689 ;  0.777521917 )   in   ( R04 ; R05 ; R06 )
                                                           V2 ( -0.101146468 ; 0.846760701 ; -0.522269766 )  in   ( R07 ; R08 ; R09 )
                                                           V3 ( -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"  >>>>    k1 = 12.90692994    ( in 2mn44s )

  and we have          k1 =  12.90692994  in R01
                               k2 =  0.185167649  in R02
                               k3 = -2.092097593  in R03

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

-The eigenvectors are not unit vectors except the first one.
-Divide V2 & V3 by  || V2 || & || V3 || 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 Oi-1.M.Oi
  and only 2 columns in M.Oi .

-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 Ok-1.A.Ok  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:        k1 =  16.97583168  ,   k2 =  6.365547530  ,   k3 = -3.301311094  ,   k4 = -5.040068160

       V1 ( 0.455772321 ;  0.346041152  ;  0.464075961  ;  0.676136537 )
       V2 ( -0.142731960 ;  0.681492880 ;  0.448494335 ; -0.560399745 )
       V3 ( 0.842568185 ;  -0.247658954 ;  0.050153515 ; -0.475634862 )
       V4 ( -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  ai j  and   aj 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 = a11    Rn+1 = a1,2   R2n    = a2,3   ..................................  R(n2+n)/2 = an-1,n
                            R02 = a22    Rn+2 = a1,3   R2n+1 = a2,4     ...................
   INPUTS          ...............         .................      ...................
                           ...............         .................    R3n-2 = a2,n
                           ................      R2n-1 = a1,n
                            Rnn = ann

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

                When the program stops:
 

                             R00 = n
                             R01 = k1       Rn+1 = V1(1)       R2n+1 = V1(2)  ....................    Rn2+1 = V1(n)
                             .............
  OUTPUTS
                             Rnn = kn        R2n =  Vn(1)         R3n =  Vn(2)    .....................   Rn2+n = Vn(n)
 

                        ( n eigenvalues  ;  1st eigenvector ;  2nd eigenvector ; .................. ;  nth 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 (3n2+n+2)/2  )
 
 
      STACK        INPUTS      OUTPUTS
           X             /            k1

                                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 )                   k1 =  12.81993499  in R01
                                                                                  k2 =  4.910741213  in R02
                                                                                  k3 = -0.730676198  in R03

and the 3 corresponding eigenvectors:  V1 (  0.351369026 ; 0.521535689 ;   0.777521917 )  in   ( R04 ; R05 ; R06 )
                                                           V2 ( -0.101146468 ; 0.846760701 ; -0.522269766 )  in   ( R07 ; R08 ; R09 )
                                                           V3 ( -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  a11  a22 ........  ann
   -then  a12   a13   a14  a15  .........  a1n
   -then  a23  a24  a25  ..........  a2n
      ............ and so on ......
   -until  an-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   ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
 where the Uk are unitary matrices. The eigenvectors are the columns of  U = U1.U2....Uk-1.Uk.... if A is Hermitian ( i-e if A equals its transconjugate )
 Actually, T is diagonal if A is Hermitian.

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

  U1  has the same elements as the Identity matrix, except that  ui i = uj j = x  and  ui j = y + i.z ,  uj i = -y + i.z

     with  y + i.z  = C.x  ,  x = ( 1 + | C |2 ) -1/2  and   C = 2.ai j / [ ai i - ai j + ( ( ai i - ai j )2 + 4 ai j aj 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 | ai 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 = [ wi j ]  by

   wi j = 0    if  i > j
   wi i = 1
   wi j = Sumk = i+1, i+2 , .... , j      ti k wk j / ( tj j - ti i )      if  i < j

-Then, the n columns of  U.W  are the eigenvectors corresponding to the n eigenvalues  t11 = k1 , ............... , tnn = kn
 

Data Registers:

                               R00 =  n
                               R01,R02 = a11     .................        R2n2-2n+1,R2n2-2n+2 = a1n
                              .................................................................................................                         ( the n2 elements of A in column order )
   INPUTS
                               R2n-1,R2n = an1  .................        R2n2-1,R2n2 = ann

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

   DURING            R01 thru R2n2 = the "infinite" product   ...Uk-1.Uk-1-1....U2-1.U1-1.A.U1.U2....Uk-1.Uk....
       THE                R2n2+1 thru R4n2 =  U1.U2........Uk.....
 ITERATIONS      R4n2+1 thru R4n2+2n: temp  ( for non-Hermitian matrices only )

  ---------

                             R00 = n
                             R01,R02 = k1          R2n+1,R2n+2 = V1(1)      ....................    R2n2+1,R2n2+2 = V1(n)
                             .......................            ........................                ....................     ..................................
  OUTPUTS
                             R2n-1,R2n = kn        R4n-1,R4n =    Vn(1)      .....................   R2n2+2n-1,R2n2+2n = Vn(n)
 

                            ( n eigenvalues     ;        1st eigenvector      ;  .................................. ;    nth 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 R2n2
475  ST+ X                                and the eigenvectors will be the columns of the matrix in registers R2n2+1 thru R4n2
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 4n2+1 if A is Hermitian / SIZE 4n2+2n+1 if A is non-Hermitian )
 
 
      STACK        INPUTS      OUTPUTS
           Y             /             y1
           X             /            x1

  where  k1 = x1 + y1  = 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

           k1 =  15.61385271 + 0.i  = ( X , Y ) = ( R01 , R02 )
           k2 =  5.230678474 + 0.i  = ( R03 , R04 )
           k3 = -6.844531162 + 0.i  = ( R05 , R06 )   and the 3 corresponding eigenvectors

    V1( 0.481585119 + 0.108201123 i , 0.393148652 + 0.527305194 i , -0.142958847 + 0.550739892 i )    =  ( R07 R08 , R09 R10 , R11 R12 )
    V2( -0.436763638 + 0.075843695 i , 0.210271354 - 0.463555612 i , -0.516799072 + 0.526598642 i )   =  ( R13 R14 , R15 R16 , R17 R18 )
    V3( -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

           k1 = 7.656606009 + 15.61073835 i  = ( X , Y ) = ( R01 , R02 )
           k2 = 1.661248138 - 1.507335315 i   = ( R03 , R04 )
           k3 = -3.317854151 - 2.103403073 i  = ( R05 , R06 )   and the 3 corresponding eigenvectors

    V1( 0.523491429 + 0.008555748 i , 0.650242685 - 0.046353000 i , 0.546288477 + 0.049882590 i )      =  ( R07 R08 , R09 R10 , R11 R12 )
    V2( -0.4777850326 - 0.196368365 i , 0.660547554 - 0.265888145 i , -0.356842635 + 0.318267519 i )  =  ( R13 R14 , R15 R16 , R17 R18 )
    V3( -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 2n2+1 )

-The eigenvalues of the 2 matrices above are now obtained in 3mn59s & 4mn39s respectively.
-Since the SIZE is now 2n2+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