C this version was is from the Pv processor ... 15 Mar 94 (sg)_ SUBROUTINE F02ABF(A, IA, N, R, V, IV, E, IFAIL) C C SIMULATES NAG DIAGONALISERS WITH EISPACK CALLS C JMH JUNE 92 C IMPLICIT DOUBLE PRECISION (A-H,O-Z) DIMENSION A(IA,N), V(IV,N), R(N), E(N) DIMENSION F(600) DATA MXF/600/ C IF(IV.NE.IA) THEN WRITE(6,100) IV,IA 100 FORMAT(' *** ERROR IN SIMULATED F02ABF: IV.NE.IA',2I5) STOP ENDIF C IF(N.GT.MXF) THEN WRITE(6,110) N,MXF 110 FORMAT(' *** ERROR IN SIMULATED F02ABF: N =',I4,' .GT.',I4) STOP ENDIF C CALL RS(IA, N, A, R, 1, V, E, F, IERR) C IF (IERR .NE. 0) THEN WRITE (6,120) IERR 120 FORMAT(' *** SIMULATED F02ABF: IERR =', I3, ' .NE. 0') STOP END IF C IFAIL=0 RETURN C ENTRY F02AAF (A, IA, N, R, E, IFAIL) C IF(N.GT.MXF) THEN WRITE(6,210) N,MXF 210 FORMAT(' *** ERROR IN SIMULATED F02AAF: N =',I4,' .GT.',I4) STOP ENDIF C CALL TRED1 (IA, N, A, R, F, E) C C GET EIGENVALUES OF TRIDIAGONAL MATRIX C CALL TQLRAT (N, R, E, IERR) C IF (IERR .NE. 0) THEN WRITE (6,220) IERR 220 FORMAT(' *** SIMULATED F02AAF: IERR =', I3, ' .NE. 0') STOP END IF C IFAIL=0 RETURN C END CSTART VAX UNIVAC UNIX-HP DEC-RISC FPS *** FROM NETLIB, TUE MAY 24 12:02:59 CDT 1988 *** * EISPACK EIGENVALUE PACKAGE SUBROUTINE RS(NM,N,A,W,MATZ,Z,FV1,FV2,IERR) C INTEGER N,NM,IERR,MATZ DOUBLE PRECISION A(NM,N),W(N),Z(NM,N),FV1(N),FV2(N) C C THIS SUBROUTINE CALLS THE RECOMMENDED SEQUENCE OF C SUBROUTINES FROM THE EIGENSYSTEM SUBROUTINE PACKAGE (EISPACK) C TO FIND THE EIGENVALUES AND EIGENVECTORS (IF DESIRED) C OF A REAL SYMMETRIC MATRIX. C C ON INPUT C C NM MUST BE SET TO THE ROW DIMENSION OF THE TWO-DIMENSIONAL C ARRAY PARAMETERS AS DECLARED IN THE CALLING PROGRAM C DIMENSION STATEMENT. C C N IS THE ORDER OF THE MATRIX A. C C A CONTAINS THE REAL SYMMETRIC MATRIX. C C MATZ IS AN INTEGER VARIABLE SET EQUAL TO ZERO IF C ONLY EIGENVALUES ARE DESIRED. OTHERWISE IT IS SET TO C ANY NON-ZERO INTEGER FOR BOTH EIGENVALUES AND EIGENVECTORS. C C ON OUTPUT C C W CONTAINS THE EIGENVALUES IN ASCENDING ORDER. C C Z CONTAINS THE EIGENVECTORS IF MATZ IS NOT ZERO. C C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR TQLRAT C AND TQL2. THE NORMAL COMPLETION CODE IS ZERO. C C FV1 AND FV2 ARE TEMPORARY STORAGE ARRAYS. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1983. C C ------------------------------------------------------------------ C IF (N .LE. NM) GO TO 10 IERR = 10 * N GO TO 50 C 10 IF (MATZ .NE. 0) GO TO 20 C .......... FIND EIGENVALUES ONLY .......... CALL TRED1(NM,N,A,W,FV1,FV2) * TQLRAT ENCOUNTERS CATASTROPHIC UNDERFLOW ON THE VAX CALL TQLRAT(N,W,FV2,IERR) * CALL TQL1(N,W,FV1,IERR) GO TO 50 C .......... FIND BOTH EIGENVALUES AND EIGENVECTORS .......... 20 CALL TRED2(NM,N,A,W,FV1,Z) CEND CSTART VAX UNIVAC UNIX-HP DEC-RISC GENERIC CALL TQL2(NM,N,W,FV1,Z,IERR) CEND CSTART FPS CTRP CALL EIS_TQL2(NM,N,W,FV1,Z,IERR) CEND CSTART VAX UNIVAX UNIX-HP DEC-RISC FPS GENERIC 50 RETURN END C**************************************************************** C Translated by IBM AIX XL FORTRAN -Pv Preprocessor/6000 C Version 02.03.0000.0023 4.01M52 on 3/14/94 at 15:06:43 C**************************************************************** C SUBROUTINE TQLRAT(N,D,E2,IERR) C INTEGER I,J,L,M,N,II,L1,MML,IERR DOUBLE PRECISION D(N),E2(N) DOUBLE PRECISION B,C,F,G,H,P,R,S,T,EPSLON,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQLRAT, C ALGORITHM 464, COMM. ACM 16, 689(1973) BY REINSCH. C C THIS SUBROUTINE FINDS THE EIGENVALUES OF A SYMMETRIC C TRIDIAGONAL MATRIX BY THE RATIONAL QL METHOD. C C ON INPUT C C N IS THE ORDER OF THE MATRIX. C C D CONTAINS THE DIAGONAL ELEMENTS OF THE INPUT MATRIX. C C E2 CONTAINS THE SQUARES OF THE SUBDIAGONAL ELEMENTS OF THE C INPUT MATRIX IN ITS LAST N-1 POSITIONS. E2(1) IS ARBITRARY. C C ON OUTPUT C C D CONTAINS THE EIGENVALUES IN ASCENDING ORDER. IF AN C ERROR EXIT IS MADE, THE EIGENVALUES ARE CORRECT AND C ORDERED FOR INDICES 1,2,...IERR-1, BUT MAY NOT BE C THE SMALLEST EIGENVALUES. C C E2 HAS BEEN DESTROYED. 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 CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATORY C C THIS VERSION DATED AUGUST 1987. C MODIFIED BY C. MOLER TO FIX UNDERFLOW/OVERFLOW DIFFICULTIES, C ESPECIALLY ON THE VAX AND OTHER MACHINES WHERE EPSLON(1.0D0)**2 C NEARLY UNDERFLOWS. SEE THE LOOP INVOLVING STATEMENT 102 AND C THE TWO STATEMENTS JUST BEFORE STATEMENT 200. C C ------------------------------------------------------------------ C INTEGER I1, I2, I3, I4 DOUBLEPRECISION PYTHAG1, P1, R1, S1, T1, U LOGICAL L2, L3, L4, L5, L6 IERR = 0 IF (N .EQ. 1) GO TO 1001 C I4 = IAND(MAX0(N - 1,0),1) IF (I4 .EQ. 1) E2(2-1) = E2(2) DO 100 I = I4 + 2, N, 2 E2(I-1) = E2(I) E2(I) = E2(I+1) 100 CONTINUE C F = 0.0D0 T = 0.0D0 E2(N) = 0.0D0 C DO 290 L = 1, N J = 0 H = DABS(D(L)) + DSQRT(E2(L)) IF (T .GT. H) GO TO 105 T = H B = EPSLON(T) C = B * B IF (C .NE. 0.0D0) GO TO 105 C SPLITING TOLERANCE UNDERFLOWED. LOOK FOR LARGER VALUE. DO 102 I = 1, N - L + 1 T = MAX(DABS(D(L+I-1))+DSQRT(E2(L+I-1)),T) 102 CONTINUE B = EPSLON(T) C = B * B C .......... LOOK FOR SMALL SQUARED SUB-DIAGONAL ELEMENT .......... 105 DO 110 M = L, N IF (E2(M) .LE. C) GO TO 120 C .......... E2(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP .......... 110 CONTINUE C 120 IF (M .EQ. L) GO TO 210 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... L1 = L + 1 S = DSQRT(E2(L)) G = D(L) P = (D(L1) - G) / (2.0D0 * S) C***** Code Expanded From Routine: PYTHAG P1 = DMAX1(DABS(P),DABS(1.0D0)) IF (P1 .EQ. 0.0D0) GO TO 77001 R1 = (DMIN1(DABS(P),DABS(1.0D0))/P1)**2 77002 CONTINUE T1 = 4.0D0 + R1 IF (T1 .EQ. 4.0D0) GO TO 77001 S1 = R1/T1 U = 1.0D0 + 2.0D0*S1 P1 = U*P1 R1 = (S1/U)**2*R1 GO TO 77002 77001 CONTINUE PYTHAG1 = P1 R = PYTHAG1 C***** End of Code Expanded From Routine: PYTHAG D(L) = S / (P + DSIGN(R,P)) H = G - D(L) C I3 = IAND(MAX0(N - L1 + 1,0),1) IF (I3 .EQ. 1) D(L1) = D(L1) - H DO 140 I = I3 + L1, N, 2 D(I) = D(I) - H D(I+1) = D(I+1) - H 140 CONTINUE C F = F + H C .......... RATIONAL QL TRANSFORMATION .......... G = D(M) IF (G .EQ. 0.0D0) G = B H = G S = 0.0D0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... C DO 200 II = 1, MML C I = M - II DO 200 I=M-1,L,-1 P = G * H R = P + E2(I) E2(I+1) = S * R S = E2(I) / R D(I+1) = H + S * (H + D(I)) G = D(I) - E2(I) / G C AVOID DIVISION BY ZERO ON NEXT PASS IF (G .EQ. 0.0D0) G = EPSLON(D(I)) H = G * (P / R) 200 CONTINUE C E2(L) = S * G D(L) = H C .......... GUARD AGAINST UNDERFLOW IN CONVERGENCE TEST .......... IF (H .EQ. 0.0D0) GO TO 210 IF (DABS(E2(L)) .LE. DABS(C/H)) GO TO 210 E2(L) = H * E2(L) IF (E2(L) .NE. 0.0D0) GO TO 130 210 P = D(L) + F C .......... ORDER EIGENVALUES .......... IF (L .EQ. 1) GO TO 250 C .......... FOR I=L STEP -1 UNTIL 2 DO -- .......... C DO 230 II = 2, L C I = L + 2 - II I1 = IAND(MAX0(L - 1,0),3) DO I = L, L + 1 - I1, -1 IF (P .GE. D(I-1)) GO TO 270 D(I) = D(I-1) END DO DO 230 I2 = L - I1, 2, -4 L6 = P .GE. D(I2-1) I = I2 IF (P .GE. D(I2-1)) GO TO 270 D(I2) = D(I2-1) L5 = P .GE. D(I2-2) I = I2 - 1 IF (P .GE. D(I2-2)) GO TO 270 D(I2-1) = D(I2-2) L4 = P .GE. D(I2-3) I = I2 - 2 IF (P .GE. D(I2-3)) GO TO 270 D(I2-2) = D(I2-3) L3 = P .GE. D(I2-4) I = I2 - 3 IF (P .GE. D(I2-4)) GO TO 270 D(I2-3) = D(I2-4) 230 CONTINUE C 250 I = 1 270 D(I) = P 290 CONTINUE C GO TO 1001 C .......... SET ERROR -- NO CONVERGENCE TO AN C EIGENVALUE AFTER 30 ITERATIONS .......... 1000 IERR = L 1001 RETURN END C**************************************************************** C Translated by IBM AIX XL FORTRAN -Pv Preprocessor/6000 C Version 02.03.0000.0023 4.01M52 on 3/14/94 at 15:06:43 C**************************************************************** C CSTART VAX UNIVAX UNIX-HP DEC-RISC GENERIC SUBROUTINE TQL2(NM,N,D,E,Z,IERR) C INTEGER I,J,K,L,M,N,II,L1,L2,NM,MML,IERR DOUBLE PRECISION D(N),E(N),Z(NM,N) DOUBLE PRECISION C,C2,C3,DL1,EL1,F,G,H,P,R,S,S2,TST1,TST2,PYTHAG C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TQL2, C NUM. MATH. 11, 293-306(1968) BY BOWDLER, MARTIN, REINSCH, AND C WILKINSON. C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 227-240(1971). C C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS C OF A SYMMETRIC TRIDIAGONAL MATRIX BY THE 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 CALLS PYTHAG FOR DSQRT(A*A + B*B) . C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATOR C C THIS VERSION DATED AUGUST 1983. C C ---------------------------------------------------------------- C INTEGER J1, I1, I2 DOUBLEPRECISION PYTHAG1, P1, R1, S1, T, U, PYTHAG2, P2, R2, S3, T1 ., U1 LOGICAL L3, L4, L5, L6, L7, L8 IERR = 0 IF (N .EQ. 1) GO TO 1001 C I2 = IAND(MAX0(N - 1,0),1) IF (I2 .EQ. 1) E(2-1) = E(2) DO 100 I = I2 + 2, N, 2 E(I-1) = E(I) E(I) = E(I+1) 100 CONTINUE C F = 0.0D0 TST1 = 0.0D0 E(N) = 0.0D0 C DO 240 L = 1, N J = 0 H = DABS(D(L)) + DABS(E(L)) IF (TST1 .LT. H) TST1 = H C .......... LOOK FOR SMALL SUB-DIAGONAL ELEMENT .......... DO 110 M = L, N TST2 = TST1 + DABS(E(M)) IF (TST2 .EQ. TST1) GO TO 120 C .......... E(N) IS ALWAYS ZERO, SO THERE IS NO EXIT C THROUGH THE BOTTOM OF THE LOOP .......... 110 CONTINUE C 120 IF (M .EQ. L) GO TO 220 130 IF (J .EQ. 30) GO TO 1000 J = J + 1 C .......... FORM SHIFT .......... L1 = L + 1 L2 = L1 + 1 G = D(L) P = (D(L1) - G) / (2.0D0 * E(L)) C***** Code Expanded From Routine: PYTHAG P1 = DMAX1(DABS(P),DABS(1.0D0)) IF (P1 .EQ. 0.0D0) GO TO 77053 R1 = (DMIN1(DABS(P),DABS(1.0D0))/P1)**2 77054 CONTINUE T = 4.0D0 + R1 IF (T .EQ. 4.0D0) GO TO 77053 S1 = R1/T U = 1.0D0 + 2.0D0*S1 P1 = U*P1 R1 = (S1/U)**2*R1 GO TO 77054 77053 CONTINUE PYTHAG1 = P1 R = PYTHAG1 C***** End of Code Expanded From Routine: PYTHAG D(L) = E(L) / (P + DSIGN(R,P)) D(L1) = E(L) * (P + DSIGN(R,P)) DL1 = D(L1) H = G - D(L) IF (L2 .GT. N) GO TO 145 C I1 = IAND(MAX0(N - L2 + 1,0),1) IF (I1 .EQ. 1) D(L2) = D(L2) - H DO 140 I = I1 + L2, N, 2 D(I) = D(I) - H D(I+1) = D(I+1) - H 140 CONTINUE C 145 F = F + H C .......... QL TRANSFORMATION .......... P = D(M) C = 1.0D0 C2 = C EL1 = E(L1) S = 0.0D0 MML = M - L C .......... FOR I=M-1 STEP -1 UNTIL L DO -- .......... C DO 200 II = 1, MML DO 200 I=M-1,L,-1 C3 = C2 C2 = C S2 = S C I = M - II G = C * E(I) H = C * P C***** Code Expanded From Routine: PYTHAG P2 = DMAX1(DABS(P),DABS(E(I))) IF (P2 .EQ. 0.0D0) GO TO 77055 R2 = (DMIN1(DABS(P),DABS(E(I)))/P2)**2 77056 CONTINUE T1 = 4.0D0 + R2 IF (T1 .EQ. 4.0D0) GO TO 77055 S3 = R2/T1 U1 = 1.0D0 + 2.0D0*S3 P2 = U1*P2 R2 = (S3/U1)**2*R2 GO TO 77056 77055 CONTINUE PYTHAG2 = P2 R = PYTHAG2 C***** End of Code Expanded From Routine: PYTHAG E(I+1) = S * R S = E(I) / R C = P / R P = C * D(I) - S * G D(I+1) = H + S * (C * G + S * D(I)) C REPLACE THIS INNER LOOP WITH BLAS CALL (MILLARD ALEXANDER 12/24/90) C .......... FORM VECTOR .......... C DO 180 K = 1, N C H = Z(K,I+1) C HH = Z(K, I) C Z(K,I+1) = S * HH + C * H C Z(K,I) = C * HH - S * H C 180 CONTINUE CEND CSTART 1/19/93 SG MODIFIED ACCORDING TO JMH SUGGESTION CALL DROT (N, Z(1,I+1), 1, Z(1,I), 1, C, S) CEND CSTART VAX UNIVAC UNIX-HP DEC-RISC GENERIC 200 CONTINUE C P = -S * S2 * C3 * EL1 * E(L) / DL1 E(L) = S * P D(L) = C * P TST2 = TST1 + DABS(E(L)) IF (TST2 .GT. TST1) GO TO 130 220 D(L) = D(L) + F 240 CONTINUE C .......... ORDER EIGENVALUES AND EIGENVECTORS .......... DO 300 II = 2, N I = II - 1 K = I P = D(I) C J1 = IAND(MAX0(N - II + 1,0),3) DO J = II, J1 - 1 + II IF (D(J) .GE. P) GO TO 77058 K = J P = D(J) 77058 CONTINUE END DO DO 260 J = J1 + II, N, 4 IF (D(J) .GE. P) GO TO 77059 K = J P = D(J) L8 = D(J+1) .GE. P 77059 CONTINUE IF (D(J+1) .GE. P) GO TO 77060 K = J + 1 P = D(J+1) L7 = D(J+2) .GE. P 77060 CONTINUE IF (D(J+2) .GE. P) GO TO 77061 K = J + 2 P = D(J+2) L6 = D(J+3) .GE. P 77061 CONTINUE IF (D(J+3) .GE. P) GO TO 77062 K = J + 3 P = D(J+3) 77062 CONTINUE 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 END C**************************************************************** C Translated by IBM AIX XL FORTRAN -Pv Preprocessor/6000 C Version 02.03.0000.0023 4.01M52 on 3/14/94 at 15:06:43 C**************************************************************** C SUBROUTINE TRED1(NM,N,A,D,E,E2) C DOUBLEPRECISION D1S,D2S INTEGER I,J,K,L,N,II,NM,JP1 DOUBLE PRECISION A(NM,N),D(N),E(N),E2(N) DOUBLE PRECISION F,G,H,SCALE C C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE TRED1, 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 C TO A SYMMETRIC TRIDIAGONAL MATRIX USING 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 A CONTAINS INFORMATION ABOUT THE ORTHOGONAL TRANS- C FORMATIONS USED IN THE REDUCTION IN ITS STRICT LOWER C TRIANGLE. THE FULL UPPER TRIANGLE OF A IS UNALTERED. 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 E2 CONTAINS THE SQUARES OF THE CORRESPONDING ELEMENTS OF E. C E2 MAY COINCIDE WITH E IF THE SQUARES ARE NOT NEEDED. C C QUESTIONS AND COMMENTS SHOULD BE DIRECTED TO BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATOR C C THIS VERSION DATED AUGUST 1983. C C ---------------------------------------------------------------- C INTEGER I1, J1, J2 I1 = 0 DO 100 I = 1, N D(I) = A(N,I) A(N,I) = A(1+I1,1) I1 = I1 + NM + 1 100 CONTINUE C .......... FOR I=N STEP -1 UNTIL 1 DO -- .......... C DO 300 II = 1, N C I = N + 1 - II DO 300 I=N,1,-1 L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 1) GO TO 130 C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... DO 120 K = 1, L 120 SCALE = SCALE + DABS(D(K)) C IF (SCALE .NE. 0.0D0) GO TO 140 C DO 125 J = 1, L D(J) = A(L,J) A(L,J) = A(I,J) A(I,J) = 0.0D0 125 CONTINUE C 130 E(I) = 0.0D0 E2(I) = 0.0D0 GO TO 300 C C 140 CONTINUE D1S = 1.D0/SCALE DO 150 K = 1, L D(K) = D(K)*D1S H = H + D(K)*D(K) 150 CONTINUE C E2(I) = SCALE * SCALE * H F = D(L) G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G D(L) = F - G IF (L .EQ. 1) GO TO 285 C .......... FORM A*U .......... J2 = IAND(MAX0(L,0),1) IF (J2 .EQ. 1) E(1) = 0.0D0 DO 170 J = J2 + 1, L, 2 E(J) = 0.0D0 E(J+1) = 0.0D0 170 CONTINUE C DO 240 J = 1, L F = D(J) G = E(J) + A(J,J) * F JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L G = G + A(K,J) * D(K) E(K) = E(K) + A(K,J) * F 200 CONTINUE C 220 E(J) = G 240 CONTINUE C .......... FORM P .......... F = 0.0D0 C D2S = 1.D0/H DO 245 J = 1, L E(J) = E(J)*D2S F = F + E(J)*D(J) 245 CONTINUE C H = F / (H + H) C .......... FORM Q .......... J1 = IAND(MAX0(L,0),1) IF (J1 .EQ. 1) E(1) = E(1) - H*D(1) DO 250 J = J1 + 1, L, 2 E(J) = E(J) - H*D(J) E(J+1) = E(J+1) - H*D(J+1) 250 CONTINUE C .......... FORM REDUCED A .......... DO 280 J = 1, L F = D(J) G = E(J) C DO 260 K = J, L 260 A(K,J) = A(K,J) - F * E(K) - G * D(K) C 280 CONTINUE C 285 DO 290 J = 1, L F = D(J) D(J) = A(L,J) A(L,J) = A(I,J) A(I,J) = F * SCALE 290 CONTINUE C 300 CONTINUE C RETURN END C**************************************************************** C Translated by IBM AIX XL FORTRAN -Pv Preprocessor/6000 C Version 02.03.0000.0023 4.01M52 on 3/14/94 at 15:06:43 C**************************************************************** C SUBROUTINE TRED2(NM,N,A,D,E,Z) C DOUBLEPRECISION D1S,D2S,D3S INTEGER I,J,K,L,N,II,NM,JP1 DOUBLE PRECISION A(NM,N),D(N),E(N),Z(NM,N) DOUBLE PRECISION F,G,H,HH,SCALE 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 BURTON S. GARBOW, C MATHEMATICS AND COMPUTER SCIENCE DIV, ARGONNE NATIONAL LABORATOR C C THIS VERSION DATED AUGUST 1983. C C ---------------------------------------------------------------- C INTEGER K1, K2, K3, K4, J1, J2, J3 DOUBLEPRECISION D1, G1, G2, G3 LOGICAL L1 DO 100 I = 1, N C J3 = IAND(MAX0(N - I + 1,0),1) IF (J3 .EQ. 1) Z(I,I) = A(I,I) DO 80 J = J3 + I, N, 2 Z(J,I) = A(J,I) Z(J+1,I) = A(J+1,I) 80 CONTINUE C D(I) = A(N,I) 100 CONTINUE C IF (N .EQ. 1) GO TO 510 C .......... FOR I=N STEP -1 UNTIL 2 DO -- .......... C DO 300 II = 2, N C I = N + 2 - II DO 300 I=N,2,-1 L = I - 1 H = 0.0D0 SCALE = 0.0D0 IF (L .LT. 2) GO TO 130 C .......... SCALE ROW (ALGOL TOL THEN NOT NEEDED) .......... DO 120 K = 1, L 120 SCALE = SCALE + DABS(D(K)) C IF (SCALE .NE. 0.0D0) GO TO 140 130 E(I) = D(L) C DO 135 J = 1, L D(J) = Z(L,J) Z(I,J) = 0.0D0 Z(J,I) = 0.0D0 135 CONTINUE C GO TO 290 C C 140 CONTINUE D1S = 1.D0/SCALE DO 150 K = 1, L D(K) = D(K)*D1S H = H + D(K)*D(K) 150 CONTINUE C F = D(L) G = -DSIGN(DSQRT(H),F) E(I) = SCALE * G H = H - F * G D(L) = F - G C .......... FORM A*U .......... J2 = IAND(MAX0(L,0),1) IF (J2 .EQ. 1) E(1) = 0.0D0 DO 170 J = J2 + 1, L, 2 E(J) = 0.0D0 E(J+1) = 0.0D0 170 CONTINUE C DO 240 J = 1, L F = D(J) Z(J,I) = F G = E(J) + Z(J,J) * F JP1 = J + 1 IF (L .LT. JP1) GO TO 220 C DO 200 K = JP1, L G = G + Z(K,J) * D(K) E(K) = E(K) + Z(K,J) * F 200 CONTINUE C 220 E(J) = G 240 CONTINUE C .......... FORM P .......... F = 0.0D0 C D2S = 1.D0/H DO 245 J = 1, L E(J) = E(J)*D2S F = F + E(J)*D(J) 245 CONTINUE C HH = F / (H + H) C .......... FORM Q .......... J1 = IAND(MAX0(L,0),1) IF (J1 .EQ. 1) E(1) = E(1) - HH*D(1) DO 250 J = J1 + 1, L, 2 E(J) = E(J) - HH*D(J) E(J+1) = E(J+1) - HH*D(J+1) 250 CONTINUE C .......... FORM REDUCED A .......... DO 280 J = 1, L F = D(J) G = E(J) C DO 260 K = J, L 260 Z(K,J) = Z(K,J) - F * E(K) - G * D(K) C D(J) = Z(L,J) Z(I,J) = 0.0D0 280 CONTINUE C 290 D(I) = H 300 CONTINUE C .......... ACCUMULATION OF TRANSFORMATION MATRICES .......... DO 500 I = 2, N L = I - 1 Z(N,L) = Z(L,L) Z(L,L) = 1.0D0 H = D(I) IF (H .EQ. 0.0D0) GO TO 380 C D3S = 1.D0/H K4 = IAND(MAX0(L,0),1) IF (K4 .EQ. 1) D(1) = Z(1,I)*D3S DO K = K4 + 1, L, 2 D(K) = Z(K,I)*D3S D(K+1) = Z(K+1,I)*D3S END DO C DO 360 J = 1, L G = 0.0D0 K3 = IAND(MAX0(L,0),3) DO K = 1, K3 G = G + Z(K,I)*Z(K,J) END DO G1 = 0.D0 G2 = 0.D0 G3 = 0.D0 DO 340 K = K3 + 1, L, 4 G = G + Z(K,I)*Z(K,J) G1 = G1 + Z(K+1,I)*Z(K+1,J) G2 = G2 + Z(K+2,I)*Z(K+2,J) G3 = G3 + Z(K+3,I)*Z(K+3,J) 340 CONTINUE G = G + G1 + G2 + G3 K2 = IAND(MAX0(L,0),1) IF (K2 .EQ. 1) Z(1,J) = Z(1,J) - G*D(1) DO K = K2 + 1, L, 2 Z(K,J) = Z(K,J) - G*D(K) Z(K+1,J) = Z(K+1,J) - G*D(K+1) END DO 360 CONTINUE C 380 CONTINUE K1 = IAND(MAX0(L,0),1) IF (K1 .EQ. 1) Z(1,I) = 0.0D0 DO 400 K = K1 + 1, L, 2 Z(K,I) = 0.0D0 Z(K+1,I) = 0.0D0 400 CONTINUE C 500 CONTINUE C 510 DO 520 I = 1, N D(I) = Z(N,I) Z(N,I) = 0.0D0 520 CONTINUE C Z(N,N) = 1.0D0 E(1) = 0.0D0 RETURN END DOUBLE PRECISION FUNCTION EPSLON (X) DOUBLE PRECISION X C C ESTIMATE UNIT ROUNDOFF IN QUANTITIES OF SIZE X. C DOUBLE PRECISION A,B,C,EPS C C THIS PROGRAM SHOULD FUNCTION PROPERLY ON ALL SYSTEMS C SATISFYING THE FOLLOWING TWO ASSUMPTIONS, C 1. THE BASE USED IN REPRESENTING FLOATING POINT C NUMBERS IS NOT A POWER OF THREE. C 2. THE QUANTITY A IN STATEMENT 10 IS REPRESENTED TO C THE ACCURACY USED IN FLOATING POINT VARIABLES C THAT ARE STORED IN MEMORY. C THE STATEMENT NUMBER 10 AND THE GO TO 10 ARE INTENDED TO C FORCE OPTIMIZING COMPILERS TO GENERATE CODE SATISFYING C ASSUMPTION 2. C UNDER THESE ASSUMPTIONS, IT SHOULD BE TRUE THAT, C A IS NOT EXACTLY EQUAL TO FOUR-THIRDS, C B HAS A ZERO FOR ITS LAST BIT OR DIGIT, C C IS NOT EXACTLY EQUAL TO ONE, C EPS MEASURES THE SEPARATION OF 1.0 FROM C THE NEXT LARGER FLOATING POINT NUMBER. C THE DEVELOPERS OF EISPACK WOULD APPRECIATE BEING INFORMED C ABOUT ANY SYSTEMS WHERE THESE ASSUMPTIONS DO NOT HOLD. C C THIS VERSION DATED 4/6/83. C A = 4.0D0/3.0D0 10 B = A - 1.0D0 C = B + B + B EPS = DABS(C-1.0D0) IF (EPS .EQ. 0.0D0) GO TO 10 EPSLON = EPS*DABS(X) RETURN END