/********************************************************************** * Calculation of eigenvalues and eigenvectors of a real square matrix * * by the QR method * * ------------------------------------------------------------------- * * SAMPLE RUN: * * * * As an example, the input file (hqr5.dat) has the form: * * 5 * * 1.0 2.0 3.0 -7.0 12.0 * * 2.0 4.0 7.0 3.0 -1.0 * * 3.0 7.0 10.0 8.0 4.0 * * -7.0 3.0 8.0 -0.75 -9.0 * * 12.0 -1.0 4.0 -9.0 10.0 * * * * The first number (5) is the dimension of the matrix followed by * * the 1st to 5th row of the 5 x 5 matrix (in our case). * * * * The output file (hqr5.lst) should contain: * * * * ------------------------------------------------------------------- * * Eigenvalues and Eigenvectors by QR algorithm * * ------------------------------------------------------------------- * * * * Dimension of the input matrix = 5 * * * * Input matrix: * * * * 1.000000 2.000000 3.000000 -7.000000 12.000000 * * 2.000000 4.000000 7.000000 3.000000 -1.000000 * * 3.000000 7.000000 10.000000 8.000000 4.000000 * * -7.000000 3.000000 8.000000 -0.750000 -9.000000 * * 12.000000 -1.000000 4.000000 -9.000000 10.000000 * * * * * * Normalized Eigenvectors (in columns): * * * * 0.467172 1.000000 0.705820 0.093053 0.208812 * * -0.007947 -0.326100 -0.012677 0.594911 1.000000 * * -0.507827 0.209620 0.131486 1.000000 -0.478027 * * 1.000000 -0.147689 -0.527499 0.455666 -0.275000 * * 0.264432 -0.815422 1.000000 0.050740 -0.216915 * * * * * * Eigenvalues: Iterations: * * * * -1.0486545165515601e+001 + 0.0000000000000000e+000 * i 0 * * -7.7745797294783081e+000 + 0.0000000000000000e+000 * i 4 * * 2.3755954765326987e+001 + 0.0000000000000000e+000 * i -4 * * 1.8291820601685462e+001 + 0.0000000000000000e+000 * i 5 * * 4.6334952798144985e-001 + 0.0000000000000000e+000 * i -5 * * * * Check sum = 2.230161e-013 (must be approximately 0) * * * * ------------------------------------------------------------------- * * Reference: "Numerical Algorithms with C by G. Engeln-Mueller and * * F. Uhlig, Springer-Verlag, 1996" * * * * C++ version with dynamic allocations * * by Jean-Pierre Moreau, Paris. * * (www.jpmoreau.fr) * **********************************************************************/ //uses modules: basis_r.cpp, vmblock.cpp, feigen.cpp, // basis.h, vmblock.h //--------------------------------------------------------------------- #include //for general declarations and macros #include //for dynamic allocations FILE *fp_in; int eigen (int, int, int, int, REAL **, REAL **, REAL *, REAL *, int *); //input data from file void Read_Mat(int n,REAL **mat) { int i,j; for (i=0; i 0 *wr, // Eigenvalues (Real part) *wi; // Eigenvalues (Imaginary parts) int n, *cnt, // Iteration counter rc, // Return Code vec = 1, // flag for eigenvectors (=0 -> none) ortho = 0, // flag for orthogonal // Hessenberg reduction ev_norm = 1; // flag for normalization of Eigenvectors register i, j, k; REAL v, w, norm; void *vmblock; FILE *fp_out; fp_in = fopen("mat5.dat","r"); fp_out = fopen("hqr5.lst","w"); WriteHead (fp_out," Eigenvalues and Eigenvectors by QR algorithm"); if (fscanf (fp_in,"%d", &n) <= 0) { LogError ("Input stream", 0, __FILE__, __LINE__); return 1; } if (n < 1) { LogError ("Dimension must be > 0", 0, __FILE__, __LINE__); return 1; } // Allocate Memory ................................................. vmblock = vminit(); mat = (REAL **)vmalloc(vmblock, MATRIX, n, n); a = (REAL **)vmalloc(vmblock, MATRIX, n, n); wr = (REAL *) vmalloc(vmblock, VEKTOR, n, 0); wi = (REAL *) vmalloc(vmblock, VEKTOR, n, 0); cnt = (int *) vmalloc(vmblock, VVEKTOR, n, sizeof(*cnt)); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return 1; } if (vec) // Only for eigenvectors ..... { ev = (REAL **)vmalloc(vmblock, MATRIX, n, n); if (! vmcomplete(vmblock)) { LogError ("No Memory", 0, __FILE__, __LINE__); return 1; } } Read_Mat(n, mat); fclose(fp_in); fprintf (fp_out,"Dimension of the input matrix = %d\n\n", n); fprintf (fp_out,"Input matrix:\n"); WriteMat (fp_out, n, n, mat); CopyMat (n, n, mat, a); fprintf(fp_out,"\n"); rc = eigen (vec, ortho, ev_norm, n, mat, ev, wr, wi, cnt); if (rc != 0) // ERROR !!! { LogError ("eigen", rc, __FILE__, __LINE__); for (i = 0; i < n; i++) printf (" %d ", cnt[i]); return 1; } // If vec != 0, print eigenvectors if (vec) { if (ev_norm) fprintf(fp_out,"Normalized Eigenvectors (in columns):\n"); else fprintf(fp_out,"not normalized Eigenvectors:\n"); WriteMat (fp_out, n, n, ev); } fprintf (fp_out,"\nEigenvalues: Iterations:\n\n"); for (i = 0; i < n; i++) { fprintf (fp_out, FORMAT_2016LE, wr[i]); fprintf (fp_out, "+ "); fprintf (fp_out, FORMAT_2016LE, wi[i]); fprintf (fp_out, "* i\t%5d\n", cnt[i]); } // Check result: sum of L1 norms of // Matrix*Eigenvector - Eigenvalue*Eigenvector if (vec) // (this must be nearly 0). { for (norm = ZERO, k = 0; k < n; k++) { if (wi[k] == ZERO) { for (i = 0; i < n; i++) { for (w = ZERO, j = 0; j < n; j++) w += a[i][j] * ev[j][k]; w -= wr[k] * ev[i][k]; norm += ABS (w); } } else { for (i = 0; i < n; i++) { for (w = ZERO, j = 0; j < n; j++) w += a[i][j] * ev[j][k]; w -= wr[k] * ev[i][k] - wi[k] * ev[i][k+1]; for (v = ZERO, j = 0; j < n; j++) v += a[i][j] * ev[j][k+1]; v -= wr[k] * ev[i][k+1] + wi[k] * ev[i][k]; norm += 2.0 * SQRT (v*v + w*w); } k++; } } fprintf (fp_out,"\nCheck sum = "); fprintf (fp_out,FORMAT_LE, norm); fprintf (fp_out,"(must be approximately 0)\n\n"); } WriteEnd (fp_out); fclose(fp_out); printf("\n Results in hqr5.lst.\n\n"); vmfree(vmblock); return (0); } // test1_hqr.cpp