/******************************************************************** * THIS PROGRAM CALCULATES THE EIGENVALUES/EIGENVECTORS OF A COMPLEX * * MATRIX C = A+I*Z BY THE JACOBI METHOD, THIS METHOD IS ALSO RECOM- * * MANDED TO DEAL WITH REAL MATRICES THE EIGENVALUES OF WHICH ARE * * COMPLEX. * * ----------------------------------------------------------------- * * SAMPLE RUN: * * Example #1: * * (Find all eigenvalues and eigenvectors of 5x5 matrix: * * 1.0 2.0 3.0 -7.0 12.0 * * 2.0 4.0 7.0 3.0 -10.0 * * 3.0 7.0 10.0 8.0 4.0 * * -7.0 3.0 8.0 -0.75 -9.0 * * 12.0 -10.0 4.0 -9.0 10.0 ) * * * * Matrix C (real parts): * * 1.000000 2.000000 3.000000 -7.000000 12.000000 * * 2.000000 4.000000 7.000000 3.000000-10.000000 * * 3.000000 7.000000 10.000000 8.000000 4.000000 * * -7.000000 3.000000 8.000000 -0.750000 -9.000000 * * 12.000000-10.000000 4.000000 -9.000000 10.000000 * * Matrix C (imaginary parts): * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * * * Eigenvalues (real parts): * * 3.349223 -10.152415 -12.544478 17.027498 26.570173 * * Eigenvalues (imaginary parts): * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * Error code (must be zero): 0 * * * * Eigenvectors in columns (real parts): * * 0.671571 0.935133 -0.563486 0.360472 0.568860 * * 1.000000 -0.231511 0.802269 0.418719 -0.516435 * * -0.432515 -0.427944 -0.567937 1.000000 -0.152070 * * -0.622247 1.000000 0.609801 0.229197 -0.576420 * * -0.290043 -0.140178 1.000000 0.295367 1.000000 * * Eigenvectors in columns (imaginary parts): * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * 0.000000 0.000000 0.000000 0.000000 0.000000 * * * * Example #2: * * Matrix C (real parts): * * 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 * * Matrix C (imaginary parts): * * 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 1.000000 1.000000 1.000000 1.000000 * * 1.000000 1.000000 1.000000 1.000000 1.000000 * * Eigenvalues (real parts): * * 0.468367 -7.775245-10.354718 18.634772 23.276824 * * Eigenvalues (imaginary parts): * * 0.040316 0.003317 0.942489 3.320560 0.693318 * * Error code (must be zero): 0 * * * * Eigenvectors in columns (real parts): * * 0.209070 1.000000 0.460630 0.199619 0.693336 * * 1.000000 -0.326405 -0.010671 0.590353 -0.101481 * * -0.481799 0.206172 -0.514461 1.000000 0.001810 * * -0.273820 -0.142932 1.000000 0.380689 -0.610896 * * -0.217450 -0.813120 0.265266 0.198153 1.000000 * * Eigenvectors in columns (imaginary parts): * * 0.001998 0.000000 -0.003799 -0.132150 0.023902 * * 0.000000 0.001431 -0.059816 0.069630 0.183347 * * -0.033383 0.016553 -0.054470 0.000000 0.276292 * * 0.009639 -0.022869 0.000000 0.218001 0.166455 * * -0.005084 -0.011084 -0.036704 -0.220804 0.000000 * * ----------------------------------------------------------------- * * REFERENCE: * * From Numath Library By Tuan Dang Trong in Fortran 77 * * [BIBLI 18]. * * * * C++ Release 1.1 By J-P Moreau, Paris. * * (www.jpmoreau.fr) * * * * Release 1.1: added example #2. *. ********************************************************************/ #include #include #define NMAX 25 #define FALSE 0 #define TRUE 1 typedef double MAT[NMAX][NMAX]; //index zero not used here. MAT AR, AI, T, U; double EN[2*NMAX]; //index zero not used here. int i,j, IER, N; char answ[2]; double ABSCD(double AR, double AI) { // ABSOLUTE VALUE OF A COMPLEX NUMBER C=AR+I*AI // ABSCD=SQRT(AR^2+AI^2) double XR,XI,W; XR=fabs(AR); XI=fabs(AI); if (XR <= XI) { W=XR; XR=XI; XI=W; } if (XI == 0.0) return XR; else return (XR*sqrt(1.0+(XI/XR)*(XI/XR))); } void DIVCD(double AR, double AI, double BR, double BI, double *ZR, double *ZI) { // COMPLEX DIVISION Z=ZR+I*ZI=(AR+I*AI)/(BR+I*BI) // DO NOT USE IF BR=BI=0. double YR,YI,W; YR=BR; YI=BI; if (fabs(YR) > 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 COMEIG (int NDIM,int N, int NN, MAT A,MAT Z,MAT T,MAT U, int *IER, double *EN) { /*------------------------------------------------------------------------------------- ! THIS SUBROUTINE CALCULATES THE EIGENVALUES/EIGENVECTORS OF A COMPLEX MATRIX ! C = A+I*Z BY THE JACOBI METHOD, THIS METHOD IS ALSO RECOMMANDED TO DEAL WITH ! REAL MATRICES THE EIGENVALUES OF WHICH ARE COMPLEX. ! DATA: ! NDIM 1ST DIMENSION OF TABLES A, Z, T, U IN MAIN PROGRAM (HERE NDIM=N) ! N REAL SIZE OF COMPLEX MATRIX C = A+I*Z ! NN MAXIMUM NUMBER OF ITERATIONS ! A TABLE STORING THE REAL PART OF GIVEN MATRIX ! Z TABLE STORING THE IMAGINARY PART OF GIVEN MATRIX ! OUTPUTS: ! A(J,J),Z(J,J),J=1,N IN MAIN DIAGONALS OF TABLES A AND Z, YOU ! HAVE NOW RESPECTIVELY THE REAL AND IMAGINARY PARTS OF EIGENVALUES. ! T,U THESE TABLES CONTAIN NOW RESPECTIVELY THE REAL AND IMAGINARY PARTS ! OF THE EIGENVECTORS MATRIX X = T+I*U (STORED IN COLUMNS). ! IER ERROR CODE ! = 0 CONVERGENCE OK ! = 1 NO CONVERGENCE AFTER NN ITERATIONS ! WORKING ZONE: ! EN TABLE OF SIZE 2*N ! NOTES: ! 1/ IN CASE OF CONVERGENCE (IER = 0), THE MATRIX EQUATION C*X = LAMDA*X ! IS VERIFIED TO THE MACHINE PRECISION. ! ! REFERENCE: ! P.J.EBERLEIN, NUMER.MATH 14, PP 232-245 (1970) !-----------------------------------------------------------------------------------*/ // Labels: e10,e20, e65,e70, e90,e95,e97; bool MARK; double B,DN,E,EMAX,EPS,EPS1,TAU,W1,W2; double BR,BI,DR,DI,ER,EI,G,HI,HJ,HR,T1,T2; double C,CA,CX,D,DE,ROOT,ROOT1,ROOT2,S,SA,SB,SH,SIG,SX,SW; double CB,CH,CN,COS2A,COTX,COT2X,ETA,SIN2A,TANH,TI,TR,TSE; double C1I,C2I,C1R,C2R,S1I,S2I,S1R,S2R,VI,VR,ZI,ZM; double AKI,AMI,ZKI,ZMI,AIK,AIM,ZIK,ZIM,TIK,TIM,UIK,UIM; int I,IM,IT,J,K,M; // CHECK MACHINE EPSILON (AROUND 1.2E-16 FOR PC) EPS = 1.0; e10: EPS = 0.5*EPS; EPS1 = EPS+1.0; if (EPS1 > 1.0) goto e10; MARK = FALSE; // INITIALIZE EIGENVECTORS for (I = 1; I<=N; I++) { T[I][I] = 1.0; U[I][I] = 0.0; for (J = I+1; J<=N; J++) { T[I][J] = 0.0; T[J][I] = 0.0; U[I][J] = 0.0; U[J][I] = 0.0; } } IT = 0; e20: IT++; // SAFETY TEST IN CASE OF NO CONVERGENCE if (IT > NN) goto e90; if (MARK) goto e95; // DEFINE CONVERGENCE CRITERIUM TAU = 0.0; for (K = 1; K<=N; K++) { W1 = 0.0; for (I = 1; I<=N; I++) if (I != K) W1 = W1+fabs(A[I][K])+fabs(Z[I][K]); TAU += W1; EN[K] = W1+fabs(A[K][K])+fabs(Z[K][K]); } // PERMUTE LINES AND COLUMNS for (K = 1; K EMAX) { EMAX = EN[J]; I = J; } if (I != K) { EN[I] = EN[K]; for (J = 1; J<=N; J++) { W2 = A[K][J]; A[K][J] = A[I][J]; A[I][J] = W2; W2 = Z[K][J]; Z[K][J] = Z[I][J]; Z[I][J] = W2; } for (J = 1; J<=N; J++) { W2 = A[J][K]; A[J][K] = A[J][I]; A[J][I] = W2; W2 = Z[J][K]; Z[J][K] = Z[J][I]; Z[J][I] = W2; W2 = T[J][K]; T[J][K] = T[J][I]; T[J][I] = W2; W2 = U[J][K]; U[J][K] = U[J][I]; U[J][I] = W2; } } } // CONVERGENCE IF TAU < 100*EPS if (TAU < 100.0*EPS) goto e95; // BEGIN ITERATIONS MARK = TRUE; for (K = 1; K= T2) { SW = 1.0; C = BR; S = EI; D = DR; DE = DI; ROOT2 = sqrt(T1); } else { SW =-1.0; C = BI; S =-ER; D = DI; DE = DR; ROOT2 = sqrt(T2); } ROOT1 = sqrt(S*S+C*C); SIG = 1.0; if (D < 0.0) SIG = -1.0; CA = 1.0; if (C < 0.0) CA = -1.0; SA = 0.0; if (ROOT1 < EPS) { SX = 0.0; SA = 0.0; CX = 1.0; CA = 1.0; if (SW > 0.0) { E = ER; B = BI; } else { E = EI; B =-BR; } DN = D*D+DE*DE; goto e65; } if (fabs(S) > EPS) { CA = C/ROOT1; SA = S/ROOT1; } COT2X = D/ROOT1; COTX = COT2X+SIG*sqrt(1.0+COT2X*COT2X); SX = SIG/sqrt(1.0+COTX*COTX); CX = SX*COTX; ETA = (ER*BR+BI*EI)/ROOT1; TSE = (BR*BI-ER*EI)/ROOT1; T1 = SIG*(TSE*D-ROOT1*DE)/ROOT2; T2 = (D*DE+ROOT1*TSE)/ROOT2; DN = ROOT2*ROOT2+T2*T2; T2 = HJ*CX*SX; COS2A = CA*CA-SA*SA; SIN2A = 2.0*CA*SA; W1 = HR*COS2A+HI*SIN2A; W2 = HI*COS2A-HR*SIN2A; HR = CX*CX*HR-SX*SX*W1-CA*T2; HI = CX*CX*HI+SX*SX*W2-SA*T2; B = SW*T1*CA+ETA*SA; E = CA*ETA-SW*T1*SA; // ROOT1 < EPS e65: S = HR-SIG*ROOT2*E; C = HI-SIG*ROOT2*B; ROOT = sqrt(C*C+S*S); if (ROOT < EPS) { CB = 1.0; CH = 1.0; SB = 0.0; SH = 0.0; goto e70; } CB = -C/ROOT; SB = S/ROOT; T2 = CB*B-E*SB; CN = T2*T2; TANH = ROOT/(G+2.0*(CN+DN)); CH = 1.0/sqrt(1.0-TANH*TANH); SH = CH*TANH; // ROOT < EPS e70: W1 = SX*SH*(SA*CB-SB*CA); C1R = CX*CH-W1; C2R = CX*CH+W1; C1I =-SX*SH*(CA*CB+SA*SB); C2I = C1I; W2 = SX*CH*CA; W1 = CX*SH*SB; S1R = W2-W1; S2R =-W2-W1; W2 = SX*CH*SA; W1 = CX*SH*CB; S1I = W2+W1; S2I = W2-W1; W1 = sqrt(S1R*S1R+S1I*S1I); W2 = sqrt(S2R*S2R+S2I*S2I); if (W1 > EPS || W2 > EPS) { // BEGIN TRANSFORMATIONS MARK = FALSE; for (I = 1; I<=N; I++) { AKI = A[K][I]; AMI = A[M][I]; ZKI = Z[K][I]; ZMI = Z[M][I]; A[K][I] = C1R*AKI-C1I*ZKI+S1R*AMI-S1I*ZMI; Z[K][I] = C1R*ZKI+C1I*AKI+S1R*ZMI+S1I*AMI; A[M][I] = S2R*AKI-S2I*ZKI+C2R*AMI-C2I*ZMI; Z[M][I] = S2R*ZKI+S2I*AKI+C2R*ZMI+C2I*AMI; } for (I = 1; I<=N; I++) { AIK = A[I][K]; AIM = A[I][M]; ZIK = Z[I][K]; ZIM = Z[I][M]; TIK = T[I][K]; TIM = T[I][M]; UIK = U[I][K]; UIM = U[I][M]; A[I][K] = C2R*AIK-C2I*ZIK-S2R*AIM+S2I*ZIM; Z[I][K] = C2R*ZIK+C2I*AIK-S2R*ZIM-S2I*AIM; A[I][M] = C1R*AIM-C1I*ZIM-S1R*AIK+S1I*ZIK; Z[I][M] = C1R*ZIM+C1I*AIM-S1R*ZIK-S1I*AIK; T[I][K] = C2R*TIK-C2I*UIK-S2R*TIM+S2I*UIM; U[I][K] = C2R*UIK+C2I*TIK-S2R*UIM-S2I*TIM; T[I][M] = C1R*TIM-C1I*UIM-S1R*TIK+S1I*UIK; U[I][M] = C1R*UIM+C1I*TIM-S1R*UIK-S1I*TIK; } } // if w1>EPS... // END TRANSFORMATIONS } // end of M loop } // end of K loop // GO TO NEXT ITERATION goto e20; // NO CONVERGENCE ! e90: *IER = 1; return; // CONVERGENCE OK. e95: *IER = 0; // SORT SOLUTIONS IN INCREASING ORDER for (J=2; J<=N; J++) { VR=A[J][J]; VI=Z[J][J]; for (K=1; K<=N; K++) { EN[K]=T[K][J]; EN[K+N]=U[K][J]; } for (I=J-1; I>0; I--) { if (ABSCD(A[I][I],Z[I][I]) <= ABSCD(VR,VI)) goto e97; A[I+1][I+1]=A[I][I]; Z[I+1][I+1]=Z[I][I]; for (K=1; K<=N; K++) { T[K][I+1]=T[K][I]; U[K][I+1]=U[K][I]; } } I=0; e97: A[I+1][I+1]=VR; Z[I+1][I+1]=VI; for (K=1; K<=N; K++) { T[K][I+1]=EN[K]; U[K][I+1]=EN[K+N]; } } // NORMALIZE VECTORS (BIGGEST COMPONENT TO UNITY) for (J=1; J<=N; J++) { ZM=0.0; for (I=1; I<=N; I++) { ZI=fabs(T[I][J]) + fabs(U[I][J]); if (ZI >= ZM) { IM=I; ZM=ZI; } } ZM=T[IM][J]; ZI=U[IM][J]; for (I=1; I<=N; I++) { DIVCD(T[I][J],U[I][J],ZM,ZI,&TR,&TI); T[I][J]=TR; U[I][J]=TI; } } } // print square matrix of type MAT with caption void PrintMat(char *caption, int N, MAT A) { int i,j; printf("%s\n", caption); for (i=1; i<=N; i++) { for (j=1; j<=N; j++) printf(" %10.6f", A[i][j]); printf("\n"); } } void main() { N=5; // size of given matrix // Define matrix A = AR + I * AI AR[1][1]= 1.0; AR[1][2]= 2.0; AR[1][3]= 3.0; AR[1][4]=-7.0; AR[1][5]= 12.0; AR[2][1]= 2.0; AR[2][2]= 4.0; AR[2][3]= 7.0; AR[2][4]= 3.0; AR[2][5]=-10.0; AR[3][1]= 3.0; AR[3][2]= 7.0; AR[3][3]=10.0; AR[3][4]= 8.0; AR[3][5]= 4.0; AR[4][1]=-7.0; AR[4][2]= 3.0; AR[4][3]= 8.0; AR[4][4]=-0.75; AR[4][5]= -9.0; AR[5][1]=12.0; AR[5][2]=-10.0; AR[5][3]= 4.0; AR[5][4]=-9.0; AR[5][5]= 10.0; for (i=1; i<=N; i++) for (j=1; j<=N; j++) AI[i][j]=0.0; // here matrix is real, symmetric. printf("\n"); PrintMat(" Matrix C (real parts):",N,AR); PrintMat(" Matrix C (imaginary parts):",N,AI); printf(" Any key + to continue... "); scanf("%s", answ); COMEIG (N,N,25,AR,AI,T,U,&IER,EN); printf("\n Eigenvalues (real parts):\n"); for (i=1; i<=N; i++) printf(" %10.6f", AR[i][i]); printf("\n Eigenvalues (imaginary parts):\n"); for (i=1; i<=N; i++) printf(" %10.6f", AI[i][i]); printf("\n Error code (must be zero): %d\n\n", IER); PrintMat(" Eigenvectors in columns (real parts):",N,T); PrintMat(" Eigenvectors in columns (imaginary parts):",N,U); printf("\n"); } // end of file tcomeig.cpp