# hp41programs

Egvl Eigenvalues & Eigenvectors for the HP-41

Overview

1°)  Characteristic Polynomial & Minimal Polynomial

a)  Program#1
b)  Program#2
c)  Minimal Polynomial

2°) The Power Method: Dominant Eigenvalue

a)  Program#1
b)  Program#2

3°) Deflation

4°) The Jacobi's Method

a)  Symmetric ( and some non-symmetric ) Matrices
b)  Symmetric Matrices only

5°) Complex Matrices

-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 )
"PMIN" ----------   minimal       ----------------  ( n < 13 )
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 - except "PMIN" - use the REGMOVE  function of the X-functions module.

1°)  Characteristic Polynomial & Minimal 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+1)     k = 0 , 1 , ....      ( the trace of a matrix is the sum of its diagonal elements )
and we have:       ck  = - trace(Mk-1.A)/k

-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

Example:   Find the characteristic polynomial of the matrix

1  2  4
A =  4  3  5
7  4  7

-First we store these 9 elements in registers R01 thru R09

R01  R04  R07         1  2  4
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" )

-Line 110 is a three-byte  GTO 01
-Do not key in lines 112-113 if you want to save the matrix A in this data file.

 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 111  LBL 07 112  "T"  113  PURFL 114  CLA 115  END

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

 STACK INPUT OUTPUT X / cn

Example1:   Find the characteristic polynomial of the matrix

1  2  4
A =  4  3  5
7  4  7

-First we store these 9 elements in registers R01 thru R09

R01  R04  R07         1  2  4
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...

c) The Minimal Polynomial

-The characteristic polynomial P verifies P(A) = 0 ( Cayley-Hamilton theorem )
-The minimal polynomial p is the unitary least degree polynomial that satisfies  p(A) = 0.
-It is a divisor of the characteristic polynomial.

-In the following program, we start with a random vector V and we compute  AV , A2V , ........... , AnV
-Then, "PMIN" seeks the smallest k such that  V , AV , A2V , ........... , AkV  are linearly dependent.

-A relation  a0 V + a1 AV + ........... + ak-1 Ak-1V + AkV = 0  is found, whence p(x)

Data Registers:             •   R00 =  n                               ( Registers R00 thru Rn2 are to be initialized before executing "PMIN" )

•   R01 = a11    ......................    •  Rn2-n+1 = a1n
...................................................................                          ( the n2 elements of A )

•   Rnn = an1    ......................    •  Rn2 = ann                            Rn2+1 to R2n2+n: temp

>>>>   When the program stops,  R01 thru Ree = the coefficients of p(x)

Flags:  /
Subroutine:  "M*"  "MCO"   ( cf "Matrix Operations for the HP-41" )
"LS3"     ( cf "Linear and Non-Linear Systems for the HP-41" )

-Lines 11-12 may be replaced by the M-code function RAND  ( cf "A few M-Code Routines fot the HP-41" )

 01  LBL "PMIN"   02  RCL 00   03  X^2   04  ST+ X   05  .1   06  %   07  +   08  RCL 00           09  ST+ Y   10  LBL 00   11  R-D   12  FRC   13  STO IND Y   14  DSE Y   15  GTO 00   16  CLX   17  RCL 00   18  .00101   19  *   20  +   21  1   22  + 23  ENTER^   24  LBL 01   25  CLX   26  .01   27  RCL 00           28  ST+ Y   29  *   30   E3   31  /   32  1   33  +   34  X<>Y   35  ENTER^   36  INT   37  RCL 00   38  -   39  XEQ "M*"   40  STO Y   41  INT   42  DSE X   43  RCL 00   44  X^2 45  -   46  X#0?   47  GTO 01   48  LASTX   49  ST+ X   50  RCL 00           51  +   52   E3   53  /   54  RCL 00   55  X^2   56  +   57  1   58  ST+ Y   59  XEQ "MCO"   60  RCL 00   61  RCL 00   62  1   63  +   64  XEQ "LS3'   65  R^   66  INT 67  STO 00   68  1   69  +   70   E5   71  /   72  RCL 00           73  X^2   74  +   75  RCL 00   76  +   77  SIGN   78  ENTER^   79  LBL 02   80  RCL IND L    81  X#Y?   82  GTO 03   83  RDN   84  ST+ Y   85  DSE L   86  GTO 02   87  ENTER^   88  LBL 03 89  LASTX   90  INT   91  R^   92  ST+ Y   93  ENTER^   94  DSE Z   95  LBL 04   96  RCL IND Z   97  CHS   98  STO IND Z    99  RDN 100  DSE Z 101  "" 102  DSE Y 103  GTO 04 104   E3 105  / 106  1 107  STO 01         108  + 109  END

( 164 bytes / SIZE 2n2+n+1 )

 STACK INPUT OUTPUT X / 1.eee

Example:   Find the minimal polynomial of the matrix

[ [  3  -1   1 ]
A = [ -1   3   1 ]
[  1   1   3 ] ]

-Store
3   -1    1                     R01  R04  R07
-1    3    1        into        R02  R05  R08        respectively     and      n = 3   STO 00
1    1    3                      R03  R06  R09

XEQ "PMIN"  >>>>   1.003     --- Execution time = 36s ---     and we get  R01 = 1   R02 = -5  R03 = 4      ( with some roundoff-errors )

-So,  p(x) = x2 - 5 x + 4       ( the characteristic polynomial is  P(x) = x3 - 9 x2 + 24 x - 16 = ( x - 4 ) p(x)  )

Notes:

