/********************************************************************* * Eigenvalues and eigenvectors of a real symmetric * * square matrix by JACOBI'S METHOD * * ------------------------------------------------------------------ * * Uses function Jacobi of module ujacobi.cpp. * * * * C++ version by J-P Moreau * * (www.jpmoreau.fr) * * ------------------------------------------------------------------ * * SAMPLE RUN: * * Input file (elpro4.dat): * * * * 4 * * 4 -2 -1 0 * * 4 0 -1 * * 4 -2 * * 4 * * * * Output file (elpro.lst): * * * * --------------------------------------------------- * * Eigenvalues and eigenvectors of a real, symmetric * * square matrix by Jacobi's iterative method * * (with dynamic allocations) * * --------------------------------------------------- * * * * Dimension of matrix: 4 * * * * Symmetric matrix is: * * * * 4.000000 -2.000000 -1.000000 0.000000 * * -2.000000 4.000000 0.000000 -1.000000 * * -1.000000 0.000000 4.000000 -2.000000 * * 0.000000 -1.000000 -2.000000 4.000000 * * * * Eigenvalues: * * 7.00000000 * * 1.00000000 * * 3.00000000 * * 5.00000000 * * * * Eigenvectors (in lines): * * 1.000000 -1.000000 -1.000000 1.000000 * * 1.000000 1.000000 1.000000 1.000000 * * 1.000000 1.000000 -1.000000 -1.000000 * * -1.000000 1.000000 -1.000000 1.000000 * * * * Number of rotations: 25 * * --------------------------------------------------- * * * ********************************************************************** Version with dynamic allocations */ #include "ujacobi.cpp" //for function Jacobi() void main() { int dimen; // Dimension of square matrix int nrot; // Number of matrix rotations int i,j; // loop variables double vmax; // normalization value REAL **Mat; // given square symmetric matrix REAL *Eigenvalues; // eigenvalues solution vector REAL **Eigenvectors; // eigenvectors solution matrix FILE *fp1, *fp2; void *vmblock = NULL; // open input/output files fp1 =fopen("elpro4.dat","r"); // read mode fp2 =fopen("elpro.lst","w"); // write mode //read size of matrix from input file fscanf(fp1,"%d", &dimen); //allocate vector and matrices vmblock = vminit(); Mat = (REAL **) vmalloc(vmblock, MATRIX, dimen+1, dimen+1); Eigenvalues = (REAL *) vmalloc(vmblock, VEKTOR, dimen+1, 0); Eigenvectors = (REAL **) vmalloc(vmblock, MATRIX, dimen+1, dimen+1); //read data from input file for (i=1; i<=dimen; i++) // read only upper half triangle for (j=i; j<=dimen; j++) { fscanf(fp1,"%lf", &Mat[i][j]); Mat[j][i]=Mat[i][j]; } fclose(fp1); //print title and input matrix fprintf(fp2,"---------------------------------------------------\n"); fprintf(fp2," Eigenvalues and eigenvectors of a real, symmetric\n"); fprintf(fp2," square matrix by Jacobi's iterative method\n"); fprintf(fp2," (with dynamic allocations)\n"); fprintf(fp2,"---------------------------------------------------\n\n"); fprintf(fp2," Dimension of matrix: %d\n\n", dimen); fprintf(fp2," Symmetric matrix is:\n"); for (i=1; i<=dimen; i++) { for (j=1; j<=dimen; j++) fprintf(fp2,"%10.6f", Mat[i][j]); fprintf(fp2,"\n"); } Jacobi(Mat,dimen,Eigenvalues,Eigenvectors,&nrot); //print results to output file fprintf(fp2,"\n\n Eigenvalues:\n"); for (i=1; i<=dimen; i++) fprintf(fp2,"%14.8f\n", Eigenvalues[i]); fprintf(fp2,"\n\n Eigenvectors (in lines):\n"); for (j=1; j<=dimen; j++) { vmax=Eigenvectors[1][j]; for (i=2; i<=dimen; i++) if(fabs(Eigenvectors[i][j])>fabs(vmax)) vmax=Eigenvectors[i][j]; for (i=1; i<=dimen; i++) fprintf(fp2,"%10.6f", Eigenvectors[i][j]/vmax); fprintf(fp2,"\n"); } fprintf(fp2,"\n\n Number of rotations: %d\n", nrot); fprintf(fp2,"---------------------------------------------------\n\n"); fclose(fp2); printf("\n Program done.\n"); printf(" Results in file elpro.lst.\n\n"); vmfree(vmblock); } // end of file tujacobi.cpp