/************************************************************************ * Eigenvalues and Eigenvectors of a general complex square matrix * * using the QR method for complex matrices * * --------------------------------------------------------------------- * * SAMPLE RUN: * * * * Example #1: * * * * Input matrix (real symmetric for comparison with Jacobi): * * 1, 2, 3, -7, 12 * * 2, 4, 7, 3, -1 * * 3, 7, 10, 8, 4 * * -7, 3, 8, -0.75,-9 * * 12, -1, 4, -9, 10 * * * * Input Matrix (real part): * * 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 * * * * Error code = 0 * * * * Eigenvalues: * * 0.463350 0.000000 * * -7.774580 0.000000 * * -10.486545 0.000000 * * 18.291821 0.000000 * * 23.755955 0.000000 * * * * Eigenvectors (real part in lines): * * 0.208812 1.000000 -0.478027 -0.275000 -0.216915 * * 1.000000 -0.326100 0.209620 -0.147689 -0.815422 * * 0.467172 -0.007947 -0.507827 1.000000 0.264432 * * 0.093053 0.594911 1.000000 0.455666 0.050740 * * 0.705820 -0.012677 0.131486 -0.527499 1.000000 * * * * Example #2: * * * * Input matrix (real part): * * 1, 5, 3, -7, 12 * * 2, 4, 7, 3, -1 * * 3, 7, 10, 8, 4 * * -7, 3, 8, -0.75,-9 * * 12, -1, 4, -9, 10 * * Input matrix (real part): * * 1, 1, 1, 1, 1 * * 1, 1, 1, 1, 1 * * 1, 1, 1, 1, 1 * * 1, 1, 1, 1, 1 * * 1, 1, 1, 1, 1 * * * * Input Matrix (real part): * * 1.000000 5.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 * * Input Matrix (imag. part): * * 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 * * * * Error code = 0 * * * * Eigenvalues: * * 0.888713 0.034707 * * -8.284268 0.024981 * * -10.363673 0.879708 * * 18.894792 3.155582 * * 23.114437 0.905021 * * * * Eigenvectors (real part in lines): * * 0.411502 1.000000 -0.532508 -0.217945 -0.417001 * * 1.000000 -0.324189 0.269140 -0.277587 -0.872694 * * 0.451233 -0.007289 -0.515266 1.000000 0.272732 * * 0.179884 0.608599 1.000000 0.456406 0.069197 * * 0.682657 -0.106663 -0.004456 -0.615170 1.000000 * * * * Eigenvectors (imag. part in lines): * * 0.000770 0.000000 -0.034750 0.009062 -0.005181 * * 0.000000 0.001676 0.051622 -0.074119 -0.037648 * * 0.043492 -0.070147 -0.050515 0.000000 -0.066367 * * -0.198515 0.075502 0.000000 0.271030 -0.314928 * * 0.044071 0.193806 0.292159 0.172780 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 function MTPRINT and Example #2 (12/04/2004). * ************************************************************************* To link with basis_r.cpp and vmblock.cpp. ----------------------------------------- */ #include #include //global variables int IERR; //error code (must be zero) int NM; //1st dimension of input complex matrix int N; //2nd dimension of input complex matrix REAL **AR, **AI; //pointers to real, imag. parts of complex matrix A REAL **ZR, **ZI; //complex matrix to store complex eigenvectors REAL *WR, *WI; //real, imag. parts of complex eigenvalues REAL *WORK; //working space, size = N boolean AVEC; //TRUE = eigenvectors asked for int I,J; //loop index void *vmblock; //pointer to a working space for dynamic allocations char c; // print a matrix with caption void MPRINT(int N, int MMAX, int M, char *Name, REAL **A) { int I,J; printf("\n%s\n",Name); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) printf("%12.6f", A[I][J]); printf("\n"); } } // print transpose matrix with caption void MTPRINT(int N, int MMAX, int M, char *Name, REAL **A) { int I,J; printf("\n%s\n",Name); for (I=1; I<=N; I++) { for (J=1; J<=N; J++) printf("%12.6f", A[J][I]); printf("\n"); } } //for debug only void VPRINT(int N, char *Name, REAL *V) { int I; printf("\n%s\n",Name); for (I=1; I<=N; I++) { printf("%12.6f", V[I]); if ((I % 5) == 0) printf("\n"); } } //Headers of functions defined below void CBALAN(int NM,int N,REAL **AR,REAL **AI,int *LOW,int *UPP,REAL *D); void CUNITH(int NM,int N,int LOW,int UPP,REAL **AR,REAL **AI,REAL *DR, REAL *DI); void COMQR2(int NM,int N,int LOW,int UPP,REAL *GR,REAL *GI, REAL **HR, REAL **HI, REAL *WR, REAL *WI, boolean AVEC, REAL **ZR, REAL **ZI, int *IERR); void CBALBK(int NM,int N,int LOW,int UPP,int M,REAL *D,REAL **Z,REAL **T); void NORMAL(int NM,int N,REAL *WR,REAL *WI,boolean AVEC,REAL **ZR,REAL **ZI,REAL *TR,REAL *TI); REAL ABSCD (REAL AR, REAL AI); void DCSQRT(REAL X,REAL Y,REAL *A,REAL *B); void DIVCD (REAL AR,REAL AI,REAL BR,REAL BI,REAL *ZR,REAL *ZI); void EXCHG1(int M,int NM,int N,REAL **A,REAL **B, REAL *D,int J,int K,int L); REAL DCABS (REAL XX,REAL YY); //Utility function int IMIN(int a,int b) { if (a=M; I--) { DR[I]=AR[I][M-1]/S; DI[I]=AI[I][M-1]/S; H += DR[I]*DR[I]+DI[I]*DI[I]; } G=SQRT(H); F=ABSCD(DR[M],DI[M]); if (F == 0.0) { DR[M]=G; AR[M][M-1]=S; } else { H=H+F*G; G=G/F; DR[M]=(1.0+G)*DR[M]; DI[M]=(1.0+G)*DI[M]; } for (J=M; J<=N; J++) { XR=0.0; XI=0.0; for (I=UPP; I>=M; I--) { XR += (DR[I]*AR[I][J]+DI[I]*AI[I][J]); XI += (DR[I]*AI[I][J]-DI[I]*AR[I][J]); } XR /= H; XI /= H; for (I=M; I<=UPP; I++) { AR[I][J] += (-XR*DR[I]+XI*DI[I]); AI[I][J] -= (XR*DI[I]+XI*DR[I]); } } for (I=1; I<=UPP; I++) { XR=0.0; XI=0.0; for (J=UPP; J>=M; J--) { XR += (DR[J]*AR[I][J]-DI[J]*AI[I][J]); XI += (DR[J]*AI[I][J]+DI[J]*AR[I][J]); } XR=XR/H; XI=XI/H; for (J=M; J<=UPP; J++) { AR[I][J] -= (XR*DR[J]+XI*DI[J]); AI[I][J] += (XR*DI[J]-XI*DR[J]); } } DR[M] *= S; DI[M] *= S; AR[M][M-1]=-G*AR[M][M-1]; AI[M][M-1]=-G*AI[M][M-1]; } // if S<>0 } //M loop } //CUNITH() void COMQR2(int NM,int N,int LOW,int UPP,REAL *GR,REAL *GI, REAL **HR, REAL **HI, REAL *WR, REAL *WI, boolean AVEC, REAL **ZR, REAL **ZI, int *IERR) { /*--------------------------------------------------------------------- QR ALGORITH APPLIED TO A COMPLEX HESSENBERG MATRIX (translated from Fortran 77). --------------------------------------------------------------------- */ // Labels: e5,e7,e10,e15,e20,e35,e45,e50,e55 REAL SR,SI,TR,TI,YR,YI,XR,XI,VR,VI,MACHEP,NORM,EPS,EPS1; int I,J,K,L,M,IN0,NA,ITS; // SEEK FLOATING POINT RELATIVE PRECISION FOR REAL EPS=1.0; do { EPS *= 0.5; EPS1=1.0+EPS; } while (EPS1 > 1.0); MACHEP=2.0*EPS; *IERR=0; if (AVEC) for (I=1; I<=N; I++) { for (J=1; J<=N; J++) { ZR[I][J]=0.0; ZI[I][J]=0.0; } ZR[I][I]=1.0; } // ACCUMULATE UNITARY TRANSFORMATIONS for (I=UPP-1; I>LOW; I--) { if (((HR[I][I-1] == 0.0) && (HI[I][I-1] == 0.0)) || ((GR[I] == 0.0) && (GI[I] == 0.0))) goto e5; NORM=HR[I][I-1]*GR[I]+HI[I][I-1]*GI[I]; for (K=I+1; K<=UPP; K++) { GR[K]=HR[K][I-1]; GI[K]=HI[K][I-1]; } for (J=I; J<=UPP; J++) { SR=0.0; SI=0.0; for (K=I; K<=UPP; K++) { SR += GR[K]*ZR[K][J]+GI[K]*ZI[K][J]; SI += GR[K]*ZI[K][J]-GI[K]*ZR[K][J]; } SR /= NORM; SI /= NORM; for (K=I; K<=UPP; K++) { ZR[K][J] += SR*GR[K]-SI*GI[K]; ZI[K][J] += SR*GI[K]+SI*GR[K]; } } e5: ;} // CREATE REAL SUBDIAGONAL L=LOW+1; for (I=L; I<=UPP; I++) { M=IMIN(I+1,UPP); if (HI[I][I-1] == 0.0) goto e7; NORM=ABSCD(HR[I][I-1],HI[I][I-1]); YR=HR[I][I-1]/NORM; YI=HI[I][I-1]/NORM; HR[I][I-1]=NORM; HI[I][I-1]=0.0; for (J=I; J<=N; J++) { SI= YR*HI[I][J]-YI*HR[I][J]; HR[I][J]=YR*HR[I][J]+YI*HI[I][J]; HI[I][J]=SI; } for (J=1; J<=M; J++) { SI= YR*HI[J][I]+YI*HR[J][I]; HR[J][I]=YR*HR[J][I]-YI*HI[J][I]; HI[J][I]=SI; } if (AVEC) for (J=LOW; J<=UPP; J++) { SI= YR*ZI[J][I]+YI*ZR[J][I]; ZR[J][I]=YR*ZR[J][I]-YI*ZI[J][I]; ZI[J][I]=SI; } e7: ;} // STORE ROOTS ALREADY ISOLATED BY FUNCTION BALANC for (I=1; I<=N; I++) if (I < LOW || I > UPP) { WR[I]=HR[I][I]; WI[I]=HI[I][I]; } IN0=UPP; TR=0.0; TI=0.0; // SEEK EIGENVALUES e10: if (IN0 < LOW) goto e50; ITS=0; NA=IN0-1; // WE LOOK AT REAL PART OF SUBDIAGONAL ELEMENT e15: for (L=IN0; L>=LOW; L--) { //here index 0 is temporarily used! SR=MACHEP*(ABS(HR[L-1][L-1])+ABS(HI[L-1][L-1]) +ABS(HR[L][L])+ABS(HI[L][L])); if (ABS(HR[L][L-1]) <= SR || L == LOW) goto e20; } // SHIFTING e20: if (L == IN0) goto e45; if (ITS == 30) goto e55; if (ITS == 10 || ITS == 20) { // EXTRA SHIFTING SR=ABS(HR[IN0][NA])+ABS(HR[NA][IN0-2]); SI=0.0; } else { SR=HR[IN0][IN0]; SI=HI[IN0][IN0]; XR=HR[NA][IN0]*HR[IN0][NA]; XI=HI[NA][IN0]*HR[IN0][NA]; if (XR != 0.0 || XI != 0.0) { YR=(HR[NA][NA]-SR)*0.5; YI=(HI[NA][NA]-SI)*0.5; DCSQRT(YR*YR-YI*YI+XR,2.0*YR*YI+XI,&VR,&VI); if (YR*VR+YI*VI < 0.0) { VR=-VR; VI=-VI; } DIVCD(XR,XI,YR+VR,YI+VI,&VR,&VI); SR=SR-VR; SI=SI-VI; } } for (I=LOW; I<=IN0; I++) { HR[I][I] -= SR; HI[I][I] -= SI; } TR += SR; TI += SI; ITS++; // TRIANGULARY DECOMPOSITION H = Q * R for (I=L+1; I<=IN0; I++) { SR=HR[I][I-1]; VR=HR[I-1][I-1]; VI=HI[I-1][I-1]; NORM=SQRT(VR*VR+VI*VI+SR*SR); XR=VR/NORM; XI=VI/NORM; WR[I-1]=XR; WI[I-1]=XI; HR[I-1][I-1]=NORM; HI[I-1][I-1]=0.0; HR[I][I-1]=0.0; HI[I][I-1]=SR/NORM; for (J=I; J<=N; J++) { YR=HR[I-1][J]; YI=HI[I-1][J]; VR=HR[I][J]; VI=HI[I][J]; HR[I-1][J]=XR*YR+XI*YI+HI[I][I-1]*VR; HI[I-1][J]=XR*YI-XI*YR+HI[I][I-1]*VI; HR[I][J] =XR*VR-XI*VI-HI[I][I-1]*YR; HI[I][J] =XR*VI+XI*VR-HI[I][I-1]*YI; } } SI=HI[IN0][IN0]; if (SI == 0.0) goto e35; NORM=ABSCD(HR[IN0][IN0],SI); SR=HR[IN0][IN0]/NORM; SI /= NORM; HI[IN0][IN0]=0.0; HR[IN0][IN0]=NORM; if (IN0 != N) for (J=IN0+1; J<=N; J++) { YR=HR[IN0][J]; YI=HI[IN0][J]; HR[IN0][J]=SR*YR+SI*YI; HI[IN0][J]=SR*YI-SI*YR; } // INVERSE OPERATION e35: for (J=L+1; J<=IN0; J++) { XR=WR[J-1]; XI=WI[J-1]; for (I=1; I<=J; I++) { YR=HR[I][J-1]; YI=0.0; VR=HR[I][J]; VI=HI[I][J]; if (I != J) { YI=HI[I][J-1]; HI[I][J-1]=XR*YI+XI*YR+HI[J][J-1]*VI; } HR[I][J-1]=XR*YR-XI*YI+HI[J][J-1]*VR; HR[I][J] =XR*VR+XI*VI-HI[J][J-1]*YR; HI[I][J] =XR*VI-XI*VR-HI[J][J-1]*YI; } if (AVEC) for (I=LOW; I<=UPP; I++) { YR=ZR[I][J-1]; YI=ZI[I][J-1]; VR=ZR[I][J]; VI=ZI[I][J]; ZR[I][J-1]=XR*YR-XI*YI+HI[J][J-1]*VR; ZI[I][J-1]=XR*YI+XI*YR+HI[J][J-1]*VI; ZR[I][J] =XR*VR+XI*VI-HI[J][J-1]*YR; ZI[I][J] =XR*VI-XI*VR-HI[J][J-1]*YI; } } if (SI == 0.0) goto e15; for (I=1; I<=IN0; I++) { YR=HR[I][IN0]; YI=HI[I][IN0]; HR[I][IN0]=SR*YR-SI*YI; HI[I][IN0]=SR*YI+SI*YR; } if (AVEC) for (I=LOW; I<=UPP; I++) { YR=ZR[I][IN0]; YI=ZI[I][IN0]; ZR[I][IN0]=SR*YR-SI*YI; ZI[I][IN0]=SR*YI+SI*YR; } goto e15; // ONE ROOT HAS BEEN FOUND e45: HR[IN0][IN0]=HR[IN0][IN0]+TR; HI[IN0][IN0]=HI[IN0][IN0]+TI; WR[IN0]=HR[IN0][IN0]; WI[IN0]=HI[IN0][IN0]; IN0=NA; goto e10; // ALL ROOTS HAVE BEEN FOUND // START CALCULATION OF EIGENVECTORS BY SUBSTITUTION e50: if (AVEC) { NORM=0.0; for (I=1; I<=N; I++) for (J=I; J<=N; J++) NORM += ABS(HR[I][J])+ABS(HI[I][J]); if (NORM == 0.0) { printf(" 1 IERR=%d\n", *IERR); return; } for (IN0=N; IN0>1; IN0--) { XR=WR[IN0]; XI=WI[IN0]; HR[IN0][IN0]=1.0; HI[IN0][IN0]=0.0; for (I=IN0-1; I>0; I--) { VR=HR[I][IN0]; VI=HI[I][IN0]; if (I != IN0-1) for (J=I+1; J UPP) for (J=I+1; J<=N; J++) { ZR[I][J]=HR[I][J]; ZI[I][J]=HI[I][J]; } } // END CALCULATION OH HESSENBERG EIGENVECTORS for (J=N; J>LOW; J--) { M=IMIN(J-1,UPP); for (I=LOW; I<=UPP; I++) { VR=ZR[I][J]; VI=ZI[I][J]; for (K=LOW; K<=M; K++) { VR += ZR[I][K]*HR[K][J]-ZI[I][K]*HI[K][J]; VI += ZR[I][K]*HI[K][J]+ZI[I][K]*HR[K][J]; } ZR[I][J]=VR; ZI[I][J]=VI; } } } //if AVEC // END CALCULATION OF EIGENVECTORS OF INPUT MATRIX return; // NO CONVERGENCE AFTER 30 ITERATIONS! e55: *IERR = IN0; } //COMQR2() void CBALAN(int NM,int N,REAL **AR,REAL **AI,int *LOW,int *UPP,REAL *D) { // BALANCE A COMPLEX MATRIX (VERSION DERIVED FROM BALANC) // No effect for a symmetrical matrix). // Labels: e10,e11,e20,e22,e30,e40,e50 REAL B2,C,F,G,R,S; boolean NOCONV; int B,J,J1,K,K1,L,L1; B=16; B2=B*B; L=1; K=N; e10: for (J=K; J>0; J--) { J1=J; K1=K; L1=L; for (I=1; I<=K; I++) { if (I != J) if (AR[J][I] != 0.0 || AI[J][I] != 0.0) goto e11; } EXCHG1(K,NM,N,AR,AI,D,J1,K1,L1); K--; goto e10; e11: ;} e20: for (J=L; J<=K; J++) { J1=J; K1=K; L1=L; for (I=L; I<=K; I++) { if (I != J) if (AR[I][J] != 0.0 || AI[I][J] != 0.0) goto e22; } EXCHG1(L,NM,N,AR,AI,D,J1,K1,L1); L++; goto e20; e22: ;} *LOW=L; *UPP=K; for (I=L; I<=K; I++) D[I]=1.0; e30: NOCONV=FALSE; for (I=L; I<=K; I++) { C=0.0; R=0.0; for (J=L; J<=K; J++) if (J != I) { C += ABS(AR[J][I])+ABS(AI[J][I]); R += ABS(AR[I][J])+ABS(AI[I][J]); } G=R/B; F=1.0; S=C+R; e40: if (C < G) { F=F*B; C=C*B2; goto e40; } G=R*B; e50: if (C >= G) { F=F/B; C=C/B2; goto e50; } if ((C+R)/F < 0.95*S) { G=1.0/F; D[I] *= F; NOCONV=TRUE; for (J=L; J<=N; J++) { AR[I][J] *= G; AI[I][J] *= G; } for (J=1; J<=K; J++) { AI[J][I] *= F; AR[J][I] *= F; } } } if (NOCONV) goto e30; } //CBALAN() void CBALBK(int NM,int N,int LOW,int UPP,int M,REAL *D,REAL **Z,REAL **T) { // This Procedure finds the eigenvectors of a complex matrix, // as output of function CBALAN. // Label: e10 REAL S; int I,J,K,L; for (I=LOW; I<=UPP; I++) { S=D[I]; for (J=1; J<=M; J++) { Z[I][J] *= S; T[I][J] *= S; } } for (L=1; L<=N; L++) { I=L; if (I >= LOW && I <= UPP) goto e10; if (I < LOW) I=LOW-L; K=int(D[I]); if (K != I) for (J=1; J<=M; J++) { S=Z[I][J]; Z[I][J]=Z[K][J]; Z[K][J]=S; S=T[I][J]; T[I][J]=T[K][J]; T[K][J]=S; } e10: ;} } void NORMAL(int NM,int N,REAL *WR,REAL *WI,boolean AVEC,REAL **ZR,REAL **ZI, REAL *TR,REAL *TI) { // Sort eigenvalues by absolute values in ascending order and // normalize. // Label: e5 REAL ZM,Z1,Z2,VR,VI; int IM,J,K; // SORT SOLUTIONS IN ASCENDING ORDER for (J=2; J<=N; J++) { VR=WR[J]; VI=WI[J]; if (AVEC) for (K=1; K<=N; K++) { TR[K]=ZR[K][J]; TI[K]=ZI[K][J]; } for (I=J-1; I>0; I--) { if (ABSCD(WR[I],WI[I]) <= ABSCD(VR,VI)) goto e5; WR[I+1]=WR[I]; WI[I+1]=WI[I]; if (AVEC) for (K=1; K<=N; K++) { ZR[K][I+1]=ZR[K][I]; ZI[K][I+1]=ZI[K][I]; } } I=0; e5: WR[I+1]=VR; WI[I+1]=VI; if (AVEC) for (K=1; K<=N; K++) { ZR[K][I+1]=TR[K]; ZI[K][I+1]=TI[K]; } } // NORMALIZE WITH RESPECT TO BIGGEST ELEMENT if (AVEC) for (J=N; J>0; J--) { ZM= 0.0; for (I=1; I<=N; I++) { VR=ZR[I][J]; VI=ZI[I][J]; Z1=ABSCD(VR,VI); if (Z1 >= ZM) { IM = I; ZM = Z1; } } Z1 = ZR[IM][J]; Z2 = ZI[IM][J]; for (I = 1; I<=N; I++) { DIVCD(ZR[I][J],ZI[I][J],Z1,Z2,&VR,&VI); ZR[I][J] = VR; ZI[I][J] = VI; } } } //NORMAL() // Utility routines for complex numbers REAL ABSCD(REAL AR,REAL AI) { // ABSOLUTE VALUE OF A COMPLEX NUMBER C=AR+I*AI // ABSCD,AR,AI REAL*8 // ABSCD=SQRT(AR**2+AI**2) REAL XR,XI,W; XR=ABS(AR); XI=ABS(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 DCSQRT(REAL X,REAL Y, REAL *A, REAL *B) { // SQUARE ROOT OF A COMPLEX NUMBER A+I*B = RAC(X+I*Y) if (X == 0.0 && Y == 0.0) { *A=0.0; *B=0.0; } else { *A=SQRT((ABS(X)+DCABS(X,Y))*0.5); if (X >= 0.0) *B=Y/(*A+*A); else if (Y < 0.0) *B=-*A; else { *B=*A; *A=Y/(*B+*B); } } } // ABSOLUTE VALUE OF A COMPLEX NUMBER X+I*Y REAL DCABS(REAL XX,REAL YY) { REAL X,Y,W; X=ABS(XX); Y=ABS(YY); if (X == 0.0) W=Y; else if (Y == 0.0) W=X; else if (X > Y) W=X*SQRT(1.0+(Y/X)*(Y/X)); else W=Y*SQRT(1.0+(X/Y)*(Y/X)); return W; } void DIVCD (REAL AR,REAL AI,REAL BR,REAL BI, REAL *ZR,REAL *ZI) { // EFFECT DIVISION Z=(ZR+I*ZI)=(AR+I*AI)/(BR+I*BI) // ---- DO NOT USE IF BR=BI=0 ---- REAL YR,YI,W; YR=BR; YI=BI; if (ABS(YR) > ABS(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; } } // Exchange elements of a complex matrix void EXCHG1(int M,int NM,int N, REAL **A,REAL **B, REAL *D, int J,int K,int L) { REAL F; D[M]=J; if (J != M) { for (I=1; I<=K; I++) { F=A[I][J]; A[I][J]=A[I][M]; A[I][M]=F; F=B[I][J]; B[I][J]=B[I][M]; B[I][M]=F; } for (I=L; I<=N; I++) { F=A[J][I]; A[J][I]=A[M][I]; A[M][I]=F; F=B[J][I]; B[J][I]=B[M][I]; B[M][I]=F; } } } void main() { NM = 5; N = 5; // size of input complex matrix // dynamic allocations vmblock = vminit(); AR = (REAL **) vmalloc(vmblock, MATRIX, NM+1, N+1); AI = (REAL **) vmalloc(vmblock, MATRIX, NM+1, N+1); ZR = (REAL **) vmalloc(vmblock, MATRIX, NM+1, N+1); ZI = (REAL **) vmalloc(vmblock, MATRIX, NM+1, N+1); WR = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); WI = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); WORK = (REAL *) vmalloc(vmblock, VEKTOR, N+1, 0); for (I=0; I<=NM; I++) for (J=0; J<=NM; J++) { AR[I][J] = 0.0; AI[I][J] = 0.0; } for (I=0; I<=N; I++) WORK[I] = 0.0; //data section #1 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] = -1.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] =-1.0; AR[5][3] = 4.0; AR[5][4] = -9.0; AR[5][5] = 10.0; /* data section #2 AR[1][1] = 1.0; AR[1][2] = 5.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] = -1.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] =-1.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] = 1.0; */ MPRINT(N,NM,N," Input Matrix (real part):", AR); /* #2 only: MPRINT(N,NM,N," Input Matrix (imag. part):", AI); printf(" Any key + Enter to continue... "); scanf("%s", c); */ //call main subroutine CEIGEN(AR,AI,NM,N,TRUE,WR,WI,ZR,ZI,WORK,&IERR); printf("\n Error code = %d\n\n", IERR); if (IERR == 0) { //no error printf(" Eigenvalues:\n"); for (I=1; I<=N; I++) printf("%12.6f%12.6f\n", WR[I], WI[I]); MTPRINT(N,NM,N," Eigenvectors (real part in lines):", ZR); /* #2 only: MTPRINT(N,NM,N," Eigenvectors (imag. part in lines):", ZI); */ printf("\n"); } //free memory vmfree(vmblock); } //end of main program // end of file tceigen.cpp