-Since we use a pseudo-random vector V, there is a small risk of getting an eigenvector, which would lead to a wrong result.
-In practice however, that seems highly improbable.
-If the coefficients of A are integers, the coefficients of p are integers too.
-Roundoff-errors may be important, especially if the system is ill-conditioned,
and E-7 is sometimes too small to identify tiny pivots ( cf "LS3" listing, line 61 ).
-If you prefer "LS2" instead of "LS3", use the listing below:

-Lines 13-14 may be replaced by the M-code function RAND  ( cf "A few M-Code Routines fot the HP-41" )

 01  LBL "PMIN"   02  RCL 00   03  X^2   04  STO Y   05  LASTX   06  +   07   E3   08  /   09  +   10  ISG X   11  RCL 00           12  LBL 00   13  R-D    14  FRC   15  STO IND Y   16  ISG Y   17  GTO 00   18  CLX   19  RCL 00   20  ST- Y   21   E5   22  /   23  + 24  ENTER^   25  LBL 01   26  CLX   27  .01   28  RCL 00   29  ST+ Y   30  *   31   E3   32  /   33  1   34  +   35  X<>Y   36  ENTER^   37  INT   38  RCL 00           39  +   40  XEQ "M*"    41  STO Y   42  INT   43  DSE X   44  RCL 00   45  X^2   46  ST+ X 47  -   48  X#0?   49  GTO 01   50  LASTX   51  RCL 00   52  +   53   E3   54  /   55  RCL 00   56  X^2   57  +   58  1   59  ST+ Y   60  XEQ "MCO"   61  RCL 00           62  RCL 00   63  1   64  +   65  XEQ "LS2'   66  X<> Z   67  INT   68  ENTER^   69  SQRT 70  STO 00           71  1   72  +   73   E2   74  /   75  +   76   E3   77  /   78  1   79  +   80  SIGN   81  ENTER^   82  LBL 02   83  RCL IND L    84  X#Y?   85  GTO 03   86  RDN   87  ST+ Y   88  ISG L   89  GTO 02   90  ENTER^   91  LBL 03   92  LASTX 93  INT   94  R^   95   E3   96  /   97  2   98  +   99  DSE Y 100  LBL 04 101  RCL IND Y 102  CHS 103  STO IND Y 104  RDN 105  DSE Y 106  "" 107  ISG X 108  GTO 04 109  FRC 110  1 111  STO 01         112  + 113  END

( 165 bytes / SIZE 2n2+n+1 )

 STACK INPUT OUTPUT X / 1.eee

-The same example gives the same results - with small differences in the last decimals.

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 74  X

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

 STACK INPUTS OUTPUTS X / k

Example:   Find the dominant eigenvalue k and the corresponding unit eigenvector U of the matrix

1  2  4
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    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

( 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 R01-R02-R03 :    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" )

-Lines 82-102 are three-byte GTO 01
-Line 142 is a three-byte GTO 12

 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   81  XY 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 145  END

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

 STACK INPUTS OUTPUTS X / k1 , k2 , ...

Example:   Find all the eigenvalues and the corresponding unit eigenvectors of the same matrix

1  2  4
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     in      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 in registers Rn+1 thru R2n
..................................................................................
the nth eigenvector is stored in registers 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:  /

-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
-Line 47 ,  the HP-41 displays the successive greatest sub-diagonal elements | aij |   ( i > j )
-Line 50 is a three-byte GTO 05
-Line 132 is a three-byte GTO 01
-Lines 177 to 269 are unuseful if the matrix is symmetric
-Line 178 is a three-byte GTO 11
-Line 268 is a three-byte GTO 07

 01  LBL "JCB"   02  RCL 00   03  X^2   04  .1   05  %   06  ST+ X   07  +   08  ISG X   09  CLRGX   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   48   E-8 49  X>Y?   50  GTO 05   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 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 178  GTO 11 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 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 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

Example1:  Let's compute all the eigenvalues and the eigenvectors of the matrix

1  2  4
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    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 )

Example2:

1  2  4
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     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
-Execution time  = 16 seconds/iteration   if n = 4
-Execution time  = 27 seconds/iteration   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!

Exercise:      Find all the eigenvalues and eigenvectors of the 4x4 matrix

1  2  4  7
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:  /

-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
-Line 56 is a three-byte GTO 10
-Line 197 is a three-byte GTO 01

 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    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   55  X>Y?   56  GTO 10   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  XY   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   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 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 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

Example:

1  2  4
A =  2  7  3
4  3  9

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

1                          R01
7  2        into        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
-Execution time  = 12 seconds/iteration   if n = 4
-Execution time  = 22 seconds/iteration   if n = 7

-For the matrix:

1  2  4  7
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:  /

-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
-Line 64 is a three-byte GTO 06
-Line 237 is a three-byte GTO 01
-Line 332 is a three-byte GTO 11
-Lines 333 thru 471 ( and lines 312 thru 328 ) are unuseful if the matrix is Hermitian

-There will be a DATA ERROR line 406 if all the eigenvalues are not distinct
-Line 471 is a three-byte  GTO 07

 01  LBL "JCBZ"   02  RCL 00   03  X^2   04  ST+ X   05  .1   06  %   07  ST+ X   08  +   09  ISG X   10  CLRGX   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    63  X>Y?   64  GTO 06   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 238  LBL 03 239  STO P 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 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 361  - 362   E3 363  / 364  + 365  STO 04         366  RCL O 367  RCL 03 368  - 369  2 E-5 370  - 371  STO P 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  /  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  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 472  LBL 11 473  RCL 00 474  X^2 475  ST+ X 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.

Example1:

1      4 -7.i   3-4.i
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
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.

Example2:

1+2.i   2+5.i   4+7.i
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
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