DECLARE SUB CDIV (AR#, AI#, BR#, BI#, ZR#, ZI#) DECLARE SUB CPRO (alpha#, c#(), c1#()) DECLARE SUB CONJ (c#(), c1#()) DECLARE SUB CMUL (c1#(), c2#(), c3#()) DECLARE SUB CDIF (c1#(), c2#(), c3#()) DECLARE SUB CADD (c1#(), c2#(), c3#()) DECLARE FUNCTION Squ# (A#) DECLARE FUNCTION CABS# (c#()) '*********************************************************************** '* 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: 1D-10 * '* Max. number of iterations: 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] * '* * '* Basic Release 1.1 By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '* * '* Release 1.1 (Nov. 30th, 2007): Added function NORMAL. * '*********************************************************************** 'Program TEPHJ DEFDBL A-H, O-Z DEFINT I-N ISIZE = 25 'integer I,it,J,M,N 'double dta DIM A(ISIZE, ISIZE, 2), VX(ISIZE, ISIZE, 2) 'Complex matrices DIM R(ISIZE) GOSUB 500 'call LME(N,A) (read complex matrix from text file) CLS PRINT PRINT " ** Compute Eigenvalues and Eigenvectors **" PRINT " of a Square Hermitian Matrix" PRINT " By Jacobi's Method" PRINT INPUT " Precision: ", dta INPUT " Max. number of iterations: ", M PRINT GOSUB 1000 'call EPHJ(dta, M, N, A, it, R, VX) GOSUB 2000 'call NORMAL(N, R, VX) F\$ = "##.######" IF it = -1 THEN PRINT " No convergence." ELSE 'display results FOR J = 1 TO N PRINT " Eigenvalue "; J; ": "; PRINT USING F\$; R(J) PRINT " Eigenvector:" FOR I = 1 TO N PRINT " "; PRINT USING F\$; VX(I, J, 1); IF VX(I, J, 2) >= 0# THEN PRINT " + "; ELSE PRINT " - "; END IF PRINT USING F\$; ABS(VX(I, J, 2)); PRINT " I" NEXT I PRINT NEXT J PRINT END IF END 'of main program 500 'Sub LME(N,A) 'Read from input file complex hermitian matrix OPEN "matc4.dat" FOR INPUT AS #1 INPUT #1, N FOR I = 1 TO N FOR J = 1 TO N INPUT #1, A(I, J, 1), A(I, J, 2) NEXT J NEXT I CLOSE #1 RETURN 1000 'Sub EPHJ(dta, M, N, A, it, R, VX) '************************************************************************ '* 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 * '************************************************************************ ' complex numbers DIM c0(2), c1(2), c2(2), c3(2), u0(2), tmp(2), u1(2), z0(2), Z1(2) z0(1) = 0#: z0(2) = 0#: Z1(1) = 1#: Z1(2) = 0# FOR I = 1 TO N FOR J = 1 TO N IF I = J THEN VX(I, J, 1) = Z1(1) VX(I, J, 2) = Z1(2) ELSE VX(I, J, 1) = z0(1) VX(I, J, 2) = z0(2) END IF NEXT J NEXT I it = -1: L = 1 WHILE L <= M AND it <> 1 s = 0# FOR I = 1 TO N FOR J = I + 1 TO N tmp(1) = A(I, J, 1): tmp(2) = A(I, J, 2) t0 = CABS(tmp()) IF t0 > s THEN s = t0: K0 = I: L0 = J END IF NEXT J NEXT I IF s = 0# THEN it = 1 ELSE tmp(1) = A(K0, L0, 1): tmp(2) = A(K0, L0, 2) delta = Squ(A(L0, L0, 1) - A(K0, K0, 1)) + 4# * Squ(CABS(tmp())) t0 = A(L0, L0, 1) - A(K0, K0, 1) + SQR(delta) t1 = A(L0, L0, 1) - A(K0, K0, 1) - SQR(delta) IF ABS(t0) >= ABS(t1) THEN w0 = t0 ELSE w0 = t1 END IF s0 = ABS(w0) / SQR(Squ(w0) + 4# * Squ(CABS(tmp()))) t0 = 2# * s0 / w0: CPRO t0, tmp(), c0() CONJ c0(), c1() FOR I = 1 TO K0 - 1 u0(1) = A(I, K0, 1): u0(2) = A(I, K0, 2) CMUL c0(), u0(), c2(): tmp(1) = A(I, L0, 1): tmp(2) = A(I, L0, 2) CPRO s0, tmp(), c3(): CADD c2(), c3(), tmp() A(I, K0, 1) = tmp(1): A(I, K0, 2) = tmp(2) tmp(1) = A(I, L0, 1): tmp(2) = A(I, L0, 2) CMUL c1(), tmp(), c2(): CPRO s0, u0(), c3(): CDIF c2(), c3(), tmp() A(I, L0, 1) = tmp(1): A(I, L0, 2) = tmp(2) NEXT I FOR K = K0 + 1 TO L0 - 1 u0(1) = A(K0, K, 1): u0(2) = A(K0, K, 2) CMUL c1(), u0(), c2(): tmp(1) = A(K, L0, 1): tmp(2) = A(K, L0, 2) CONJ tmp(), u1(): CPRO s0, u1(), c3() CADD c2(), c3(), tmp(): A(K0, K, 1) = tmp(1): A(K0, K, 2) = tmp(2) tmp(1) = A(K, L0, 1): tmp(2) = A(K, L0, 2) CMUL c1(), tmp(), c2(): CONJ u0(), u1(): CPRO s0, u1(), c3() CDIF c2(), c3(), tmp() A(K, L0, 1) = tmp(1): A(K, L0, 2) = tmp(2) NEXT K FOR J = L0 + 1 TO N u0(1) = A(K0, J, 1): u0(2) = A(K0, J, 2) CMUL c1(), u0(), c2(): tmp(1) = A(L0, J, 1): tmp(2) = A(L0, J, 2) CPRO s0, tmp(), c3(): CADD c2(), c3(), tmp() A(K0, J, 1) = tmp(1): A(K0, J, 2) = tmp(2) tmp(1) = A(L0, J, 1): tmp(2) = A(L0, J, 2) CMUL c0(), tmp(), c2(): CPRO s0, u0(), c3(): CDIF c2(), c3(), tmp() A(L0, J, 1) = tmp(1): A(L0, J, 2) = tmp(2) NEXT J t0 = A(K0, K0, 1) tmp(1) = A(K0, L0, 1): tmp(2) = A(K0, L0, 2) t1 = 4# * Squ(s0 * CABS(tmp())) / w0 A(K0, K0, 1) = Squ(CABS(c0())) * t0 + t1 + Squ(s0) * A(L0, L0, 1) A(L0, L0, 1) = Squ(s0) * t0 - t1 + Squ(CABS(c0())) * A(L0, L0, 1) A(K0, L0, 1) = z0(1): A(K0, L0, 2) = z0(2) FOR I = 1 TO N u0(1) = VX(I, K0, 1): u0(2) = VX(I, K0, 2) CMUL c0(), u0(), c2(): tmp(1) = VX(I, L0, 1): tmp(2) = VX(I, L0, 2) CPRO s0, tmp(), c3(): CADD c2(), c3(), tmp() VX(I, K0, 1) = tmp(1): VX(I, K0, 2) = tmp(2) tmp(1) = VX(I, L0, 1): tmp(2) = VX(I, L0, 2) CMUL c1(), tmp(), c2(): CPRO s0, u0(), c3(): CDIF c2(), c3(), tmp() VX(I, L0, 1) = tmp(1): VX(I, L0, 2) = tmp(2) NEXT I t0 = 0# FOR I = 1 TO N - 1 FOR J = I + 1 TO N tmp(1) = A(I, J, 1): tmp(2) = A(I, J, 2) t0 = t0 + Squ(CABS(tmp())) NEXT J NEXT I s = 2# * t0 IF (s < dta) THEN it = 1 ELSE L = L + 1 END IF END IF WEND IF it = 1 THEN FOR I = 1 TO N R(I) = A(I, I, 1) NEXT I END IF RETURN 2000 'Sub NORMAL(N, R, VX) '-------------------------------------------------------------- ' Sort eigenvalues by fabsolute values in ascending order and ' normalize. '-------------------------------------------------------------- DIM Tr(ISIZE), Ti(ISIZE) 'double ZM,Z1,Z2,VR DIM V(2) 'Complex number ' SORT SOLUTIONS IN ASCENDING ORDER FOR J = 2 TO N VR = R(J) FOR K = 1 TO N Tr(K) = VX(K, J, 1) Ti(K) = VX(K, J, 2) NEXT K FOR I = J - 1 TO 1 STEP -1 IF ABS(R(I)) <= ABS(VR) THEN GOTO 5 R(I + 1) = R(I) FOR K = 1 TO N VX(K, I + 1, 1) = VX(K, I, 1) VX(K, I + 1, 2) = VX(K, I, 2) NEXT K NEXT I I = 0 5 R(I + 1) = VR FOR K = 1 TO N VX(K, I + 1, 1) = Tr(K) VX(K, I + 1, 2) = Ti(K) NEXT K NEXT J ' NORMALIVXE WITH RESPECT TO BIGGEST ELEMENT FOR J = N TO 1 STEP -1 ZM = 0# FOR I = 1 TO N V(1) = VX(I, J, 1) V(2) = VX(I, J, 2) Z1 = CABS(V()) IF Z1 >= ZM THEN IM = I ZM = Z1 END IF NEXT I Z1 = VX(IM, J, 1) Z2 = VX(IM, J, 2) FOR I = 1 TO N CDIV VX(I, J, 1), VX(I, J, 2), Z1, Z2, V(1), V(2) VX(I, J, 1) = V(1) VX(I, J, 2) = V(2) NEXT I NEXT J RETURN FUNCTION CABS (c()) 'Absolute value of a complex number CABS = SQR(Squ(c(1)) + Squ(c(2))) END FUNCTION SUB CADD (c1(), c2(), c3()) 'Add two complex numbers c3(1) = c1(1) + c2(1): c3(2) = c1(2) + c2(2) END SUB SUB CDIF (c1(), c2(), c3()) 'Substract two complex numbers c3(1) = c1(1) - c2(1): c3(2) = c1(2) - c2(2) END SUB SUB CDIV (AR, AI, BR, BI, ZR, ZI) ' EFFECT 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 ABS(YR) > ABS(YI) THEN 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 END IF END SUB SUB CMUL (c1(), c2(), c3()) 'Multiply two complex numbers c3(1) = c1(1) * c2(1) - c1(2) * c2(2) c3(2) = c1(1) * c2(2) + c1(2) * c2(1) END SUB SUB CONJ (c(), c1()) 'Return conjugate of a complex number c1(1) = c(1): c1(2) = -c(2) END SUB SUB CPRO (alpha, c(), c1()) 'Multiply a complex number by a real number c1(1) = alpha * c(1): c1(2) = alpha * c(2) END SUB FUNCTION Squ (A) Squ = A * A END FUNCTION