/*************************************************** * Calculate the eigenvalues and the eigenvectors * * of a real symmetric tridiagonal matrix * * ------------------------------------------------ * * Ref.: "NUMERICAL RECIPES, Cambridge University * * Press, 1986" [BIBLI 08]. * * * * SAMPLE RUN: * * (find eigenvalues and eigenvectors of matrix: * * 1 2 0 0 0 * * 2 4 7 0 0 * * A = 0 7 10 8 0 * * 0 0 8 -0.75 -9 * * 0 0 0 -9 10 ) * * * * Eigenvalues 1: -0.78071442 * * Eigenvector: * * 1.000000 -0.890357 0.322363 0.344649 0.287721 * * * * Eigenvalues 2: 2.53046878 * * Eigenvector: * * 1.000000 0.765234 -0.446362 -0.252815 -0.304616 * * * * Eigenvalues 3: -9.15659229 * * Eigenvector: * *-0.056408 0.286458 -0.522285 1.000000 0.469812 * * * * Eigenvalues 4: 12.53989066 * * Eigenvector: * * 0.097161 0.560616 0.656182 -0.282210 1.000000 * * * * Eigenvalues 5: 19.11694726 * * Eigenvector: * * 0.051876 0.469920 1.000000 0.728439 -0.719095 * * * * C++ version by J-P Moreau, Paris * * (www.jpmoreau.fr) * * (with dynamic allocations using module vmblock) * ***************************************************/ #include #include LONG_REAL *D, *E; LONG_REAL **Z; int i,j,n; LONG_REAL max; void *vmblock = NULL; LONG_REAL Sign(LONG_REAL a,LONG_REAL b) { if (b>=0) return fabs(a); else return -fabs(a); } //read vectors D and E from screen void Read_data(int n, LONG_REAL *D, LONG_REAL *E) { int i; printf("\n Input %d elements of subdiagonal:\n",n-1); E[1]=1.0; //arbitrary value for (i=1; i1) { for (i=2; i<=n; i++) E[i-1]=E[i]; E[n]=0.0; for (l=1; l<=n; l++) { iter=0; l1: for (m=l; m>n; m--) { dd=fabs(D[m])+fabs(D[m+1]); if (fabs(E[m]-dd)<1e-8) goto l2; } m=n; l2: if (m!=l) { if (iter==30) { printf("\n 30 iterations.\n"); return; } iter++; g=(D[l+1]-D[l])/(2.0*E[l]); r=sqrt(g*g+1.0); g=D[m]-D[l]+E[l]/(g+Sign(r,g)); s=1.0; c=1.0; p=0.0; for (i=m-1; i>l-1; i--) { f=s*E[i]; b=c*E[i]; if (fabs(f)>=fabs(g)) { c=g/f; r=sqrt(c*c+1.0); E[i+1]=f*r; s=1.0/r; c=c*s; } else { s=f/g; r=sqrt(s*s+1.0); E[i+1]=g*r; c=1.0/r; s=s*c; } g=D[i+1]-p; r=(D[i]-g)*s+2.0*c*b; p=s*r; D[i+1]=g+p; g=c*r-b; for (k=1; k<=n; k++) { f=Z[k][i+1]; Z[k][i+1]=s*Z[k][i]+c*f; Z[k][i]=c*Z[k][i]-s*f; } } D[l]=D[l]-p; E[l]=g; E[m]=0.0; goto l1; } } //of l loop if (iter<30) printf(" %d iterations.\n", iter); } } int main() { //read size of matrix printf("\n Input size of matrix: "); scanf("%d", &n); //allocate vectors and matrix (index 0 not used here) //see file vmblock.cpp vmblock = vminit(); D = (LONG_REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); E = (LONG_REAL *) vmalloc(vmblock, VEKTOR, n+1, 0); Z = (LONG_REAL **)vmalloc(vmblock, MATRIX, n+1, n+1); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return 1; } Read_data(n,D,E); //initialize Z to identity matrix for (i=1; i<=n; i++) for (j=1; j<=n; j++) if (i!=j) Z[i][j] = 0.0; else Z[i][j] = 1.0; //calculate eigenvalues and eigenvectors TQLI(D,E,n,Z); //print results for (j=1; j<=n; j++) { printf("\n Eigenvalue %d: %f\n",j, D[j]); printf(" Eigenvector:\n"); //normalize eigenvector before printing max=fabs(Z[1][j]); for (i=2; i<=n; i++) if (fabs(Z[i][j])>fabs(max)) max=Z[i][j]; for (i=1; i<=n; i++) printf(" %f", Z[i][j]/max); printf("\n"); } printf("\n\n"); free(vmblock); return 0; } // end of file elprotd.cpp