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
c Use a wrapper for C99 hypot, which is guaranteed to handle special values
c such as NaN.
DOUBLE PRECISION FUNCTION PYTHAG(A,B)
DOUBLE PRECISION A,B,P
CALL HYPOT(A,B,P)
PYTHAG = P
RETURN
END
C$$$ DOUBLE PRECISION FUNCTION PYTHAG(A,B)
C$$$ DOUBLE PRECISION A,B
C$$$C
C$$$C FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
C$$$C
C$$$ DOUBLE PRECISION P,R,S,T,U
C$$$ P = DMAX1(DABS(A),DABS(B))
C$$$c 'nan' did lead to infinite loop before 2018-02:
C$$$ IF ((risfinite(P) .ne. 0) .and. (P .NE. 0.0D0)) THEN
C$$$ R = (DMIN1(DABS(A),DABS(B))/P)**2
C$$$ 10 CONTINUE
C$$$ T = 4.0D0 + R
C$$$ IF (T .EQ. 4.0D0) GO TO 20
C$$$ S = R/T
C$$$ U = 1.0D0 + 2.0D0*S
C$$$ P = U*P
C$$$ R = (S/U)**2 * R
C$$$ GO TO 10
C$$$ END IF
C$$$ 20 PYTHAG = P
C$$$ RETURN
C$$$ END
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)
CALL TQLRAT(N,W,FV2,IERR)
GO TO 50
C .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
20 CALL TRED2(NM,N,A,W,FV1,Z)
CALL TQL2(NM,N,W,FV1,Z,IERR)
50 RETURN
END
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 LABORATORY
C
C THIS VERSION DATED AUGUST 1983.
C
C ------------------------------------------------------------------
c
c unnecessary initialization of C3 and S2 to keep g77 -Wall happy
c
C3 = 0.0D0
S2 = 0.0D0
C
IERR = 0
IF (N .EQ. 1) GO TO 1001
C
DO I = 2, N
E(I-1) = E(I)
end do
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 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 ..........
end do
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))
R = PYTHAG(P,1.0D0)
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
DO I = L2, N
D(I) = D(I) - H
end do
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 -- ..........
DO 200 II = 1, MML
C3 = C2
C2 = C
S2 = S
I = M - II
G = C * E(I)
H = C * P
R = PYTHAG(P,E(I))
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 .......... FORM VECTOR ..........
DO 180 K = 1, N
H = Z(K,I+1)
Z(K,I+1) = S * Z(K,I) + C * H
Z(K,I) = C * Z(K,I) - S * H
180 CONTINUE
C
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
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
END
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 1983.
C
C ------------------------------------------------------------------
c
c unnecessary initialization of B and C to keep g77 -Wall happy
c
B = 0.0D0
C = 0.0D0
C
IERR = 0
IF (N .EQ. 1) GO TO 1001
C
DO I = 2, N
E2(I-1) = E2(I)
end do
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
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)
R = PYTHAG(P,1.0D0)
D(L) = S / (P + DSIGN(R,P))
H = G - D(L)
C
DO I = L1, N
D(I) = D(I) - H
end do
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 -- ..........
DO 200 II = 1, MML
I = M - II
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
IF (G .EQ. 0.0D0) G = B
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 -- ..........
DO 230 II = 2, L
I = L + 2 - II
IF (P .GE. D(I-1)) GO TO 270
D(I) = D(I-1)
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
SUBROUTINE TRED1(NM,N,A,D,E,E2)
C
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 LABORATORY
C
C THIS VERSION DATED AUGUST 1983.
C
C ------------------------------------------------------------------
C
DO 100 I = 1, N
D(I) = A(N,I)
A(N,I) = A(I,I)
100 CONTINUE
C .......... FOR I=N STEP -1 UNTIL 1 DO -- ..........
DO 300 II = 1, N
I = N + 1 - II
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 K = 1, L
SCALE = SCALE + DABS(D(K))
end do
C
IF (SCALE .NE. 0.0D0) GO TO 140
C
DO J = 1, L
D(J) = A(L,J)
A(L,J) = A(I,J)
A(I,J) = 0.0D0
end do
C
130 E(I) = 0.0D0
E2(I) = 0.0D0
GO TO 300
C
140 DO K = 1, L
D(K) = D(K) / SCALE
H = H + D(K) * D(K)
end do
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 ..........
DO J = 1, L
E(J) = 0.0D0
end do
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 K = JP1, L
G = G + A(K,J) * D(K)
E(K) = E(K) + A(K,J) * F
end do
C
220 E(J) = G
240 CONTINUE
C .......... FORM P ..........
F = 0.0D0
C
DO J = 1, L
E(J) = E(J) / H
F = F + E(J) * D(J)
end do
C
H = F / (H + H)
C .......... FORM Q ..........
DO J = 1, L
E(J) = E(J) - H * D(J)
end do
C .......... FORM REDUCED A ..........
DO J = 1, L
F = D(J)
G = E(J)
C
DO K = J, L
A(K,J) = A(K,J) - F * E(K) - G * D(K)
end do
C
end do
C
285 DO J = 1, L
F = D(J)
D(J) = A(L,J)
A(L,J) = A(I,J)
A(I,J) = F * SCALE
end do
C
300 CONTINUE
C
RETURN
END
SUBROUTINE TRED2(NM,N,A,D,E,Z)
C
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 LABORATORY
C
C THIS VERSION DATED AUGUST 1983.
C
C ------------------------------------------------------------------
C
DO I = 1, N
DO J = I, N
Z(J,I) = A(J,I)
end do
D(I) = A(N,I)
end do
C
IF (N .EQ. 1) GO TO 510
C .......... FOR I=N STEP -1 UNTIL 2 DO -- ..........
DO 300 II = 2, N
I = N + 2 - II
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 K = 1, L
SCALE = SCALE + DABS(D(K))
end do
C
IF (SCALE .NE. 0.0D0) GO TO 140
130 E(I) = D(L)
C
DO J = 1, L
D(J) = Z(L,J)
Z(I,J) = 0.0D0
Z(J,I) = 0.0D0
end do
C
GO TO 290
C
140 DO K = 1, L
D(K) = D(K) / SCALE
H = H + D(K) * D(K)
end do
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 ..........
DO J = 1, L
E(J) = 0.0D0
end do
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 K = JP1, L
G = G + Z(K,J) * D(K)
E(K) = E(K) + Z(K,J) * F
end do
C
220 E(J) = G
240 CONTINUE
C .......... FORM P ..........
F = 0.0D0
C
DO J = 1, L
E(J) = E(J) / H
F = F + E(J) * D(J)
end do
C
HH = F / (H + H)
C .......... FORM Q ..........
DO J = 1, L
E(J) = E(J) - HH * D(J)
end do
C .......... FORM REDUCED A ..........
DO 280 J = 1, L
F = D(J)
G = E(J)
C
DO K = J, L
Z(K,J) = Z(K,J) - F * E(K) - G * D(K)
end do
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 .ne. 0.0D0) then
DO K = 1, L
D(K) = Z(K,I) / H
end do
C
DO J = 1, L
G = 0.0D0
DO K = 1, L
G = G + Z(K,I) * Z(K,J)
end do
C
DO K = 1, L
Z(K,J) = Z(K,J) - G * D(K)
end do
end do
end if
C 380
DO K = 1, L
Z(K,I) = 0.0D0
end do
C
500 CONTINUE
C
510 DO I = 1, N
D(I) = Z(N,I)
Z(N,I) = 0.0D0
end do
C
Z(N,N) = 1.0D0
E(1) = 0.0D0
RETURN
END