C ALGORITHM 598, COLLECTED ALGORITHMS FROM ACM.
C ALGORITHM APPEARED IN ACM-TRANS. MATH. SOFTWARE, VOL.9, NO. 2,
C JUN., 1983, P. 246-254.
COMPLEX A(50,50), B(50,50), C(50,50), S(50,50) MAN 10
COMPLEX WORK(17550) MAN 20
COMPLEX CMPLX MAN 30
REAL REAL MAN 40
REAL FXNORM, TOL MAN 50
INTEGER NM, N, I, ICASE, IERR, IGUESS, J, MAXITS, NW MAN 60
NM = 50 MAN 70
NW = 17550 MAN 80
C MAN 90
C CASE 1: DENNIS, TRAUB AND WEBER. MAN 100
C MAN 110
ICASE = 1 MAN 120
NOUT = 6 MAN 130
WRITE (NOUT,99999) ICASE MAN 140
N = 2 MAN 150
A(1,1) = CMPLX(1.0,0.0) MAN 160
A(1,2) = CMPLX(0.0,0.0) MAN 170
A(2,1) = CMPLX(0.0,0.0) MAN 180
A(2,2) = CMPLX(1.0,0.0) MAN 190
B(1,1) = CMPLX(7.0,0.0) MAN 200
B(1,2) = CMPLX(8.0,0.0) MAN 210
B(2,1) = CMPLX(8.0,0.0) MAN 220
B(2,2) = CMPLX(10.0,0.0) MAN 230
C(1,1) = CMPLX(9.0,0.0) MAN 240
C(1,2) = CMPLX(3.0,0.0) MAN 250
C(2,1) = CMPLX(4.0,0.0) MAN 260
C(2,2) = CMPLX(4.0,0.0) MAN 270
IGUESS = 0 MAN 280
TOL = 0.0 MAN 290
MAXITS = 0 MAN 300
CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 310
* IERR) MAN 320
IF (IERR.NE.0) GO TO 10 MAN 330
FXNORM = REAL(WORK(1)) MAN 340
WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 350
CALL PRINT(NM, N, S) MAN 360
GO TO 20 MAN 370
10 WRITE (NOUT,99998) ICASE, IERR MAN 380
C MAN 390
C CASE 2: ONE-BY-ONE TEST. MAN 400
C MAN 410
20 ICASE = 2 MAN 420
WRITE (NOUT,99999) ICASE MAN 430
N = 1 MAN 440
A(1,1) = CMPLX(1.0,0.0) MAN 450
B(1,1) = CMPLX(-3.0,0.0) MAN 460
C(1,1) = CMPLX(2.0,0.0) MAN 470
IGUESS = 0 MAN 480
TOL = 0.0 MAN 490
MAXITS = 0 MAN 500
CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 510
* IERR) MAN 520
IF (IERR.NE.0) GO TO 30 MAN 530
FXNORM = REAL(WORK(1)) MAN 540
WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 550
CALL PRINT(NM, N, S) MAN 560
GO TO 40 MAN 570
30 WRITE (NOUT,99998) ICASE, IERR MAN 580
C MAN 590
C CASE 3: COMPLEX SOLVENT. MAN 600
C MAN 610
40 ICASE = 3 MAN 620
WRITE (NOUT,99999) ICASE MAN 630
N = 3 MAN 640
A(1,1) = CMPLX(1.0,0.0) MAN 650
A(1,2) = CMPLX(2.0,0.0) MAN 660
A(1,3) = CMPLX(-1.0,0.0) MAN 670
A(2,1) = CMPLX(3.0,0.0) MAN 680
A(2,2) = CMPLX(4.0,0.0) MAN 690
A(2,3) = CMPLX(0.0,0.0) MAN 700
A(3,1) = CMPLX(-1.0,0.0) MAN 710
A(3,2) = CMPLX(2.0,0.0) MAN 720
A(3,3) = CMPLX(1.0,0.0) MAN 730
B(1,1) = CMPLX(-2.0,0.0) MAN 740
B(1,2) = CMPLX(0.0,1.0) MAN 750
B(1,3) = CMPLX(0.0,0.0) MAN 760
B(2,1) = CMPLX(0.0,1.0) MAN 770
B(2,2) = CMPLX(2.0,0.0) MAN 780
B(2,3) = CMPLX(0.0,0.0) MAN 790
B(3,1) = CMPLX(0.0,0.0) MAN 800
B(3,2) = CMPLX(0.0,0.0) MAN 810
B(3,3) = CMPLX(0.0,1.0) MAN 820
C(1,1) = CMPLX(2.0,0.0) MAN 830
C(1,2) = CMPLX(1.0,-7.0) MAN 840
C(1,3) = CMPLX(4.0,0.0) MAN 850
C(2,1) = CMPLX(1.0,-7.0) MAN 860
C(2,2) = CMPLX(-8.0,-17.0) MAN 870
C(2,3) = CMPLX(0.0,0.0) MAN 880
C(3,1) = CMPLX(0.0,2.0) MAN 890
C(3,2) = CMPLX(2.0,-2.0) MAN 900
C(3,3) = CMPLX(-4.0,-2.0) MAN 910
IGUESS = 0 MAN 920
TOL = 0.001 MAN 930
MAXITS = 0 MAN 940
CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 950
* IERR) MAN 960
IF (IERR.NE.0) GO TO 50 MAN 970
FXNORM = REAL(WORK(1)) MAN 980
WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 990
CALL PRINT(NM, N, S) MAN 1000
GO TO 60 MAN 1010
50 WRITE (NOUT,99998) ICASE, IERR MAN 1020
C MAN 1030
C CASE 4: NO SOLVENT. MAN 1040
C MAN 1050
60 ICASE = 4 MAN 1060
WRITE (NOUT,99999) ICASE MAN 1070
N = 2 MAN 1080
A(1,1) = CMPLX(1.0,0.0) MAN 1090
A(1,2) = CMPLX(0.0,0.0) MAN 1100
A(2,1) = CMPLX(0.0,0.0) MAN 1110
A(2,2) = CMPLX(0.0,0.0) MAN 1120
B(1,1) = CMPLX(0.0,0.0) MAN 1130
B(1,2) = CMPLX(0.0,0.0) MAN 1140
B(2,1) = CMPLX(0.0,0.0) MAN 1150
B(2,2) = CMPLX(0.0,0.0) MAN 1160
C(1,1) = CMPLX(1.0,0.0) MAN 1170
C(1,2) = CMPLX(0.0,0.0) MAN 1180
C(2,1) = CMPLX(0.0,0.0) MAN 1190
C(2,2) = CMPLX(1.0,0.0) MAN 1200
IGUESS = 0 MAN 1210
TOL = 0.0 MAN 1220
MAXITS = 0 MAN 1230
CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 1240
* IERR) MAN 1250
IF (IERR.NE.0) GO TO 70 MAN 1260
FXNORM = REAL(WORK(1)) MAN 1270
WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 1280
CALL PRINT(NM, N, S) MAN 1290
GO TO 80 MAN 1300
70 WRITE (NOUT,99998) ICASE, IERR MAN 1310
C MAN 1320
C CASE 5: EXAMPLE CITED IN PAPER. MAN 1330
C MAN 1340
80 ICASE = 5 MAN 1350
WRITE (NOUT,99999) ICASE MAN 1360
N = 2 MAN 1370
A(1,1) = CMPLX(1.0,0.0) MAN 1380
A(1,2) = CMPLX(0.0,0.0) MAN 1390
A(2,1) = CMPLX(0.0,0.0) MAN 1400
A(2,2) = CMPLX(1.0,0.0) MAN 1410
B(1,1) = CMPLX(-1.0,0.0) MAN 1420
B(1,2) = CMPLX(-1.0,0.0) MAN 1430
B(2,1) = CMPLX(1.0,0.0) MAN 1440
B(2,2) = CMPLX(-1.0,0.0) MAN 1450
C(1,1) = CMPLX(0.0,0.0) MAN 1460
C(1,2) = CMPLX(1.0,0.0) MAN 1470
C(2,1) = CMPLX(-1.0,0.0) MAN 1480
C(2,2) = CMPLX(0.0,0.0) MAN 1490
IGUESS = 0 MAN 1500
TOL = 0.0 MAN 1510
MAXITS = 0 MAN 1520
CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 1530
* IERR) MAN 1540
IF (IERR.NE.0) GO TO 90 MAN 1550
FXNORM = REAL(WORK(1)) MAN 1560
WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 1570
CALL PRINT(NM, N, S) MAN 1580
GO TO 100 MAN 1590
90 WRITE (NOUT,99998) ICASE, IERR MAN 1600
C MAN 1610
C CASE 6: COMPLEX GUESS FOR CASE 5. MAN 1620
C MAN 1630
100 ICASE = 6 MAN 1640
WRITE (NOUT,99999) ICASE MAN 1650
N = 2 MAN 1660
A(1,1) = CMPLX(1.0,0.0) MAN 1670
A(1,2) = CMPLX(0.0,0.0) MAN 1680
A(2,1) = CMPLX(0.0,0.0) MAN 1690
A(2,2) = CMPLX(1.0,0.0) MAN 1700
B(1,1) = CMPLX(-1.0,0.0) MAN 1710
B(1,2) = CMPLX(-1.0,0.0) MAN 1720
B(2,1) = CMPLX(1.0,0.0) MAN 1730
B(2,2) = CMPLX(-1.0,0.0) MAN 1740
C(1,1) = CMPLX(0.0,0.0) MAN 1750
C(1,2) = CMPLX(1.0,0.0) MAN 1760
C(2,1) = CMPLX(-1.0,0.0) MAN 1770
C(2,2) = CMPLX(0.0,0.0) MAN 1780
IGUESS = 1 MAN 1790
DO 120 I=1,N MAN 1800
DO 110 J=1,N MAN 1810
S(I,J) = CMPLX(0.0,0.0) MAN 1820
110 CONTINUE MAN 1830
S(I,I) = CMPLX(0.0,1.0) MAN 1840
120 CONTINUE MAN 1850
TOL = 0.0 MAN 1860
MAXITS = 0 MAN 1870
CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 1880
* IERR) MAN 1890
IF (IERR.NE.0) GO TO 130 MAN 1900
FXNORM = REAL(WORK(1)) MAN 1910
WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 1920
CALL PRINT(NM, N, S) MAN 1930
GO TO 140 MAN 1940
130 WRITE (NOUT,99998) ICASE, IERR MAN 1950
C MAN 1960
C CASE 7: BAD GUESS FOR CASE 5. MAN 1970
C MAN 1980
140 ICASE = 7 MAN 1990
WRITE (NOUT,99999) ICASE MAN 2000
N = 2 MAN 2010
A(1,1) = CMPLX(1.0,0.0) MAN 2020
A(1,2) = CMPLX(0.0,0.0) MAN 2030
A(2,1) = CMPLX(0.0,0.0) MAN 2040
A(2,2) = CMPLX(1.0,0.0) MAN 2050
B(1,1) = CMPLX(-1.0,0.0) MAN 2060
B(1,2) = CMPLX(-1.0,0.0) MAN 2070
B(2,1) = CMPLX(1.0,0.0) MAN 2080
B(2,2) = CMPLX(-1.0,0.0) MAN 2090
C(1,1) = CMPLX(0.0,0.0) MAN 2100
C(1,2) = CMPLX(1.0,0.0) MAN 2110
C(2,1) = CMPLX(-1.0,0.0) MAN 2120
C(2,2) = CMPLX(0.0,0.0) MAN 2130
IGUESS = 1 MAN 2140
DO 160 I=1,N MAN 2150
DO 150 J=1,N MAN 2160
S(I,J) = CMPLX(500.0,0.0) MAN 2170
150 CONTINUE MAN 2180
160 CONTINUE MAN 2190
TOL = 0.0 MAN 2200
MAXITS = 90 MAN 2210
CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 2220
* IERR) MAN 2230
IF (IERR.NE.0) GO TO 170 MAN 2240
FXNORM = REAL(WORK(1)) MAN 2250
WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 2260
CALL PRINT(NM, N, S) MAN 2270
GO TO 180 MAN 2280
170 WRITE (NOUT,99998) ICASE, IERR MAN 2290
C MAN 2300
C CASE 8: SINGULAR A,C,SOLVENT AND GUESS. MAN 2310
C MAN 2320
180 ICASE = 8 MAN 2330
WRITE (NOUT,99999) ICASE MAN 2340
N = 2 MAN 2350
A(1,1) = CMPLX(1.0E-8,0.0) MAN 2360
A(1,2) = CMPLX(0.0,0.0) MAN 2370
A(2,1) = CMPLX(0.0,0.0) MAN 2380
A(2,2) = CMPLX(1.0,0.0) MAN 2390
B(1,1) = CMPLX(1.0,0.0) MAN 2400
B(1,2) = CMPLX(0.0,0.0) MAN 2410
B(2,1) = CMPLX(0.0,0.0) MAN 2420
B(2,2) = CMPLX(1.0,0.0) MAN 2430
C(1,1) = CMPLX(0.0,0.0) MAN 2440
C(1,2) = CMPLX(0.0,0.0) MAN 2450
C(2,1) = CMPLX(0.0,0.0) MAN 2460
C(2,2) = CMPLX(-2.0,0.0) MAN 2470
IGUESS = 1 MAN 2480
S(1,1) = CMPLX(1.0,0.0) MAN 2490
S(1,2) = CMPLX(1.0,0.0) MAN 2500
S(2,1) = CMPLX(1.0,0.0) MAN 2510
S(2,2) = CMPLX(1.0,0.0) MAN 2520
TOL = 0.0 MAN 2530
MAXITS = 0 MAN 2540
CALL SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, MAXITS, MAN 2550
* IERR) MAN 2560
IF (IERR.NE.0) GO TO 190 MAN 2570
FXNORM = REAL(WORK(1)) MAN 2580
WRITE (NOUT,99997) IGUESS, IGUESS, FXNORM MAN 2590
CALL PRINT(NM, N, S) MAN 2600
GO TO 200 MAN 2610
190 WRITE (NOUT,99998) ICASE, IERR MAN 2620
C MAN 2630
200 CONTINUE MAN 2640
99999 FORMAT (//6H CASE , I2, 18H FOLLOWS..........) MAN 2650
99998 FORMAT (5H CASE, I2, 8H IERR IS, I3, 1H.) MAN 2660
99997 FORMAT (15H CONVERGENCE IN, I3, 23H ITERATIONS. NORM(F(X(, I2, MAN 2670
* 6H))) = , E15.6) MAN 2680
STOP MAN 2690
END MAN 2700
SUBROUTINE PRINT(NM, N, A) PRI 10
C
C
C SUBROUTINE PRINT PRINTS OUT AN INPUT MATRIX.
C
C
C ON ENTRY,
C
C NM IS THE LEADING DIMENSION OF A IN THE MAIN PROGRAM.
C
C N IS THE ORDER OF THE MATRIX A.
C
C A IS THE MATRIX TO BE PRINTED.
C
C
C ON RETURN,
C
C NM, N AND A ARE UNALTERED.
C
C
INTEGER I, J, JJ, JP1, N, NM
COMPLEX A(NM,N)
DO 20 J=1,N,2
JP1 = J + 1
IF (JP1.GT.N) JP1 = N
NOUT = 6
WRITE (NOUT,99999) J, JP1
DO 10 I=1,N
WRITE (NOUT,99998) (A(I,JJ),JJ=J,JP1)
10 CONTINUE
IF (JP1.EQ.N) RETURN
20 CONTINUE
99999 FORMAT (//8H COLUMNS, I3, 8H THROUGH, I3, 11H ..........)
99998 FORMAT (/3(E15.6, 2H +, E15.6, 3H *I, 2X))
RETURN
END
SUBROUTINE SQUINT(NM, N, A, B, C, IGUESS, S, WORK, NW, TOL, SQU 10
* MAXITS, IERR)
C
C
C SUBROUTINE SQUINT BREAKS DOWN THE WORK ARRAY 'WORK' INTO
C SMALLER PIECES. THE ACTUAL SOLUTION TO AX**2 + BX + C = 0 IS
C DONE IN SUBROUTINE SQUIN2. THIS SUBROUTINE MERELY RELIEVES
C THE USER FROM A LONG CALLING SEQUENCE.
C
C
C ON ENTRY,
C
C NM IS THE LEADING DIMENSION OF ALL THE MATRICES
C IN THE CALLING PROGRAM.
C
C N IS THE ORDER OF THE MATRICES A, B AND C.
C
C A IS THE MATRIX COEFFICIENT OF X**2.
C
C B IS THE MATRIX COEFFICIENT OF X.
C
C C IS THE CONSTANT MATRIX.
C
C IGUESS IS AN INTEGER. IF IGUESS.NE.0, THE USER SUPPLIES AN
C INITIAL GUESS AT A SOLVENT. THIS GUESS IS STORED IN
C ARRAY S. IF IGUESS.EQ.0, THE SUBROUTINE PROVIDES ITS
C OWN INITIAL GUESS.
C
C S CONTAINS THE USERS INITIAL GUESS AT A SOLVENT, IF IGUESS
C HAS BEEN SET TO A NONZERO QUANTITY. OTHERWISE THE INPUT
C CONTENTS IN S ARE IGNORED.
C
C WORK IS A WORK VECTOR. IT MUST BE DIMENSIONED AT LEAST
C (7N**2 + N), WHERE N IS THE ORDER OF A, B, C AND S.
C
C NW IS THE DIMENSION OF THE ARRAY WORK IN THE CALLING
C PROGRAM.
C
C TOL IS A USER-SUPPLIED ACCURACY TOLERANCE. SETTING TOL = 0.0
C CAUSES ITERATION TO PROCEED UNTIL FULL MACHINE PRECISION
C IS ATTAINED. OTHERWISE, EXECUTION TERMINATES WHEN
C NORM(AS**2+BS+C).LT.TOL.
C
C MAXITS IS AN INTEGER. IF MAXITS.NE.0, THE USER SPECIFIES THE
C MOST INTERATIONS THE ALGORITHM IS TO TAKE. IF MAXITS.
C .LE.0, IT IS RESET TO 30.
C
C
C ON RETURN,
C
C A,B,C ARE DESTROYED.
C
C IGUESS CONTAINS THE NUMBER OF ITERATIONS PERFORMED TO
C COMPUTE S.
C
C S CONTAINS THE RIGHT SOLVENT.
C
C WORK(1) IS A COMPLEX NUMBER WITH REAL PART EQUAL TO THE NORM
C OF AS**2+BS+C.
C
C IERR IS AN INTEGER ERROR RETURN.
C
C IERR = 0 FOR A NORMAL RETURN.
C
C IERR = 1 INDICATES FAILURE OF SQUINT TO CONVERGE TO
C A SOLVENT IN THE MAXIMUM NUMBER OF ITERATIONS.
C
C IERR = 2 INDICATES FAILURE IN THE UPPER REDUCTION IN
C CQZIT.
C
C IERR = 3 INDICATES FAILURE IN THE LOWER REDUCTION IN
C CQZIT.
C
C IERR = 10 + N INDICATES AN ERROR RETURN FROM TRISLV
C ON ITERATION N, DESIGNATING INCONSISTENCY OF
C THE TRIANGULAR SYSTEM.
C
C IERR = 999 INDICATES IMPROPER DIMENSIONING. THE
C CONDITIONS NM.GE.N.GT.0 AND
C NW.GE.(7*N*N + N) MUST HOLD.
C
INTEGER NM, N, IGUESS, NW, MAXITS, IERR, I1, I2, I3, I4, I5, I6,
* I7
COMPLEX WORK(NW)
COMPLEX A(NM,N), B(NM,N), C(NM,N), S(NM,N)
REAL TOL
I1 = N*N + 1
I2 = N*N + I1
I3 = N*N + I2
I4 = N*N + I3
I5 = N*N + I4
I6 = N*N + I5
I7 = N*N + I6
CALL SQUIN2(NM, N, A, B, C, IGUESS, S, WORK(1), WORK(I1),
* WORK(I2), WORK(I3), WORK(I4), WORK(I5), WORK(I6), WORK(I7), NW,
* TOL, MAXITS, IERR)
RETURN
END
C SQU 10
SUBROUTINE SQUIN2(NM, N, A, B, C, IGUESS, S, L, U, V, Z, R, XOLD, SQU 20
* EYE, TEMP, NW, TOL, MAXITS, IERR)
C
C SUBROUTINE SQUIN2 FINDS A RIGHT SOLVENT OF THE MATRIX
C EQUATION AX**2 + BX + C = 0.
C
C ON ENTRY,
C
C NM IS THE LEADING DIMENSION OF ALL THE MATRICES IN
C THE CALLING PROGRAM.
C
C N IS THE ORDER OF THE MATRICES A, B AND C.
C
C A IS THE MATRIX COEFFICIENT OF X**2.
C
C B IS THE MATRIX COEFFICIENT OF X.
C
C C IS THE CONSTANT MATRIX.
C
C IGUESS IS AN INTEGER SET IN THE CALL TO SQUINT.
C
C S IS A MATRIX SET IN THE CALL TO SQUINT.
C
C NW, TOL, AND MAXITS ARE INTEGER AND REAL PARAMETERS SET IN
C THE CALL TO SQUINT.
C
C
C THE FOLLOWING ARE INTERNAL VARIABLES:
C
C
C L IS A MATRIX CONTAINING THE ITERATE X(I) FOR
C REDUCTION TO LOWER TRIANGULAR FORM.
C
C U IS A MATRIX CONTAINING AX(I) + B FOR REDUCTION
C TO UPPER TRIANGULAR FORM.
C
C V IS A MATRIX CONTAINING A FOR REDUCTION TO UPPER
C TRIANGULAR FORM.
C
C Z,R ARE MATRICES CONTAINING THE HISTORY OF THE
C TRANSFORMATIONS IN THE REDUCTIONS.
C
C XOLD IS A MATRIX HOLDING THE CURRENT ITERATE X(I).
C
C EYE CONTAINS AN IDENTITY MATRIX FOR THE LOWER REDUCTION
C STEP.
C
C TEMP IS A WORK VECTOR.
C
C
C ON RETURN,
C
C A, B, C, IGUESS, S AND IERR HAVE THE SAME PROPERTIES AS
C DESCRIBED IN THE RETURN FROM SUBROUTINE SQUINT.
C
C L(1,1) IS A COMPLEX NUMBER WITH REAL PART EQUAL TO THE NORM
C OF AS**2+BS+C.
C
C
INTEGER NM, N, IGUESS, NW, MAXITS, IERR
INTEGER ITS, MATS, I, J, K
COMPLEX A(NM,N), B(NM,N), C(NM,N), S(NM,N), L(NM,N), U(NM,N)
COMPLEX V(NM,N), Z(NM,N), R(NM,N), XOLD(NM,N), EYE(NM,N), TEMP(N)
COMPLEX CMPLX, CONJG
REAL ANORM, ANI, BNORM, BNI, CNORM, CNI, XNORM, XNI
REAL FXNORM, FXNI, GNORM, TNORM, T, TOL
REAL SQRT, CABS, FLOAT
K = 7*N*N + N
IF (NW.LT.K) GO TO 460
IF (NM.LT.N) GO TO 460
IF (N.LE.0) GO TO 460
IF (MAXITS.LE.0) MAXITS = 30
C ********** INITIALIZE ARRAYS **********
DO 20 I=1,N
DO 10 J=1,N
L(I,J) = CMPLX(0.0,0.0)
U(I,J) = CMPLX(0.0,0.0)
V(I,J) = CMPLX(0.0,0.0)
Z(I,J) = CMPLX(0.0,0.0)
XOLD(I,J) = CMPLX(0.0,0.0)
EYE(I,J) = CMPLX(0.0,0.0)
10 CONTINUE
TEMP(I) = CMPLX(0.0,0.0)
20 CONTINUE
C ********** SET INITIAL GUESS(ES) **********
ANORM = 0.0
BNORM = 0.0
CNORM = 0.0
DO 40 I=1,N
ANI = 0.0
BNI = 0.0
CNI = 0.0
DO 30 J=1,N
ANI = ANI + CABS(A(I,J))
BNI = BNI + CABS(B(I,J))
CNI = CNI + CABS(C(I,J))
30 CONTINUE
IF (ANI.GT.ANORM) ANORM = ANI
IF (BNI.GT.BNORM) BNORM = BNI
IF (CNI.GT.CNORM) CNORM = CNI
40 CONTINUE
GNORM = (BNORM+SQRT(BNORM**2+4.0*ANORM*CNORM))/(2.0*ANORM)
IF (IGUESS.EQ.0) GO TO 70
DO 60 I=1,N
DO 50 J=1,N
XOLD(I,J) = S(I,J)
50 CONTINUE
60 CONTINUE
GO TO 100
C
70 DO 90 I=1,N
DO 80 J=1,N
XOLD(I,J) = CMPLX(0.0,0.0)
80 CONTINUE
XOLD(I,I) = CMPLX(GNORM,0.0)
90 CONTINUE
C
100 DO 360 ITS=1,MAXITS
IF (ITS.NE.31) GO TO 130
DO 120 I=1,N
DO 110 J=1,N
XOLD(I,J) = CMPLX(0.0,0.0)
110 CONTINUE
XOLD(I,I) = CMPLX(0.0,GNORM)
120 CONTINUE
C
130 IF (ITS.NE.61) GO TO 160
DO 150 I=1,N
DO 140 J=1,N
XOLD(I,J) = C(I,J)
140 CONTINUE
150 CONTINUE
C ********** SET UP U AND RIGHT HAND SIDE **********
160 CALL MULT(NM, N, A, XOLD, U)
DO 180 I=1,N
DO 170 J=1,N
U(I,J) = U(I,J) + B(I,J)
170 CONTINUE
180 CONTINUE
CALL MULT(NM, N, U, XOLD, S)
DO 200 I=1,N
DO 190 J=1,N
S(I,J) = S(I,J) + C(I,J)
190 CONTINUE
200 CONTINUE
C ********** CHECK FOR CONVERGENCE **********
XNORM = 0.0
FXNORM = 0.0
DO 220 I=1,N
XNI = 0.0
FXNI = 0.0
DO 210 J=1,N
XNI = XNI + CABS(XOLD(I,J))
FXNI = FXNI + CABS(S(I,J))
210 CONTINUE
IF (XNI.GT.XNORM) XNORM = XNI
IF (FXNI.GT.FXNORM) FXNORM = FXNI
220 CONTINUE
IF (TOL.LE.0.0) GO TO 230
IF (FXNORM.LT.TOL) GO TO 370
230 TNORM = 8.0*FLOAT(N)*ANORM*XNORM**2 + 5.0*FLOAT(N)*BNORM*XNORM
* + CNORM
T = 1.0 + FXNORM/TNORM
IF (T.EQ.1.0) GO TO 370
IF (ITS.GE.MAXITS) GO TO 400
C ********** UPPER TRIANGULARIZATION **********
IF (ITS.NE.1) GO TO 240
MATS = 1
CALL CQZHES(NM, N, U, A, MATS, Z, S, B, C)
240 DO 260 I=1,N
DO 250 J=1,N
V(I,J) = A(I,J)
L(I,J) = CONJG(XOLD(J,I))
EYE(I,J) = CMPLX(0.0,0.0)
250 CONTINUE
EYE(I,I) = CMPLX(1.0,0.0)
260 CONTINUE
IF (ITS.EQ.1) GO TO 270
MATS = 2
CALL CQZHES(NM, N, U, V, MATS, Z, S, B, C)
270 CALL CQZIT(NM, N, U, V, 0.0, MATS, Z, S, IERR)
IF (IERR.NE.0) GO TO 430
C ********** LOWER TRIANGULARIZATION **********
MATS = 3
CALL CQZHES(NM, N, L, EYE, MATS, R, S, B, C)
CALL CQZIT(NM, N, L, EYE, 0.0, MATS, R, S, IERR)
IF (IERR.NE.0) GO TO 440
CALL CTRANS(NM, N, L)
C ********** UPDATE S WITH R **********
DO 310 I=1,N
DO 280 J=1,N
TEMP(J) = S(I,J)
S(I,J) = CMPLX(0.0,0.0)
280 CONTINUE
DO 300 J=1,N
DO 290 K=1,N
S(I,J) = S(I,J) + TEMP(K)*R(K,J)
290 CONTINUE
300 CONTINUE
310 CONTINUE
DO 330 J=1,N
DO 320 I=1,N
L(I,J) = L(I,J)*EYE(J,J)
320 CONTINUE
EYE(J,J) = CMPLX(1.0,0.0)
330 CONTINUE
C ********** BACKSOLVE THE TRANSFORMED SYSTEM **********
CALL TRISLV(NM, N, U, V, L, S, TEMP, IERR)
IF (IERR.NE.0) GO TO 450
C ********** TRANSLATE BACK TO THE SOLUTION **********
CALL MULT(NM, N, Z, S, L)
CALL CTRANS(NM, N, R)
CALL MULT(NM, N, L, R, S)
DO 350 I=1,N
DO 340 J=1,N
XOLD(I,J) = XOLD(I,J) - S(I,J)
340 CONTINUE
350 CONTINUE
360 CONTINUE
C ********** CONVERGENCE **********
370 IGUESS = ITS - 1
L(1,1) = CMPLX(FXNORM,0.0)
DO 390 I=1,N
DO 380 J=1,N
S(I,J) = XOLD(I,J)
380 CONTINUE
390 CONTINUE
IERR = 0
RETURN
C ********** ERROR RETURNS **********
400 IERR = 1
L(1,1) = CMPLX(FXNORM,0.0)
DO 420 I=1,N
DO 410 J=1,N
S(I,J) = XOLD(I,J)
410 CONTINUE
420 CONTINUE
RETURN
430 IERR = 2
RETURN
440 IERR = 3
RETURN
450 IERR = ITS + 10
RETURN
460 IERR = 999
RETURN
END
SUBROUTINE CTRANS(NM, N, A) CTR 10
C
C
C SUBROUTINE CTRANS FINDS THE COMPLEX CONJUGATE OF AN INPUT
C MATRIX.
C
C
C ON ENTRY,
C
C NM IS THE LEADING DIMENSION OF MATRIX A IN THE MAIN PROGRAM.
C
C N IS THE ORDER OF MATRIX A.
C
C A IS THE INPUT MATRIX.
C
C
C ON RETURN,
C
C A CONTAINS ITS CONJUGATE TRANSPOSE.
C
C
INTEGER I, J, N, NM
COMPLEX A(NM,N), TEMP
COMPLEX CONJG
DO 20 I=1,N
DO 10 J=I,N
TEMP = A(I,J)
A(I,J) = CONJG(A(J,I))
A(J,I) = CONJG(TEMP)
10 CONTINUE
20 CONTINUE
RETURN
END
SUBROUTINE MULT(NM, N, A, B, C) MUL 10
C
C
C SUBROUTINE MULT MULTIPLIES TWO MATRICES TOGETHER.
C
C
C ON ENTRY,
C
C NM IS THE LEADING DIMENSION OF MATRICES A, B AND C IN THE
C MAIN PROGRAM.
C
C N IS THE ORDER OF THE MATRICES.
C
C A,B ARE THE MATRICES TO BE MULTIPLIED.
C
C C IS A WORK MATRIX WHICH IS INTERNALLY INITIALIZED TO ZERO.
C
C
C ON RETURN,
C
C A,B ARE UNALTERED.
C
C C CONTAINS THE MATRIX PRODUCT A*B.
C
C
INTEGER I, J, K, N, NM
COMPLEX A(NM,N), B(NM,N), C(NM,N)
COMPLEX CMPLX
DO 30 I=1,N
DO 20 J=1,N
C(I,J) = CMPLX(0.0,0.0)
DO 10 K=1,N
C(I,J) = C(I,J) + A(I,K)*B(K,J)
10 CONTINUE
20 CONTINUE
30 CONTINUE
RETURN
END
SUBROUTINE CQZHES(NM, N, A, B, MATS, Z, F, G, H) CQZ 10
C
C
C SUBROUTINE CQZHES IS A MODIFICATION OF THE EISPACK SUBROUTINE
C QZHES. ALL OPERATIONS ARE PERFORMED IN COMPLEX ARITHMETIC,
C AND THE LEFT TRANSFORMATIONS MAY ALSO BE APPLIED TO AUXILIARY
C MATRICES F, G AND H.
C
C
C ON ENTRY,
C
C NM IS THE LEADING DIMENSION OF THE MATRICES A AND B IN
C THE MAIN PROGRAM.
C
C N IS THE ORDER OF THE MATRICES A AND B.
C
C A CONTAINS THE MATRIX TO BE REDUCED TO UPPER HESSENBERG
C FORM.
C
C B CONTAINS THE MATRIX TO BE REDUCED TO UPPER TRAINGULAR
C FORM.
C
C MATS IS AN INTEGER INPUT VARIABLE.
C
C IF MATS = 0, THE ACCUMULATION OF THE TRANSFORMATIONS
C IS NOT DESIRED.
C
C IF MATS = ANY OTHER NUMBER BUT 0, THE TRANSFORMATIONS
C ARE ACCUMULATED.
C
C IF MATS = 1, THE AUXILIARY MATRICES G AND H ARE UPDATED
C WITH THE UNITARY MATRIX Q.
C
C IF MATS = 2, MATRIX B IS ASSUMED UPPER TRIANGULAR.
C
C IF MATS = 3, THE AUXILIARY MATRIX F IS NOT UPDATED
C WITH THE UNITARY MATRIX Q.
C
C F, G AND H ARE AUXILIARY MATRICES.
C
C
C ON RETURN,
C
C A IS UPPER HESSENBERG.
C
C B IS UPPER TRIANGULAR.
C
C Z CONTAINS THE HISTORY OF THE TRANSFORMATIONS, IF DESIRED.
C
C F, G AND H ARE UPDATED, IF DESIRED.
C
C
INTEGER I, J, K, L, N, LB, L1, NM, NK1, NM1, NM2
COMPLEX A(NM,N), B(NM,N), Z(NM,N)
COMPLEX RR, T, U1, U2, V1, V2, RHO
COMPLEX F(NM,N), TF, G(NM,N), TG, H(NM,N), TH
COMPLEX CMPLX, CONJG
REAL R, S
REAL SQRT, CABS
INTEGER MATS
IF (MATS.EQ.0) GO TO 30
DO 20 I=1,N
DO 10 J=1,N
Z(I,J) = CMPLX(0.0,0.0)
10 CONTINUE
Z(I,I) = CMPLX(1.0,0.0)
20 CONTINUE
C ********** REDUCE B TO UPPER TRIANGULAR FORM **********
30 IF (N.LE.1) GO TO 260
NM1 = N - 1
IF (MATS.EQ.2) GO TO 140
DO 130 L=1,NM1
L1 = L + 1
S = 0.0
DO 40 I=L1,N
S = S + CABS(B(I,L))
40 CONTINUE
IF (S.EQ.0.0) GO TO 130
S = S + CABS(B(L,L))
R = 0.0
DO 50 I=L,N
B(I,L) = B(I,L)/CMPLX(S,0.0)
R = R + CABS(B(I,L))**2
50 CONTINUE
R = SQRT(R)
RR = CMPLX(R,0.0)
IF (CABS(B(L,L)).NE.0.0) RR = (B(L,L)/CABS(B(L,L)))*RR
B(L,L) = B(L,L) + RR
RHO = CONJG(RR)*B(L,L)
DO 80 J=L1,N
T = CMPLX(0.0,0.0)
DO 60 I=L,N
T = T + CONJG(B(I,L))*B(I,J)
60 CONTINUE
T = -T/RHO
DO 70 I=L,N
B(I,J) = B(I,J) + T*B(I,L)
70 CONTINUE
80 CONTINUE
DO 110 J=1,N
T = CMPLX(0.0,0.0)
TF = CMPLX(0.0,0.0)
TG = CMPLX(0.0,0.0)
TH = CMPLX(0.0,0.0)
DO 90 I=L,N
T = T + CONJG(B(I,L))*A(I,J)
IF (MATS.EQ.3) GO TO 90
TF = TF + CONJG(B(I,L))*F(I,J)
IF (MATS.NE.1) GO TO 90
TG = TG + CONJG(B(I,L))*G(I,J)
TH = TH + CONJG(B(I,L))*H(I,J)
90 CONTINUE
T = -T/RHO
TF = -TF/RHO
TG = -TG/RHO
TH = -TH/RHO
DO 100 I=L,N
A(I,J) = A(I,J) + T*B(I,L)
IF (MATS.EQ.3) GO TO 100
F(I,J) = F(I,J) + TF*B(I,L)
IF (MATS.NE.1) GO TO 100
G(I,J) = G(I,J) + TG*B(I,L)
H(I,J) = H(I,J) + TH*B(I,L)
100 CONTINUE
110 CONTINUE
B(L,L) = -CMPLX(S,0.0)*RR
DO 120 I=L1,N
B(I,L) = CMPLX(0.0,0.0)
120 CONTINUE
130 CONTINUE
C ********** REDUCE A TO UPPER HESSENBERG FORM, WHILE
C KEEPING B TRIANGULAR **********
140 IF (N.EQ.2) GO TO 260
NM2 = N - 2
DO 250 K=1,NM2
NK1 = NM1 - K
DO 240 LB=1,NK1
L = N - LB
L1 = L + 1
C ********** ZERO A(L+1,K) **********
S = CABS(A(L,K)) + CABS(A(L1,K))
IF (S.EQ.0.0) GO TO 240
U1 = A(L,K)/CMPLX(S,0.0)
U2 = A(L1,K)/CMPLX(S,0.0)
R = SQRT(CABS(U1)**2+CABS(U2)**2)
RR = CMPLX(R,0.0)
IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR
V1 = -(U1+RR)/RR
V2 = -U2/RR
U2 = V2/V1
DO 150 J=K,N
T = A(L,J) + CONJG(U2)*A(L1,J)
A(L,J) = A(L,J) + T*V1
A(L1,J) = A(L1,J) + T*V2
150 CONTINUE
A(L1,K) = CMPLX(0.0,0.0)
DO 160 J=L,N
T = B(L,J) + CONJG(U2)*B(L1,J)
B(L,J) = B(L,J) + T*V1
B(L1,J) = B(L1,J) + T*V2
160 CONTINUE
IF (MATS.EQ.3) GO TO 180
DO 170 J=1,N
TF = F(L,J) + CONJG(U2)*F(L1,J)
F(L,J) = F(L,J) + TF*V1
F(L1,J) = F(L1,J) + TF*V2
170 CONTINUE
180 IF (MATS.NE.1) GO TO 200
DO 190 J=1,N
TG = G(L,J) + CONJG(U2) + G(L1,J)
TH = H(L,J) + CONJG(U2) + H(L1,J)
G(L,J) = G(L,J) + TG*V1
H(L,J) = H(L,J) + TH*V1
G(L1,J) = G(L1,J) + TG*V2
H(L1,J) = H(L1,J) + TH*V2
190 CONTINUE
C ********** ZERO B(L+1,L) **********
200 S = CABS(B(L1,L1)) + CABS(B(L1,L))
IF (S.EQ.0.0) GO TO 240
U1 = B(L1,L1)/CMPLX(S,0.0)
U2 = B(L1,L)/CMPLX(S,0.0)
R = SQRT(CABS(U1)**2+CABS(U2)**2)
RR = CMPLX(R,0.0)
IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR
V1 = -(U1+RR)/RR
V2 = -U2/RR
U2 = V2/V1
DO 210 I=1,L1
T = B(I,L1) + CONJG(U2)*B(I,L)
B(I,L1) = B(I,L1) + T*V1
B(I,L) = B(I,L) + T*V2
210 CONTINUE
B(L1,L) = CMPLX(0.0,0.0)
DO 220 I=1,N
T = A(I,L1) + CONJG(U2)*A(I,L)
A(I,L1) = A(I,L1) + T*V1
A(I,L) = A(I,L) + T*V2
220 CONTINUE
IF (MATS.EQ.0) GO TO 240
DO 230 I=1,N
T = Z(I,L1) + CONJG(U2)*Z(I,L)
Z(I,L1) = Z(I,L1) + T*V1
Z(I,L) = Z(I,L) + T*V2
230 CONTINUE
240 CONTINUE
250 CONTINUE
260 RETURN
END
SUBROUTINE CQZIT(NM, N, A, B, EPS1, MATS, Z, F, IERR) CQZ 10
C
C
C SUBROUTINE CQZIT IS A MODIFICATION OF THE EISPACK SUBROUTINE
C QZIT. ALL OPERATIONS ARE PERFORMED IN COMPLEX ARITHMETIC,
C AND THE LEFT TRANSFORMATIONS MAY ALSO BE APPLIED TO AN AUXILIARY
C MATRIX F.
C
C
C ON ENTRY,
C
C NM IS THE LEADING DIMENSION OF THE MATRICES A AND B IN
C THE MAIN PROGRAM.
C
C N IS THE ORDER OF THE MATRICES A AND B.
C
C A CONTAINS AN UPPER HESSENBERG MATRIX FROM CQZHES.
C
C B CONTAINS AN UPPER TRIANGULAR MATRIX FROM CQZHES.
C
C EPS1 IS A REAL NUMBER DEFINING THE TOLERANCE USED TO DETERMINE
C NEGLIGIBLE ELEMENTS OF A AND B IN THE COURSE OF THE ALG-
C ORITHM. AN ELEMENT OF EITHER MATRIX WILL BE CONSIDERED
C NEGLIGIBLE AND RESET TO ZERO IF IT IS NOT LARGER THAN THE
C PRODUCT OF EPS1 AND THE NORM OF THE MATRIX. IF EPS1.LE.0,
C RELATIVE MACHINE PRECISION WILL BE COMPUTED AND
C USED INSTEAD.
C
C MATS IS AN INTEGER INPUT VARIABLE. IT IS SET PRIOR
C TO THE CALL TO CQZHES.
C
C F CONTAINS AN AUXILIARY MATRIX.
C
C
C ON RETURN,
C
C A IS UPPER TRIANGULAR.
C
C B IS UPPER TRIANGULAR.
C
C Z CONTAINS THE HISTORY OF THE TRANSFORMATIONS, IF DESIRED.
C
C F CONTAINS THE AUXILIARY MATRIX, UPDATED IF DESIRED.
C
C IERR IS AN INTEGER ERROR RETURN WHICH INDICATES FAILURE
C OF THE QZ ALGORITHM TO REDUCE A SUBDIAGONAL ELEMENT
C TO ZERO AFTER 50 ITERATIONS.
C
INTEGER I, J, K, L, N, EN, JJ, K1, K2, LD, LL, L1, NA, NM, ISH,
* ITS, KM1, LM1
INTEGER ENM2, IERR, LOR1, ENORN
COMPLEX A(NM,N), B(NM,N), Z(NM,N)
COMPLEX A11, A21, A33, A34, A43, A44, B11, B22, B33, B34, B44
COMPLEX A1, A2, U1, U2, V1, V2, T, RR, SH, SS
COMPLEX F(NM,N), TF
COMPLEX CMPLX, CONJG, CSQRT
REAL EPS1, EPSA, EPSB, ANORM, BNORM, ANI, BNI, SRELPR, R, S
REAL SQRT, CABS
INTEGER MAX0, MIN0, MATS
IERR = 0
C ********** COMPUTE EPSA, EPSB **********
ANORM = 0.0
BNORM = 0.0
DO 20 I=1,N
ANI = 0.0
IF (I.NE.1) ANI = CABS(A(I,I-1))
BNI = 0.0
DO 10 J=I,N
ANI = ANI + CABS(A(I,J))
BNI = BNI + CABS(B(I,J))
10 CONTINUE
IF (ANI.GT.ANORM) ANORM = ANI
IF (BNI.GT.BNORM) BNORM = BNI
20 CONTINUE
IF (ANORM.EQ.0.0) ANORM = 1.0
IF (BNORM.EQ.0.0) BNORM = 1.0
SRELPR = EPS1
IF (SRELPR.GT.0.0) GO TO 40
C ********** COMPUTE ROUNDOFF LEVEL IF EPS1 IS ZERO **********
SRELPR = 1.0
30 SRELPR = SRELPR/2.0
IF (1.0+SRELPR.GT.1.0) GO TO 30
SRELPR = SRELPR + SRELPR
40 EPSA = SRELPR*ANORM
EPSB = SRELPR*BNORM
C ********** REDUCE A TO TRIANGULAR FORM, WHILE
C KEEPING B TRIANGULAR **********
LOR1 = 1
ENORN = N
EN = N
C ********** BEGIN QZ STEP **********
50 IF (EN.LE.1) GO TO 220
IF (MATS.EQ.0) ENORN = EN
ITS = 0
NA = EN - 1
ENM2 = NA - 1
60 ISH = 1
C ********** CHECK FOR CONVERGENCE OR REDUCIBILITY **********
DO 70 LL=1,EN
LM1 = EN - LL
L = LM1 + 1
IF (L.EQ.1) GO TO 90
IF (CABS(A(L,LM1)).LE.EPSA) GO TO 80
70 CONTINUE
80 A(L,LM1) = CMPLX(0.0,0.0)
IF (L.LT.NA) GO TO 90
C ********** 1-BY-1 BLOCK ISOLATED **********
EN = LM1
GO TO 50
C ********** CHECK FOR SMALL TOP OF B **********
90 LD = L
L1 = L + 1
B11 = B(L,L)
IF (CABS(B11).GT.EPSB) GO TO 120
B(L,L) = CMPLX(0.0,0.0)
S = CABS(A(L,L)) + CABS(A(L1,L))
U1 = A(L,L)/CMPLX(S,0.0)
U2 = A(L1,L)/CMPLX(S,0.0)
R = SQRT(CABS(U1)**2+CABS(U2)**2)
RR = CMPLX(R,0.0)
IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR
V1 = -(U1+RR)/RR
V2 = -U2/RR
U2 = V2/V1
DO 110 J=L,ENORN
T = A(L,J) + CONJG(U2)*A(L1,J)
A(L,J) = A(L,J) + T*V1
A(L1,J) = A(L1,J) + T*V2
T = B(L,J) + CONJG(U2)*B(L1,J)
B(L,J) = B(L,J) + T*V1
B(L1,J) = B(L1,J) + T*V2
IF (MATS.EQ.3) GO TO 110
DO 100 JJ=1,N
TF = F(L,JJ) + CONJG(U2)*F(L1,JJ)
F(L,JJ) = F(L,JJ) + TF*V1
F(L1,JJ) = F(L1,JJ) + TF*V2
100 CONTINUE
110 CONTINUE
IF (L.NE.1) A(L,LM1) = -A(L,LM1)
LM1 = L
L = L1
GO TO 80
120 A11 = A(L,L)/B11
A21 = A(L1,L)/B11
C ********** ITERATION STRATEGY **********
IF (ITS.EQ.50) GO TO 210
C ********** DETERMINE SHIFT **********
B22 = B(L1,L1)
IF (CABS(B22).LT.EPSB) B22 = CMPLX(EPSB,0.0)
B33 = B(NA,NA)
IF (CABS(B33).LT.EPSB) B33 = CMPLX(EPSB,0.0)
B44 = B(EN,EN)
IF (CABS(B44).LT.EPSB) B44 = CMPLX(EPSB,0.0)
A33 = A(NA,NA)/B33
A34 = A(NA,EN)/B44
A43 = A(EN,NA)/B33
A44 = A(EN,EN)/B44
B34 = B(NA,EN)/B44
T = CMPLX(0.5,0.0)*(A43*B34-A33-A44)
RR = T*T + A34*A43 - A33*A44
C ********** DETERMINE SINGLE SHIFT ZEROTH COLUMN OF A **********
RR = CSQRT(RR)
SH = -T + RR
SS = -T - RR
IF (CABS(SS-A44).LT.CABS(SH-A44)) SH = SS
A1 = A11 - SH
A2 = A21
IF (L.NE.LD) A(L,LM1) = -A(L,LM1)
IF (ITS.NE.10) GO TO 130
A1 = CMPLX(1.0,0.0)
A2 = CMPLX(1.1605,0.0)
130 ITS = ITS + 1
IF (MATS.EQ.0) LOR1 = LD
C ********** MAIN LOOP **********
DO 200 K=L,NA
K1 = K + 1
K2 = K + 2
KM1 = MAX0(K-1,L)
LL = MIN0(EN,K1+ISH)
C ********** ZERO A(K+1,K-1) **********
IF (K.EQ.L) GO TO 140
A1 = A(K,KM1)
A2 = A(K1,KM1)
140 S = CABS(A1) + CABS(A2)
IF (S.EQ.0.0) GO TO 60
U1 = A1/CMPLX(S,0.0)
U2 = A2/CMPLX(S,0.0)
R = SQRT(CABS(U1)**2+CABS(U2)**2)
RR = CMPLX(R,0.0)
IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR
V1 = -(U1+RR)/RR
V2 = -U2/RR
U2 = V2/V1
DO 150 J=KM1,ENORN
T = A(K,J) + CONJG(U2)*A(K1,J)
A(K,J) = A(K,J) + T*V1
A(K1,J) = A(K1,J) + T*V2
T = B(K,J) + CONJG(U2)*B(K1,J)
B(K,J) = B(K,J) + T*V1
B(K1,J) = B(K1,J) + T*V2
150 CONTINUE
IF (K.NE.L) A(K1,KM1) = CMPLX(0.0,0.0)
IF (MATS.EQ.3) GO TO 170
DO 160 J=1,N
TF = F(K,J) + CONJG(U2)*F(K1,J)
F(K,J) = F(K,J) + TF*V1
F(K1,J) = F(K1,J) + TF*V2
160 CONTINUE
C ********** ZERO B(K+1,K) **********
170 S = CABS(B(K1,K1)) + CABS(B(K1,K))
IF (S.EQ.0.0) GO TO 200
U1 = B(K1,K1)/CMPLX(S,0.0)
U2 = B(K1,K)/CMPLX(S,0.0)
R = SQRT(CABS(U1)**2+CABS(U2)**2)
RR = CMPLX(R,0.0)
IF (CABS(U1).NE.0.0) RR = (U1/CABS(U1))*RR
V1 = -(U1+RR)/RR
V2 = -U2/RR
U2 = V2/V1
DO 180 I=LOR1,LL
T = A(I,K1) + CONJG(U2)*A(I,K)
A(I,K1) = A(I,K1) + T*V1
A(I,K) = A(I,K) + T*V2
T = B(I,K1) + CONJG(U2)*B(I,K)
B(I,K1) = B(I,K1) + T*V1
B(I,K) = B(I,K) + T*V2
180 CONTINUE
B(K1,K) = CMPLX(0.0,0.0)
IF (MATS.EQ.0) GO TO 200
DO 190 I=1,N
T = Z(I,K1) + CONJG(U2)*Z(I,K)
Z(I,K1) = Z(I,K1) + T*V1
Z(I,K) = Z(I,K) + T*V2
190 CONTINUE
200 CONTINUE
C ********** END QZ STEP **********
GO TO 60
C ********** SET ERROR -- NEITHER BOTTOM SUBDIAGONAL ELEMENT
C HAS BECOME NEGLIGIBLE AFTER 50 ITERATIONS **********
210 IERR = EN
C ********** SAVE EPSB FOR USE BY CQZVEC **********
220 IF (N.GT.1) B(N,1) = CMPLX(EPSB,0.0)
RETURN
END
SUBROUTINE TRISLV(NM, N, U, V, L, F, TEMP, IERR) TRI 10
C
C
C SUBROUTINE TRISLV BACKSOLVES A SYSTEM OF THE FORM
C UY + VYL = F, WHERE U AND V ARE UPPER TRIANGULAR, AND
C L IS LOWER TRIANGULAR.
C
C
C ON ENTRY,
C
C NM IS THE LEADING DIMENSION OF THE MATRICES U, V, L AND F IN
C THE MAIN PROGRAM.
C
C N IS THE ORDER OF THE MATRICES U, V, L AND F.
C
C U CONTAINS AN UPPER TRIANGULAR MATRIX. IT IS THE LEFT
C COEFFICIENT OF Y IN THE FIRST TERM IN UY + VYL.
C
C V CONTAINS AN UPPER TRIANGULAR MATRIX. IT IS THE LEFT
C COEFFICIENT OF Y IN THE SECOND TERM OF UY + VYL.
C
C L CONTAINS A LOWER TRIANGULAR MATRIX. IT IS THE RIGHT
C COEFFICIENT OF Y IN THE SECOND TERM OF UY + VYL.
C
C F CONTAINS THE RIGHT HAND SIDE OF UY + VYL = F.
C
C TEMP CONTAINS A WORK VECTOR OF LENGTH AT LEAST N.
C
C ON RETURN,
C
C F CONTAINS THE SOLUTION Y.
C
C IERR IS AN ERROR RETURN DESIGNATING INCONSISTENCY OF THE
C ORIGINAL SYSTEM.
C IERR.EQ.0 FOR A NORMAL RETURN.
C IERR.EQ.1 IF THE TRIANGULAR SYSTEM IS INCONSISTENT.
C
C
INTEGER I, IERR, J, JP1, K, KK, KM1, M, N, NM, NM1
COMPLEX U(NM,N), V(NM,N), L(NM,N), F(NM,N)
COMPLEX TEMP(N)
COMPLEX CMPLX, DENOM, SUM
REAL CABS, S, T
IERR = 0
NM1 = N - 1
DO 120 KK=1,N
C ********** BACKSUBSTITUTE FOR ROW K. **********
K = N - KK + 1
IF (CABS(F(K,N)).NE.0.0) GO TO 10
F(K,N) = CMPLX(0.0,0.0)
GO TO 30
10 DENOM = U(K,K) + V(K,K)*L(N,N)
S = CABS(DENOM)
T = 1.0 + S/CABS(F(K,N))
IF (T.GT.1.0) GO TO 20
IERR = 1
RETURN
20 F(K,N) = F(K,N)/DENOM
IF (N.EQ.1) RETURN
30 DO 70 I=1,NM1
J = N - I
JP1 = J + 1
SUM = CMPLX(0.0,0.0)
DO 40 M=JP1,N
SUM = SUM + F(K,M)*L(M,J)
40 CONTINUE
SUM = F(K,J) - V(K,K)*SUM
IF (CABS(SUM).NE.0.0) GO TO 50
F(K,J) = CMPLX(0.0,0.0)
GO TO 70
50 DENOM = U(K,K) + V(K,K)*L(J,J)
S = CABS(DENOM)
T = 1.0 + S/CABS(SUM)
IF (T.GT.1.0) GO TO 60
IERR = 1
RETURN
60 F(K,J) = SUM/DENOM
70 CONTINUE
C ********** FORM TEMP = YK-TRANS*L. **********
IF (K.EQ.1) RETURN
KM1 = K - 1
DO 90 I=1,N
TEMP(I) = CMPLX(0.0,0.0)
DO 80 J=1,N
TEMP(I) = TEMP(I) + F(K,J)*L(J,I)
80 CONTINUE
90 CONTINUE
C ********** PREPARE F' WHICH IS (K-1) BY N. **********
DO 110 I=1,KM1
DO 100 J=1,N
F(I,J) = F(I,J) - U(I,K)*F(K,J) - V(I,K)*TEMP(J)
100 CONTINUE
110 CONTINUE
120 CONTINUE
RETURN
END