'********************************************************* '* 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