SUBROUTINE SIMITZ(N, P, KM, EPS, IP, OP, INF, EM, X, MN, D, WK) 10 C*********************************************************************** C IDENTIFICATION C SIMITZ - ITERATIVE COMPUTATION OF EIGENVALUES LARGEST IN MAGNI- C TUDE AND CORRESPONDING EIGENVECTORS OF A REAL GENERAL- C IZED SYMMETRIC MATRIX C FORTRAN SUBROUTINE SUBPROGRAM C US AIR FORCE FLIGHT DYNAMICS LABORATORY C WRIGHT-PATTERSON AFB, OHIO 45433 C PURPOSE C A REAL N-SQUARE MATRIX C IS B-SYMMETRIC RELATIVE TO AN N-SQUARE C POSITIVE DEFINITE MATRIX B IN CASE BC = (C-TRANSPOSED)B. C GIVEN AS OPTIONAL INPUT A SET OF P INITIAL APPROXIMATE C EIGENVECTORS OF A REAL N-SQUARE B-SYMMETRIC MATRIX C CORRES- C PONDING TO P EIGENVALUES OF C LARGEST IN MAGNITUDE, SIMITZ COM- C PUTES EM EIGENVALUES AND EM CORRESPONDING EIGENVECTORS TO A C PRECISION DEPENDENT ON THE STRUCTURE OF B AND C AND ON A GIVEN C TOLERANCE EPS. THE MATRIX B IS PRESENTED TO SIMITZ AS AN ALGO- C RITHM FOR CALCULATING THE STANDARD INNER PRODUCT (W, BZ) = C (W-TRANSPOSED)BZ GIVEN COLUMN N-VECTORS W AND Z IMPLEMENTED AS C A FORTRAN COMPATIBLE REAL FUNCTION SUBPROGRAM. THE MATRIX C IS C PRESENTED AS A SUBROUTINE SUBPROGRAM WHICH GIVEN A COLUMN C N-VECTOR Z CALCULATES ITS IMAGE W = CZ UNDER THE MATRIX C. C DEPENDING ON THE CHOICE OF B AND C, SIMITZ APPLIES TO A WIDE C VARIETY OF SYMMETRIC EIGENPROBLEMS. C CONTROL C C DIMENSION X(MN,P), D(P), WK(K) C INTEGER P, EM C REAL IP C EXTERNAL IP, INF, OP C . C . C . C CALL SIMITZ(N, P, KM, EPS, IP, OP, INF, EM, X, MN, D, WK) C C WHERE C N IS AN INTEGER INPUT VARIABLE, THE ORDER OF THE MATRIX C. C P IS AN INTEGER INPUT VARIABLE, THE NUMBER OF SIMULTANEOUS C ITERATION VECTORS. C KM AS AN INTEGER INPUT VARIABLE IS IN MAGNITUDE THE MAXIMUM C NUMBER OF ITERATION STEPS TO BE EXECUTED. IF KM IDENTIFIES C A NEGATIVE VALUE THEN P INITIAL APPROXIMATE EIGENVECTORS C ARE ASSUMED TO BE PRESENT IN THE ARRAY X. OTHERWISE SIMITZ C SUPPLIES RANDOM INITIAL EIGENVECTORS. C KM AS AN INTEGER OUTPUT VARIABLE IDENTIFIES THE NUMBER KS OF C ITERATION STEPS FINALLY USED IN THE CALCULATION OF EM C EIGENVECTORS. C EPS IS A REAL INPUT VARIABLE, THE TOLERANCE FOR ACCEPTING C EIGENVECTORS. AS SOON AS SUCCESSIVE ITERATES OF THE RITZ C VALUES ABS(D(H+1)) DIFFER BY LESS THAN ABS(D(H+1))*EPS/10.0 C THEN D(H+1) IS ACCEPTED AS AN EIGENVALUE AND H, THE NUMBER C OF PREVIOUSLY ACCEPTED EIGENVALUES, IS INCREASED BY 1. AS C SOON AS THE ERROR QUANTITIES F(I), NORMALIZED RESIDUALS, C SATISFY D(I)*F(I)/(D(I) - D(P)) .LT. EPS, THEN G, THE NUM- C BER OF ALREADY ACCEPTED RITZ VECTORS, IS INCREASED TO C G + L, I = G + 1, ..., L. THE F(I) ARE DISCOUNTED WITH C SUCCESSIVE ITERATIONS TO FORCE CONVERGENCE IN CASE OF UN- C FORTUNATE CHOICE OF PARAMETERS. IF M SIGNIFICANT DIGITS C OF ACCURACY ARE REQUIRED OF THE EIGENVALUES, THEN SET C EPS EQUAL TO 10.0**(-M) AS A GENERAL RULE. C IP IS AN EXTERNAL INPUT VARIABLE, THE NAME OF A FORTRAN COM- C PATIBLE REAL FUNCTION SUBPROGRAM OF THE FORM IP(N, Z, W) C WHICH MUST RETURN THE INNER PRODUCT (W, BZ) = C (W-TRANSPOSED)BZ OF THE VECTORS IDENTIFIED BY THE N-ARRAYS C Z AND W RELATIVE TO THE POSITIVE DEFINITE MATRIX B. C OP IS AN EXTERNAL INPUT VARIABLE, THE NAME OF A FORTRAN COM- C PATIBLE SUBROUTINE SUBPROGRAM OF THE FORM OP(N, Z, W) C WHICH MUST CALCULATE THE IMAGE W OF THE VECTOR IDENTIFIED C BY THE N-ARRAY Z UNDER THE N-SQUARE MATRIX C WITHOUT OVER- C WRITING Z. C INF IS AN EXTERNAL INPUT VARIABLE, THE NAME OF A FORTRAN COM- C PATIBLE SUBROUTINE SUBPROGRAM WHICH MAY BE USED FOR C OBTAINING INFORMATION OR TO EXERT CONTROL DURING EXECUTION C OF SIMITZ. INF HAS THE FORM INF(KS, G, H, F) WHERE C KS IS AN INTEGER OUTPUT VARIABLE, THE NUMBER OF THE NEXT C ITERATION STEP. C G IS AN INTEGER OUTPUT VARIABLE, THE NUMBER OF ALREADY C ACCEPTED EIGENVECTORS. C H IS AN INTEGER OUTPUT VARIABLE, THE NUMBER OF ALREADY C ACCEPTED EIGENVALUES. C F IS A REAL OUTPUT VARIABLE P-ARRAY, ERROR QUANTITIES C MEASURING RESPECTIVELY THE STATE OF CONVERGENCE OF C THE P SIMULTANEOUS ITERATION VECTORS. IN ADDITION, C IF CONVERGENCE FAILS IN SUBROUTINE IMTQL2 AFTER G C EIGENVECTORS HAVE BEEN ACCEPTED, THEN F(G+1) IS RE- C PLACED BY 1000.*FLOAT(IERR) WHERE IERR IS THE ERROR C INDICATOR OUTPUT BY IMTQL2. EACH ELEMENT OF THE ARRAY C F IS INITIALLY SET BY SIMITZ TO THE VALUE 4.0. C EM AS AN INTEGER INPUT VARIABLE IS THE NUMBER OF EIGENVALUES C TO BE COMPUTED, 0 .LT. EM .LT. P .LE. N .LE. MN. C EM AS AN INTEGER OUTPUT VARIABLE IS THE NUMBER OF EIGENVECTORS C COMPUTED THROUGH KM ITERATION STEPS. C X AS A REAL N-BY-P INPUT ARRAY IS A SET OF P OPTIONAL INITIAL C APPROXIMATE EIGENVECTORS X(I,1), ..., X(I,P), I = 1, ..., C N, INTERPRETED BY SIMITZ IF KM IS NEGATIVE. C X AS A REAL N-BY-P OUTPUT ARRAY IS A SET OF EM EIGENVECTORS C X(I,1), ..., X(I,EM), I = 1, ..., N, COMPUTED THROUGH C IABS(KM) ITERATION STEPS WITH THE REMAINDER OF X CONSISTING C OF P - EM APPROXIMATE EIGENVECTORS. THE N-BY-P MATRIX X C SATISFIES (X-TRANSPOSED)BX = I, THAT IS, THE EIGENVECTORS C OF C ARE B-ORTHONORMAL. C MN IS AN INTEGER INPUT VARIABLE WHICH IDENTIFIES THE LEADING C DIMENSION IN THE CALLING PROGRAM OF THE ARRAY X. C D IS A REAL OUTPUT P-ARRAY OF WHICH D(1), ..., D(EM) ARE THE C EIGENVALUES OF C LARGEST IN MAGNITUDE IN DECREASING ORDER C CORRESPONDING TO THE COMPUTED EIGENVECTORS X(I,1), ..., C X(I,EM), I = 1, ..., N. D(EM+1), ..., D(P-1) CONTAIN C APPROXIMATIONS TO PROGRESSIVELY SMALLER SUCH EIGENVALUES. C D(P) CONTAINS THE MOST RECENTLY COMPUTED VALUE OF E, WHERE C THE INTERVAL (-E, E) IS THE INTERVAL OVER WHICH THE C CHEBYSHEV ACCELERATION WAS PERFORMED. C WK THE INITIAL LOCATION OF AT LEAST P**2 + 3*P + 2*N = K C CONSECUTIVE STORAGE LOCATIONS WHICH MAY NOT BE OVER- C WRITTEN WHILE SIMITZ IS IN EXECUTION. C OTHER PROGRAMMING INFORMATION C SIMITZ EMPLOYS A DATA STATEMENT TO ASSIGN TO A MACHINE DEPEN- C DENT REAL VARIABLE MT THE QUOTIENT OF THE SMALLEST POSITIVE C REAL VALUE REPRESENTABLE BY FORTRAN AND THE SMALLEST REAL VALUE C WHOSE SUM WITH 1.0 EXCEEDS 1.0. C C THE PERFORMANCE OF SIMITZ IS STRONGLY DEPENDENT UPON THE CHOICE C OF INPUT PARAMETERS AND UPON THE CAREFUL PREPARATION OF THE C SUBPROGRAMS IP AND OP. THE USER SHOULD CONSIDER USING HIS OWN C ACTIVE SUBROUTINE INF TO MONITOR PROGRESS OF SIMITZ RELATIVE TO C HIS CHOICE OF INPUT PARAMETERS IF NO INFORMATION IS OTHERWISE C AVAILABLE CONCERNING THE LOCATIONS OF THE RELEVANT EIGENVALUES. C RECALL THAT SIMITZ MAY BE REENTERED WITH KM .LT. 0 WITHOUT LOSS C OF INFORMATION TO PERMIT CONSERVATIVE INITIAL CHOICES OF C ABS(KM), EPS AND P. C OTHER PROGRAMS REQUIRED C FUNCTION RANF C RETURNS UNIFORMLY DISTRIBUTED RANDOM NUMBERS ON THE OPEN C INTERVAL (0, 1) GIVEN ANY ONE ARGUMENT OF ANY TYPE. C SUBROUTINE TRED2 C IS THE EISPACK (4) PROGRAM WHICH COMPUTES A HOUSEHOLDER C TRIDIAGONAL FORM OF A REAL SYMMETRIC MATRIX. C SUBROUTINE IMTQL2 C IS THE EISPACK PROGRAM WHICH COMPUTES THE EIGENVALUES AND C ORTHONORMAL EIGENVECTORS OF A SYMMETRIC TRIDIAGONAL MATRIX. C FUNCTION IP C IS DESCRIBED ABOVE. C SUBROUTINE OP C IS DESCRIBED ABOVE. C SUBROUTINE INF C IS DESCRIBED ABOVE. C METHOD C SIMITZ REPRESENTS RESULTS OF EXTENSIVE MODIFICATIONS AND TESTS C OF SUBROUTINE RITZIT (1), AN ANSI FORTRAN TRANSLATION OF THE C ALGOL 60 PROCEDURE OF THE SAME NAME (3). THE BASIC RUTISHAUSER C -REINSCH ALGORITHM IS PRESERVED. C REFERENCES C (1) PAUL J. NIKOLAI AND NAI-KUAN TSAO, THE ARL LINEAR ALGEBRA C LIBRARY HANDBOOK, ARL TR 74-0106, AEROSPACE RESEARCH LABOR- C ATORIES, WRIGHT-PATTERSON AFB, OHIO, 1974. C (2) HEINZ RUTISHAUSER, COMPUTATIONAL ASPECTS OF F.L. BAUER S C SIMULTANEOUS ITERATION METHOD, NUMER. MATH. 13(1969), 4-13. C (3) -----------------, SIMULTANEOUS ITERATION METHOD FOR SYM- C METRIC MATRICES, NUMER. MATH. 16(1970), 205-223. C (4) B.T. SMITH ET AL, MATRIX EIGENSYSTEM ROUTINES-EISPACK C GUIDE, LECTURE NOTES IN COMPUTER SCIENCE 6, SPRINGER-VERLAG C NEW YORK, 1974. C C PLEASE REFER QUESTIONS, COMMENTS OR SUGGESTIONS TO C C PAUL J. NIKOLAI C AFFDL/FBR C WRIGHT-PATTERSON AFB, OH 45433 C (513)-255-5350 C C*********************************************************************** EXTERNAL INF, IP, OP INTEGER EM, G, H, I, IG, IK, J, * JK, JP, K, KM, KS, L, LF, * L1, M, MN, M1, N, P, Z1, * Z2 LOGICAL ORIG REAL D, E, EPS, E1, E2, IP, MT, * S, T, WK, X DIMENSION X(MN,1), D(1), WK(P,P,1) C DATA MT /.220360641585062E-279/ C C THE LOCAL VARIABLE ARRAYS FROM (3) ARE ASSIGNED TO THE C VARIABLE ARRAY WK AS FOLLOWS C C WK(I,J,1) = B(I,J), I, J = 1, ..., P. C WK(I,1,2) = CX(I), I = 1, ..., P. C WK(I,2,2) = F(I), I = 1, ..., P. C WK(I,3,2) = RQ(I), I = 1, ..., P. C WK(I,4,2) = U(I), I = 1, ..., N. C WK(I+N,4,2) = W(I), I = 1, ..., N. C NOT NEEDED ARE V(I), I = 1, ..., N, R(I), I = 1, ..., P, AND C Q(I,J), I, J = 1, ..., P. C C THE NEXT STATEMENT IS START. E = .0E+00 G = 0 IG = 1 H = 0 Z1 = 0 Z2 = 0 KS = 0 M = 1 WK(P,1,2) = .0E+00 DO 10 L = 1, P WK(L,2,2) = .4E+01 WK(L,3,2) = .0E+00 10 CONTINUE IF (KM) 50, 50, 20 20 DO 40 L = 1, P DO 30 J = 1, N X(J,L) = .2E+01*RANF(.2E+01) - .1E+01 30 CONTINUE 40 CONTINUE 50 KM = IABS(KM) ASSIGN 60 TO IK LF = IG L1 = P GO TO 990 C RAYLEIGH-RITZ STEP C STATEMENT 60 IS LOOP. 60 DO 80 K = IG, P CALL OP(N, X(1,K), WK(1,4,2)) DO 70 J = 1, N X(J,K) = WK(J,4,2) 70 CONTINUE 80 CONTINUE ASSIGN 90 TO IK LF = IG L1 = P GO TO 990 90 IF (KS) 150, 100, 150 C MEASURES AGAINST UNHAPPY CHOICE OF INITIAL VECTORS 100 DO 130 K = 1, P IF (WK(K,K,1)) 130, 110, 130 110 DO 120 I = 1, N X(I,K) = .2E+01*RANF(.2E+01) - .1E+01 120 CONTINUE KS = 1 130 CONTINUE IF (KS - 1) 150, 140, 150 140 ASSIGN 60 TO IK LF = 1 L1 = P GO TO 990 150 DO 180 K = IG, P DO 170 L = K, P S = .0E+00 DO 160 I = L, P S = S + WK(I,K,1)*WK(I,L,1) 160 CONTINUE WK(L,K,1) = -S 170 CONTINUE 180 CONTINUE CALL TRED2(P, P - G, WK(IG,IG,1), D(IG), WK(1,4,2), WK(IG,IG,1)) CALL IMTQL2(P, P - G, D(IG), WK(1,4,2), WK(IG,IG,1), L) WK(IG,2,2) = AMAX1(WK(IG,2,2), .1E+04*FLOAT(L)) DO 190 K = IG, P D(K) = SQRT(AMAX1(-D(K), .0E+00)) 190 CONTINUE C REORDERING EIGENVALUES AND EIGENVECTORS ACCORDING TO SIZE OF C THE FORMER IS ACCOMPLISHED IN SUBROUTINE IMTQL2. DO 230 J = 1, N DO 210 K = IG, P S = .0E+00 DO 200 L = IG, P S = S + X(J,L)*WK(L,K,1) 200 CONTINUE WK(K,4,2) = S 210 CONTINUE DO 220 K = IG, P X(J,K) = WK(K,4,2) 220 CONTINUE 230 CONTINUE KS = KS + 1 E = AMAX1(D(P), E) C RANDOMIZATION IF (3 - Z1) 260, 240, 240 240 DO 250 J = 1, N X(J,P) = .2E+01*RANF(.2E+01) - .1E+01 250 CONTINUE JP = P - 1 ASSIGN 260 TO IK LF = P L1 = P GO TO 990 C COMPUTE CONTROL QUANTITIES CX(I). 260 DO 310 K = IG, JP S = (D(K) - E)*(D(K) + E) IF (S) 270, 270, 280 270 WK(K,1,2) = .0E+00 GO TO 310 280 IF (E) 300, 290, 300 290 WK(K,1,2) = .1E+04 + ALOG(D(K)) GO TO 310 300 WK(K,1,2) = ALOG((D(K) + SQRT(S))/E) 310 CONTINUE C ACCEPTANCE TEST FOR EIGENVALUES INCLUDING ADJUSTMENT OF EM AND C H SUCH THAT D(EM) .GT. E, D(H) .GT. E AND D(EM) DOES NOT C OSCILLATE STRONGLY I = Z1 - 1 K = G 320 K = K + 1 IF (EM - K) 370, 330, 330 330 IF (D(K) - E) 360, 360, 340 340 IF (I) 320, 320, 350 350 IF (D(K) - .999E+00*WK(K,3,2)) 360, 360, 320 360 CONTINUE EM = K - 1 C STATEMENT 370 IS EX4. 370 IF (EM) 380, 1130, 380 380 K = H S = .1E+01 + .1E+00*EPS 390 K = K + 1 IF (D(K)) 400, 410, 400 400 IF (D(K) - S*WK(K,3,2)) 390, 390, 410 410 CONTINUE H = K - 1 K = EM 420 K = K + 1 IF (K - H) 430, 430, 450 430 IF (D(K) - E) 440, 440, 420 440 CONTINUE H = K - 1 C ACCEPTANCE TEST FOR EIGENVECTORS 450 L = G E2 = .0E+00 DO 590 K = IG, JP IF (K - (L + 1)) 510, 460, 510 C CHECK FOR NESTED EIGENVALUES 460 L = K L1 = K S = .5E+00/FLOAT(KS) T = .1E+01/FLOAT(KS*M) 470 L = L + 1 IF (L - JP) 480, 480, 490 480 IF (WK(L,1,2)*(WK(L,1,2) + S) + T - WK(L-1,1,2)*WK(L-1,1,2)) * 490, 490, 470 490 CONTINUE L = L - 1 C THE NEXT STATEMENT IS EX5. IF (L - H) 510, 510, 500 500 L = L1 - 1 GO TO 600 510 CALL OP(N, X(1,K), WK(1,4,2)) S = .0E+00 DO 540 J = 1, L IF (ABS(D(J) - D(K)) - .1E-01*D(K)) 520, 540, 540 520 T = IP(N, WK(1,4,2), X(1,J)) DO 530 I = 1, N WK(I,4,2) = WK(I,4,2) - T*X(I,J) 530 CONTINUE S = S + T*T 540 CONTINUE T = .1E+01 IF (S .NE. .0E+00) T = IP(N, WK(1,4,2), WK(1,4,2)) E2 = AMAX1(E2, SQRT(T/(S + T))) IF (K - L) 590, 550, 590 C TEST FOR ACCEPTANCE OF GROUP OF EIGENVECTORS 550 IF (L .GE. EM .AND. D(EM)*WK(EM,2,2) .LT. EPS*(D(EM) - E)) * G = EM IF (E2 - WK(L,2,2)) 560, 580, 580 560 DO 570 J = L1, L WK(J,2,2) = E2 570 CONTINUE 580 IF (L .LE. EM .AND. D(L)*WK(L,2,2) .LT. EPS*(D(L) - E)) G = L 590 CONTINUE C ADJUST M. C STATEMENT 600 IS EX6. 600 IG = G + 1 IF (E - .4E-01*D(1)) 610, 610, 620 610 M = 1 K = 1 GO TO 630 620 E2 = .2E+01/E E1 = .51E+00*E2 K = 2*INT(.4E+01/AMIN1(WK(1,1,2), .4E+01)) M = MIN0(M, K) C REDUCE EM IF CONVERGENCE WOULD BE TOO SLOW. 630 IF (WK(EM,2,2)) 640, 690, 640 640 IF (FLOAT(KS) - .9E+00*FLOAT(KM)) 650, 690, 690 650 S = FLOAT(K)*WK(EM,1,2) IF (S - .5E-01) 660, 670, 670 660 T = .5E+00*S*WK(EM,1,2) GO TO 680 670 T = WK(EM,1,2) + ALOG(.5E+00 + .5E+00*EXP(-.2E+01*S))/FLOAT(K) 680 S = ALOG(D(EM)*WK(EM,2,2)/(EPS*(D(EM) - E))) IF (S*FLOAT(KS) .GT. T*FLOAT((KM - KS)*KM)) EM = EM - 1 C STATEMENT 690 IS EX2. 690 DO 700 K = IG, JP WK(K,3,2) = D(K) 700 CONTINUE CALL INF(KS, G, H, WK(1,2,2)) IF (G .GE. EM .OR. KS .GE. KM) GO TO 1130 C STATEMENT 710 IS EX1. 710 IF (KS + M - KM) 730, 730, 720 720 Z2 = -1 IF (M .GT. 1) M = 2*((KM - KS + 1)/2) 730 M1 = M C SHORTCUT LAST INTERMEDIATE BLOCK IF ALL F(I) ARE SUFFICIENTLY C SMALL. IF (L - EM) 780, 740, 740 740 S = D(EM)*WK(EM,2,2)/(EPS*(D(EM) - E)) T = S*S - .1E+01 IF (T) 60, 60, 750 750 S = ALOG(S + SQRT(T))/(WK(EM,1,2) - WK(H+1,1,2)) M1 = 2*INT(.5E+00*S + .101E+01) IF (M1 - M) 770, 770, 760 760 M1 = M GO TO 780 770 Z2 = -1 C CHEBYSHEV ITERATION 780 IF (M - 1) 900, 790, 820 790 DO 810 K = IG, P CALL OP(N, X(1,K), WK(1,4,2)) DO 800 I = 1, N X(I,K) = WK(I,4,2) 800 CONTINUE 810 CONTINUE GO TO 900 820 L1 = M1 - 4 DO 890 K = IG, P CALL OP(N, X(1,K), WK(1,4,2)) DO 830 I = 1, N IK = I + N WK(IK,4,2) = E1*WK(I,4,2) 830 CONTINUE CALL OP(N, WK(N+1,4,2), WK(1,4,2)) DO 840 I = 1, N X(I,K) = E2*WK(I,4,2) - X(I,K) 840 CONTINUE IF (L1) 890, 850, 850 850 DO 880 J = 4, M1, 2 CALL OP(N, X(1,K), WK(1,4,2)) DO 860 I = 1, N IK = I + N WK(IK,4,2) = E2*WK(I,4,2) - WK(IK,4,2) 860 CONTINUE CALL OP(N, WK(N+1,4,2), WK(1,4,2)) DO 870 I = 1, N X(I,K) = E2*WK(I,4,2) - X(I,K) 870 CONTINUE 880 CONTINUE 890 CONTINUE 900 ASSIGN 910 TO IK LF = IG L1 = P GO TO 990 C DISCOUNTING THE ERROR QUANTITIES F 910 IF (G - H) 920, 970, 970 920 IF (M - 1) 950, 930, 950 930 DO 940 K = IG, H WK(K,2,2) = WK(K,2,2)*(D(H+1)/D(K)) 940 CONTINUE GO TO 970 950 T = EXP(-FLOAT(M1)*WK(H+1,1,2)) DO 960 K = IG, H S = EXP(-FLOAT(M1)*(WK(K,1,2) - WK(H+1,1,2))) WK(K,2,2) = S*WK(K,2,2)*(.1E+01 + T*T)/(.1E+01 + (S*T)*(S*T)) 960 CONTINUE 970 KS = KS + M1 Z2 = Z2 - M1 C POSSIBLE REPETITION OF INTERMEDIATE STEPS IF (Z2) 980, 710, 710 980 Z1 = Z1 + 1 Z2 = 2*Z1 M = M + M GO TO 60 C PERFORMS ORTHONORMALIZATION OF COLUMNS 1 THROUGH L1 OF ARRAY C X ASSUMING THAT COLUMNS 1 THROUGH LF - 1 ARE ALREADY ORTHO- C NORMAL 990 DO 1120 K = LF, L1 ORIG = .TRUE. 1000 T = .0E+00 JK = K - 1 IF (JK) 1040, 1040, 1010 1010 DO 1030 I = 1, JK S = IP(N, X(1,I), X(1,K)) IF (ORIG) WK(K,I,1) = S T = T + S*S DO 1020 J = 1, N X(J,K) = X(J,K) - S*X(J,I) 1020 CONTINUE 1030 CONTINUE 1040 S = IP(N, X(1,K), X(1,K)) T = S + T IF (S - .1E-01*T) 1060, 1060, 1050 1050 IF (T - MT) 1060, 1060, 1080 1060 ORIG = .FALSE. IF (S - MT) 1070, 1070, 1000 1070 S = .0E+00 1080 S = SQRT(S) WK(K,K,1) = S IF (S) 1090, 1100, 1090 1090 S = .1E+01/S 1100 DO 1110 J = 1, N X(J,K) = S*X(J,K) 1110 CONTINUE 1120 CONTINUE GO TO IK, (60, 90, 260, 910, 1140) C STATEMENT 1130 IS EX. 1130 EM = G C SOLVE EIGENVALUE PROBLEM OF PROJECTION OF MATRIX C. ASSIGN 1140 TO IK LF = 1 L1 = JP GO TO 990 1140 DO 1160 K = 1, JP CALL OP(N, X(1,K), X(1,P)) DO 1150 I = 1, K WK(K,I,1) = -IP(N, X(1,I), X(1,P)) 1150 CONTINUE 1160 CONTINUE CALL TRED2(P, JP, WK, D, WK(1,4,2), WK) CALL IMTQL2(P, JP, D, WK(1,4,2), WK, L) WK(IG,2,2) = AMAX1(WK(IG,2,2), .1E+04*FLOAT(L)) C ARRANGE EIGENVALUES IN ORDER OF DECREASING ABSOLUTE VALUE. DO 1210 J = 1, JP K = J DO 1170 I = J, JP IF (ABS(D(I)) .GT. ABS(D(K))) K = I 1170 CONTINUE IF (K - J) 1200, 1200, 1180 1180 T = D(K) D(K) = D(J) D(J) = T DO 1190 I = 1, JP T = WK(I,K,1) WK(I,K,1) = WK(I,J,1) WK(I,J,1) = T 1190 CONTINUE 1200 D(J) = -D(J) 1210 CONTINUE DO 1250 J = 1, N DO 1230 I = 1, JP S = .0E+00 DO 1220 K = 1, JP S = S + X(J,K)*WK(K,I,1) 1220 CONTINUE WK(I,4,2) = S 1230 CONTINUE DO 1240 I = 1, JP X(J,I) = WK(I,4,2) 1240 CONTINUE 1250 CONTINUE KM = KS D(P) = E RETURN END C PROGRAM TESTB(INPUT, OUTPUT, TAPE5 = INPUT, TAPE6 = OUTPUT) 5510 DIMENSION A(50,6), T(50,6), Y(50) 5520 COMMON NT, T, NA, A, ND, Y 5530 INTEGER Y 5540 DIMENSION X(50,10), D(10), WK(115,2) 5550 EXTERNAL INTROS, ABAND, BAND 5560 C 5570 C THE VARIABLE ND MUST IDENTIFY A VALUE AT LEAST EQUAL TO N. 5580 C THE VALUE ASSIGNED TO THE VARIABLE ND MUST AGREE WITH THE LEADING 5590 C DIMENSION OF THE VARIABLES A, T, X, AND Y. THE LEADING DIMENSION OF 5600 C THE VARIABLE WK MUST EQUAL OR EXCEED ND + 65. THE SECOND DIMENSION 5610 C OF A AND T MUST EXCEED NA AND NT, RESPECTIVELY. OTHERWISE RECOMPILE. 5620 C NOTE THAT DIMENSIONS OF THE VARIABLES A, T, AND Y MUST AGREE WITH 5630 C THEIR COUNTERPARTS IN FUNCTION BAND AND SUBROUTINE ABAND. 5640 C 5650 C N IS AN INTEGER INPUT VARIABLE, THE ORDER OF THE MATRICES A AND 5660 C B. IF N .LE. 0, PROCESSING CEASES. 5670 C NA IS AN INTEGER INPUT VARIABLE, THE NUMBER OF ADJACENT NON-ZERO 5680 C DIAGONALS STRICTLY BELOW THE MAIN DIAGONAL OF THE MATRIX A. 5690 C NT IS AN INTEGER INPUT VARIABLE, THE NUMBER OF ADJACENT NON-ZERO 5700 C DIAGONALS STRICTLY BELOW THE MAIN DIAGONAL OF THE MATRIX T. 5710 C KM IS AN INTEGER INPUT VARIABLE, THE MAXIMUM NUMBER OF ITERATION 5720 C STEPS TO BE EXECUTED BY SIMITZ ON EACH CALL. IF KM IDENTIFIES 5730 C A NEGATIVE VALUE, SIMITZ WILL USE PREVIOUSLY COMPUTED EIGENVEC 5740 C TORS AS INITIAL APPROXIMATIONS FOR A GIVEN VALUE OF P WHILE EM 5750 C IS INCREASED. 5760 C EPS IS A REAL INPUT VARIABLE, THE TOLERANCE FOR ACCEPTING EIGEN- 5770 C VECTORS. 5780 C SEED IS A REAL INPUT VARIABLE TO INITIALIZE THE RANDOM NUMBER GEN- 5790 C ERATOR RANF USING THE INITIALIZING SUBROUTINE RANSET. IF SEED 5800 C IDENTIFIES A NEGATIVE VALUE, NO PRINTING OF THE NON-ZERO 5810 C DIAGONALS WILL BE DONE. 5820 C TRANSA IS A REAL INPUT VARIABLE, THE LEFT ENDPOINT OF THE INTERVAL 5830 C FROM WHICH THE NON-ZERO ENTRIES OF A AND T ARE SELECTED. 5840 C DIALA IS A REAL INPUT VARIABLE, A POSITIVE NUMBER SELECTED SO THAT 5850 C TRANSA + DIALA IS THE RIGHT ENDPOINT OF THE INTERVAL OVER 5860 C WHICH THE NON-ZERO ENTRIES OF A AND T ARE SELECTED. 5870 C 5880 ND = 50 5890 10 READ (5, 20) N, NA, NT, KM, EPS, SEED, TRANSA, DIALA 5900 20 FORMAT(4I5,1P4E10.2) 5910 IF (N .LE. 0) STOP 5920 CALL RANSET(ABS(SEED)) 5930 MP1 = NT + 1 5940 DO 40 I = 1, N 5950 MAX = MAX0(1, -I + (MP1 + 1)) 5960 DO 30 J = MAX, MP1 5970 T(I,J) = AINT(TRANSA + DIALA*RANF(1)) 5980 30 CONTINUE 5990 T(I,MP1) = AMAX1(ABS(T(I,MP1)), ABS(AINT(TRANSA + DIALA)) - 6000 * ABS(T(I,MP1))) 6010 40 CONTINUE 6020 WRITE (6, 50) N, NA, NT, KM, EPS, SEED, TRANSA, DIALA 6030 50 FORMAT(1H1,6X,1HN,11X,2HNA,11X,2HNB,11X,2HKM,18X,3HEPS,17X,4HSEED, 6040 * 10X,11HTRANSLATION,13X,8HDILATION/1H /I8,3I13,1P4E21.2/1H0) 6050 IF (SEED) 110, 60, 60 6060 60 WRITE (6, 70) 6070 70 FORMAT(1H0,36X,59HADJACENT NON-ZERO LOWER DIAGONALS OF T, T(T-TRAN 6080 *SPOSED) = B) 6090 MAX = MP1 6100 DO 90 J = 1, MP1 6110 DO 80 I = J, N 6120 Y(I) = IFIX(T(I,MAX)) 6130 80 CONTINUE 6140 WRITE (6, 100) (Y(I), I = J, N) 6150 MAX = MAX - 1 6160 90 CONTINUE 6170 100 FORMAT(1H /(16X,15I7)) 6180 110 MP1 = NA + 1 6190 DO 130 I = 1, N 6200 MAX = MAX0(1, -I + (MP1 + 1)) 6210 DO 120 J = MAX, MP1 6220 A(I,J) = AINT(TRANSA + DIALA*RANF(1)) 6230 120 CONTINUE 6240 130 CONTINUE 6250 IF (SEED) 180, 140, 140 6260 140 WRITE (6, 150) 6270 150 FORMAT(1H0/39X,53HADJACENT NON-ZERO LOWER DIAGONALS OF A = A-TRANS 6280 *POSED) 6290 MAX = MP1 6300 DO 170 J = 1, MP1 6310 DO 160 I = J, N 6320 Y(I) = IFIX(A(I,MAX)) 6330 160 CONTINUE 6340 WRITE (6, 100) (Y(I), I = J, N) 6350 MAX = MAX - 1 6360 170 CONTINUE 6370 180 WRITE (6, 190) 6380 190 FORMAT(1H0/1H0,21X,1HP,5X,6HEM(IN),4X,7HEM(OUT),9X,2HKS,10X,1HK, 6390 *6X,12HMAX RESIDUAL,5X,13HMEAN RESIDUAL,8X,10HTIME(SECS)/1H ) 6400 MP = MIN0(N/5, 10) 6410 IF (KM) 192, 198, 198 6420 192 DO 196 J = 1, MP 6430 DO 194 I = 1, N 6440 X(I,J) = 2.0*RANF(2.0) - 1.0 6450 194 CONTINUE 6460 196 CONTINUE 6470 198 DO 280 IP = 2, MP 6480 MAXM = IP - 1 6490 DO 270 M = 1, MAXM 6500 MM = M 6510 KS = KM 6520 T1 = SECOND(S) 6530 CALL SIMITZ(N, IP, KS, EPS, BAND, ABAND, INTROS, MM, X, ND, D, 6540 * WK) 6550 T1 = SECOND(S) - T1 6560 IF (MM) 200, 200, 210 6570 200 L = 0 6580 ANMAX = 1.0 6590 ANMEAN = 1.0 6600 GO TO 250 6610 210 PROD = 1.0 6620 ANMAX = 0.0 6630 L = 1 6640 220 DO 240 J = 1, MM 6650 CALL BANDOP(N, NA, A, ND, X(1,J), WK(1,1)) 6660 CALL TANDOP(N, NT, T, ND, X(1,J), WK(1,2)) 6670 S = 0.0 6680 DO 230 I = 1, N 6690 S = S + (WK(I,1) - D(J)*WK(I,2))*(WK(I,1) - D(J)*WK(I,2)) 6700 230 CONTINUE 6710 PROD = PROD*S 6720 ANMAX = AMAX1(ANMAX, S) 6730 L = MAX1(FLOAT(L), SIGN(FLOAT(J), S - ANMAX)) 6740 240 CONTINUE 6750 S = ABS(D(1)) - D(IP) 6760 ANMAX = SQRT(ANMAX)/S 6770 ANMEAN = SQRT(PROD**(1.0/FLOAT(MM)))/S 6780 250 WRITE (6,260) IP, M, MM, KS, L, ANMAX, ANMEAN, T1 6790 260 FORMAT(I23,4I11,1P2E18.2,0PF18.2) 6800 270 CONTINUE 6810 280 CONTINUE 6820 WRITE (6,290) (D(J), J = 1, MAXM) 6830 290 FORMAT(1H0/1H0,57X,17HFINAL EIGENVALUES/1H /(1PE37.12,4E21.12)) 6840 GO TO 10 6850 END 6860 SUBROUTINE INTROS(KS, G, H, F) 6870 C INTROS IS A DUMMY SUBROUTINE INF. INTEGER G, H REAL F(1) RETURN END SUBROUTINE BANDOP(N, NA, A, ND, V, W) 6930 C CALCULATE THE IMAGE W OF THE N-VECTOR V UNDER C THE MATRIX A WHERE A IS SYMMETRIC AND BANDED OF BAND WIDTH 2*NA + 1. DIMENSION A(ND,1), V(1), W(1) INTEGER P, Q, R MP1 = NA + 1 DO 70 I = 1, N P = MAX0(1, -I + (MP1 + 1)) Q = MAX0(1, I + (MP1 - N)) T = A(I,MP1)*V(I) IF (NA - P) 30, 10, 10 10 R = I - (MP1 - P) DO 20 K = P, NA T = T + A(I,K)*V(R) R = R + 1 20 CONTINUE 30 IF (NA - Q) 60, 40, 40 40 R = I + (MP1 - Q) DO 50 K = Q, NA T = T + A(R,K)*V(R) R = R - 1 50 CONTINUE 60 W(I) = T 70 CONTINUE RETURN END SUBROUTINE TANDOP(N, NA, A, ND, V, W) 7190 C COMPUTE THE IMAGE W OF THE N-VECTOR V UNDER THE N-SQUARE MATRIX C A(A-TRANSPOSED), A LOWER TRIANGULAR OF BANDWIDTH NA + 1. DIMENSION A(ND,1), V(1), W(1) INTEGER P, Q, R MP1 = NA + 1 DO 40 I = 1, N P = MAX0(1, I + (MP1 - N)) T = A(I,MP1)*V(I) IF (NA - P) 30, 10, 10 10 R = I + (MP1 - P) DO 20 K = P, NA T = T + A(R,K)*V(R) R = R - 1 20 CONTINUE 30 W(I) = T 40 CONTINUE I = N DO 80 Q = 1, N P = MAX0(1, -I + (MP1 + 1)) T = A(I,MP1)*W(I) IF (NA - P) 70, 50, 50 50 R = I - (MP1 - P) DO 60 K = P, NA T = T + A(I,K)*W(R) R = R + 1 60 CONTINUE 70 W(I) = T I = I - 1 80 CONTINUE RETURN END FUNCTION BAND(N, Z, W) 7510 C CALCULATE THE INNER PRODUCT (W-TRANSPOSED)BZ WHERE B = C T(T-TRANSPOSED). DIMENSION A(50,6), T(50,6), Y(50) COMMON NT, T, NA, A, ND, Y DIMENSION W(1), Z(1) 10 CALL TANDOP(N, NT, T, ND, Z, Y) S = 0.0 DO 15 I = 1, N S = S + W(I)*Y(I) 15 CONTINUE BAND = S RETURN END SUBROUTINE ABAND(N, Z, W) 7650 C CALCULATE THE SOLUTION W OF THE LINEAR SYSTEM BW = AZ GIVEN THE C N-VECTOR Z, THE BANDED SYMMETRIC N-SQUARE MATRIX A OF BAND WIDTH C 2*NA + 1, AND THE BANDED TRIANGULAR FACTOR T OF BAND WIDTH NT + 1 C OF THE POSITIVE DEFINITE MATRIX B, T(T-TRANSPOSED) = B. DIMENSION W(1), Z(1) DIMENSION A(50,6), T(50,6), Y(50) COMMON NT, T, NA, A, ND, Y CALL BANDOP(N, NA, A, ND, Z, W) CALL SOLVE(N, NT, T, ND, W) RETURN END SUBROUTINE SOLVE(N, M, A, N1, B) 7770 C SOLVE A LINEAR SYSTEM GIVEN A TRIANGULAR FACTORIZATION OF A BANDED C POSITIVE DEFINITE COEFFICIENT MATRIX OF BAND WIDTH 2*M + 1. DIMENSION A(N1,1), B(N) INTEGER P, Q, R MP1 = M + 1 COMMENT SOLUTION OF LY = B DO 130 I = 1, N P = MAX0(1, MP1 - I + 1) T = B(I) IF (M - P) 120, 100, 100 100 Q = I - MP1 + P DO 110 K = P, M T = T - A(I,K)*B(Q) Q = Q + 1 110 CONTINUE 120 B(I) = T/A(I,MP1) 130 CONTINUE COMMENT SOLUTION OF UX = Y I = N DO 170 R = 1, N P = MAX0(1, MP1 - N + I) T = B(I) IF (M - P) 160, 140, 140 140 Q = I + MP1 - P DO 150 K = P, M T = T - A(Q,K)*B(Q) Q = Q - 1 150 CONTINUE 160 B(I) = T/A(I,MP1) I = I - 1 170 CONTINUE RETURN END FUNCTION RANF(X) 8110 C RETURNS A RANDOM NUMBER ON THE OPEN UNIT INTERVAL GIVEN ANY ARGUMENT. COMMON /TSEED/ K K = MOD(3125*K, 65536) RANF = ABS(FLOAT(K - 32768)/32768.0) RETURN END SUBROUTINE RANSET(SEED) 8180 C INITIALIZES THE RANDOM NUMBER SEED K USING THE REAL INPUT C VARIABLE SEED. COMMON /TSEED/ K T = SEED IF (ABS(SEED) .GT. 1.0) T = 1.0/T K = 2*IFIX(16384.0*T) - IFIX(SIGN(1.0, T)) RETURN END FUNCTION SECOND(T) 8270 C DUMMY FUNCTION SECOND RETURNS THE VALUES SECOND = T = O.O. C C SECOND MAY BE MODIFIED TO CALL AN OPERATING SYSTEM DEPENDENT PROGRAM C WHICH RETURNS CENTRAL PROCESSOR TIME IN SECONDS RELATIVE TO C AN ARBITRARY STARTING POINT AS A REAL VALUE T. C T = 0.0 SECOND = 0.0 RETURN END C 8380 C ------------------------------------------------------------------ 8390 C 8400 SUBROUTINE IMTQL2(NM,N,D,E,Z,IERR) 8410 C INTEGER I,J,K,L,M,N,II,NM,MML,IERR REAL D(N),E(N),Z(NM,N) REAL B,C,F,G,P,R,S,MACHEP C REAL SQRT,ABS,SIGN C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE IMTQL2, C NUM. MATH. 12, 377-383(1968) BY MARTIN AND WILKINSON, C AS MODIFIED IN NUM. MATH. 15, 450(1970) BY DUBRULLE. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 241-248(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE IMPLICIT QL METHOD. C THE EIGENVECTORS OF A FULL SYMMETRIC MATRIX CAN ALSO C BE FOUND IF TRED2 HAS BEEN USED TO REDUCE THIS C FULL MATRIX TO TRIDIAGONAL FORM. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX, C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX C IN ITS LAST N-1 POSITIONS. E(1) IS ARBITRARY, C C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE C REDUCTION BY TRED2, IF PERFORMED. IF THE EIGENVECTORS C OF THE TRIDIAGONAL MATRIX ARE DESIRED, Z MUST CONTAIN C THE IDENTITY MATRIX. C C ON OUTPUT- C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT BUT C UNORDERED FOR INDICES 1,2,...,IERR-1, C C E HAS BEEN DESTROYED, C C Z CONTAINS ORTHONORMAL EIGENVECTORS OF THE SYMMETRIC C TRIDIAGONAL (OR FULL) MATRIX. IF AN ERROR EXIT IS MADE, C Z CONTAINS THE EIGENVECTORS ASSOCIATED WITH THE STORED C EIGENVALUES, C C IERR IS SET TO C ZERO FOR NORMAL RETURN, C J IF THE J-TH EIGENVALUE HAS NOT BEEN C DETERMINED AFTER 30 ITERATIONS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C C ********** MACHEP IS A MACHINE DEPENDENT PARAMETER SPECIFYING C THE RELATIVE PRECISION OF FLOATING POINT ARITHMETIC. C C ********** MACHEP = 2.**(-47) C IERR = 0 IF (N .EQ. 1) GO TO 1001 C DO 100 I = 2, N 100 E(I-1) = E(I) C E(N) = 0.0 C DO 240 L = 1, N J = 0 C ********** LOOK FOR SMALL SUB-DIAGONAL ELEMENT ********** 105 DO 110 M = L, N IF (M .EQ. N) GO TO 120 IF (ABS(E(M)) .LE. MACHEP * (ABS(D(M)) + ABS(D(M+1)))) X GO TO 120 110 CONTINUE C 120 P = D(L) IF (M .EQ. L) GO TO 240 IF (J .EQ. 30) GO TO 1000 J = J + 1 C ********** FORM SHIFT ********** G = (D(L+1) - P) / (2.0 * E(L)) R = SQRT(G*G+1.0) G = D(M) - P + E(L) / (G + SIGN(R,G)) S = 1.0 C = 1.0 P = 0.0 MML = M - L C ********** FOR I=M-1 STEP -1 UNTIL L DO -- ********** DO 200 II = 1, MML I = M - II F = S * E(I) B = C * E(I) IF (ABS(F) .LT. ABS(G)) GO TO 150 C = G / F R = SQRT(C*C+1.0) E(I+1) = F * R S = 1.0 / R C = C * S GO TO 160 150 S = F / G R = SQRT(S*S+1.0) E(I+1) = G * R C = 1.0 / R S = S * C 160 G = D(I+1) - P R = (D(I) - G) * S + 2.0 * C * B P = S * R D(I+1) = G + P G = C * R - B C ********** FORM VECTOR ********** DO 180 K = 1, N F = Z(K,I+1) Z(K,I+1) = S * Z(K,I) + C * F Z(K,I) = C * Z(K,I) - S * F 180 CONTINUE C 200 CONTINUE C D(L) = D(L) - P E(L) = G E(M) = 0.0 GO TO 105 240 CONTINUE C ********** ORDER EIGENVALUES AND EIGENVECTORS ********** DO 300 II = 2, N I = II - 1 K = I P = D(I) C DO 260 J = II, N IF (D(J) .GE. P) GO TO 260 K = J P = D(J) 260 CONTINUE C IF (K .EQ. I) GO TO 300 D(K) = D(I) D(I) = P C DO 280 J = 1, N P = Z(J,I) Z(J,I) = Z(J,K) Z(J,K) = P 280 CONTINUE C 300 CONTINUE C GO TO 1001 C ********** SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS ********** 1000 IERR = L 1001 RETURN C ********** LAST CARD OF IMTQL2 ********** END C 10020 C ------------------------------------------------------------------ 10030 C 10040 SUBROUTINE TRED2(NM,N,A,D,E,Z) 10050 C INTEGER I,J,K,L,N,II,NM,JP1 REAL A(NM,N),D(N),E(N),Z(NM,N) REAL F,G,H,HH,SCALE C REAL SQRT,ABS,SIGN C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED2, C NUM. MATH. 11, 181-195(1968) BY MARTIN, REINSCH, AND WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971). C C THIS SUBROUTINE REDUCES A REAL SYMMETRIC MATRIX TO A C SYMMETRIC TRIDIAGONAL MATRIX USING AND ACCUMULATING C ORTHOGONAL SIMILARITY TRANSFORMATIONS. C C ON INPUT- C C NM MUST BE SET TO THE ROW DIMENSION OF TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT, C C N IS THE ORDER OF THE MATRIX, C C A CONTAINS THE REAL SYMMETRIC INPUT MATRIX. ONLY THE C LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED. C C ON OUTPUT- C C D CONTAINS THE DIAGONAL ELEMENTS OF THE TRIDIAGONAL MATRIX, C C E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE TRIDIAGONAL C MATRIX IN ITS LAST N-1 POSITIONS. E(1) IS SET TO ZERO, C C Z CONTAINS THE ORTHOGONAL TRANSFORMATION MATRIX C PRODUCED IN THE REDUCTION, C C A AND Z MAY COINCIDE. IF DISTINCT, A IS UNALTERED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO B. S. GARBOW, C APPLIED MATHEMATICS DIVISION, ARGONNE NATIONAL LABORATORY C C ------------------------------------------------------------------ C DO 100 I = 1, N C DO 100 J = 1, I Z(I,J) = A(I,J) 100 CONTINUE C IF (N .EQ. 1) GO TO 320 C ********** FOR I=N STEP -1 UNTIL 2 DO -- ********** DO 300 II = 2, N I = N + 2 - II L = I - 1 H = 0.0 SCALE = 0.0 IF (L .LT. 2) GO TO 130 C ********** SCALE ROW (ALGOL TOL THEN NOT NEEDED) ********** DO 120 K = 1, L 120 SCALE = SCALE + ABS(Z(I,K)) C IF (SCALE .NE. 0.0) GO TO 140 130 E(I) = Z(I,L) GO TO 290 C 140 DO 150 K = 1, L Z(I,K) = Z(I,K) / SCALE H = H + Z(I,K) * Z(I,K) 150 CONTINUE C F = Z(I,L) G = -SIGN(SQRT(H),F) E(I) = SCALE * G H = H - F * G Z(I,L) = F - G F = 0.0 C DO 240 J = 1, L Z(J,I) = Z(I,J) / H G = 0.0 C ********** FORM ELEMENT OF A*U ********** DO 180 K = 1, J 180 G = G + Z(J,K) * Z(I,K) C JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L 200 G = G + Z(K,J) * Z(I,K) C ********** FORM ELEMENT OF P ********** 220 E(J) = G / H F = F + E(J) * Z(I,J) 240 CONTINUE C HH = F / (H + H) C ********** FORM REDUCED A ********** DO 260 J = 1, L F = Z(I,J) G = E(J) - HH * F E(J) = G C DO 260 K = 1, J Z(J,K) = Z(J,K) - F * E(K) - G * Z(I,K) 260 CONTINUE C 290 D(I) = H 300 CONTINUE C 320 D(1) = 0.0 E(1) = 0.0 C ********** ACCUMULATION OF TRANSFORMATION MATRICES ********** DO 500 I = 1, N L = I - 1 IF (D(I) .EQ. 0.0) GO TO 380 C DO 360 J = 1, L G = 0.0 C DO 340 K = 1, L 340 G = G + Z(I,K) * Z(K,J) C DO 360 K = 1, L Z(K,J) = Z(K,J) - G * Z(K,I) 360 CONTINUE C 380 D(I) = Z(I,I) Z(I,I) = 1.0 IF (L .LT. 1) GO TO 500 C DO 400 J = 1, L Z(I,J) = 0.0 Z(J,I) = 0.0 400 CONTINUE C 500 CONTINUE C RETURN C ********** LAST CARD OF TRED2 ********** END 30 5 4 200 1.00E-10 1.00E-01 -1.00E+01 2.00E+01 11440 30 5 4 -200 1.00E-10 -1.00E-01 -1.00E+01 2.00E+01 11450 40 4 3 200 1.00E-10 2.00E-01 -1.00E+02 2.00E+02 11460 40 4 3 -200 1.00E-10 -2.00E-01 -1.00E+02 2.00E+02 11470 50 2 1 200 1.00E-10 3.00E-01 -1.00E+02 2.00E+02 11480 50 2 1 -200 1.00E-10 -3.00E-01 -1.00E+02 2.00E+02 11490 0 11500 1 N NA NB KM EPS 11510 11520 30 5 4 200 1.00E-10 11530 0 11540 0 ADJACENT NON-ZERO LOWER DI 11550 11560 7 7 7 5 10 9 10 6 11570 5 9 7 7 5 8 7 6 11580 11590 8 5 -5 7 -5 3 -4 8 11600 5 6 -4 -5 -9 -7 2 5 11610 11620 4 5 3 3 -6 9 -6 -6 11630 0 -2 2 7 6 -1 4 6 11640 11650 0 -7 -2 0 -9 4 4 -9 11660 -8 -7 0 4 6 -8 -5 -3 11670 11680 -6 9 6 -4 3 7 2 1 11690 1 -7 -8 8 9 -5 3 3 11700 0 11710 ADJACENT NON-ZERO LOWER DI 11720 11730 0 7 6 -4 -6 5 8 3 11740 -8 1 -4 6 -9 8 5 3 11750 11760 -7 -1 0 4 8 1 -5 0 11770 -4 0 4 4 5 -2 2 -7 11780 11790 -2 -1 -8 7 1 0 -6 6 11800 2 1 9 5 -6 -6 -6 0 11810 11820 0 0 -8 0 -9 0 0 7 11830 6 5 -6 9 -8 5 -4 -9 11840 11850 -7 0 -5 8 0 -6 -9 5 11860 6 -5 2 0 6 -4 -5 -2 11870 11880 -8 0 8 3 9 4 1 9 11890 3 3 -8 0 2 -6 -4 0 11900 0 11910 0 P EM(IN) EM(OUT) KS K 11920 11930 2 1 1 29 1 11940 3 1 1 15 1 11950 3 2 2 24 1 11960 4 1 1 15 1 11970 4 2 2 15 1 11980 4 3 3 51 1 11990 5 1 1 13 1 12000 5 2 2 13 1 12010 5 3 3 21 1 12020 5 4 4 21 1 12030 6 1 1 13 1 12040 6 2 2 13 1 12050 6 3 3 21 1 12060 6 4 4 13 1 12070 6 5 5 40 1 12080 0 12090 0 FINAL EIGENV 12100 12110 4.550109524972E+01 2.427103576410E+01 -7.628352951 12120 1 N NA NB KM EPS 12130 12140 30 5 4 -200 1.00E-10 12150 0 12160 0 12170 0 P EM(IN) EM(OUT) KS K 12180 12190 2 1 1 30 1 12200 3 1 1 3 1 12210 3 2 2 15 1 12220 4 1 1 3 1 12230 4 2 2 3 1 12240 4 3 3 31 1 12250 5 1 1 3 1 12260 5 2 2 3 1 12270 5 3 3 3 1 12280 5 4 4 3 1 12290 6 1 1 3 1 12300 6 2 2 3 1 12310 6 3 3 3 1 12320 6 4 4 3 1 12330 6 5 5 18 1 12340 0 12350 0 FINAL EIGENV 12360 12370 4.550109524972E+01 2.427103576410E+01 -7.628352951 12380 1 N NA NB KM EPS 12390 12400 40 4 3 200 1.00E-10 12410 0 12420 0 ADJACENT NON-ZERO LOWER DI 12430 12440 50 62 67 76 78 67 67 98 12450 82 66 64 84 66 89 76 50 12460 58 61 99 56 100 91 55 83 12470 12480 72 -73 -61 70 -94 -36 -80 29 12490 -3 76 25 -85 78 21 37 3 12500 -15 52 -96 76 -91 26 52 -39 12510 12520 50 16 35 20 73 -50 76 -75 12530 -89 95 -96 60 90 -81 -71 -71 12540 -19 80 40 23 -44 39 35 -91 12550 12560 -76 97 -47 47 36 74 -44 -18 12570 24 34 -61 -86 66 -25 10 16 12580 -3 10 55 25 -74 -19 14 12590 0 12600 ADJACENT NON-ZERO LOWER DI 12610 12620 -53 73 -15 -63 70 39 12 -48 12630 13 -21 78 -34 -17 62 64 -17 12640 85 -13 6 12 -20 -76 -20 -57 12650 12660 66 77 -46 66 -49 -86 40 -14 12670 8 -93 -93 -98 28 51 40 -33 12680 73 27 -14 96 -31 -15 -71 6 12690 12700 64 2 49 46 -60 24 63 -50 12710 54 99 -5 51 -9 -43 -34 -46 12720 15 -7 -2 33 77 41 72 -81 12730 12740 -66 -40 42 -45 10 -97 -57 -14 12750 -48 12 93 0 21 -79 7 3 12760 62 -5 16 -52 82 -28 43 12770 12780 -77 -79 -4 -63 61 31 -3 55 12790 83 8 -44 81 -68 -44 9 -8 12800 -64 8 19 -35 42 99 12810 0 12820 0 P EM(IN) EM(OUT) KS K 12830 12840 2 1 1 13 1 12850 3 1 1 13 1 12860 3 2 2 13 1 12870 4 1 1 7 1 12880 4 2 2 13 1 12890 4 3 3 54 1 12900 5 1 1 7 1 12910 5 2 2 13 1 12920 5 3 3 34 1 12930 5 4 4 170 1 12940 6 1 1 7 1 12950 6 2 2 13 1 12960 6 3 3 30 1 12970 6 4 4 37 1 12980 6 5 5 26 1 12990 7 1 1 7 1 13000 7 2 2 13 1 13010 7 3 3 25 1 13020 7 4 4 27 1 13030 7 5 5 27 1 13040 7 6 5 91 1 13050 8 1 1 7 1 13060 8 2 2 13 1 13070 8 3 3 30 1 13080 8 4 4 37 1 13090 8 5 5 28 1 13100 8 6 6 180 1 13110 8 7 6 177 1 13120 0 13130 0 FINAL EIGENV 13140 13150 -1.059611108016E+03 -2.114191261132E+01 2.825002039 13160 -6.809394216930E-02 -6.243506342913E-02 13170 1 N NA NB KM EPS 13180 13190 40 4 3 -200 1.00E-10 13200 0 13210 0 13220 0 P EM(IN) EM(OUT) KS K 13230 13240 2 1 1 13 1 13250 3 1 1 3 1 13260 3 2 2 3 1 13270 4 1 1 3 1 13280 4 2 2 7 1 13290 4 3 3 40 1 13300 5 1 1 3 1 13310 5 2 2 3 1 13320 5 3 3 18 1 13330 5 4 4 100 1 13340 6 1 1 3 1 13350 6 2 2 3 1 13360 6 3 3 11 1 13370 6 4 4 13 1 13380 6 5 5 13 1 13390 7 1 1 3 1 13400 7 2 2 3 1 13410 7 3 3 7 1 13420 7 4 4 19 1 13430 7 5 5 13 1 13440 7 6 6 187 1 13450 8 1 1 3 1 13460 8 2 2 3 1 13470 8 3 3 7 1 13480 8 4 4 10 1 13490 8 5 5 12 1 13500 8 6 6 55 1 13510 8 7 6 157 1 13520 0 13530 0 FINAL EIGENV 13540 13550 -1.059611108016E+03 -2.114191261132E+01 2.825002039 13560 -6.809394215568E-02 -6.243506344448E-02 13570 1 N NA NB KM EPS 13580 13590 50 2 1 200 1.00E-10 13600 0 13610 0 ADJACENT NON-ZERO LOWER DI 13620 13630 74 79 100 89 95 93 60 59 13640 71 61 74 78 70 95 78 64 13650 51 96 98 71 55 72 84 51 13660 72 95 62 68 82 13670 13680 46 11 -81 16 -40 50 -96 12 13690 59 90 -54 17 88 -97 -18 79 13700 -20 -20 73 2 -16 94 64 86 13710 92 -5 1 -63 13720 0 13730 ADJACENT NON-ZERO LOWER DI 13740 13750 -40 84 -74 54 -91 80 -90 -60 13760 90 74 -97 -46 -43 -95 37 -12 13770 61 -57 -55 13 44 41 -57 -79 13780 -27 77 -98 62 49 13790 13800 -67 0 -19 68 -12 50 53 -44 13810 -6 -63 -4 -86 -42 2 -26 -3 13820 -38 73 35 -83 32 -58 -32 28 13830 73 26 -66 -18 13840 13850 13 -74 -67 -20 -77 53 92 -74 13860 57 61 56 -91 -62 -90 -14 45 13870 62 30 -37 88 25 -70 -33 -99 13880 -34 -47 -78 13890 0 13900 0 P EM(IN) EM(OUT) KS K 13910 13920 2 1 1 46 1 13930 3 1 1 43 1 13940 3 2 2 60 2 13950 4 1 1 30 1 13960 4 2 2 43 1 13970 4 3 3 63 3 13980 5 1 1 26 1 13990 5 2 2 39 2 14000 5 3 3 52 1 14010 5 4 4 61 1 14020 6 1 1 26 1 14030 6 2 2 39 1 14040 6 3 3 39 3 14050 6 4 4 50 2 14060 6 5 5 107 1 14070 7 1 1 24 1 14080 7 2 2 24 1 14090 7 3 3 24 3 14100 7 4 4 33 1 14110 7 5 5 35 3 14120 7 6 6 35 3 14130 8 1 1 24 1 14140 8 2 2 24 1 14150 8 3 3 24 1 14160 8 4 4 24 4 14170 8 5 5 35 4 14180 8 6 6 35 4 14190 8 7 7 114 4 14200 9 1 1 24 1 14210 9 2 2 24 1 14220 9 3 3 24 1 14230 9 4 4 24 1 14240 9 5 5 35 4 14250 9 6 6 24 6 14260 9 7 7 63 5 14270 9 8 7 63 7 14280 10 1 1 24 1 14290 10 2 2 24 1 14300 10 3 3 24 1 14310 10 4 4 24 1 14320 10 5 5 24 1 14330 10 6 6 24 6 14340 10 7 7 48 1 14350 10 8 8 48 8 14360 10 9 9 48 9 14370 0 14380 0 FINAL EIGENV 14390 14400 -2.536352449033E-01 2.038265111579E-01 -1.695641085 14410 -1.194676112328E-01 7.012023880100E-02 6.100076826 14420 1 N NA NB KM EPS 14430 14440 50 2 1 -200 1.00E-10 14450 0 14460 0 14470 0 P EM(IN) EM(OUT) KS K 14480 14490 2 1 1 44 1 14500 3 1 1 3 1 14510 3 2 2 43 2 14520 4 1 1 3 1 14530 4 2 2 3 1 14540 4 3 3 56 3 14550 5 1 1 3 1 14560 5 2 2 3 1 14570 5 3 3 3 1 14580 5 4 4 44 4 14590 6 1 1 3 1 14600 6 2 2 3 1 14610 6 3 3 3 1 14620 6 4 4 3 1 14630 6 5 5 107 5 14640 7 1 1 3 1 14650 7 2 2 3 1 14660 7 3 3 3 1 14670 7 4 4 3 1 14680 7 5 5 8 1 14690 7 6 6 6 6 14700 8 1 1 3 1 14710 8 2 2 3 1 14720 8 3 3 3 1 14730 8 4 4 3 1 14740 8 5 5 4 1 14750 8 6 6 3 1 14760 8 7 7 80 7 14770 9 1 1 3 1 14780 9 2 2 3 1 14790 9 3 3 3 1 14800 9 4 4 3 1 14810 9 5 5 3 1 14820 9 6 6 3 1 14830 9 7 7 3 1 14840 9 8 8 99 1 14850 10 1 1 3 1 14860 10 2 2 3 1 14870 10 3 3 3 1 14880 10 4 4 3 1 14890 10 5 5 3 1 14900 10 6 6 3 1 14910 10 7 7 3 1 14920 10 8 8 18 1 14930 10 9 9 3 1 14940 0 14950 0 FINAL EIGENV 14960 14970 -2.536352449033E-01 2.038265111579E-01 -1.695641085 14980 -1.194676112328E-01 7.012023880101E-02 6.100076826 14990