SUBROUTINE BALANC(NM,N,A,LOW,IGH,SCALE)
C
INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC
DOUBLE PRECISION A(NM,N),SCALE(N)
DOUBLE PRECISION C,F,G,R,S,B2,RADIX
LOGICAL NOCONV
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALANCE,
C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
C
C THIS SUBROUTINE BALANCES A REAL MATRIX AND ISOLATES
C EIGENVALUES WHENEVER POSSIBLE.
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 INPUT MATRIX TO BE BALANCED.
C
C ON OUTPUT
C
C A CONTAINS THE BALANCED MATRIX.
C
C LOW AND IGH ARE TWO INTEGERS SUCH THAT A(I,J)
C IS EQUAL TO ZERO IF
C (1) I IS GREATER THAN J AND
C (2) J=1,...,LOW-1 OR I=IGH+1,...,N.
C
C SCALE CONTAINS INFORMATION DETERMINING THE
C PERMUTATIONS AND SCALING FACTORS USED.
C
C SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH
C HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED
C WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS
C OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN
C SCALE(J) = P(J), FOR J = 1,...,LOW-1
C = D(J,J), J = LOW,...,IGH
C = P(J) J = IGH+1,...,N.
C THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1,
C THEN 1 TO LOW-1.
C
C NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY.
C
C THE ALGOL PROCEDURE EXC CONTAINED IN BALANCE APPEARS IN
C BALANC IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS
C K,L HAVE BEEN REVERSED.)
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
RADIX = 16.0D0
C
B2 = RADIX * RADIX
K = 1
L = N
GO TO 100
C .......... IN-LINE PROCEDURE FOR ROW AND
C COLUMN EXCHANGE ..........
20 SCALE(M) = J
IF (J .EQ. M) GO TO 50
C
DO 30 I = 1, L
F = A(I,J)
A(I,J) = A(I,M)
A(I,M) = F
30 CONTINUE
C
DO 40 I = K, N
F = A(J,I)
A(J,I) = A(M,I)
A(M,I) = F
40 CONTINUE
C
50 GO TO (80,130), IEXC
C .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE
C AND PUSH THEM DOWN ..........
80 IF (L .EQ. 1) GO TO 280
L = L - 1
C .......... FOR J=L STEP -1 UNTIL 1 DO -- ..........
100 DO 120 JJ = 1, L
J = L + 1 - JJ
C
DO 110 I = 1, L
IF (I .EQ. J) GO TO 110
IF (A(J,I) .NE. 0.0D0) GO TO 120
110 CONTINUE
C
M = L
IEXC = 1
GO TO 20
120 CONTINUE
C
GO TO 140
C .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE
C AND PUSH THEM LEFT ..........
130 K = K + 1
C
140 DO 170 J = K, L
C
DO 150 I = K, L
IF (I .EQ. J) GO TO 150
IF (A(I,J) .NE. 0.0D0) GO TO 170
150 CONTINUE
C
M = K
IEXC = 2
GO TO 20
170 CONTINUE
C .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L ..........
DO 180 I = K, L
180 SCALE(I) = 1.0D0
C .......... ITERATIVE LOOP FOR NORM REDUCTION ..........
190 NOCONV = .FALSE.
C
DO 270 I = K, L
C = 0.0D0
R = 0.0D0
C
DO 200 J = K, L
IF (J .EQ. I) GO TO 200
C = C + DABS(A(J,I))
R = R + DABS(A(I,J))
200 CONTINUE
C .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ..........
IF (C .EQ. 0.0D0 .OR. R .EQ. 0.0D0) GO TO 270
G = R / RADIX
F = 1.0D0
S = C + R
210 IF (C .GE. G) GO TO 220
F = F * RADIX
C = C * B2
GO TO 210
220 G = R * RADIX
230 IF (C .LT. G) GO TO 240
F = F / RADIX
C = C / B2
GO TO 230
C .......... NOW BALANCE ..........
240 IF ((C + R) / F .GE. 0.95D0 * S) GO TO 270
G = 1.0D0 / F
SCALE(I) = SCALE(I) * F
NOCONV = .TRUE.
C
DO 250 J = K, N
250 A(I,J) = A(I,J) * G
C
DO 260 J = 1, L
260 A(J,I) = A(J,I) * F
C
270 CONTINUE
C
IF (NOCONV) GO TO 190
C
280 LOW = K
IGH = L
RETURN
END
SUBROUTINE BALBAK(NM,N,LOW,IGH,SCALE,M,Z)
C
INTEGER I,J,K,M,N,II,NM,IGH,LOW
DOUBLE PRECISION SCALE(N),Z(NM,M)
DOUBLE PRECISION S
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE BALBAK,
C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
C
C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A REAL GENERAL
C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
C BALANCED MATRIX DETERMINED BY BALANC.
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 LOW AND IGH ARE INTEGERS DETERMINED BY BALANC.
C
C SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS
C AND SCALING FACTORS USED BY BALANC.
C
C M IS THE NUMBER OF COLUMNS OF Z TO BE BACK TRANSFORMED.
C
C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGEN-
C VECTORS TO BE BACK TRANSFORMED IN ITS FIRST M COLUMNS.
C
C ON OUTPUT
C
C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE
C TRANSFORMED EIGENVECTORS IN ITS FIRST M COLUMNS.
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 (M .EQ. 0) GO TO 200
IF (IGH .EQ. LOW) GO TO 120
C
DO 110 I = LOW, IGH
S = SCALE(I)
C .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED
C IF THE FOREGOING STATEMENT IS REPLACED BY
C S=1.0D0/SCALE(I). ..........
DO 100 J = 1, M
100 Z(I,J) = Z(I,J) * S
C
110 CONTINUE
C ......... FOR I=LOW-1 STEP -1 UNTIL 1,
C IGH+1 STEP 1 UNTIL N DO -- ..........
120 DO 140 II = 1, N
I = II
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 140
IF (I .LT. LOW) I = LOW - II
K = SCALE(I)
IF (K .EQ. I) GO TO 140
C
DO 130 J = 1, M
S = Z(I,J)
Z(I,J) = Z(K,J)
Z(K,J) = S
130 CONTINUE
C
140 CONTINUE
C
200 RETURN
END
SUBROUTINE CBABK2(NM,N,LOW,IGH,SCALE,M,ZR,ZI)
C
INTEGER I,J,K,M,N,II,NM,IGH,LOW
DOUBLE PRECISION SCALE(N),ZR(NM,M),ZI(NM,M)
DOUBLE PRECISION S
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE
C CBABK2, WHICH IS A COMPLEX VERSION OF BALBAK,
C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
C
C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX GENERAL
C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
C BALANCED MATRIX DETERMINED BY CBAL.
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 LOW AND IGH ARE INTEGERS DETERMINED BY CBAL.
C
C SCALE CONTAINS INFORMATION DETERMINING THE PERMUTATIONS
C AND SCALING FACTORS USED BY CBAL.
C
C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
C
C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVECTORS TO BE
C BACK TRANSFORMED IN THEIR FIRST M COLUMNS.
C
C ON OUTPUT
C
C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
C IN THEIR FIRST M COLUMNS.
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 (M .EQ. 0) GO TO 200
IF (IGH .EQ. LOW) GO TO 120
C
DO 110 I = LOW, IGH
S = SCALE(I)
C .......... LEFT HAND EIGENVECTORS ARE BACK TRANSFORMED
C IF THE FOREGOING STATEMENT IS REPLACED BY
C S=1.0D0/SCALE(I). ..........
DO 100 J = 1, M
ZR(I,J) = ZR(I,J) * S
ZI(I,J) = ZI(I,J) * S
100 CONTINUE
C
110 CONTINUE
C .......... FOR I=LOW-1 STEP -1 UNTIL 1,
C IGH+1 STEP 1 UNTIL N DO -- ..........
120 DO 140 II = 1, N
I = II
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 140
IF (I .LT. LOW) I = LOW - II
K = SCALE(I)
IF (K .EQ. I) GO TO 140
C
DO 130 J = 1, M
S = ZR(I,J)
ZR(I,J) = ZR(K,J)
ZR(K,J) = S
S = ZI(I,J)
ZI(I,J) = ZI(K,J)
ZI(K,J) = S
130 CONTINUE
C
140 CONTINUE
C
200 RETURN
END
SUBROUTINE CBAL(NM,N,AR,AI,LOW,IGH,SCALE)
C
INTEGER I,J,K,L,M,N,JJ,NM,IGH,LOW,IEXC
DOUBLE PRECISION AR(NM,N),AI(NM,N),SCALE(N)
DOUBLE PRECISION C,F,G,R,S,B2,RADIX
LOGICAL NOCONV
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE
C CBALANCE, WHICH IS A COMPLEX VERSION OF BALANCE,
C NUM. MATH. 13, 293-304(1969) BY PARLETT AND REINSCH.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 315-326(1971).
C
C THIS SUBROUTINE BALANCES A COMPLEX MATRIX AND ISOLATES
C EIGENVALUES WHENEVER POSSIBLE.
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 AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX MATRIX TO BE BALANCED.
C
C ON OUTPUT
C
C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE BALANCED MATRIX.
C
C LOW AND IGH ARE TWO INTEGERS SUCH THAT AR(I,J) AND AI(I,J)
C ARE EQUAL TO ZERO IF
C (1) I IS GREATER THAN J AND
C (2) J=1,...,LOW-1 OR I=IGH+1,...,N.
C
C SCALE CONTAINS INFORMATION DETERMINING THE
C PERMUTATIONS AND SCALING FACTORS USED.
C
C SUPPOSE THAT THE PRINCIPAL SUBMATRIX IN ROWS LOW THROUGH IGH
C HAS BEEN BALANCED, THAT P(J) DENOTES THE INDEX INTERCHANGED
C WITH J DURING THE PERMUTATION STEP, AND THAT THE ELEMENTS
C OF THE DIAGONAL MATRIX USED ARE DENOTED BY D(I,J). THEN
C SCALE(J) = P(J), FOR J = 1,...,LOW-1
C = D(J,J) J = LOW,...,IGH
C = P(J) J = IGH+1,...,N.
C THE ORDER IN WHICH THE INTERCHANGES ARE MADE IS N TO IGH+1,
C THEN 1 TO LOW-1.
C
C NOTE THAT 1 IS RETURNED FOR IGH IF IGH IS ZERO FORMALLY.
C
C THE ALGOL PROCEDURE EXC CONTAINED IN CBALANCE APPEARS IN
C CBAL IN LINE. (NOTE THAT THE ALGOL ROLES OF IDENTIFIERS
C K,L HAVE BEEN REVERSED.)
C
C ARITHMETIC IS REAL THROUGHOUT.
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
RADIX = 16.0D0
C
B2 = RADIX * RADIX
K = 1
L = N
GO TO 100
C .......... IN-LINE PROCEDURE FOR ROW AND
C COLUMN EXCHANGE ..........
20 SCALE(M) = J
IF (J .EQ. M) GO TO 50
C
DO 30 I = 1, L
F = AR(I,J)
AR(I,J) = AR(I,M)
AR(I,M) = F
F = AI(I,J)
AI(I,J) = AI(I,M)
AI(I,M) = F
30 CONTINUE
C
DO 40 I = K, N
F = AR(J,I)
AR(J,I) = AR(M,I)
AR(M,I) = F
F = AI(J,I)
AI(J,I) = AI(M,I)
AI(M,I) = F
40 CONTINUE
C
50 GO TO (80,130), IEXC
C .......... SEARCH FOR ROWS ISOLATING AN EIGENVALUE
C AND PUSH THEM DOWN ..........
80 IF (L .EQ. 1) GO TO 280
L = L - 1
C .......... FOR J=L STEP -1 UNTIL 1 DO -- ..........
100 DO 120 JJ = 1, L
J = L + 1 - JJ
C
DO 110 I = 1, L
IF (I .EQ. J) GO TO 110
IF (AR(J,I) .NE. 0.0D0 .OR. AI(J,I) .NE. 0.0D0) GO TO 120
110 CONTINUE
C
M = L
IEXC = 1
GO TO 20
120 CONTINUE
C
GO TO 140
C .......... SEARCH FOR COLUMNS ISOLATING AN EIGENVALUE
C AND PUSH THEM LEFT ..........
130 K = K + 1
C
140 DO 170 J = K, L
C
DO 150 I = K, L
IF (I .EQ. J) GO TO 150
IF (AR(I,J) .NE. 0.0D0 .OR. AI(I,J) .NE. 0.0D0) GO TO 170
150 CONTINUE
C
M = K
IEXC = 2
GO TO 20
170 CONTINUE
C .......... NOW BALANCE THE SUBMATRIX IN ROWS K TO L ..........
DO 180 I = K, L
180 SCALE(I) = 1.0D0
C .......... ITERATIVE LOOP FOR NORM REDUCTION ..........
190 NOCONV = .FALSE.
C
DO 270 I = K, L
C = 0.0D0
R = 0.0D0
C
DO 200 J = K, L
IF (J .EQ. I) GO TO 200
C = C + DABS(AR(J,I)) + DABS(AI(J,I))
R = R + DABS(AR(I,J)) + DABS(AI(I,J))
200 CONTINUE
C .......... GUARD AGAINST ZERO C OR R DUE TO UNDERFLOW ..........
IF (C .EQ. 0.0D0 .OR. R .EQ. 0.0D0) GO TO 270
G = R / RADIX
F = 1.0D0
S = C + R
210 IF (C .GE. G) GO TO 220
F = F * RADIX
C = C * B2
GO TO 210
220 G = R * RADIX
230 IF (C .LT. G) GO TO 240
F = F / RADIX
C = C / B2
GO TO 230
C .......... NOW BALANCE ..........
240 IF ((C + R) / F .GE. 0.95D0 * S) GO TO 270
G = 1.0D0 / F
SCALE(I) = SCALE(I) * F
NOCONV = .TRUE.
C
DO 250 J = K, N
AR(I,J) = AR(I,J) * G
AI(I,J) = AI(I,J) * G
250 CONTINUE
C
DO 260 J = 1, L
AR(J,I) = AR(J,I) * F
AI(J,I) = AI(J,I) * F
260 CONTINUE
C
270 CONTINUE
C
IF (NOCONV) GO TO 190
C
280 LOW = K
IGH = L
RETURN
END
SUBROUTINE CDIV(AR,AI,BR,BI,CR,CI)
DOUBLE PRECISION AR,AI,BR,BI,CR,CI
C
C COMPLEX DIVISION, (CR,CI) = (AR,AI)/(BR,BI)
C
DOUBLE PRECISION S,ARS,AIS,BRS,BIS
S = DABS(BR) + DABS(BI)
ARS = AR/S
AIS = AI/S
BRS = BR/S
BIS = BI/S
S = BRS**2 + BIS**2
CR = (ARS*BRS + AIS*BIS)/S
CI = (AIS*BRS - ARS*BIS)/S
RETURN
END
SUBROUTINE COMQR(NM,N,LOW,IGH,HR,HI,WR,WI,IERR)
C
INTEGER I,J,L,N,EN,LL,NM,IGH,ITN,ITS,LOW,LP1,ENM1,IERR
DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N)
DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2,
X PYTHAG
C
C THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
C ALGOL PROCEDURE COMLR, NUM. MATH. 12, 369-376(1968) BY MARTIN
C AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 396-403(1971).
C THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
C (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.
C
C THIS SUBROUTINE FINDS THE EIGENVALUES OF A COMPLEX
C UPPER HESSENBERG MATRIX BY THE QR METHOD.
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 LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED,
C SET LOW=1, IGH=N.
C
C HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
C THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN
C INFORMATION ABOUT THE UNITARY TRANSFORMATIONS USED IN
C THE REDUCTION BY CORTH, IF PERFORMED.
C
C ON OUTPUT
C
C THE UPPER HESSENBERG PORTIONS OF HR AND HI HAVE BEEN
C DESTROYED. THEREFORE, THEY MUST BE SAVED BEFORE
C CALLING COMQR IF SUBSEQUENT CALCULATION OF
C EIGENVECTORS IS TO BE PERFORMED.
C
C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR
C EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
C FOR INDICES IERR+1,...,N.
C
C IERR IS SET TO
C ZERO FOR NORMAL RETURN,
C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
C
C CALLS CDIV FOR COMPLEX DIVISION.
C CALLS CSROOT FOR COMPLEX SQUARE ROOT.
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 L to keep g77 -Wall happy
c
L = 0
C
IERR = 0
IF (LOW .EQ. IGH) GO TO 180
C .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
L = LOW + 1
C
DO 170 I = L, IGH
LL = MIN0(I+1,IGH)
IF (HI(I,I-1) .EQ. 0.0D0) GO TO 170
NORM = PYTHAG(HR(I,I-1),HI(I,I-1))
YR = HR(I,I-1) / NORM
YI = HI(I,I-1) / NORM
HR(I,I-1) = NORM
HI(I,I-1) = 0.0D0
C
DO 155 J = I, IGH
SI = YR * HI(I,J) - YI * HR(I,J)
HR(I,J) = YR * HR(I,J) + YI * HI(I,J)
HI(I,J) = SI
155 CONTINUE
C
DO 160 J = LOW, LL
SI = YR * HI(J,I) + YI * HR(J,I)
HR(J,I) = YR * HR(J,I) - YI * HI(J,I)
HI(J,I) = SI
160 CONTINUE
C
170 CONTINUE
C .......... STORE ROOTS ISOLATED BY CBAL ..........
180 DO 200 I = 1, N
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
WR(I) = HR(I,I)
WI(I) = HI(I,I)
200 CONTINUE
C
EN = IGH
TR = 0.0D0
TI = 0.0D0
ITN = 30*N
C .......... SEARCH FOR NEXT EIGENVALUE ..........
220 IF (EN .LT. LOW) GO TO 1001
ITS = 0
ENM1 = EN - 1
C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
C FOR L=EN STEP -1 UNTIL LOW D0 -- ..........
240 DO 260 LL = LOW, EN
L = EN + LOW - LL
IF (L .EQ. LOW) GO TO 300
TST1 = DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))
X + DABS(HR(L,L)) + DABS(HI(L,L))
TST2 = TST1 + DABS(HR(L,L-1))
IF (TST2 .EQ. TST1) GO TO 300
260 CONTINUE
C .......... FORM SHIFT ..........
300 IF (L .EQ. EN) GO TO 660
IF (ITN .EQ. 0) GO TO 1000
IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
SR = HR(EN,EN)
SI = HI(EN,EN)
XR = HR(ENM1,EN) * HR(EN,ENM1)
XI = HI(ENM1,EN) * HR(EN,ENM1)
IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 340
YR = (HR(ENM1,ENM1) - SR) / 2.0D0
YI = (HI(ENM1,ENM1) - SI) / 2.0D0
CALL CSROOT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)
IF (YR * ZZR + YI * ZZI .GE. 0.0D0) GO TO 310
ZZR = -ZZR
ZZI = -ZZI
310 CALL CDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
SR = SR - XR
SI = SI - XI
GO TO 340
C .......... FORM EXCEPTIONAL SHIFT ..........
320 SR = DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2))
SI = 0.0D0
C
340 DO 360 I = LOW, EN
HR(I,I) = HR(I,I) - SR
HI(I,I) = HI(I,I) - SI
360 CONTINUE
C
TR = TR + SR
TI = TI + SI
ITS = ITS + 1
ITN = ITN - 1
C .......... REDUCE TO TRIANGLE (ROWS) ..........
LP1 = L + 1
C
DO 500 I = LP1, EN
SR = HR(I,I-1)
HR(I,I-1) = 0.0D0
NORM = PYTHAG(PYTHAG(HR(I-1,I-1),HI(I-1,I-1)),SR)
XR = HR(I-1,I-1) / NORM
WR(I-1) = XR
XI = HI(I-1,I-1) / NORM
WI(I-1) = XI
HR(I-1,I-1) = NORM
HI(I-1,I-1) = 0.0D0
HI(I,I-1) = SR / NORM
C
DO 490 J = I, EN
YR = HR(I-1,J)
YI = HI(I-1,J)
ZZR = HR(I,J)
ZZI = HI(I,J)
HR(I-1,J) = XR * YR + XI * YI + HI(I,I-1) * ZZR
HI(I-1,J) = XR * YI - XI * YR + HI(I,I-1) * ZZI
HR(I,J) = XR * ZZR - XI * ZZI - HI(I,I-1) * YR
HI(I,J) = XR * ZZI + XI * ZZR - HI(I,I-1) * YI
490 CONTINUE
C
500 CONTINUE
C
SI = HI(EN,EN)
IF (SI .EQ. 0.0D0) GO TO 540
NORM = PYTHAG(HR(EN,EN),SI)
SR = HR(EN,EN) / NORM
SI = SI / NORM
HR(EN,EN) = NORM
HI(EN,EN) = 0.0D0
C .......... INVERSE OPERATION (COLUMNS) ..........
540 DO 600 J = LP1, EN
XR = WR(J-1)
XI = WI(J-1)
C
DO 580 I = L, J
YR = HR(I,J-1)
YI = 0.0D0
ZZR = HR(I,J)
ZZI = HI(I,J)
IF (I .EQ. J) GO TO 560
YI = HI(I,J-1)
HI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
560 HR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
HR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
HI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
580 CONTINUE
C
600 CONTINUE
C
IF (SI .EQ. 0.0D0) GO TO 240
C
DO 630 I = L, EN
YR = HR(I,EN)
YI = HI(I,EN)
HR(I,EN) = SR * YR - SI * YI
HI(I,EN) = SR * YI + SI * YR
630 CONTINUE
C
GO TO 240
C .......... A ROOT FOUND ..........
660 WR(EN) = HR(EN,EN) + TR
WI(EN) = HI(EN,EN) + TI
EN = ENM1
GO TO 220
C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
C CONVERGED AFTER 30*N ITERATIONS ..........
1000 IERR = EN
1001 RETURN
END
SUBROUTINE COMQR2(NM,N,LOW,IGH,ORTR,ORTI,HR,HI,WR,WI,ZR,ZI,IERR)
C
INTEGER I,J,K,L,M,N,EN,II,JJ,LL,NM,NN,IGH,IP1,
X ITN,ITS,LOW,LP1,ENM1,IEND,IERR
DOUBLE PRECISION HR(NM,N),HI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N),
X ORTR(IGH),ORTI(IGH)
DOUBLE PRECISION SI,SR,TI,TR,XI,XR,YI,YR,ZZI,ZZR,NORM,TST1,TST2,
X PYTHAG
C
C THIS SUBROUTINE IS A TRANSLATION OF A UNITARY ANALOGUE OF THE
C ALGOL PROCEDURE COMLR2, NUM. MATH. 16, 181-204(1970) BY PETERS
C AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
C THE UNITARY ANALOGUE SUBSTITUTES THE QR ALGORITHM OF FRANCIS
C (COMP. JOUR. 4, 332-345(1962)) FOR THE LR ALGORITHM.
C
C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C OF A COMPLEX UPPER HESSENBERG MATRIX BY THE QR
C METHOD. THE EIGENVECTORS OF A COMPLEX GENERAL MATRIX
C CAN ALSO BE FOUND IF CORTH HAS BEEN USED TO REDUCE
C THIS GENERAL MATRIX TO HESSENBERG 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 LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED,
C SET LOW=1, IGH=N.
C
C ORTR AND ORTI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
C FORMATIONS USED IN THE REDUCTION BY CORTH, IF PERFORMED.
C ONLY ELEMENTS LOW THROUGH IGH ARE USED. IF THE EIGENVECTORS
C OF THE HESSENBERG MATRIX ARE DESIRED, SET ORTR(J) AND
C ORTI(J) TO 0.0D0 FOR THESE ELEMENTS.
C
C HR AND HI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX UPPER HESSENBERG MATRIX.
C THEIR LOWER TRIANGLES BELOW THE SUBDIAGONAL CONTAIN FURTHER
C INFORMATION ABOUT THE TRANSFORMATIONS WHICH WERE USED IN THE
C REDUCTION BY CORTH, IF PERFORMED. IF THE EIGENVECTORS OF
C THE HESSENBERG MATRIX ARE DESIRED, THESE ELEMENTS MAY BE
C ARBITRARY.
C
C ON OUTPUT
C
C ORTR, ORTI, AND THE UPPER HESSENBERG PORTIONS OF HR AND HI
C HAVE BEEN DESTROYED.
C
C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVALUES. IF AN ERROR
C EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
C FOR INDICES IERR+1,...,N.
C
C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVECTORS. THE EIGENVECTORS
C ARE UNNORMALIZED. IF AN ERROR EXIT IS MADE, NONE OF
C THE EIGENVECTORS HAS BEEN FOUND.
C
C IERR IS SET TO
C ZERO FOR NORMAL RETURN,
C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
C
C CALLS CDIV FOR COMPLEX DIVISION.
C CALLS CSROOT FOR COMPLEX SQUARE ROOT.
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 L to keep g77 -Wall happy
c
L = 0
C
IERR = 0
C .......... INITIALIZE EIGENVECTOR MATRIX ..........
DO 101 J = 1, N
C
DO 100 I = 1, N
ZR(I,J) = 0.0D0
ZI(I,J) = 0.0D0
100 CONTINUE
ZR(J,J) = 1.0D0
101 CONTINUE
C .......... FORM THE MATRIX OF ACCUMULATED TRANSFORMATIONS
C FROM THE INFORMATION LEFT BY CORTH ..........
IEND = IGH - LOW - 1
C IF (IEND) 180, 150, 105
IF (IEND .GT. 0) GOTO 180
IF (IEND .EQ. 0) GOTO 150
C .......... FOR I=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
105 DO 140 II = 1, IEND
I = IGH - II
IF (ORTR(I) .EQ. 0.0D0 .AND. ORTI(I) .EQ. 0.0D0) GO TO 140
IF (HR(I,I-1) .EQ. 0.0D0 .AND. HI(I,I-1) .EQ. 0.0D0) GO TO 140
C .......... NORM BELOW IS NEGATIVE OF H FORMED IN CORTH ..........
NORM = HR(I,I-1) * ORTR(I) + HI(I,I-1) * ORTI(I)
IP1 = I + 1
C
DO 110 K = IP1, IGH
ORTR(K) = HR(K,I-1)
ORTI(K) = HI(K,I-1)
110 CONTINUE
C
DO 130 J = I, IGH
SR = 0.0D0
SI = 0.0D0
C
DO 115 K = I, IGH
SR = SR + ORTR(K) * ZR(K,J) + ORTI(K) * ZI(K,J)
SI = SI + ORTR(K) * ZI(K,J) - ORTI(K) * ZR(K,J)
115 CONTINUE
C
SR = SR / NORM
SI = SI / NORM
C
DO 120 K = I, IGH
ZR(K,J) = ZR(K,J) + SR * ORTR(K) - SI * ORTI(K)
ZI(K,J) = ZI(K,J) + SR * ORTI(K) + SI * ORTR(K)
120 CONTINUE
C
130 CONTINUE
C
140 CONTINUE
C .......... CREATE REAL SUBDIAGONAL ELEMENTS ..........
150 L = LOW + 1
C
DO 170 I = L, IGH
LL = MIN0(I+1,IGH)
IF (HI(I,I-1) .EQ. 0.0D0) GO TO 170
NORM = PYTHAG(HR(I,I-1),HI(I,I-1))
YR = HR(I,I-1) / NORM
YI = HI(I,I-1) / NORM
HR(I,I-1) = NORM
HI(I,I-1) = 0.0D0
C
DO 155 J = I, N
SI = YR * HI(I,J) - YI * HR(I,J)
HR(I,J) = YR * HR(I,J) + YI * HI(I,J)
HI(I,J) = SI
155 CONTINUE
C
DO 160 J = 1, LL
SI = YR * HI(J,I) + YI * HR(J,I)
HR(J,I) = YR * HR(J,I) - YI * HI(J,I)
HI(J,I) = SI
160 CONTINUE
C
DO 165 J = LOW, IGH
SI = YR * ZI(J,I) + YI * ZR(J,I)
ZR(J,I) = YR * ZR(J,I) - YI * ZI(J,I)
ZI(J,I) = SI
165 CONTINUE
C
170 CONTINUE
C .......... STORE ROOTS ISOLATED BY CBAL ..........
180 DO 200 I = 1, N
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 200
WR(I) = HR(I,I)
WI(I) = HI(I,I)
200 CONTINUE
C
EN = IGH
TR = 0.0D0
TI = 0.0D0
ITN = 30*N
C .......... SEARCH FOR NEXT EIGENVALUE ..........
220 IF (EN .LT. LOW) GO TO 680
ITS = 0
ENM1 = EN - 1
C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
C FOR L=EN STEP -1 UNTIL LOW DO -- ..........
240 DO 260 LL = LOW, EN
L = EN + LOW - LL
IF (L .EQ. LOW) GO TO 300
TST1 = DABS(HR(L-1,L-1)) + DABS(HI(L-1,L-1))
X + DABS(HR(L,L)) + DABS(HI(L,L))
TST2 = TST1 + DABS(HR(L,L-1))
IF (TST2 .EQ. TST1) GO TO 300
260 CONTINUE
C .......... FORM SHIFT ..........
300 IF (L .EQ. EN) GO TO 660
IF (ITN .EQ. 0) GO TO 1000
IF (ITS .EQ. 10 .OR. ITS .EQ. 20) GO TO 320
SR = HR(EN,EN)
SI = HI(EN,EN)
XR = HR(ENM1,EN) * HR(EN,ENM1)
XI = HI(ENM1,EN) * HR(EN,ENM1)
IF (XR .EQ. 0.0D0 .AND. XI .EQ. 0.0D0) GO TO 340
YR = (HR(ENM1,ENM1) - SR) / 2.0D0
YI = (HI(ENM1,ENM1) - SI) / 2.0D0
CALL CSROOT(YR**2-YI**2+XR,2.0D0*YR*YI+XI,ZZR,ZZI)
IF (YR * ZZR + YI * ZZI .GE. 0.0D0) GO TO 310
ZZR = -ZZR
ZZI = -ZZI
310 CALL CDIV(XR,XI,YR+ZZR,YI+ZZI,XR,XI)
SR = SR - XR
SI = SI - XI
GO TO 340
C .......... FORM EXCEPTIONAL SHIFT ..........
320 SR = DABS(HR(EN,ENM1)) + DABS(HR(ENM1,EN-2))
SI = 0.0D0
C
340 DO 360 I = LOW, EN
HR(I,I) = HR(I,I) - SR
HI(I,I) = HI(I,I) - SI
360 CONTINUE
C
TR = TR + SR
TI = TI + SI
ITS = ITS + 1
ITN = ITN - 1
C .......... REDUCE TO TRIANGLE (ROWS) ..........
LP1 = L + 1
C
DO 500 I = LP1, EN
SR = HR(I,I-1)
HR(I,I-1) = 0.0D0
NORM = PYTHAG(PYTHAG(HR(I-1,I-1),HI(I-1,I-1)),SR)
XR = HR(I-1,I-1) / NORM
WR(I-1) = XR
XI = HI(I-1,I-1) / NORM
WI(I-1) = XI
HR(I-1,I-1) = NORM
HI(I-1,I-1) = 0.0D0
HI(I,I-1) = SR / NORM
C
DO 490 J = I, N
YR = HR(I-1,J)
YI = HI(I-1,J)
ZZR = HR(I,J)
ZZI = HI(I,J)
HR(I-1,J) = XR * YR + XI * YI + HI(I,I-1) * ZZR
HI(I-1,J) = XR * YI - XI * YR + HI(I,I-1) * ZZI
HR(I,J) = XR * ZZR - XI * ZZI - HI(I,I-1) * YR
HI(I,J) = XR * ZZI + XI * ZZR - HI(I,I-1) * YI
490 CONTINUE
C
500 CONTINUE
C
SI = HI(EN,EN)
IF (SI .EQ. 0.0D0) GO TO 540
NORM = PYTHAG(HR(EN,EN),SI)
SR = HR(EN,EN) / NORM
SI = SI / NORM
HR(EN,EN) = NORM
HI(EN,EN) = 0.0D0
IF (EN .EQ. N) GO TO 540
IP1 = EN + 1
C
DO 520 J = IP1, N
YR = HR(EN,J)
YI = HI(EN,J)
HR(EN,J) = SR * YR + SI * YI
HI(EN,J) = SR * YI - SI * YR
520 CONTINUE
C .......... INVERSE OPERATION (COLUMNS) ..........
540 DO 600 J = LP1, EN
XR = WR(J-1)
XI = WI(J-1)
C
DO 580 I = 1, J
YR = HR(I,J-1)
YI = 0.0D0
ZZR = HR(I,J)
ZZI = HI(I,J)
IF (I .EQ. J) GO TO 560
YI = HI(I,J-1)
HI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
560 HR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
HR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
HI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
580 CONTINUE
C
DO 590 I = LOW, IGH
YR = ZR(I,J-1)
YI = ZI(I,J-1)
ZZR = ZR(I,J)
ZZI = ZI(I,J)
ZR(I,J-1) = XR * YR - XI * YI + HI(J,J-1) * ZZR
ZI(I,J-1) = XR * YI + XI * YR + HI(J,J-1) * ZZI
ZR(I,J) = XR * ZZR + XI * ZZI - HI(J,J-1) * YR
ZI(I,J) = XR * ZZI - XI * ZZR - HI(J,J-1) * YI
590 CONTINUE
C
600 CONTINUE
C
IF (SI .EQ. 0.0D0) GO TO 240
C
DO 630 I = 1, EN
YR = HR(I,EN)
YI = HI(I,EN)
HR(I,EN) = SR * YR - SI * YI
HI(I,EN) = SR * YI + SI * YR
630 CONTINUE
C
DO 640 I = LOW, IGH
YR = ZR(I,EN)
YI = ZI(I,EN)
ZR(I,EN) = SR * YR - SI * YI
ZI(I,EN) = SR * YI + SI * YR
640 CONTINUE
C
GO TO 240
C .......... A ROOT FOUND ..........
660 HR(EN,EN) = HR(EN,EN) + TR
WR(EN) = HR(EN,EN)
HI(EN,EN) = HI(EN,EN) + TI
WI(EN) = HI(EN,EN)
EN = ENM1
GO TO 220
C .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND
C VECTORS OF UPPER TRIANGULAR FORM ..........
680 NORM = 0.0D0
C
DO 720 I = 1, N
C
DO 720 J = I, N
TR = DABS(HR(I,J)) + DABS(HI(I,J))
IF (TR .GT. NORM) NORM = TR
720 CONTINUE
C
IF (N .EQ. 1 .OR. NORM .EQ. 0.0D0) GO TO 1001
C .......... FOR EN=N STEP -1 UNTIL 2 DO -- ..........
DO 800 NN = 2, N
EN = N + 2 - NN
XR = WR(EN)
XI = WI(EN)
HR(EN,EN) = 1.0D0
HI(EN,EN) = 0.0D0
ENM1 = EN - 1
C .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
DO 780 II = 1, ENM1
I = EN - II
ZZR = 0.0D0
ZZI = 0.0D0
IP1 = I + 1
C
DO 740 J = IP1, EN
ZZR = ZZR + HR(I,J) * HR(J,EN) - HI(I,J) * HI(J,EN)
ZZI = ZZI + HR(I,J) * HI(J,EN) + HI(I,J) * HR(J,EN)
740 CONTINUE
C
YR = XR - WR(I)
YI = XI - WI(I)
IF (YR .NE. 0.0D0 .OR. YI .NE. 0.0D0) GO TO 765
TST1 = NORM
YR = TST1
760 YR = 0.01D0 * YR
TST2 = NORM + YR
IF (TST2 .GT. TST1) GO TO 760
765 CONTINUE
CALL CDIV(ZZR,ZZI,YR,YI,HR(I,EN),HI(I,EN))
C .......... OVERFLOW CONTROL ..........
TR = DABS(HR(I,EN)) + DABS(HI(I,EN))
IF (TR .EQ. 0.0D0) GO TO 780
TST1 = TR
TST2 = TST1 + 1.0D0/TST1
IF (TST2 .GT. TST1) GO TO 780
DO 770 J = I, EN
HR(J,EN) = HR(J,EN)/TR
HI(J,EN) = HI(J,EN)/TR
770 CONTINUE
C
780 CONTINUE
C
800 CONTINUE
C .......... END BACKSUBSTITUTION ..........
ENM1 = N - 1
C .......... VECTORS OF ISOLATED ROOTS ..........
DO 840 I = 1, ENM1
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
IP1 = I + 1
C
DO 820 J = IP1, N
ZR(I,J) = HR(I,J)
ZI(I,J) = HI(I,J)
820 CONTINUE
C
840 CONTINUE
C .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
C VECTORS OF ORIGINAL FULL MATRIX.
C FOR J=N STEP -1 UNTIL LOW+1 DO -- ..........
DO 880 JJ = LOW, ENM1
J = N + LOW - JJ
M = MIN0(J,IGH)
C
DO 880 I = LOW, IGH
ZZR = 0.0D0
ZZI = 0.0D0
C
DO 860 K = LOW, M
ZZR = ZZR + ZR(I,K) * HR(K,J) - ZI(I,K) * HI(K,J)
ZZI = ZZI + ZR(I,K) * HI(K,J) + ZI(I,K) * HR(K,J)
860 CONTINUE
C
ZR(I,J) = ZZR
ZI(I,J) = ZZI
880 CONTINUE
C
GO TO 1001
C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
C CONVERGED AFTER 30*N ITERATIONS ..........
1000 IERR = EN
1001 RETURN
END
SUBROUTINE CORTH(NM,N,LOW,IGH,AR,AI,ORTR,ORTI)
C
INTEGER I,J,M,N,II,JJ,LA,MP,NM,IGH,KP1,LOW
DOUBLE PRECISION AR(NM,N),AI(NM,N),ORTR(IGH),ORTI(IGH)
DOUBLE PRECISION F,G,H,FI,FR,SCALE,PYTHAG
C
C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
C THE ALGOL PROCEDURE ORTHES, NUM. MATH. 12, 349-368(1968)
C BY MARTIN AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
C
C GIVEN A COMPLEX GENERAL MATRIX, THIS SUBROUTINE
C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
C UNITARY 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 LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C SUBROUTINE CBAL. IF CBAL HAS NOT BEEN USED,
C SET LOW=1, IGH=N.
C
C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX INPUT MATRIX.
C
C ON OUTPUT
C
C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE HESSENBERG MATRIX. INFORMATION
C ABOUT THE UNITARY TRANSFORMATIONS USED IN THE REDUCTION
C IS STORED IN THE REMAINING TRIANGLES UNDER THE
C HESSENBERG MATRIX.
C
C ORTR AND ORTI CONTAIN FURTHER INFORMATION ABOUT THE
C TRANSFORMATIONS. ONLY ELEMENTS LOW THROUGH IGH ARE USED.
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
LA = IGH - 1
KP1 = LOW + 1
IF (LA .LT. KP1) GO TO 200
C
DO 180 M = KP1, LA
H = 0.0D0
ORTR(M) = 0.0D0
ORTI(M) = 0.0D0
SCALE = 0.0D0
C .......... SCALE COLUMN (ALGOL TOL THEN NOT NEEDED) ..........
DO 90 I = M, IGH
90 SCALE = SCALE + DABS(AR(I,M-1)) + DABS(AI(I,M-1))
C
IF (SCALE .EQ. 0.0D0) GO TO 180
MP = M + IGH
C .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
DO 100 II = M, IGH
I = MP - II
ORTR(I) = AR(I,M-1) / SCALE
ORTI(I) = AI(I,M-1) / SCALE
H = H + ORTR(I) * ORTR(I) + ORTI(I) * ORTI(I)
100 CONTINUE
C
G = DSQRT(H)
F = PYTHAG(ORTR(M),ORTI(M))
IF (F .EQ. 0.0D0) GO TO 103
H = H + F * G
G = G / F
ORTR(M) = (1.0D0 + G) * ORTR(M)
ORTI(M) = (1.0D0 + G) * ORTI(M)
GO TO 105
C
103 ORTR(M) = G
AR(M,M-1) = SCALE
C .......... FORM (I-(U*UT)/H) * A ..........
105 DO 130 J = M, N
FR = 0.0D0
FI = 0.0D0
C .......... FOR I=IGH STEP -1 UNTIL M DO -- ..........
DO 110 II = M, IGH
I = MP - II
FR = FR + ORTR(I) * AR(I,J) + ORTI(I) * AI(I,J)
FI = FI + ORTR(I) * AI(I,J) - ORTI(I) * AR(I,J)
110 CONTINUE
C
FR = FR / H
FI = FI / H
C
DO 120 I = M, IGH
AR(I,J) = AR(I,J) - FR * ORTR(I) + FI * ORTI(I)
AI(I,J) = AI(I,J) - FR * ORTI(I) - FI * ORTR(I)
120 CONTINUE
C
130 CONTINUE
C .......... FORM (I-(U*UT)/H)*A*(I-(U*UT)/H) ..........
DO 160 I = 1, IGH
FR = 0.0D0
FI = 0.0D0
C .......... FOR J=IGH STEP -1 UNTIL M DO -- ..........
DO 140 JJ = M, IGH
J = MP - JJ
FR = FR + ORTR(J) * AR(I,J) - ORTI(J) * AI(I,J)
FI = FI + ORTR(J) * AI(I,J) + ORTI(J) * AR(I,J)
140 CONTINUE
C
FR = FR / H
FI = FI / H
C
DO 150 J = M, IGH
AR(I,J) = AR(I,J) - FR * ORTR(J) - FI * ORTI(J)
AI(I,J) = AI(I,J) + FR * ORTI(J) - FI * ORTR(J)
150 CONTINUE
C
160 CONTINUE
C
ORTR(M) = SCALE * ORTR(M)
ORTI(M) = SCALE * ORTI(M)
AR(M,M-1) = -G * AR(M,M-1)
AI(M,M-1) = -G * AI(M,M-1)
180 CONTINUE
C
200 RETURN
END
SUBROUTINE CSROOT(XR,XI,YR,YI)
DOUBLE PRECISION XR,XI,YR,YI
C
C (YR,YI) = COMPLEX DSQRT(XR,XI)
C BRANCH CHOSEN SO THAT YR .GE. 0.0 AND SIGN(YI) .EQ. SIGN(XI)
C
DOUBLE PRECISION S,TR,TI,PYTHAG
TR = XR
TI = XI
S = DSQRT(0.5D0*(PYTHAG(TR,TI) + DABS(TR)))
IF (TR .GE. 0.0D0) YR = S
IF (TI .LT. 0.0D0) S = -S
IF (TR .LE. 0.0D0) YI = S
IF (TR .LT. 0.0D0) YR = 0.5D0*(TI/YI)
IF (TR .GT. 0.0D0) YI = 0.5D0*(TI/YR)
RETURN
END
SUBROUTINE ELMHES(NM,N,LOW,IGH,A,INT)
C
INTEGER I,J,M,N,LA,NM,IGH,KP1,LOW,MM1,MP1
DOUBLE PRECISION A(NM,N)
DOUBLE PRECISION X,Y
INTEGER INT(IGH)
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMHES,
C NUM. MATH. 12, 349-368(1968) BY MARTIN AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 339-358(1971).
C
C GIVEN A REAL GENERAL MATRIX, THIS SUBROUTINE
C REDUCES A SUBMATRIX SITUATED IN ROWS AND COLUMNS
C LOW THROUGH IGH TO UPPER HESSENBERG FORM BY
C STABILIZED ELEMENTARY 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 LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED,
C SET LOW=1, IGH=N.
C
C A CONTAINS THE INPUT MATRIX.
C
C ON OUTPUT
C
C A CONTAINS THE HESSENBERG MATRIX. THE MULTIPLIERS
C WHICH WERE USED IN THE REDUCTION ARE STORED IN THE
C REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX.
C
C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
C INTERCHANGED IN THE REDUCTION.
C ONLY ELEMENTS LOW THROUGH IGH ARE USED.
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
LA = IGH - 1
KP1 = LOW + 1
IF (LA .LT. KP1) GO TO 200
C
DO 180 M = KP1, LA
MM1 = M - 1
X = 0.0D0
I = M
C
DO 100 J = M, IGH
IF (DABS(A(J,MM1)) .LE. DABS(X)) GO TO 100
X = A(J,MM1)
I = J
100 CONTINUE
C
INT(M) = I
IF (I .EQ. M) GO TO 130
C .......... INTERCHANGE ROWS AND COLUMNS OF A ..........
DO 110 J = MM1, N
Y = A(I,J)
A(I,J) = A(M,J)
A(M,J) = Y
110 CONTINUE
C
DO 120 J = 1, IGH
Y = A(J,I)
A(J,I) = A(J,M)
A(J,M) = Y
120 CONTINUE
C .......... END INTERCHANGE ..........
130 IF (X .EQ. 0.0D0) GO TO 180
MP1 = M + 1
C
DO 160 I = MP1, IGH
Y = A(I,MM1)
IF (Y .EQ. 0.0D0) GO TO 160
Y = Y / X
A(I,MM1) = Y
C
DO 140 J = M, N
140 A(I,J) = A(I,J) - Y * A(M,J)
C
DO 150 J = 1, IGH
150 A(J,M) = A(J,M) + Y * A(J,I)
C
160 CONTINUE
C
180 CONTINUE
C
200 RETURN
END
SUBROUTINE ELTRAN(NM,N,LOW,IGH,A,INT,Z)
C
INTEGER I,J,N,KL,MM,MP,NM,IGH,LOW,MP1
DOUBLE PRECISION A(NM,IGH),Z(NM,N)
INTEGER INT(IGH)
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE ELMTRANS,
C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
C
C THIS SUBROUTINE ACCUMULATES THE STABILIZED ELEMENTARY
C SIMILARITY TRANSFORMATIONS USED IN THE REDUCTION OF A
C REAL GENERAL MATRIX TO UPPER HESSENBERG FORM BY ELMHES.
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 LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED,
C SET LOW=1, IGH=N.
C
C A CONTAINS THE MULTIPLIERS WHICH WERE USED IN THE
C REDUCTION BY ELMHES IN ITS LOWER TRIANGLE
C BELOW THE SUBDIAGONAL.
C
C INT CONTAINS INFORMATION ON THE ROWS AND COLUMNS
C INTERCHANGED IN THE REDUCTION BY ELMHES.
C ONLY ELEMENTS LOW THROUGH IGH ARE USED.
C
C ON OUTPUT
C
C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED IN THE
C REDUCTION BY ELMHES.
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 .......... INITIALIZE Z TO IDENTITY MATRIX ..........
DO 80 J = 1, N
C
DO 60 I = 1, N
60 Z(I,J) = 0.0D0
C
Z(J,J) = 1.0D0
80 CONTINUE
C
KL = IGH - LOW - 1
IF (KL .LT. 1) GO TO 200
C .......... FOR MP=IGH-1 STEP -1 UNTIL LOW+1 DO -- ..........
DO 140 MM = 1, KL
MP = IGH - MM
MP1 = MP + 1
C
DO 100 I = MP1, IGH
100 Z(I,MP) = A(I,MP-1)
C
I = INT(MP)
IF (I .EQ. MP) GO TO 140
C
DO 130 J = MP, IGH
Z(MP,J) = Z(I,J)
Z(I,J) = 0.0D0
130 CONTINUE
C
Z(I,MP) = 1.0D0
140 CONTINUE
C
200 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
SUBROUTINE HQR(NM,N,LOW,IGH,H,WR,WI,IERR)
C
INTEGER I,J,K,L,M,N,EN,LL,MM,NA,NM,IGH,ITN,ITS,LOW,MP2,ENM2,IERR
DOUBLE PRECISION H(NM,N),WR(N),WI(N)
DOUBLE PRECISION P,Q,R,S,T,W,X,Y,ZZ,NORM,TST1,TST2
LOGICAL NOTLAS
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR,
C NUM. MATH. 14, 219-231(1970) BY MARTIN, PETERS, AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 359-371(1971).
C
C THIS SUBROUTINE FINDS THE EIGENVALUES OF A REAL
C UPPER HESSENBERG MATRIX BY THE QR METHOD.
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 LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED,
C SET LOW=1, IGH=N.
C
C H CONTAINS THE UPPER HESSENBERG MATRIX. INFORMATION ABOUT
C THE TRANSFORMATIONS USED IN THE REDUCTION TO HESSENBERG
C FORM BY ELMHES OR ORTHES, IF PERFORMED, IS STORED
C IN THE REMAINING TRIANGLE UNDER THE HESSENBERG MATRIX.
C
C ON OUTPUT
C
C H HAS BEEN DESTROYED. THEREFORE, IT MUST BE SAVED
C BEFORE CALLING HQR IF SUBSEQUENT CALCULATION AND
C BACK TRANSFORMATION OF EIGENVECTORS IS TO BE PERFORMED.
C
C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES
C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS
C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE
C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN
C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
C FOR INDICES IERR+1,...,N.
C
C IERR IS SET TO
C ZERO FOR NORMAL RETURN,
C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
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 L M P Q R to keep g77 -Wall happy
c
L = 0
M = 0
P = 0.0D0
Q = 0.0D0
R = 0.0D0
C
IERR = 0
NORM = 0.0D0
K = 1
C .......... STORE ROOTS ISOLATED BY BALANC
C AND COMPUTE MATRIX NORM ..........
DO 50 I = 1, N
C
DO 40 J = K, N
40 NORM = NORM + DABS(H(I,J))
C
K = I
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
WR(I) = H(I,I)
WI(I) = 0.0D0
50 CONTINUE
C
EN = IGH
T = 0.0D0
ITN = 30*N
C .......... SEARCH FOR NEXT EIGENVALUES ..........
60 IF (EN .LT. LOW) GO TO 1001
ITS = 0
NA = EN - 1
ENM2 = NA - 1
C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
C FOR L=EN STEP -1 UNTIL LOW DO -- ..........
70 DO 80 LL = LOW, EN
L = EN + LOW - LL
IF (L .EQ. LOW) GO TO 100
S = DABS(H(L-1,L-1)) + DABS(H(L,L))
IF (S .EQ. 0.0D0) S = NORM
TST1 = S
TST2 = TST1 + DABS(H(L,L-1))
IF (TST2 .EQ. TST1) GO TO 100
80 CONTINUE
C .......... FORM SHIFT ..........
100 X = H(EN,EN)
IF (L .EQ. EN) GO TO 270
Y = H(NA,NA)
W = H(EN,NA) * H(NA,EN)
IF (L .EQ. NA) GO TO 280
IF (ITN .EQ. 0) GO TO 1000
IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
C .......... FORM EXCEPTIONAL SHIFT ..........
T = T + X
C
DO 120 I = LOW, EN
120 H(I,I) = H(I,I) - X
C
S = DABS(H(EN,NA)) + DABS(H(NA,ENM2))
X = 0.75D0 * S
Y = X
W = -0.4375D0 * S * S
130 ITS = ITS + 1
ITN = ITN - 1
C .......... LOOK FOR TWO CONSECUTIVE SMALL
C SUB-DIAGONAL ELEMENTS.
C FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
DO 140 MM = L, ENM2
M = ENM2 + L - MM
ZZ = H(M,M)
R = X - ZZ
S = Y - ZZ
P = (R * S - W) / H(M+1,M) + H(M,M+1)
Q = H(M+1,M+1) - ZZ - R - S
R = H(M+2,M+1)
S = DABS(P) + DABS(Q) + DABS(R)
P = P / S
Q = Q / S
R = R / S
IF (M .EQ. L) GO TO 150
TST1 = DABS(P)*(DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1)))
TST2 = TST1 + DABS(H(M,M-1))*(DABS(Q) + DABS(R))
IF (TST2 .EQ. TST1) GO TO 150
140 CONTINUE
C
150 MP2 = M + 2
C
DO 160 I = MP2, EN
H(I,I-2) = 0.0D0
IF (I .EQ. MP2) GO TO 160
H(I,I-3) = 0.0D0
160 CONTINUE
C .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
C COLUMNS M TO EN ..........
DO 260 K = M, NA
NOTLAS = K .NE. NA
IF (K .EQ. M) GO TO 170
P = H(K,K-1)
Q = H(K+1,K-1)
R = 0.0D0
IF (NOTLAS) R = H(K+2,K-1)
X = DABS(P) + DABS(Q) + DABS(R)
IF (X .EQ. 0.0D0) GO TO 260
P = P / X
Q = Q / X
R = R / X
170 S = DSIGN(DSQRT(P*P+Q*Q+R*R),P)
IF (K .EQ. M) GO TO 180
H(K,K-1) = -S * X
GO TO 190
180 IF (L .NE. M) H(K,K-1) = -H(K,K-1)
190 P = P + S
X = P / S
Y = Q / S
ZZ = R / S
Q = Q / P
R = R / P
IF (NOTLAS) GO TO 225
C .......... ROW MODIFICATION ..........
DO 200 J = K, N
P = H(K,J) + Q * H(K+1,J)
H(K,J) = H(K,J) - P * X
H(K+1,J) = H(K+1,J) - P * Y
200 CONTINUE
C
J = MIN0(EN,K+3)
C .......... COLUMN MODIFICATION ..........
DO 210 I = 1, J
P = X * H(I,K) + Y * H(I,K+1)
H(I,K) = H(I,K) - P
H(I,K+1) = H(I,K+1) - P * Q
210 CONTINUE
GO TO 255
225 CONTINUE
C .......... ROW MODIFICATION ..........
DO 230 J = K, N
P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J)
H(K,J) = H(K,J) - P * X
H(K+1,J) = H(K+1,J) - P * Y
H(K+2,J) = H(K+2,J) - P * ZZ
230 CONTINUE
C
J = MIN0(EN,K+3)
C .......... COLUMN MODIFICATION ..........
DO 240 I = 1, J
P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2)
H(I,K) = H(I,K) - P
H(I,K+1) = H(I,K+1) - P * Q
H(I,K+2) = H(I,K+2) - P * R
240 CONTINUE
255 CONTINUE
C
260 CONTINUE
C
GO TO 70
C .......... ONE ROOT FOUND ..........
270 WR(EN) = X + T
WI(EN) = 0.0D0
EN = NA
GO TO 60
C .......... TWO ROOTS FOUND ..........
280 P = (Y - X) / 2.0D0
Q = P * P + W
ZZ = DSQRT(DABS(Q))
X = X + T
IF (Q .LT. 0.0D0) GO TO 320
C .......... REAL PAIR ..........
ZZ = P + DSIGN(ZZ,P)
WR(NA) = X + ZZ
WR(EN) = WR(NA)
IF (ZZ .NE. 0.0D0) WR(EN) = X - W / ZZ
WI(NA) = 0.0D0
WI(EN) = 0.0D0
GO TO 330
C .......... COMPLEX PAIR ..........
320 WR(NA) = X + P
WR(EN) = X + P
WI(NA) = ZZ
WI(EN) = -ZZ
330 EN = ENM2
GO TO 60
C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
C CONVERGED AFTER 30*N ITERATIONS ..........
1000 IERR = EN
1001 RETURN
END
SUBROUTINE HQR2(NM,N,LOW,IGH,H,WR,WI,Z,IERR)
C
INTEGER I,J,K,L,M,N,EN,II,JJ,LL,MM,NA,NM,NN,
X IGH,ITN,ITS,LOW,MP2,ENM2,IERR
DOUBLE PRECISION H(NM,N),WR(N),WI(N),Z(NM,N)
DOUBLE PRECISION P,Q,R,S,T,W,X,Y,RA,SA,VI,VR,ZZ,NORM,TST1,TST2
LOGICAL NOTLAS
C
C THIS SUBROUTINE IS A TRANSLATION OF THE ALGOL PROCEDURE HQR2,
C NUM. MATH. 16, 181-204(1970) BY PETERS AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 372-395(1971).
C
C THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS
C OF A REAL UPPER HESSENBERG MATRIX BY THE QR METHOD. THE
C EIGENVECTORS OF A REAL GENERAL MATRIX CAN ALSO BE FOUND
C IF ELMHES AND ELTRAN OR ORTHES AND ORTRAN HAVE
C BEEN USED TO REDUCE THIS GENERAL MATRIX TO HESSENBERG FORM
C AND TO ACCUMULATE THE 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 LOW AND IGH ARE INTEGERS DETERMINED BY THE BALANCING
C SUBROUTINE BALANC. IF BALANC HAS NOT BEEN USED,
C SET LOW=1, IGH=N.
C
C H CONTAINS THE UPPER HESSENBERG MATRIX.
C
C Z CONTAINS THE TRANSFORMATION MATRIX PRODUCED BY ELTRAN
C AFTER THE REDUCTION BY ELMHES, OR BY ORTRAN AFTER THE
C REDUCTION BY ORTHES, IF PERFORMED. IF THE EIGENVECTORS
C OF THE HESSENBERG MATRIX ARE DESIRED, Z MUST CONTAIN THE
C IDENTITY MATRIX.
C
C ON OUTPUT
C
C H HAS BEEN DESTROYED.
C
C WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVALUES. THE EIGENVALUES
C ARE UNORDERED EXCEPT THAT COMPLEX CONJUGATE PAIRS
C OF VALUES APPEAR CONSECUTIVELY WITH THE EIGENVALUE
C HAVING THE POSITIVE IMAGINARY PART FIRST. IF AN
C ERROR EXIT IS MADE, THE EIGENVALUES SHOULD BE CORRECT
C FOR INDICES IERR+1,...,N.
C
C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS.
C IF THE I-TH EIGENVALUE IS REAL, THE I-TH COLUMN OF Z
C CONTAINS ITS EIGENVECTOR. IF THE I-TH EIGENVALUE IS COMPLEX
C WITH POSITIVE IMAGINARY PART, THE I-TH AND (I+1)-TH
C COLUMNS OF Z CONTAIN THE REAL AND IMAGINARY PARTS OF ITS
C EIGENVECTOR. THE EIGENVECTORS ARE UNNORMALIZED. IF AN
C ERROR EXIT IS MADE, NONE OF THE EIGENVECTORS HAS BEEN FOUND.
C
C IERR IS SET TO
C ZERO FOR NORMAL RETURN,
C J IF THE LIMIT OF 30*N ITERATIONS IS EXHAUSTED
C WHILE THE J-TH EIGENVALUE IS BEING SOUGHT.
C
C CALLS CDIV FOR COMPLEX DIVISION.
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 L M P R S to keep g77 -Wall happy
c
L = 0
M = 0
P = 0.0D0
R = 0.0D0
S = 0.0D0
C
IERR = 0
NORM = 0.0D0
K = 1
C .......... STORE ROOTS ISOLATED BY BALANC
C AND COMPUTE MATRIX NORM ..........
DO 50 I = 1, N
C
DO 40 J = K, N
40 NORM = NORM + DABS(H(I,J))
C
K = I
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 50
WR(I) = H(I,I)
WI(I) = 0.0D0
50 CONTINUE
C
EN = IGH
T = 0.0D0
ITN = 30*N
C .......... SEARCH FOR NEXT EIGENVALUES ..........
60 IF (EN .LT. LOW) GO TO 340
ITS = 0
NA = EN - 1
ENM2 = NA - 1
C .......... LOOK FOR SINGLE SMALL SUB-DIAGONAL ELEMENT
C FOR L=EN STEP -1 UNTIL LOW DO -- ..........
70 DO 80 LL = LOW, EN
L = EN + LOW - LL
IF (L .EQ. LOW) GO TO 100
S = DABS(H(L-1,L-1)) + DABS(H(L,L))
IF (S .EQ. 0.0D0) S = NORM
TST1 = S
TST2 = TST1 + DABS(H(L,L-1))
IF (TST2 .EQ. TST1) GO TO 100
80 CONTINUE
C .......... FORM SHIFT ..........
100 X = H(EN,EN)
IF (L .EQ. EN) GO TO 270
Y = H(NA,NA)
W = H(EN,NA) * H(NA,EN)
IF (L .EQ. NA) GO TO 280
IF (ITN .EQ. 0) GO TO 1000
IF (ITS .NE. 10 .AND. ITS .NE. 20) GO TO 130
C .......... FORM EXCEPTIONAL SHIFT ..........
T = T + X
C
DO 120 I = LOW, EN
120 H(I,I) = H(I,I) - X
C
S = DABS(H(EN,NA)) + DABS(H(NA,ENM2))
X = 0.75D0 * S
Y = X
W = -0.4375D0 * S * S
130 ITS = ITS + 1
ITN = ITN - 1
C .......... LOOK FOR TWO CONSECUTIVE SMALL
C SUB-DIAGONAL ELEMENTS.
C FOR M=EN-2 STEP -1 UNTIL L DO -- ..........
DO 140 MM = L, ENM2
M = ENM2 + L - MM
ZZ = H(M,M)
R = X - ZZ
S = Y - ZZ
P = (R * S - W) / H(M+1,M) + H(M,M+1)
Q = H(M+1,M+1) - ZZ - R - S
R = H(M+2,M+1)
S = DABS(P) + DABS(Q) + DABS(R)
P = P / S
Q = Q / S
R = R / S
IF (M .EQ. L) GO TO 150
TST1 = DABS(P)*(DABS(H(M-1,M-1)) + DABS(ZZ) + DABS(H(M+1,M+1)))
TST2 = TST1 + DABS(H(M,M-1))*(DABS(Q) + DABS(R))
IF (TST2 .EQ. TST1) GO TO 150
140 CONTINUE
C
150 MP2 = M + 2
C
DO 160 I = MP2, EN
H(I,I-2) = 0.0D0
IF (I .EQ. MP2) GO TO 160
H(I,I-3) = 0.0D0
160 CONTINUE
C .......... DOUBLE QR STEP INVOLVING ROWS L TO EN AND
C COLUMNS M TO EN ..........
DO 260 K = M, NA
NOTLAS = K .NE. NA
IF (K .EQ. M) GO TO 170
P = H(K,K-1)
Q = H(K+1,K-1)
R = 0.0D0
IF (NOTLAS) R = H(K+2,K-1)
X = DABS(P) + DABS(Q) + DABS(R)
IF (X .EQ. 0.0D0) GO TO 260
P = P / X
Q = Q / X
R = R / X
170 S = DSIGN(DSQRT(P*P+Q*Q+R*R),P)
IF (K .EQ. M) GO TO 180
H(K,K-1) = -S * X
GO TO 190
180 IF (L .NE. M) H(K,K-1) = -H(K,K-1)
190 P = P + S
X = P / S
Y = Q / S
ZZ = R / S
Q = Q / P
R = R / P
IF (NOTLAS) GO TO 225
C .......... ROW MODIFICATION ..........
DO 200 J = K, N
P = H(K,J) + Q * H(K+1,J)
H(K,J) = H(K,J) - P * X
H(K+1,J) = H(K+1,J) - P * Y
200 CONTINUE
C
J = MIN0(EN,K+3)
C .......... COLUMN MODIFICATION ..........
DO 210 I = 1, J
P = X * H(I,K) + Y * H(I,K+1)
H(I,K) = H(I,K) - P
H(I,K+1) = H(I,K+1) - P * Q
210 CONTINUE
C .......... ACCUMULATE TRANSFORMATIONS ..........
DO 220 I = LOW, IGH
P = X * Z(I,K) + Y * Z(I,K+1)
Z(I,K) = Z(I,K) - P
Z(I,K+1) = Z(I,K+1) - P * Q
220 CONTINUE
GO TO 255
225 CONTINUE
C .......... ROW MODIFICATION ..........
DO 230 J = K, N
P = H(K,J) + Q * H(K+1,J) + R * H(K+2,J)
H(K,J) = H(K,J) - P * X
H(K+1,J) = H(K+1,J) - P * Y
H(K+2,J) = H(K+2,J) - P * ZZ
230 CONTINUE
C
J = MIN0(EN,K+3)
C .......... COLUMN MODIFICATION ..........
DO 240 I = 1, J
P = X * H(I,K) + Y * H(I,K+1) + ZZ * H(I,K+2)
H(I,K) = H(I,K) - P
H(I,K+1) = H(I,K+1) - P * Q
H(I,K+2) = H(I,K+2) - P * R
240 CONTINUE
C .......... ACCUMULATE TRANSFORMATIONS ..........
DO 250 I = LOW, IGH
P = X * Z(I,K) + Y * Z(I,K+1) + ZZ * Z(I,K+2)
Z(I,K) = Z(I,K) - P
Z(I,K+1) = Z(I,K+1) - P * Q
Z(I,K+2) = Z(I,K+2) - P * R
250 CONTINUE
255 CONTINUE
C
260 CONTINUE
C
GO TO 70
C .......... ONE ROOT FOUND ..........
270 H(EN,EN) = X + T
WR(EN) = H(EN,EN)
WI(EN) = 0.0D0
EN = NA
GO TO 60
C .......... TWO ROOTS FOUND ..........
280 P = (Y - X) / 2.0D0
Q = P * P + W
ZZ = DSQRT(DABS(Q))
H(EN,EN) = X + T
X = H(EN,EN)
H(NA,NA) = Y + T
IF (Q .LT. 0.0D0) GO TO 320
C .......... REAL PAIR ..........
ZZ = P + DSIGN(ZZ,P)
WR(NA) = X + ZZ
WR(EN) = WR(NA)
IF (ZZ .NE. 0.0D0) WR(EN) = X - W / ZZ
WI(NA) = 0.0D0
WI(EN) = 0.0D0
X = H(EN,NA)
S = DABS(X) + DABS(ZZ)
P = X / S
Q = ZZ / S
R = DSQRT(P*P+Q*Q)
P = P / R
Q = Q / R
C .......... ROW MODIFICATION ..........
DO 290 J = NA, N
ZZ = H(NA,J)
H(NA,J) = Q * ZZ + P * H(EN,J)
H(EN,J) = Q * H(EN,J) - P * ZZ
290 CONTINUE
C .......... COLUMN MODIFICATION ..........
DO 300 I = 1, EN
ZZ = H(I,NA)
H(I,NA) = Q * ZZ + P * H(I,EN)
H(I,EN) = Q * H(I,EN) - P * ZZ
300 CONTINUE
C .......... ACCUMULATE TRANSFORMATIONS ..........
DO 310 I = LOW, IGH
ZZ = Z(I,NA)
Z(I,NA) = Q * ZZ + P * Z(I,EN)
Z(I,EN) = Q * Z(I,EN) - P * ZZ
310 CONTINUE
C
GO TO 330
C .......... COMPLEX PAIR ..........
320 WR(NA) = X + P
WR(EN) = X + P
WI(NA) = ZZ
WI(EN) = -ZZ
330 EN = ENM2
GO TO 60
C .......... ALL ROOTS FOUND. BACKSUBSTITUTE TO FIND
C VECTORS OF UPPER TRIANGULAR FORM ..........
340 IF (NORM .EQ. 0.0D0) GO TO 1001
C .......... FOR EN=N STEP -1 UNTIL 1 DO -- ..........
DO 800 NN = 1, N
EN = N + 1 - NN
P = WR(EN)
Q = WI(EN)
NA = EN - 1
C IF (Q) 710, 600, 800
IF (Q .GT. 0.0D0) GOTO 710
IF (Q .LT. 0.0D0) GOTO 800
C .......... REAL VECTOR ..........
600 M = EN
H(EN,EN) = 1.0D0
IF (NA .EQ. 0) GO TO 800
C .......... FOR I=EN-1 STEP -1 UNTIL 1 DO -- ..........
DO 700 II = 1, NA
I = EN - II
W = H(I,I) - P
R = 0.0D0
C
DO 610 J = M, EN
610 R = R + H(I,J) * H(J,EN)
C
IF (WI(I) .GE. 0.0D0) GO TO 630
ZZ = W
S = R
GO TO 700
630 M = I
IF (WI(I) .NE. 0.0D0) GO TO 640
T = W
IF (T .NE. 0.0D0) GO TO 635
TST1 = NORM
T = TST1
632 T = 0.01D0 * T
TST2 = NORM + T
IF (TST2 .GT. TST1) GO TO 632
635 H(I,EN) = -R / T
GO TO 680
C .......... SOLVE REAL EQUATIONS ..........
640 X = H(I,I+1)
Y = H(I+1,I)
Q = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I)
T = (X * S - ZZ * R) / Q
H(I,EN) = T
IF (DABS(X) .LE. DABS(ZZ)) GO TO 650
H(I+1,EN) = (-R - W * T) / X
GO TO 680
650 H(I+1,EN) = (-S - Y * T) / ZZ
C
C .......... OVERFLOW CONTROL ..........
680 T = DABS(H(I,EN))
IF (T .EQ. 0.0D0) GO TO 700
TST1 = T
TST2 = TST1 + 1.0D0/TST1
IF (TST2 .GT. TST1) GO TO 700
DO 690 J = I, EN
H(J,EN) = H(J,EN)/T
690 CONTINUE
C
700 CONTINUE
C .......... END REAL VECTOR ..........
GO TO 800
C .......... COMPLEX VECTOR ..........
710 M = NA
C .......... LAST VECTOR COMPONENT CHOSEN IMAGINARY SO THAT
C EIGENVECTOR MATRIX IS TRIANGULAR ..........
IF (DABS(H(EN,NA)) .LE. DABS(H(NA,EN))) GO TO 720
H(NA,NA) = Q / H(EN,NA)
H(NA,EN) = -(H(EN,EN) - P) / H(EN,NA)
GO TO 730
720 CALL CDIV(0.0D0,-H(NA,EN),H(NA,NA)-P,Q,H(NA,NA),H(NA,EN))
730 H(EN,NA) = 0.0D0
H(EN,EN) = 1.0D0
ENM2 = NA - 1
IF (ENM2 .EQ. 0) GO TO 800
C .......... FOR I=EN-2 STEP -1 UNTIL 1 DO -- ..........
DO 795 II = 1, ENM2
I = NA - II
W = H(I,I) - P
RA = 0.0D0
SA = 0.0D0
C
DO 760 J = M, EN
RA = RA + H(I,J) * H(J,NA)
SA = SA + H(I,J) * H(J,EN)
760 CONTINUE
C
IF (WI(I) .GE. 0.0D0) GO TO 770
ZZ = W
R = RA
S = SA
GO TO 795
770 M = I
IF (WI(I) .NE. 0.0D0) GO TO 780
CALL CDIV(-RA,-SA,W,Q,H(I,NA),H(I,EN))
GO TO 790
C .......... SOLVE COMPLEX EQUATIONS ..........
780 X = H(I,I+1)
Y = H(I+1,I)
VR = (WR(I) - P) * (WR(I) - P) + WI(I) * WI(I) - Q * Q
VI = (WR(I) - P) * 2.0D0 * Q
IF (VR .NE. 0.0D0 .OR. VI .NE. 0.0D0) GO TO 784
TST1 = NORM * (DABS(W) + DABS(Q) + DABS(X)
X + DABS(Y) + DABS(ZZ))
VR = TST1
783 VR = 0.01D0 * VR
TST2 = TST1 + VR
IF (TST2 .GT. TST1) GO TO 783
784 CALL CDIV(X*R-ZZ*RA+Q*SA,X*S-ZZ*SA-Q*RA,VR,VI,
X H(I,NA),H(I,EN))
IF (DABS(X) .LE. DABS(ZZ) + DABS(Q)) GO TO 785
H(I+1,NA) = (-RA - W * H(I,NA) + Q * H(I,EN)) / X
H(I+1,EN) = (-SA - W * H(I,EN) - Q * H(I,NA)) / X
GO TO 790
785 CALL CDIV(-R-Y*H(I,NA),-S-Y*H(I,EN),ZZ,Q,
X H(I+1,NA),H(I+1,EN))
C
C .......... OVERFLOW CONTROL ..........
790 T = DMAX1(DABS(H(I,NA)), DABS(H(I,EN)))
IF (T .EQ. 0.0D0) GO TO 795
TST1 = T
TST2 = TST1 + 1.0D0/TST1
IF (TST2 .GT. TST1) GO TO 795
DO 792 J = I, EN
H(J,NA) = H(J,NA)/T
H(J,EN) = H(J,EN)/T
792 CONTINUE
C
795 CONTINUE
C .......... END COMPLEX VECTOR ..........
800 CONTINUE
C .......... END BACK SUBSTITUTION.
C VECTORS OF ISOLATED ROOTS ..........
DO 840 I = 1, N
IF (I .GE. LOW .AND. I .LE. IGH) GO TO 840
C
DO 820 J = I, N
820 Z(I,J) = H(I,J)
C
840 CONTINUE
C .......... MULTIPLY BY TRANSFORMATION MATRIX TO GIVE
C VECTORS OF ORIGINAL FULL MATRIX.
C FOR J=N STEP -1 UNTIL LOW DO -- ..........
DO 880 JJ = LOW, N
J = N + LOW - JJ
M = MIN0(J,IGH)
C
DO 880 I = LOW, IGH
ZZ = 0.0D0
C
DO 860 K = LOW, M
860 ZZ = ZZ + Z(I,K) * H(K,J)
C
Z(I,J) = ZZ
880 CONTINUE
C
GO TO 1001
C .......... SET ERROR -- ALL EIGENVALUES HAVE NOT
C CONVERGED AFTER 30*N ITERATIONS ..........
1000 IERR = EN
1001 RETURN
END
SUBROUTINE HTRIBK(NM,N,AR,AI,TAU,M,ZR,ZI)
C
INTEGER I,J,K,L,M,N,NM
DOUBLE PRECISION AR(NM,N),AI(NM,N),TAU(2,N),ZR(NM,M),ZI(NM,M)
DOUBLE PRECISION H,S,SI
C
C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
C THE ALGOL PROCEDURE TRBAK1, NUM. MATH. 11, 181-195(1968)
C BY MARTIN, REINSCH, AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C THIS SUBROUTINE FORMS THE EIGENVECTORS OF A COMPLEX HERMITIAN
C MATRIX BY BACK TRANSFORMING THOSE OF THE CORRESPONDING
C REAL SYMMETRIC TRIDIAGONAL MATRIX DETERMINED BY HTRIDI.
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 AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
C FORMATIONS USED IN THE REDUCTION BY HTRIDI IN THEIR
C FULL LOWER TRIANGLES EXCEPT FOR THE DIAGONAL OF AR.
C
C TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
C
C M IS THE NUMBER OF EIGENVECTORS TO BE BACK TRANSFORMED.
C
C ZR CONTAINS THE EIGENVECTORS TO BE BACK TRANSFORMED
C IN ITS FIRST M COLUMNS.
C
C ON OUTPUT
C
C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE TRANSFORMED EIGENVECTORS
C IN THEIR FIRST M COLUMNS.
C
C NOTE THAT THE LAST COMPONENT OF EACH RETURNED VECTOR
C IS REAL AND THAT VECTOR EUCLIDEAN NORMS ARE PRESERVED.
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 (M .EQ. 0) GO TO 200
C .......... TRANSFORM THE EIGENVECTORS OF THE REAL SYMMETRIC
C TRIDIAGONAL MATRIX TO THOSE OF THE HERMITIAN
C TRIDIAGONAL MATRIX. ..........
DO 50 K = 1, N
C
DO 50 J = 1, M
ZI(K,J) = -ZR(K,J) * TAU(2,K)
ZR(K,J) = ZR(K,J) * TAU(1,K)
50 CONTINUE
C
IF (N .EQ. 1) GO TO 200
C .......... RECOVER AND APPLY THE HOUSEHOLDER MATRICES ..........
DO 140 I = 2, N
L = I - 1
H = AI(I,I)
IF (H .EQ. 0.0D0) GO TO 140
C
DO 130 J = 1, M
S = 0.0D0
SI = 0.0D0
C
DO 110 K = 1, L
S = S + AR(I,K) * ZR(K,J) - AI(I,K) * ZI(K,J)
SI = SI + AR(I,K) * ZI(K,J) + AI(I,K) * ZR(K,J)
110 CONTINUE
C .......... DOUBLE DIVISIONS AVOID POSSIBLE UNDERFLOW ..........
S = (S / H) / H
SI = (SI / H) / H
C
DO 120 K = 1, L
ZR(K,J) = ZR(K,J) - S * AR(I,K) - SI * AI(I,K)
ZI(K,J) = ZI(K,J) - SI * AR(I,K) + S * AI(I,K)
120 CONTINUE
C
130 CONTINUE
C
140 CONTINUE
C
200 RETURN
END
SUBROUTINE HTRIDI(NM,N,AR,AI,D,E,E2,TAU)
C
INTEGER I,J,K,L,N,II,NM,JP1
DOUBLE PRECISION AR(NM,N),AI(NM,N),D(N),E(N),E2(N),TAU(2,N)
DOUBLE PRECISION F,G,H,FI,GI,HH,SI,SCALE,PYTHAG
C
C THIS SUBROUTINE IS A TRANSLATION OF A COMPLEX ANALOGUE OF
C THE ALGOL PROCEDURE TRED1, NUM. MATH. 11, 181-195(1968)
C BY MARTIN, REINSCH, AND WILKINSON.
C HANDBOOK FOR AUTO. COMP., VOL.II-LINEAR ALGEBRA, 212-226(1971).
C
C THIS SUBROUTINE REDUCES A COMPLEX HERMITIAN MATRIX
C TO A REAL SYMMETRIC TRIDIAGONAL MATRIX USING
C UNITARY 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 AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX HERMITIAN INPUT MATRIX.
C ONLY THE LOWER TRIANGLE OF THE MATRIX NEED BE SUPPLIED.
C
C ON OUTPUT
C
C AR AND AI CONTAIN INFORMATION ABOUT THE UNITARY TRANS-
C FORMATIONS USED IN THE REDUCTION IN THEIR FULL LOWER
C TRIANGLES. THEIR STRICT UPPER TRIANGLES AND THE
C DIAGONAL OF AR ARE UNALTERED.
C
C D CONTAINS THE DIAGONAL ELEMENTS OF THE 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 TAU CONTAINS FURTHER INFORMATION ABOUT THE TRANSFORMATIONS.
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
TAU(1,N) = 1.0D0
TAU(2,N) = 0.0D0
C
DO 100 I = 1, N
100 D(I) = AR(I,I)
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 120 K = 1, L
120 SCALE = SCALE + DABS(AR(I,K)) + DABS(AI(I,K))
C
IF (SCALE .NE. 0.0D0) GO TO 140
TAU(1,L) = 1.0D0
TAU(2,L) = 0.0D0
130 E(I) = 0.0D0
E2(I) = 0.0D0
GO TO 290
C
140 DO 150 K = 1, L
AR(I,K) = AR(I,K) / SCALE
AI(I,K) = AI(I,K) / SCALE
H = H + AR(I,K) * AR(I,K) + AI(I,K) * AI(I,K)
150 CONTINUE
C
E2(I) = SCALE * SCALE * H
G = DSQRT(H)
E(I) = SCALE * G
F = PYTHAG(AR(I,L),AI(I,L))
C .......... FORM NEXT DIAGONAL ELEMENT OF MATRIX T ..........
IF (F .EQ. 0.0D0) GO TO 160
TAU(1,L) = (AI(I,L) * TAU(2,I) - AR(I,L) * TAU(1,I)) / F
SI = (AR(I,L) * TAU(2,I) + AI(I,L) * TAU(1,I)) / F
H = H + F * G
G = 1.0D0 + G / F
AR(I,L) = G * AR(I,L)
AI(I,L) = G * AI(I,L)
IF (L .EQ. 1) GO TO 270
GO TO 170
160 TAU(1,L) = -TAU(1,I)
SI = TAU(2,I)
AR(I,L) = G
170 F = 0.0D0
C
DO 240 J = 1, L
G = 0.0D0
GI = 0.0D0
C .......... FORM ELEMENT OF A*U ..........
DO 180 K = 1, J
G = G + AR(J,K) * AR(I,K) + AI(J,K) * AI(I,K)
GI = GI - AR(J,K) * AI(I,K) + AI(J,K) * AR(I,K)
180 CONTINUE
C
JP1 = J + 1
IF (L .LT. JP1) GO TO 220
C
DO 200 K = JP1, L
G = G + AR(K,J) * AR(I,K) - AI(K,J) * AI(I,K)
GI = GI - AR(K,J) * AI(I,K) - AI(K,J) * AR(I,K)
200 CONTINUE
C .......... FORM ELEMENT OF P ..........
220 E(J) = G / H
TAU(2,J) = GI / H
F = F + E(J) * AR(I,J) - TAU(2,J) * AI(I,J)
240 CONTINUE
C
HH = F / (H + H)
C .......... FORM REDUCED A ..........
DO 260 J = 1, L
F = AR(I,J)
G = E(J) - HH * F
E(J) = G
FI = -AI(I,J)
GI = TAU(2,J) - HH * FI
TAU(2,J) = -GI
C
DO 260 K = 1, J
AR(J,K) = AR(J,K) - F * E(K) - G * AR(I,K)
X + FI * TAU(2,K) + GI * AI(I,K)
AI(J,K) = AI(J,K) - F * TAU(2,K) - G * AI(I,K)
X - FI * E(K) - GI * AR(I,K)
260 CONTINUE
C
270 DO 280 K = 1, L
AR(I,K) = SCALE * AR(I,K)
AI(I,K) = SCALE * AI(I,K)
280 CONTINUE
C
TAU(2,L) = -SI
290 HH = D(I)
D(I) = AR(I,I)
AR(I,I) = HH
AI(I,I) = SCALE * DSQRT(H)
300 CONTINUE
C
RETURN
END
DOUBLE PRECISION FUNCTION PYTHAG(A,B)
DOUBLE PRECISION A,B
C
C FINDS DSQRT(A**2+B**2) WITHOUT OVERFLOW OR DESTRUCTIVE UNDERFLOW
C
DOUBLE PRECISION P,R,S,T,U
P = DMAX1(DABS(A),DABS(B))
IF (P .EQ. 0.0D0) GO TO 20
R = (DMIN1(DABS(A),DABS(B))/P)**2
10 CONTINUE
T = 4.0D0 + R
IF (T .EQ. 4.0D0) GO TO 20
S = R/T
U = 1.0D0 + 2.0D0*S
P = U*P
R = (S/U)**2 * R
GO TO 10
20 PYTHAG = P
RETURN
END
SUBROUTINE TQL1(N,D,E,IERR)
C
INTEGER I,J,L,M,N,II,L1,L2,MML,IERR
DOUBLE PRECISION D(N),E(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 TQL1,
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 OF A SYMMETRIC
C TRIDIAGONAL MATRIX BY THE 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 E CONTAINS THE SUBDIAGONAL ELEMENTS OF THE INPUT MATRIX
C IN ITS LAST N-1 POSITIONS. E(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 E 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 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 100 I = 2, N
100 E(I-1) = E(I)
C
F = 0.0D0
TST1 = 0.0D0
E(N) = 0.0D0
C
DO 290 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 210
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 140 I = L2, N
140 D(I) = D(I) - H
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))
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
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 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 100 I = 2, N
100 E(I-1) = E(I)
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))
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 140 I = L2, N
140 D(I) = D(I) - H
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 100 I = 2, N
100 E2(I-1) = E2(I)
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 140 I = L1, N
140 D(I) = D(I) - H
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 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
140 DO 150 K = 1, L
D(K) = D(K) / SCALE
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 ..........
DO 170 J = 1, L
170 E(J) = 0.0D0
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
DO 245 J = 1, L
E(J) = E(J) / H
F = F + E(J) * D(J)
245 CONTINUE
C
H = F / (H + H)
C .......... FORM Q ..........
DO 250 J = 1, L
250 E(J) = E(J) - H * D(J)
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
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 100 I = 1, N
C
DO 80 J = I, N
80 Z(J,I) = A(J,I)
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 -- ..........
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 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
140 DO 150 K = 1, L
D(K) = D(K) / SCALE
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 ..........
DO 170 J = 1, L
170 E(J) = 0.0D0
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
DO 245 J = 1, L
E(J) = E(J) / H
F = F + E(J) * D(J)
245 CONTINUE
C
HH = F / (H + H)
C .......... FORM Q ..........
DO 250 J = 1, L
250 E(J) = E(J) - HH * D(J)
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
DO 330 K = 1, L
330 D(K) = Z(K,I) / H
C
DO 360 J = 1, L
G = 0.0D0
C
DO 340 K = 1, L
340 G = G + Z(K,I) * Z(K,J)
C
DO 360 K = 1, L
Z(K,J) = Z(K,J) - G * D(K)
360 CONTINUE
C
380 DO 400 K = 1, L
400 Z(K,I) = 0.0D0
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
SUBROUTINE RG(NM,N,A,WR,WI,MATZ,Z,IV1,FV1,IERR)
C
INTEGER N,NM,IS1,IS2,IERR,MATZ
DOUBLE PRECISION A(NM,N),WR(N),WI(N),Z(NM,N),FV1(N)
INTEGER IV1(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 GENERAL 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 GENERAL 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 WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVALUES. COMPLEX CONJUGATE
C PAIRS OF EIGENVALUES APPEAR CONSECUTIVELY WITH THE
C EIGENVALUE HAVING THE POSITIVE IMAGINARY PART FIRST.
C
C Z CONTAINS THE REAL AND IMAGINARY PARTS OF THE EIGENVECTORS
C IF MATZ IS NOT ZERO. IF THE J-TH EIGENVALUE IS REAL, THE
C J-TH COLUMN OF Z CONTAINS ITS EIGENVECTOR. IF THE J-TH
C EIGENVALUE IS COMPLEX WITH POSITIVE IMAGINARY PART, THE
C J-TH AND (J+1)-TH COLUMNS OF Z CONTAIN THE REAL AND
C IMAGINARY PARTS OF ITS EIGENVECTOR. THE CONJUGATE OF THIS
C VECTOR IS THE EIGENVECTOR FOR THE CONJUGATE EIGENVALUE.
C
C IERR IS AN INTEGER OUTPUT VARIABLE SET EQUAL TO AN ERROR
C COMPLETION CODE DESCRIBED IN THE DOCUMENTATION FOR HQR
C AND HQR2. THE NORMAL COMPLETION CODE IS ZERO.
C
C IV1 AND FV1 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 CALL BALANC(NM,N,A,IS1,IS2,FV1)
CALL ELMHES(NM,N,IS1,IS2,A,IV1)
IF (MATZ .NE. 0) GO TO 20
C .......... FIND EIGENVALUES ONLY ..........
CALL HQR(NM,N,IS1,IS2,A,WR,WI,IERR)
GO TO 50
C .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
20 CALL ELTRAN(NM,N,IS1,IS2,A,IV1,Z)
CALL HQR2(NM,N,IS1,IS2,A,WR,WI,Z,IERR)
IF (IERR .NE. 0) GO TO 50
CALL BALBAK(NM,N,IS1,IS2,FV1,N,Z)
50 RETURN
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 CG(NM,N,AR,AI,WR,WI,MATZ,ZR,ZI,FV1,FV2,FV3,IERR)
C
INTEGER N,NM,IS1,IS2,IERR,MATZ
DOUBLE PRECISION AR(NM,N),AI(NM,N),WR(N),WI(N),ZR(NM,N),ZI(NM,N),
X FV1(N),FV2(N),FV3(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 COMPLEX GENERAL 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=(AR,AI).
C
C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX GENERAL 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 WR AND WI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE EIGENVALUES.
C
C ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF 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 COMQR
C AND COMQR2. THE NORMAL COMPLETION CODE IS ZERO.
C
C FV1, FV2, AND FV3 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 CALL CBAL(NM,N,AR,AI,IS1,IS2,FV1)
CALL CORTH(NM,N,IS1,IS2,AR,AI,FV2,FV3)
IF (MATZ .NE. 0) GO TO 20
C .......... FIND EIGENVALUES ONLY ..........
CALL COMQR(NM,N,IS1,IS2,AR,AI,WR,WI,IERR)
GO TO 50
C .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
20 CALL COMQR2(NM,N,IS1,IS2,FV2,FV3,AR,AI,WR,WI,ZR,ZI,IERR)
IF (IERR .NE. 0) GO TO 50
CALL CBABK2(NM,N,IS1,IS2,FV1,N,ZR,ZI)
50 RETURN
END
SUBROUTINE CH(NM,N,AR,AI,W,MATZ,ZR,ZI,FV1,FV2,FM1,IERR)
C
INTEGER I,J,N,NM,IERR,MATZ
DOUBLE PRECISION AR(NM,N),AI(NM,N),W(N),ZR(NM,N),ZI(NM,N),
X FV1(N),FV2(N),FM1(2,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 COMPLEX HERMITIAN 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=(AR,AI).
C
C AR AND AI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF THE COMPLEX HERMITIAN 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 ZR AND ZI CONTAIN THE REAL AND IMAGINARY PARTS,
C RESPECTIVELY, OF 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, FV2, AND FM1 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 CALL HTRIDI(NM,N,AR,AI,W,FV1,FV2,FM1)
IF (MATZ .NE. 0) GO TO 20
C .......... FIND EIGENVALUES ONLY ..........
CALL TQLRAT(N,W,FV2,IERR)
GO TO 50
C .......... FIND BOTH EIGENVALUES AND EIGENVECTORS ..........
20 DO 40 I = 1, N
C
DO 30 J = 1, N
ZR(J,I) = 0.0D0
30 CONTINUE
C
ZR(I,I) = 1.0D0
40 CONTINUE
C
CALL TQL2(NM,N,W,FV1,ZR,IERR)
IF (IERR .NE. 0) GO TO 50
CALL HTRIBK(NM,N,AR,AI,FM1,N,ZR,ZI)
50 RETURN
END