'*********************************************************
'*   Calculate the coefficients of the characteristic    *
'*   polynomial P(l)=A-l*I of a square symmetric real    *
'*         matrix A(i,j) by Lanczos's method             *
'* ----------------------------------------------------- *
'* Ref.: "Algèbre, Algorithmes et programmes en Pascal   *
'*        By Jean-Louis Jardrin, DUNOD Paris, 1988" [10].*
'* ----------------------------------------------------- *
'* SAMPLE RUN:                                           *
'* (Find the caracteristic polynomial of matrix:         *
'*          10  7   8   7                                *
'*      A =  7  5   6   5                                *
'*           8  6  10   9                                *
'*           7  5   9  10                                *
'*                                                       *
'* Input matrix size (max 50): 4                         *
'*  Line 1:                                              *
'*   Element 1: 10                                       *
'*   Element 2: 7                                        *
'*   Element 3: 8                                        *
'*   Element 4: 7                                        *
'*  Line 2:                                              *
'*   Element 2: 5                                        *
'*   Element 3: 6                                        *
'*   Element 4: 5                                        *
'*  Line 3:                                              *
'*   Element 3: 10                                       *
'*   Element 4: 9                                        *
'*  Line 4:                                              *
'*   Element 4: 10                                       *
'*                                                       *
'* Chracteristic polynomial P(l):                        *
'* Coefficient for degree 4:     1.000000                *
'* Coefficient for degree 3:   -35.000000                *
'* Coefficient for degree 2:   146.000000                *
'* Coefficient for degree 1:  -100.000000                *
'* Coefficient for degree 0:     1.000000                *
'*                                                       *
'* (The characteristic polynomial of matrix A is:        *
'*   P(l)=l^4-35*l^3+146*l^2-100*l+1 )                   *
'*                                                       *
'*                Basic version by J-P Moreau, Paris     *
'*                         (www.jpmoreau.fr)             *
'*********************************************************
'Program CARPOL3
defdbl a-h,o-z
defint i-n

NMAX=20

dim A(NMAX,NMAX), Tau(NMAX,NMAX)
dim D(NMAX), E(NMAX), P(NMAX+1)

  cls
' read size n of matrix A(n,n)
  print
  input " Input matrix size: ", n
' read A matrix
  gosub 1000  'call Read_data
' fix precision
  eps=1e-10
' tridiagonalize symmetric A matrix by Lanczos's method
  gosub 2000  'call TDSL(eps,n,A,it,D,E,TAU)
' if success, use method for a tridiagonal matrix
  if it=1 then gosub 3000  'call PCTR(n,E,D,P)
' print results
  print
  if it=0 then
    print " Method failed."
  else
    print " Characteristic polynomial P(l):"
    for k=1 to n+1
      print " Coefficient for degree "; n-k+1; ": ";
      print using "#####.######"; P(k)
    next k
  end if
  print
  print

END


1000 'read symmetric matrix A(n,n)
  for i=1 to n
    print
    print " Line "; i
    for j=i to n
      print "  Element "; j; ": "; : input "", A(i,j)
      if j<>i then  A(j,i)=A(i,j)
    next j
  next i
return

'*********************************************************
'* The TDSL procedure tridiagonalizes the symmetric real *
'* matrix A(n,n) by Lanczos's method. Only the main dia- *
'* gonal terms and the supradiagonal terms are calculated*
'* ----------------------------------------------------- *
'* Inputs:                                               *
'*          eps  precision (double)                      *
'*          n    size of square A matrix (integer)       *
'*          A    square real symmetric matrix (n,n)      *
'* Outputs:                                              *
'*          it   flag=0 if method failed, else =1        *
'*          D    vector of size n+1 storing the n        *
'*               elements of the main diagonal.          *
'*          E    vector of size n+1 storing the n-1      *
'*               elements of the supradiagonal.          *
'*          Tau  Utility matrix used by the method.      *
'*                                                       *
'*********************************************************
2000 'Subroutine TDSL(eps,n,A,it,D,E,Tau)
dim W(NMAX)
  Tau(1,1)=1#
  for i=2 to n
    Tau(i,1)=0#
  next i
  u0=1# : it=1 : k=1
10  for i=1 to n
      s=0#
      for j=1 to n
        s=s+A(i,j)*Tau(j,k)
      next j
      W(i)=s
    next i
    s=0#
    for i=1 to n
      s=s+Tau(i,k)*W(i)
    next i
    D(k)=s/u0
    if k<>n then
      for i=1 to n
        if k=1 then
          Tau(i,k+1)=W(i)-D(k)*Tau(i,k)
        else
          Tau(i,k+1)=W(i)-D(k)*Tau(i,k)-E(k-1)*Tau(i,k-1)
        end if
      next i
      v0=0#
      for i=1 to n
        v0=v0+Tau(i,k+1)*Tau(i,k+1)
      next i
      if  v0<eps then
        it=0
      else
        E(k)=v0/u0 : u0=v0
      end if
    end if
    k=k+1
  if k<=n and it=1 then goto 10
  return
End

'*****************************************************
'* The procedure PCTR calculates the coefficients of *
'* the characteristic polynomial P(l)=A-l*I of a real*
'* square tridiagonal matrix A(i,j).                 *
'*                                                   *
'* Note: the roots of P(l) are the eigenvalues of    *
'*       matrix A(i,j).                              *
'* ------------------------------------------------- *
'* INPUTS:                                           *
'*          n: size of matrix (integer)              *
'*          B: vector of codiagonal terms products   *
'*             (see main program).                   *
'*          D: vector of main diagonal terms         *
'* OUTPUT:                                           *
'*          P: vector of coefficients of P(l)        *
'* ------------------------------------------------- *
'* (For demo of procedure PCTR alone, see program    *
'*  carpol.bas).                                     *
'*****************************************************
3000 'Subroutine PCTR(n,E,D,P)
dim T(NMAX)
  P(2)=D(1) : T(1)=1# : P(1)=-T(1)
  for k=2 to n
    P(k+1)=D(k)*P(k)-E(k-1)*T(k-1)
    for i=k to 3 step -1
      T(i)=P(i)
      P(i)=-T(i)+D(k)*P(i-1)-E(k-1)*T(i-2)
    next i
    T(2)=P(2) : P(2)=-T(2)+D(k)*P(1)
    T(1)=P(1) : P(1)=-T(1)
  next k
return

'end of file carpol3.bas