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