CLS
REM ==CYCLIC JACOBI METHOD==
DIM A(10,10):DIM B(10):DIM U(10,10)
REM ==INPUT DATA==
GOSUB [SEZE]
REM ==INITIALIZE==
GOSUB [SISIZE]
GOSUB [NIEIZE]
REM ==OUTPUT RESULTS==
GOSUB [SEFIZE]
GOTO [THZEZEZE]
REM ==INPUT DATA==
[SEZE]
PRINT "EIGENVALUES & EIGENVECTORS"
PRINT "USING THE CYCLIC JACOBI METHOD"
PRINT
PRINT "ENTER SIZE OF PROBLEM (MATRIX)"
INPUT "N= ? ";N
PRINT
INPUT "NUMBER OF SIGNIFICANT FIGURES? ";S1
PRINT
PRINT "ENTER UPPER TRIANGLE"
PRINT "OF MATRIX A , BY COLUMNS"
PRINT
FOR J = 1 TO N
PRINT
PRINT "ENTER UPPER PART OF COLUMN ";J
FOR I = 1 TO J
PRINT "A("; I;",";J;") =";
INPUT AAA
A(I,J)=AAA
NEXT I
NEXT J
FOR I = 1 TO N
I1 = I -1
FOR J = 1 TO I1
A(I,J) = A(J,I)
NEXT J
NEXT I
PRINT
PRINT "ENTER ELEMENTS OF "
PRINT "DIAGONAL MATRIX B"
FOR I = 1 TO N
PRINT "B(";I;") = ";
INPUT BBB
B(I) = BBB
NEXT I
PRINT
RETURN
REM ==INITIALIZE TOLERANCE AND MAX NO. OF ROTATIONS==
[SISIZE]
Z = 2 * S1
T1 = 1/ (10^Z)
PRINT
PRINT "TOLERANCE= ";T1
R = 5 * N * N
R1 = 0
T2 = 0.1
N1 = N -1
RETURN
REM ==OUTPUT EIGENSOLUTION==
[SEFIZE]
PRINT
PRINT
PRINT "SOLUTION"
PRINT "********"
PRINT "NO. OF ROTATIONS REQUIRED = ";R1
PRINT
FOR J= 1 TO N
PRINT
PRINT "EIGENVALUE ";J;" IS ";B(J)
PRINT "ITS EIGENVECTOR IS "
PRINT
FOR I = 1 TO N
PRINT U(I,J)
NEXT I
PRINT
PRINT "PRESS ENTER TO CONTINUE"
INPUT C$
PRINT
NEXT J
RETURN
REM ==EIGENPROBLEM SOLUTION==
[NIEIZE]
GOSUB [ONONTHZE]
REM ==PERFORM ONE CYCLE OF ROTATIONS==
[NINIZE]
GOSUB [ONTWNIZE]
REM ==CHECK TOLERANCE==
IF X1 < T1 THEN [ONONONZE]
REM CHECK NO. OF ROTATIONS
IF R1 > R THEN [ONZEFOZE]
T2 = 0.1 * X1
GOTO [NINIZE]
[ONZEFOZE]
PRINT
PRINT" *** ERROR ***"
PRINT "NO CONVERGENCE ATTAINED"
PRINT "WITH ";R1;" ROTATIONS"
END
[ONONONZE]
GOSUB [ONEITHZE]
RETURN
[ONONTHZE]
FOR I = 1 TO N
FOR J = 1 TO N
U(I,J) = 0
NEXT J
U(I, I) = 1
NEXT I
FOR I = 1 TO N
B1 =(B(I))^.5
B(I) = 1 / B1
NEXT I
FOR I = 1 TO N
FOR J= 1 TO N
A(I,J) = B(I)* A(I,J) * B(J)
NEXT J
NEXT I
RETURN
[ONTWNIZE]
X1 = 0
FOR K = 1 TO N1
K1 = K + 1
FOR L = K1 TO N
A1 = A(K,K)
A2 = A(K,L)
A3 = A(L,L)
X = A2 * A2 / (A1 * A3)
REM ==CHECK IF ROTATION IS NEEDED==
IF X > X1 THEN [ONTHNIZE]
GOTO [ONFOZEZE]
[ONTHNIZE]
X1 = X
[ONFOZEZE]
IF X < T2 THEN [ONEIZEZE]
R1 = R1 + 1
REM ==COMPUTE ANGLE==
IF A1 = A3 THEN [ONFOSEZE]
Z = 0.5 * (A1 -A3) / A2
Z1 = 1 + 1 / (Z * Z)
T = -1*Z*(1 +Z1^.5)
GOTO [ONFOEIZE]
[ONFOSEZE]
T = 1
[ONFOEIZE]
C = 1 / (1 + T * T)^.5
S = C * T
S2 = S * S
C2 = C * C
A(K,L) = 0
REM ==TRANSFORM DIAGONALELEMENTS==
A0 = 2 * A2 * C * S
A(K,K) = A1 * C2 + A0 + A3 * S2
A(L,L) = A1 * S2 -A0 + A3 * C2
REM ==TRANSFORM OFFDIAGONALELEMENTS==
FOR I = 1 TO N
IF I < K THEN [ONSIZEZE]
IF I > K THEN [ONSIFOZE]
GOTO [ONSEFOZE]
[ONSIZEZE]
A0 = A(I,K)
A(I,K) = C * A0 + S * A(I,L)
A(I,L) = -1*S * A0 + C * A(I,L)
GOTO [ONSEFOZE]
[ONSIFOZE]
IF I < L THEN [ONSISEZE]
IF I > L THEN [ONSEONZE]
GOTO [ONSEFOZE]
[ONSISEZE]
A0 = A(K,I)
A(K,I) = C * A0 + S * A(I,L)
A(I,L) = -1*S * A0 + C * A(I,L)
GOTO [ONSEFOZE]
[ONSEONZE]
A0 = A(K,I)
A(K,I) = C * A0 + S * A(L,I)
A(L,I) = -1*S * A0 + C * A(L,I)
[ONSEFOZE]
NEXT I
REM ==TRANSFORM MATRIX U TO GENERATE EIGENVECTORS==
FOR I = 1 TO N
U0 = U(I,K)
U(I,K) = C * U0 + S * U(I,L)
U(I,L)= -1*S * U0 + C * U(I,L)
NEXT I
[ONEIZEZE]
NEXT L
NEXT K
RETURN
REM ==NORMALIZE EIGENVECTORS==
[ONEITHZE]
FOR I = 1 TO N
FOR J = 1 TO N
U(I,J) = U(I,J) * B(I)
NEXT J
NEXT I
FOR I = 1 TO N
B(I) = A(I, I)
NEXT I
REM ==ORDER EIGENSOLUTION==
FOR I = 1 TO N1
I1 = I+1
Z = B(I)
M = I
FOR J = I1 TO N
IF Z < B(J) THEN [ONNINIZE]
Z = B(J)
M = J
[ONNINIZE]
NEXT I1
B(M) = B(I)
B(I) = Z
FOR J = 1 TO N
Z = U(J,I)
U(J,I) = U(J,M)
U(J,M) = Z
NEXT J
NEXT I
RETURN
[THZEZEZE]
PRINT "END OF PROGRAM"
END