/******************************************************** * 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" * * [BIBLI 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 ) * * * * C++ version by J-P Moreau, Paris * * (with dynamically allocated vectors and matrices). * * (www.jpmoreau.fr) * ********************************************************/ #include #include int it,k,n; LONG_REAL eps; LONG_REAL **A, **Tau; LONG_REAL *D,*E,*P; void *vmblock = NULL; //read symmetric matrix A from screen void Read_data(int n) { int i,j; for (i=1; i<=n; i++) { printf("\n Line %d\n", i); for (j=i; j<=n; j++) { printf(" Element %d: ", j); scanf("%lf",&A[i][j]); A[j][i] = A[i][j]; } } } /******************************************************** * 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. * * * * The pointers to a real vector or a real matrix must * * be declared and allocated in the calling program by * * using the functions of module VMBLOCK.CPP. * ********************************************************/ void TDSL(double eps,int n, LONG_REAL **A, int *it,LONG_REAL *D, LONG_REAL *E, LONG_REAL **Tau) { int i,j,k; LONG_REAL s,u0,v0; LONG_REAL *W; void *vmblock1 = NULL; //initialize local vector W(n+1) vmblock1 = vminit(); W = (LONG_REAL *) vmalloc(vmblock1, VEKTOR, n+1, 0); Tau[1][1]=1.0; for (i=2; i<=n; i++) Tau[i][1]=0.0; u0=1.0; *it=1; k=1; do { for (i=1; i<=n; i++) { s=0.0; for (j=1; j<=n; j++) s += A[i][j]*Tau[j][k]; W[i]=s; } s=0.0; for (i=1; i<=n; i++) s += Tau[i][k]*W[i]; D[k]=s/u0; if (k!=n) { for (i=1; i<=n; i++) { if (k==1) 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]; printf(" i=%d k+1=%d Tau=%f\n", i, k+1, Tau[i][k+1]); } v0=0.0; for (i=1; i<=n; i++) v0 += (Tau[i][k+1]*Tau[i][k+1]); if (v02; i--) { T[i]=P[i]; P[i]=-T[i]+D[k]*P[i-1]-B[k-1]*T[i-2]; } T[2]=P[2]; P[2]=-T[2]+D[k]*P[1]; T[1]=P[1]; P[1]=-T[1]; } free(vmblock1); } void main() { //read size n of matrix A(n,n) printf("\n Input matrix size: "); scanf("%d", &n); //allocate matrices and vectors vmblock = vminit(); A = (LONG_REAL **)vmalloc(vmblock, MATRIX, n+1, n+1); Tau = (LONG_REAL **)vmalloc(vmblock, MATRIX, n+1, n+1); D = (LONG_REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); E = (LONG_REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); P = (LONG_REAL *) vmalloc(vmblock, VEKTOR, n+2, 0); //read A matrix from screen Read_data(n); //fix precision eps=1e-10; //tridiagonalize symmetric A matrix by Lanczos's method TDSL(eps,n,A,&it,D,E,Tau); //if success, use method for a tridiagonal matrix if (it==1) PCTR(n,E,D,P); //print results switch(it) { case 0: printf("\n Method failed.\n"); break; case 1: printf("\n Characteristic polynomial P(l):\n"); for (k=1; k