hp41programs

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 k_i and the eigenvectors, provided that | k₁ | >  | k₂ | > ....... >  | k_n | > 0   ( n < 17 )
  4°) "JCB" uses the Jacobi's method to produce all the n real eigenvalues and eigenvectors of a symmetric matrix of order n ( n < 13 )
        -However, it also works if the matrix is not symmetric provided all the eigenvalues are real and distinct. ( n < 13 )
        "JCB2" runs faster and uses less data registers but it works for symmetric matrices only ( n < 15 )
  5°) "JCBZ" uses a generalization of the Jacobi's method to compute all the eigenvalues and eigenvectors of a Hermitian matrix. ( n < 9 )
        -It also works for a general complex matrix provided all the eigenvalues are distinct. ( n < 9 )
 

-All these programs - 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) = xⁿ + c₁.x^n-1 + c₂.x^n-2 + ......... + c_n-2.x² + c_n-1.x + c_n is defined by p(x) = det ( x.I - A )
    where I is the n x n unit matrix.

-Thus, its roots are actually the n eigenvalues of A.

-"CP" uses the following formulae: we define a sequence of matrices M_k by:

      M₀ = I   and  M_k+1 = M_k.A - trace(M_k.A) I / (k+1)     k = 0 , 1 , ....      ( the trace of a matrix is the sum of its diagonal elements )
      and we have:       c_k  = - trace(M_k-1.A)/k

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

Data Registers:

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

   ---------

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

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

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

 
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"   >>>>  c₃ = 5 = R03   ( in 31 seconds )
   and  R01 = -11  ,  R02 = -25  ,  R03 = 5     Thus p(x) = x³ - 11 x² - 25 x + 5
 

Notes:

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

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

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

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

     b) Program#2  ( HP-41CX )
 

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

Data Registers:

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

   ---------

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

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

-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.
 
 

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

 
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"   >>>>  c₃ = 5 = R03   ( in 37 seconds )
   and  R01 = -11  ,  R02 = -25  ,  R03 = 5     whence    p(x) = x³ - 11 x² - 25 x + 5
 

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

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

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

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

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

     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 , A²V , ........... , AⁿV
-Then, "PMIN" seeks the smallest k such that  V , AV , A²V , ........... , A^kV  are linearly dependent.

-A relation  a₀ V + a₁ AV + ........... + a_k-1 A^k-1V + A^kV = 0  is found, whence p(x)
 

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

                                        •   R01 = a₁₁    ......................    •  Rn²-n+1 = a_1n
                                           ...................................................................                          ( the n² elements of A )

                                        •   Rnn = a_n1    ......................    •  Rn² = a_nn                            Rn²+1 to R2n²+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" )
 
 

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

 
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) = x² - 5 x + 4       ( the characteristic polynomial is  P(x) = x³ - 9 x² + 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" )
 
 

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

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

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

Data Registers:

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

   ---------

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

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

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

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

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

     b) Program#2
 

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

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

 
-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  Rn²+n+1 ...... Rn²+2n  ( in this example,  R13  R14  R15 )
-Of course, the last decimals may differ slightly.
 

3°) Deflation
 

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

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

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

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

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

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

   ... etc ...
 

Data Registers:

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

   ---------

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

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

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

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

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

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

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

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

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

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

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

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

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

Warning:

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

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

4°) The Jacobi's method
 

     a) Symmetric ( and some non-symmetric ) Matrices
 

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

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

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

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

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

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

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

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

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

*** When the following program stops:

   the n eigenvalues are stored in registers R01  thru Rnn
  the first eigenvector is stored in registers Rn+1 thru R2n
   ..................................................................................
   the n^th eigenvector is stored in registers Rn²+1 thru Rn²+n  as shown below.
 
 

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

   ---------
 

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

  ---------

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

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

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

-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 | a_ij |   ( 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
 
 

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

 
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 )

                 k₁ =  12.81993499  in R01
                 k₂ =  4.910741214  in R02
                 k₃ = -0.730676199  in R03

and the 3 corresponding eigenvectors:

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

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"  >>>>    k₁ = 12.90692994    ( in 2mn44s )  and we have

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

and the 3 corresponding eigenvectors:

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

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

Notes:
 

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

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

-Execution time  = 13 seconds/iteration   if n = 3
-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 O_k^-1.A.O_k  if the matrix is non-symmetric!
-We could use the formula of paragraph 5 below. Unfortunately, it involves complex arithmetic if the radicand is negative.
-But we can replace this radicand by 0 in this case!

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

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

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

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:        k₁ =  16.97583168  ,   k₂ =  6.365547530  ,   k₃ = -3.301311094  ,   k₄ = -5.040068160

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

-These results are obtained in 4mn22s.
 

     b) Symmetric  Matrices only
 

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

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

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

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

            When the program stops:
 

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

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

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
 
 

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

 
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 )

                  k₁ =  12.81993499  in R01
                  k₂ =  4.910741213  in R02
                  k₃ = -0.730676198  in R03

and the 3 corresponding eigenvectors:

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

Notes:

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

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

-Execution time  = 10 seconds/iteration   if n = 3
-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   ...U_k^-1.U_k-1^-1....U₂^-1.U₁^-1.A.U₁.U₂....U_k-1.U_k....
 where the U_k are unitary matrices. The eigenvectors are the columns of  U = U₁.U₂....U_k-1.U_k.... if A is Hermitian ( i-e if A equals its transconjugate )
 Actually, T is diagonal if A is Hermitian.

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

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

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

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

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

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

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

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

Data Registers:

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

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

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

  ---------

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

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

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

-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
 
 

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

  where  k₁ = x₁ + y₁  = 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

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

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

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

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

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

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

Notes:

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

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

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

References:

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