/********************************************************************** * Calculate eigenvalues and eigenvectors of a Square Hermitian * * Matrix By Jacobi's Method * * ------------------------------------------------------------------- * * SAMPLE RUN: * * File Matc4.dat contains: * * 4 * * 12 0 -1 -1 2 0 3 3 * * -1 1 12 0 1 -1 -2 0 * * 2 0 1 1 8 0 -1 -1 * * 3 -3 -2 0 -1 1 8 0 * * * * Precision: 1e-10 * * Max. number of oterations: 25 * * * * Eigenvalue 1: 4.000000 * * Eigenvector: * * -0.500000 - 0.500000 I * * 0.000000 + 0.000000 I * * 0.500000 + 0.500000 I * * 1.000000 + 0.000000 I * * Eigenvalue 2: 8.000000 * * Eigenvector: * * -0.000000 + 0.000000 I * * -0.500000 + 0.500000 I * * 1.000000 + 0.000000 I * * -0.500000 + 0.500000 I * * Eigenvalue 3: 12.000000 * * Eigenvector: * * 0.500000 + 0.500000 I * * 1.000000 + 0.000000 I * * 0.500000 + 0.500000 I * * -0.000000 + 0.000000 I * * Eigenvalue 4: 16.000000 * * Eigenvector: * * 1.000000 + 0.000000 I * * -0.500000 + 0.500000 I * * -0.000000 + 0.000000 I * * 0.500000 - 0.500000 I * * * * ------------------------------------------------------------------- * * Reference: "ALGEBRE Algorithmes et programmes en Pascal * * By Jean-Louis Jardrin - Dunod BO-PRE 1988" [BIBLI 10] * * * * C++ Release 1.1 By J-P Moreau, Paris. * * (www.jpmoreau.fr) * * * * Release 1.1 (Nov. 30th, 2007): Added function NORMAL. * **********************************************************************/ #include #include #define SIZE 25 //complex number typedef struct { double R,I; //algebraic form } COMPLEX; typedef COMPLEX CMat[SIZE][SIZE]; typedef double Vec[SIZE]; int I,it,J,M,N; double dta; CMat A,VX; Vec R; FILE *fp; // input text file double Sqr(double a) { return a*a; } //fabsolute value of a complex number double CABS(COMPLEX C) { return sqrt(Sqr(C.R)+Sqr(C.I)); } //Add two complex numbers void CADD(COMPLEX c1, COMPLEX c2, COMPLEX *c3) { c3->R=c1.R+c2.R; c3->I=c1.I+c2.I; } //Substract two complex numbers void CDIF(COMPLEX c1, COMPLEX c2, COMPLEX *c3) { c3->R=c1.R-c2.R; c3->I=c1.I-c2.I; } //Multiply two complex numbers void CMUL(COMPLEX c1, COMPLEX c2, COMPLEX *c3) { c3->R=c1.R*c2.R - c1.I*c2.I; c3->I=c1.R*c2.I + c1.I*c2.R; } //Return conjugate of a complex number void CONJ(COMPLEX c, COMPLEX *c1) { c1->R=c.R; c1->I=-c.I; } //Multiply a complex number by a real number void CPRO(double alpha, COMPLEX c, COMPLEX *c1) { c1->R=alpha*c.R; c1->I=alpha*c.I; } /*********************************************************************** * Compute all eigenvalues/eigenvectors of a square hermitian matrix * * using Jacobi's method. * * -------------------------------------------------------------------- * * Inputs: * * dta: required precision (double) * * M : maximum number of iterations (integer) * * N : size of matrix * * A : pointer to hermitian (complex) matrix * * Outputs: * * it : flag for convergence (integer) =-1, no convergence * * R : vector(1..N) of eigenvalues (double) * * VX : matrix storing complex eigenvectors in columns * ***********************************************************************/ void EPHJ(double dta, int M, int N, CMat A, int *it, Vec R, CMat VX) { int I,J,K,K0,L,L0; double delta,s,s0,t0,t1,w0; COMPLEX c0,c1,c2,c3,u0,u1,z0,z1; z0.R=0.0; z0.I=0.0; z1.R=1.0; z1.I=0.0; for (I=1; I<=N; I++) for (J=1; J<=N; J++) if (I==J) VX[I][J]=z1; else VX[I][J]=z0; *it=-1; L=1; while (L<=M && *it!=1) { s=0.0; for (I=1; I s) { s=t0; K0=I; L0=J; } } if (s==0.0) *it=1; else { delta=Sqr(A[L0][L0].R-A[K0][K0].R)+4.0*Sqr(CABS(A[K0][L0])); t0=A[L0][L0].R-A[K0][K0].R + sqrt(delta); t1=A[L0][L0].R-A[K0][K0].R - sqrt(delta); if (fabs(t0) >= fabs(t1)) w0=t0; else w0=t1; s0=fabs(w0)/sqrt(Sqr(w0)+4.0*Sqr(CABS(A[K0][L0]))); t0=2.0*s0/w0; CPRO(t0,A[K0][L0],&c0); CONJ(c0,&c1); for (I=1; I fabs(YI)) { W=YI/YR; YR=W*YI+YR; *ZR=(AR+W*AI)/YR; *ZI=(AI-W*AR)/YR; } else { W=YR/YI; YI=W*YR+YI; *ZR=(W*AR+AI)/YI; *ZI=(W*AI-AR)/YI; } } void NORMAL(int N, Vec R, CMat Z) { /*------------------------------------------------------------- Sort eigenvalues by fabsolute values in ascending order and normalize. -------------------------------------------------------------*/ // Label e5 Vec Tr,Ti; double ZM,Z1,Z2,VR; COMPLEX V; int IM,J,K; // SORT SOLUTIONS IN ASCENDING ORDER for (J=2; J<=N; J++) { VR=R[J]; for (K=1; K<=N; K++) { Tr[K]=Z[K][J].R; Ti[K]=Z[K][J].I; } for (I=J-1; I>0; I--) { if (fabs(R[I]) <= fabs(VR)) goto e5; R[I+1]=R[I]; for (K=1; K<=N; K++) { Z[K][I+1].R=Z[K][I].R; Z[K][I+1].I=Z[K][I].I; } } I=0; e5: R[I+1]=VR; for (K=1; K<=N; K++) { Z[K][I+1].R=Tr[K]; Z[K][I+1].I=Ti[K]; } } // NORMALIZE WITH RESPECT TO BIGGEST ELEMENT for (J=N; J>0; J--) { ZM = 0.0; for (I=1; I<=N; I++) { V.R=Z[I][J].R; V.I=Z[I][J].I; Z1=CABS(V); if (Z1 >= ZM) { IM = I; ZM = Z1; } } Z1 = Z[IM][J].R; Z2 = Z[IM][J].I; for (I=1; I<=N; I++) { CDIV(Z[I][J].R,Z[I][J].I,Z1,Z2,&V.R,&V.I); Z[I][J].R = V.R; Z[I][J].I = V.I; } } } // NORMAL() // Read from input file complex hermitian matrix void LME() { int I,J; fp = fopen("matc4.dat","r"); fscanf(fp,"%d", &N); for (I=1; I<=N; I++) for (J=1; J<=N; J++) fscanf(fp, "%lf %lf", &A[I][J].R, &A[I][J].I); fclose(fp); } void main() { LME(); // read complex matrix from text file printf("\n ** Compute Eigenvalues and Eigenvectors **\n"); printf(" of a Square Hermitian Matrix\n"); printf(" By Jacobi's Method\n\n"); printf(" Precision: "); scanf("%lf", &dta); printf(" Max. number of iterations: "); scanf("%d", &M); printf("\n"); EPHJ(dta, M, N, A, &it, R, VX); NORMAL(N, R, VX); if (it==-1) printf(" No convergence.\n"); else // display results for (J=1; J<=N; J++) { printf(" Eigenvalue %d: %f\n", J, R[J]); printf(" Eigenvector:\n"); for (I=1; I<=N; I++) { printf(" %f ", VX[I][J].R); if (VX[I][J].I >= 0.0) printf(" + "); else printf(" - "); printf("%f I\n", fabs(VX[I][J].I)); } } } // end of file Tephj.cpp