      SUBROUTINE F01BCE(N,TOL,Z,IZ,W,IW,D,E,C,S)
C
C     TRECX2
C     F01BCE REDUCES A COMPLEX HERMITIAN MATRIX TO REAL
C     TRIDIAGONAL FORM FROM WHICH THE EIGENVALUES AND EIGENVECTORS
C     CAN BE FOUND USING SUBROUTINE F02AYE,(CXTQL2). THE HERMITIAN
C     MATRIX A=A(1) IS REDUCED TO THE TRIDIAGONAL MATRIX A(N-1) BY
C     N-2 UNITARY TRANSFORMATIONS. THE HOUSEHOLDER REDUCTION ITSELF
C     DOES NOT GIVE A REAL TRIDIAGONAL MATRIX, THE OFF-DIAGONAL
C     ELEMENTS ARE COMPLEX. THEY ARE SUBSEQUENTLY MADE REAL BY A
C     DIAGONAL TRANSFORMATION.
C
C     REVISED TO INTRODUCE SCALING INTO
C     THE GENERATION OF HOUSEHOLDER MATRICES AS PROPOSED BY
C     G.W. STEWART, INTRODUCTION TO MATRIX COMPUTATIONS, CHAPTER 7.
C     TOL IS NOW A DUMMY PARAMETER.
C
C     .. Scalar Arguments ..
      INTEGER           IW, IZ, N
      REAL              TOL
C     .. Array Arguments ..
      DIMENSION C(N), D(N), E(N), S(N), W(IW,N), Z(IZ,N)
C     .. Local Scalars ..
      REAL              CO, F, FI, FR, G, GI, GR, H, HH, R, SCALE, SI
      INTEGER           I, II, J, K, KD, KE, L
C     .. External Functions ..
      REAL              DDOT
      INTEGER           IDAMAX
      EXTERNAL          DDOT, IDAMAX
C     .. External Subroutines ..
      EXTERNAL          F01BCY, F01BCZ, DSCAL
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, SQRT
C     .. Executable Statements ..
      DO 20 I = 1, N
         D(I) = Z(N,I)
         E(I) = -W(N,I)
   20 CONTINUE
      IF (N.EQ.1) GO TO 500
      DO 320 II = 2, N
         I = N - II + 2
         L = I - 1
         FR = D(I-1)
         FI = E(I-1)
         IF (ABS(FR)+ABS(FI).NE.0.0E0) GO TO 40
         R = 0.0E0
         CO = 1.0E0
         C(I) = 1.0E0
         SI = 0.0E0
         S(I) = 0.0E0
         GO TO 100
   40    IF (ABS(FR).LT.ABS(FI)) GO TO 60
         R = ABS(FR)*SQRT(1.0E0+(FI/FR)**2)
         GO TO 80
   60    R = ABS(FI)*SQRT(1.0E0+(FR/FI)**2)
   80    SI = FI/R
         S(I) = -SI
         CO = FR/R
         C(I) = CO
  100    IF (L.EQ.1) GO TO 240
C        FIND THE ELEMENTS OF LARGEST ABSOLUTE VALUE IN D AND E
         KD = IDAMAX(L,D,1)
         KE = IDAMAX(L,E,1)
         SCALE = MAX(ABS(D(KD)),ABS(E(KE)))
C        IF (D,E) IS A NULL VECTOR THEN SKIP THE TRANSFORMATION
         IF (SCALE.EQ.0.0E0) GO TO 240
         CALL DSCAL(L,1.0E0/SCALE,D,1)
         CALL DSCAL(L,1.0E0/SCALE,E,1)
         H = DDOT(L,D,1,D,1) + DDOT(L,E,1,E,1)
         G = -SQRT(H)
         E(I) = G*SCALE
C        E(I) HAS ITS FINAL REAL VALUE
         R = R/SCALE
         H = H - R*G
C        S*S + SR
         D(I-1) = (R-G)*CO
         E(I-1) = (R-G)*SI
         DO 120 J = 1, L
            Z(J,I) = D(J)
            W(J,I) = E(J)
  120    CONTINUE
         CALL F01BCZ(Z,IZ,W,IW,L,D,E,C,S)
C        FORM P
         DO 140 J = 1, L
            C(J) = C(J)/H
            S(J) = S(J)/H
  140    CONTINUE
         FR = 0.0D0
         DO 160 J = 1, L
            FR = FR + C(J)*D(J) + S(J)*E(J)
  160    CONTINUE
C        FORM K
         HH = FR/(H+H)
C        FORM Q
         DO 180 J = 1, L
            C(J) = C(J) - HH*D(J)
            S(J) = S(J) - HH*E(J)
  180    CONTINUE
C        NOW FORM REDUCED A
         DO 220 J = 1, L
            FR = D(J)
            FI = E(J)
            GR = C(J)
            GI = S(J)
            DO 200 K = J, L
               Z(K,J) = (((Z(K,J)-GR*D(K))-GI*E(K))-FR*C(K)) - FI*S(K)
               W(K,J) = (((W(K,J)-GR*E(K))+GI*D(K))-FR*S(K)) + FI*C(K)
  200       CONTINUE
            D(J) = Z(L,J)
            Z(I,J) = 0.0E0
            E(J) = -W(L,J)
            W(I,J) = 0.0E0
            W(J,J) = 0.0E0
  220    CONTINUE
         GO TO 300
  240    E(I) = R
         H = 0.0E0
         DO 260 J = 1, L
            Z(J,I) = D(J)
            W(J,I) = E(J)
  260    CONTINUE
         DO 280 J = 1, L
            Z(I,J) = 0.0E0
            D(J) = Z(I-1,J)
            W(I,J) = 0.0E0
            E(J) = -W(I-1,J)
  280    CONTINUE
  300    D(I) = H
  320 CONTINUE
C     WE NOW FORM THE PRODUCT OF THE
C     HOUSEHOLDER MATRICES, OVERWRITING
C     ON Z AND W
      DO 460 I = 2, N
         L = I - 1
         Z(N,L) = Z(L,L)
         Z(L,L) = 1.0E0
         W(N,L) = E(L)
         W(L,L) = 0.0E0
         H = D(I)
         IF (H.EQ.0.0E0) GO TO 420
         DO 340 K = 1, L
            D(K) = 0.0E0
            E(K) = 0.0E0
  340    CONTINUE
         CALL F01BCY(Z,IZ,W,IW,L,L,Z(1,I),W(1,I),D,E)
         DO 360 K = 1, L
            D(K) = D(K)/H
            E(K) = -E(K)/H
  360    CONTINUE
         DO 400 J = 1, L
            DO 380 K = 1, L
               Z(K,J) = Z(K,J) - Z(K,I)*D(J) + W(K,I)*E(J)
               W(K,J) = W(K,J) - Z(K,I)*E(J) - W(K,I)*D(J)
  380       CONTINUE
  400    CONTINUE
  420    DO 440 J = 1, L
            Z(J,I) = 0.0E0
            W(J,I) = 0.0E0
  440    CONTINUE
  460 CONTINUE
      W(N,N) = E(N)
      DO 480 I = 1, N
         D(I) = Z(N,I)
         Z(N,I) = 0.0E0
         E(I) = W(N,I)
         W(N,I) = 0.0E0
  480 CONTINUE
  500 Z(N,N) = 1.0E0
      W(N,N) = 0.0E0
      E(1) = 0.0E0
C     NOW WE MULTIPLY BY THE
C     COSTHETA + I SINTHETA COLUMN
C     FACTORS
      CO = 1.0E0
      SI = 0.0E0
      IF (N.EQ.1) RETURN
      DO 540 I = 2, N
         F = CO*C(I) - SI*S(I)
         SI = CO*S(I) + SI*C(I)
         CO = F
         DO 520 J = 1, N
            F = Z(J,I)*CO - W(J,I)*SI
            W(J,I) = Z(J,I)*SI + W(J,I)*CO
            Z(J,I) = F
  520    CONTINUE
  540 CONTINUE
      RETURN
      END
C
      SUBROUTINE F01BCY(AR,IAR,AI,IAI,M,N,BR,BI,CR,CI)
C
C     COMPUTES  C = C +  (A**H)*B  (COMPLEX) WHERE
C     A IS RECTANGULAR M BY N.
C     C MUST BE DISTINCT FROM B.
C
C
C     .. Scalar Arguments ..
      INTEGER           IAI, IAR, M, N
C     .. Array Arguments ..
      DIMENSION AI(IAI,N), AR(IAR,N), BI(M), BR(M), CI(N), CR(N)
C     .. Local Scalars ..
      REAL XI, XR
      INTEGER           I, J
C     .. Executable Statements ..
      DO 40 I = 1, N
         XR = CR(I)
         XI = CI(I)
         DO 20 J = 1, M
            XR = XR + AR(J,I)*BR(J) + AI(J,I)*BI(J)
            XI = XI + AR(J,I)*BI(J) - AI(J,I)*BR(J)
   20    CONTINUE
         CR(I) = XR
         CI(I) = XI
   40 CONTINUE
      RETURN
      END
C
      SUBROUTINE F01BCZ(AR,IAR,AI,IAI,N,BR,BI,CR,CI)
C
C     COMPUTES  C = A*B  (COMPLEX) WHERE
C     A IS A HERMITIAN N-BY-N MATRIX,
C     WHOSE LOWER TRIANGLE IS STORED IN A.
C     C MUST BE DISTINCT FROM B.
C
C
C     .. Scalar Arguments ..
      INTEGER           IAI, IAR, N
C     .. Array Arguments ..
      DIMENSION AI(IAI,N), AR(IAR,N), BI(N), BR(N), CI(N), CR(N)
C     .. Local Scalars ..
      REAL              YI, YR
      INTEGER           I, IP1, J, NM1
C     .. Executable Statements ..
      DO 20 I = 1, N
         CR(I) = 0.0E0
         CI(I) = 0.0E0
   20 CONTINUE
      IF (N.EQ.1) GO TO 100
      NM1 = N - 1
      DO 80 I = 1, NM1
         DO 40 J = I, N
            CR(J) = CR(J) + AR(J,I)*BR(I) - AI(J,I)*BI(I)
            CI(J) = CI(J) + AR(J,I)*BI(I) + AI(J,I)*BR(I)
   40    CONTINUE
         YR = CR(I)
         YI = CI(I)
         IP1 = I + 1
         DO 60 J = IP1, N
            YR = YR + AR(J,I)*BR(J) + AI(J,I)*BI(J)
            YI = YI + AR(J,I)*BI(J) - AI(J,I)*BR(J)
   60    CONTINUE
         CR(I) = YR
         CI(I) = YI
   80 CONTINUE
  100 CR(N) = CR(N) + AR(N,N)*BR(N) - AI(N,N)*BI(N)
      CI(N) = CI(N) + AR(N,N)*BI(N) + AI(N,N)*BR(N)
      RETURN
      END
C
      SUBROUTINE F06QXF( SIDE, PIVOT, DIRECT, M, N, K1, K2, C, S, A,
     $                   LDA )
C
C     .. Scalar Arguments ..
      INTEGER            K1, K2, LDA, M, N
      CHARACTER*1        DIRECT, PIVOT, SIDE
C     .. Array Arguments ..
      DIMENSION A( LDA, * ), C( * ), S( * )
C     ..
C
C  F06QXF  performs the transformation
C
C     A := P*A,   when   SIDE = 'L' or 'l'  (  Left-hand side )
C
C     A := A*P',  when   SIDE = 'R' or 'r'  ( Right-hand side )
C
C  where A is an m by n matrix and P is an orthogonal matrix, consisting
C  of a  sequence  of  plane  rotations,  applied  in  planes  k1 to k2,
C  determined by the parameters PIVOT and DIRECT as follows:
C
C     When  PIVOT  = 'V' or 'v'  ( Variable pivot )
C     and   DIRECT = 'F' or 'f'  ( Forward sequence ) then
C
C        P is given as a sequence of plane rotation matrices
C
C           P = P( k2 - 1 )*...*P( k1 + 1 )*P( k1 ),
C
C        where  P( k )  is a plane rotation matrix for the  ( k, k + 1 )
C        plane.
C
C     When  PIVOT  = 'V' or 'v'  ( Variable pivot )
C     and   DIRECT = 'B' or 'b'  ( Backward sequence ) then
C
C        P is given as a sequence of plane rotation matrices
C
C           P = P( k1 )*P( k1 + 1 )*...*P( k2 - 1 ),
C
C        where  P( k )  is a plane rotation matrix for the  ( k, k + 1 )
C        plane.
C
C     When  PIVOT  = 'T' or 't'  ( Top pivot )
C     and   DIRECT = 'F' or 'f'  ( Forward sequence ) then
C
C        P is given as a sequence of plane rotation matrices
C
C           P = P( k2 - 1 )*P( k2 - 2 )*...*P( k1 ),
C
C        where  P( k )  is a plane rotation matrix for the ( k1, k + 1 )
C        plane.
C
C     When  PIVOT  = 'T' or 't'  ( Top pivot )
C     and   DIRECT = 'B' or 'b'  ( Backward sequence ) then
C
C        P is given as a sequence of plane rotation matrices
C
C           P = P( k1 )*P( k1 + 1 )*...*P( k2 - 1 ),
C
C        where  P( k )  is a plane rotation matrix for the ( k1, k + 1 )
C        plane.
C
C     When  PIVOT  = 'B' or 'b'  ( Bottom pivot )
C     and   DIRECT = 'F' or 'f'  ( Forward sequence ) then
C
C        P is given as a sequence of plane rotation matrices
C
C           P = P( k2 - 1 )*P( k2 - 2 )*...*P( k1 ),
C
C        where  P( k )  is a  plane rotation  matrix  for the  ( k, k2 )
C        plane.
C
C     When  PIVOT  = 'B' or 'b'  ( Bottom pivot )
C     and   DIRECT = 'B' or 'b'  ( Backward sequence ) then
C
C        P is given as a sequence of plane rotation matrices
C
C           P = P( k1 )*P( k1 + 1 )*...*P( k2 - 1 ),
C
C        where  P( k )  is a  plane rotation  matrix  for the  ( k, k2 )
C        plane.
C
C  c( k ) and s( k )  must contain the  cosine and sine  that define the
C  matrix  P( k ).  The  two by two  plane rotation  part of the  matrix
C  P( k ), R( k ), is assumed to be of the form
C
C     R( k ) = (  c( k )  s( k ) ).
C              ( -s( k )  c( k ) )
C
C  If m, n or k1 are less than unity,  or k2 is not greater than k1,  or
C  SIDE = 'L' or 'l'  and  k2  is greater than  m, or  SIDE = 'R' or 'r'
C  and  k2  is greater than  n,  then an  immediate return  is effected.
C
C     .. Parameters ..
      REAL ONE, ZERO
      PARAMETER          ( ONE = 1.0E+0, ZERO = 0.0E+0 )
C     .. Local Scalars ..
      REAL AIJ, CTEMP, STEMP, TEMP
      INTEGER            I, J
      LOGICAL            LEFT, RIGHT
C     .. Intrinsic Functions ..
      INTRINSIC          MIN
C     ..
C     .. Executable Statements ..
      LEFT = ( SIDE.EQ.'L' ).OR.( SIDE.EQ.'l' )
      RIGHT = ( SIDE.EQ.'R' ).OR.( SIDE.EQ.'r' )
      IF( ( MIN( M, N, K1 ).LT.1 ).OR.( K2.LE.K1 ).OR.
     $    ( ( LEFT ).AND.( K2.GT.M ) ).OR.
     $    ( ( RIGHT ).AND.( K2.GT.N ) ) )RETURN
      IF( LEFT )THEN
         IF( ( PIVOT.EQ.'V' ).OR.( PIVOT.EQ.'v' ) )THEN
            IF( ( DIRECT.EQ.'F' ).OR.( DIRECT.EQ.'f' ) )THEN
               DO 20 J = 1, N
                  AIJ = A( K1, J )
                  DO 10 I = K1, K2 - 1
                     TEMP = A( I + 1, J )
                     A( I, J ) = S( I )*TEMP + C( I )*AIJ
                     AIJ = C( I )*TEMP - S( I )*AIJ
   10             CONTINUE
                  A( K2, J ) = AIJ
   20          CONTINUE
            ELSE IF( ( DIRECT.EQ.'B' ).OR.( DIRECT.EQ.'b' ) )THEN
               DO 40 J = 1, N
                  AIJ = A( K2, J )
                  DO 30 I = K2 - 1, K1, -1
                     TEMP = A( I, J )
                     A( I + 1, J ) = C( I )*AIJ - S( I )*TEMP
                     AIJ = S( I )*AIJ + C( I )*TEMP
   30             CONTINUE
                  A( K1, J ) = AIJ
   40          CONTINUE
            END IF
         ELSE IF( ( PIVOT.EQ.'T' ).OR.( PIVOT.EQ.'t' ) )THEN
            IF( ( DIRECT.EQ.'F' ).OR.( DIRECT.EQ.'f' ) )THEN
               DO 60 J = 1, N
                  TEMP = A( K1, J )
                  DO 50 I = K1, K2 - 1
                     AIJ = A( I + 1, J )
                     A( I + 1, J ) = C( I )*AIJ - S( I )*TEMP
                     TEMP = S( I )*AIJ + C( I )*TEMP
   50             CONTINUE
                  A( K1, J ) = TEMP
   60          CONTINUE
            ELSE IF( ( DIRECT.EQ.'B' ).OR.( DIRECT.EQ.'b' ) )THEN
               DO 80 J = 1, N
                  TEMP = A( K1, J )
                  DO 70 I = K2 - 1, K1, -1
                     AIJ = A( I + 1, J )
                     A( I + 1, J ) = C( I )*AIJ - S( I )*TEMP
                     TEMP = S( I )*AIJ + C( I )*TEMP
   70             CONTINUE
                  A( K1, J ) = TEMP
   80          CONTINUE
            END IF
         ELSE IF( ( PIVOT.EQ.'B' ).OR.( PIVOT.EQ.'b' ) )THEN
            IF( ( DIRECT.EQ.'F' ).OR.( DIRECT.EQ.'f' ) )THEN
               DO 100 J = 1, N
                  TEMP = A( K2, J )
                  DO 90 I = K1, K2 - 1
                     AIJ = A( I, J )
                     A( I, J ) = S( I )*TEMP + C( I )*AIJ
                     TEMP = C( I )*TEMP - S( I )*AIJ
   90             CONTINUE
                  A( K2, J ) = TEMP
  100          CONTINUE
            ELSE IF( ( DIRECT.EQ.'B' ).OR.( DIRECT.EQ.'b' ) )THEN
               DO 120 J = 1, N
                  TEMP = A( K2, J )
                  DO 110 I = K2 - 1, K1, -1
                     AIJ = A( I, J )
                     A( I, J ) = S( I )*TEMP + C( I )*AIJ
                     TEMP = C( I )*TEMP - S( I )*AIJ
  110             CONTINUE
                  A( K2, J ) = TEMP
  120          CONTINUE
            END IF
         END IF
      ELSE IF( RIGHT )THEN
         IF( ( PIVOT.EQ.'V' ).OR.( PIVOT.EQ.'v' ) )THEN
            IF( ( DIRECT.EQ.'F' ).OR.( DIRECT.EQ.'f' ) )THEN
               DO 140 J = K1, K2 - 1
                  IF( ( C( J ).NE.ONE ).OR.( S( J ).NE.ZERO ) )THEN
                     CTEMP = C( J )
                     STEMP = S( J )
                     DO 130 I = 1, M
                        TEMP = A( I, J + 1 )
                        A( I, J + 1 ) = CTEMP*TEMP - STEMP*A( I, J )
                        A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
  130                CONTINUE
                  END IF
  140          CONTINUE
            ELSE IF( ( DIRECT.EQ.'B' ).OR.( DIRECT.EQ.'b' ) )THEN
               DO 160 J = K2 - 1, K1, -1
                  IF( ( C( J ).NE.ONE ).OR.( S( J ).NE.ZERO ) )THEN
                     CTEMP = C( J )
                     STEMP = S( J )
                     DO 150 I = M, 1, -1
                        TEMP = A( I, J + 1 )
                        A( I, J + 1 ) = CTEMP*TEMP - STEMP*A( I, J )
                        A( I, J ) = STEMP*TEMP + CTEMP*A( I, J )
  150                CONTINUE
                  END IF
  160          CONTINUE
            END IF
         ELSE IF( ( PIVOT.EQ.'T' ).OR.( PIVOT.EQ.'t' ) )THEN
            IF( ( DIRECT.EQ.'F' ).OR.( DIRECT.EQ.'f' ) )THEN
               DO 180 J = K1 + 1, K2
                  CTEMP = C( J - 1 )
                  STEMP = S( J - 1 )
                  IF( ( CTEMP.NE.ONE ).OR.( STEMP.NE.ZERO ) )THEN
                     DO 170 I = 1, M
                        TEMP = A( I, J )
                        A( I, J ) = CTEMP*TEMP - STEMP*A( I, K1 )
                        A( I, K1 ) = STEMP*TEMP + CTEMP*A( I, K1 )
  170                CONTINUE
                  END IF
  180          CONTINUE
            ELSE IF( ( DIRECT.EQ.'B' ).OR.( DIRECT.EQ.'b' ) )THEN
               DO 200 J = K2, K1 + 1, -1
                  CTEMP = C( J - 1 )
                  STEMP = S( J - 1 )
                  IF( ( CTEMP.NE.ONE ).OR.( STEMP.NE.ZERO ) )THEN
                     DO 190 I = M, 1, -1
                        TEMP = A( I, J )
                        A( I, J ) = CTEMP*TEMP - STEMP*A( I, K1 )
                        A( I, K1 ) = STEMP*TEMP + CTEMP*A( I, K1 )
  190                CONTINUE
                  END IF
  200          CONTINUE
            END IF
         ELSE IF( ( PIVOT.EQ.'B' ).OR.( PIVOT.EQ.'b' ) )THEN
            IF( ( DIRECT.EQ.'F' ).OR.( DIRECT.EQ.'f' ) )THEN
               DO 220 J = K1, K2 - 1
                  IF( ( C( J ).NE.ONE ).OR.( S( J ).NE.ZERO ) )THEN
                     CTEMP = C( J )
                     STEMP = S( J )
                     DO 210 I = 1, M
                        TEMP = A( I, J )
                        A( I, J ) = STEMP*A( I, K2 ) + CTEMP*TEMP
                        A( I, K2 ) = CTEMP*A( I, K2 ) - STEMP*TEMP
  210                CONTINUE
                  END IF
  220          CONTINUE
            ELSE IF( ( DIRECT.EQ.'B' ).OR.( DIRECT.EQ.'b' ) )THEN
               DO 240 J = K2 - 1, K1, -1
                  IF( ( C( J ).NE.ONE ).OR.( S( J ).NE.ZERO ) )THEN
                     CTEMP = C( J )
                     STEMP = S( J )
                     DO 230 I = M, 1, -1
                        TEMP = A( I, J )
                        A( I, J ) = STEMP*A( I, K2 ) + CTEMP*TEMP
                        A( I, K2 ) = CTEMP*A( I, K2 ) - STEMP*TEMP
  230                CONTINUE
                  END IF
  240          CONTINUE
            END IF
         END IF
      END IF
C
      RETURN
      END
C
      SUBROUTINE P01ABZ
C
C     Terminates execution when a hard failure occurs.
C
C     ******************** IMPLEMENTATION NOTE ********************
C     The following STOP statement may be replaced by a call to an
C     implementation-dependent routine to display a message and/or
C     to abort the program.
C     *************************************************************
C     .. Executable Statements ..
      STOP
      END
C
      SUBROUTINE X04AAF(I,NERR)
C
C     IF I = 0, SETS NERR TO CURRENT ERROR MESSAGE UNIT NUMBER
C     (STORED IN NERR1).
C     IF I = 1, CHANGES CURRENT ERROR MESSAGE UNIT NUMBER TO
C     VALUE SPECIFIED BY NERR.
C
C     .. Scalar Arguments ..
      INTEGER           I, NERR
C     .. Local Scalars ..
      INTEGER           NERR1
C     .. Save statement ..
      SAVE              NERR1
C     .. Data statements ..
      DATA              NERR1/6/
C     .. Executable Statements ..
      IF (I.EQ.0) NERR = NERR1
      IF (I.EQ.1) NERR1 = NERR
      RETURN
      END
C
      SUBROUTINE X04BAF(NOUT,REC)
C
C     X04BAF writes the contents of REC to the unit defined by NOUT.
C
C     Trailing blanks are not output, except that if REC is entirely
C     blank, a single blank character is output.
C     If NOUT.lt.0, i.e. if NOUT is not a valid Fortran unit identifier,
C     then no output occurs.
C
C     .. Scalar Arguments ..
      INTEGER           NOUT
      CHARACTER*(*)     REC
C     .. Local Scalars ..
      INTEGER           I
C     .. Intrinsic Functions ..
      INTRINSIC         LEN
C     .. Executable Statements ..
      IF (NOUT.GE.0) THEN
C        Remove trailing blanks
         DO 20 I = LEN(REC), 2, -1
            IF (REC(I:I).NE.' ') GO TO 40
   20    CONTINUE
C        Write record to external file
   40    WRITE (NOUT,FMT=99999) REC(1:I)
      END IF
      RETURN
C
99999 FORMAT (A)
      END
C
      FUNCTION F06EAF( N, X, INCX, Y, INCY )
C
C     .. Entry Points ..
      REAL                      DDOT
      ENTRY                     DDOT  ( N, X, INCX, Y, INCY )
C     .. Scalar Arguments ..
      INTEGER                           INCX, INCY, N
C     .. Array Arguments ..
      DIMENSION                         X( * ), Y( * )
C     ..
C
C  F06EAF returns the value
C
C     F06EAF = x'y
C
C     .. Parameters ..
      REAL                  ZERO
      PARAMETER           ( ZERO = 0.0E+0 )
C     .. Local Scalars ..
      REAL                  SUM
      INTEGER               I, IX, IY
C     ..
C     .. Executable Statements ..
      SUM = ZERO
      IF( N.GT.0 )THEN
         IF( ( INCX.EQ.INCY ).AND.( INCX.GT.0 ) )THEN
            DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX
               SUM = SUM + X( IX )*Y( IX )
   10       CONTINUE
         ELSE
            IF( INCY.GE.0 )THEN
               IY = 1
            ELSE
               IY = 1 - ( N - 1 )*INCY
            END IF
            IF( INCX.GT.0 )THEN
               DO 20, IX = 1, 1 + ( N - 1 )*INCX, INCX
                  SUM = SUM + X( IX )*Y( IY )
                  IY  = IY  + INCY
   20          CONTINUE
            ELSE
               IX = 1 - ( N - 1 )*INCX
               DO 30, I = 1, N
                  SUM = SUM + X( IX )*Y( IY )
                  IX  = IX  + INCX
                  IY  = IY  + INCY
   30          CONTINUE
            END IF
         END IF
      END IF
C
      F06EAF = SUM
      RETURN
      END
C
      INTEGER FUNCTION IDAMAX( N, X, INCX )
C     .. Scalar Arguments ..
      INTEGER                  INCX, N
C     .. Array Arguments ..
      DIMENSION                X( * )
C     ..
C
C  F06JLF returns the smallest value of i such that
C
C     abs( x( i ) ) = max( abs( x( j ) ) )
C                      j
C
C     .. Local Scalars ..
      REAL                     XMAX
      INTEGER                  I, IMAX, IX
C     .. Intrinsic Functions ..
      INTRINSIC                ABS
C     ..
C     .. Executable Statements ..
      IF( N.GT.0 )THEN
         IMAX = 1
         IF( N.GT.1 )THEN
            XMAX = ABS( X( 1 ) )
            IX   = 1
            DO 10, I = 2, N
               IX = IX + INCX
               IF( XMAX.LT.ABS( X( IX ) ) )THEN
                  XMAX = ABS( X( IX ) )
                  IMAX = I
               END IF
   10       CONTINUE
         END IF
      ELSE
         IMAX = 0
      END IF
C
      IDAMAX = IMAX
      RETURN
      END
C
      INTEGER FUNCTION P01ABF(IFAIL,IERROR,SRNAME,NREC,REC)
C
C     P01ABF either returns the value of IERROR through the routine
C     name (soft failure), or terminates execution of the program
C     (hard failure). Diagnostic messages may be output.
C
C     If IERROR = 0 (successful exit from the calling routine),
C     the value 0 is returned through the routine name, and no
C     message is output
C
C     If IERROR is non-zero (abnormal exit from the calling routine),
C     the action taken depends on the value of IFAIL.
C
C     IFAIL =  1: soft failure, silent exit (i.e. no messages are
C                 output)
C     IFAIL = -1: soft failure, noisy exit (i.e. messages are output)
C     IFAIL =-13: soft failure, noisy exit but standard messages from
C                 P01ABF are suppressed
C     IFAIL =  0: hard failure, noisy exit
C
C     P01ABF also allows an alternative specification of IFAIL in which
C     it is regarded as a decimal integer with least significant digits
C     cba. Then
C
C     a = 0: hard failure  a = 1: soft failure
C     b = 0: silent exit   b = 1: noisy exit
C
C     except that hard failure now always implies a noisy exit.
C
C     .. Scalar Arguments ..
      INTEGER                 IERROR, IFAIL, NREC
      CHARACTER*(*)           SRNAME
C     .. Array Arguments ..
      CHARACTER*(*)           REC(*)
C     .. Local Scalars ..
      INTEGER                 I, NERR
      CHARACTER*72            MESS
C     .. External Subroutines ..
      EXTERNAL                P01ABZ, X04AAF, X04BAF
C     .. Intrinsic Functions ..
      INTRINSIC               ABS, MOD
C     .. Executable Statements ..
      IF (IERROR.NE.0) THEN
C        Abnormal exit from calling routine
         IF (IFAIL.EQ.-1 .OR. IFAIL.EQ.0 .OR. IFAIL.EQ.-13 .OR.
     *       (IFAIL.GT.0 .AND. MOD(IFAIL/10,10).NE.0)) THEN
C           Noisy exit
            CALL X04AAF(0,NERR)
            DO 20 I = 1, NREC
               CALL X04BAF(NERR,REC(I))
   20       CONTINUE
            IF (IFAIL.NE.-13) THEN
               WRITE (MESS,FMT=99999) SRNAME, IERROR
               CALL X04BAF(NERR,MESS)
               IF (ABS(MOD(IFAIL,10)).NE.1) THEN
C                 Hard failure
                  CALL X04BAF(NERR,
     *                     ' ** NAG hard failure - execution terminated'
     *                        )
                  CALL P01ABZ
               ELSE
C                 Soft failure
                  CALL X04BAF(NERR,
     *                        ' ** NAG soft failure - control returned')
               END IF
            END IF
         END IF
      END IF
      P01ABF = IERROR
      RETURN
C
99999 FORMAT (' ** ABNORMAL EXIT from Library routine ',A,': IFAIL',
     *  ' =',I6)
      END
C
      FUNCTION X02AJE()
C
C     RETURNS  (1/2)*B**(1-P)  IF ROUNDS IS .TRUE.
C     RETURNS  B**(1-P)  OTHERWISE
C
C     .. Constants ..
      DOUBLE PRECISION Z
      DATA Z/1.1102230246251568D-16/
C
C     Z has been derived from (2**(-53))+2**(-105)
C
C     .. Executable Statements ..
      X02AJE = Z   
      RETURN
      END
C
      SUBROUTINE F02AYE(N,EPS,D,E,Z,IZ,W,IW,IFAIL)
C
C     CXTQL2
C     THIS SUBROUTINE FINDS THE EIGENVALUES AND EIGENVECTORS OF A
C     HERMITIAN MATRIX, WHICH HAS BEEN REDUCED TO A REAL
C     TRIDIAGONAL MATRIX, T, GIVEN WITH ITS DIAGONAL ELEMENTS IN
C     THE ARRAY D(N) AND ITS SUB-DIAGONAL ELEMENTS IN THE LAST N
C     - 1 STORES OF THE ARRAY E(N), USING QL TRANSFORMATIONS. THE
C     EIGENVALUES ARE OVERWRITTEN ON THE DIAGONAL ELEMENTS IN THE
C     ARRAY D IN ASCENDING ORDER. THE REAL AND IMAGINARY PARTS OF
C     THE EIGENVECTORS ARE FORMED IN THE ARRAYS Z,W(N,N)
C     RESPECTIVELY, OVERWRITING THE ACCUMULATED TRANSFORMATIONS AS
C     SUPPLIED BY THE SUBROUTINE F01BCE. THE SUBROUTINE WILL FAIL
C     IF ALL EIGENVALUES TAKE MORE THAN 30*N ITERATIONS
C
C     .. Parameters ..
      INTEGER           VLEN
      PARAMETER         (VLEN=128)
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='F02AYE')
C     .. Scalar Arguments ..
      REAL              EPS
      INTEGER           IFAIL, IW, IZ, N
C     .. Array Arguments ..
      DIMENSION         D(N), E(N), W(IW,N), Z(IZ,N)
C     .. Local Scalars ..
      REAL              B, C, F, G, H, P, R, S
      INTEGER           I, I1, IPOS, ISAVE, ISEG, J, K, L, M
C     .. Local Arrays ..
      DIMENSION         CC(VLEN), SS(VLEN)
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. External Subroutines ..
      EXTERNAL          F06QXF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, SQRT
C     .. Executable Statements ..
      ISAVE = IFAIL
      IF (N.EQ.1) GO TO 40
      DO 20 I = 2, N
         E(I-1) = E(I)
   20 CONTINUE
   40 E(N) = 0.0E0
      B = 0.0E0
      F = 0.0E0
      J = 30*N
      DO 300 L = 1, N
         H = EPS*(ABS(D(L))+ABS(E(L)))
         IF (B.LT.H) B = H
C        LOOK FOR SMALL SUB-DIAG ELEMENT
         DO 60 M = L, N
            IF (ABS(E(M)).LE.B) GO TO 80
   60    CONTINUE
   80    IF (M.EQ.L) GO TO 280
  100    IF (J.LE.0) GO TO 400
         J = J - 1
C        FORM SHIFT
         G = D(L)
         H = D(L+1) - G
         IF (ABS(H).GE.ABS(E(L))) GO TO 120
         P = H*0.5E0/E(L)
         R = SQRT(P*P+1.0E0)
         H = P + R
         IF (P.LT.0.0E0) H = P - R
         D(L) = E(L)/H
         GO TO 140
  120    P = 2.0E0*E(L)/H
         R = SQRT(P*P+1.0E0)
         D(L) = E(L)*P/(1.0E0+R)
  140    H = G - D(L)
         I1 = L + 1
         IF (I1.GT.N) GO TO 180
         DO 160 I = I1, N
            D(I) = D(I) - H
  160    CONTINUE
  180    F = F + H
C        QL TRANSFORMATION
         P = D(M)
         C = 1.0E0
         S = 0.0E0
         DO 260 K = M - 1, L, -VLEN
            ISEG = MAX(K-VLEN+1,L)
            DO 240 I = K, ISEG, -1
               G = C*E(I)
               H = C*P
               IF (ABS(P).LT.ABS(E(I))) GO TO 200
               C = E(I)/P
               R = SQRT(C*C+1.0E0)
               E(I+1) = S*P*R
               S = C/R
               C = 1.0E0/R
               GO TO 220
  200          C = P/E(I)
               R = SQRT(C*C+1.0E0)
               E(I+1) = S*E(I)*R
               S = 1.0E0/R
               C = C/R
  220          P = C*D(I) - S*G
               D(I+1) = H + S*(C*G+S*D(I))
C           STORE ROTATIONS
               CC(VLEN-K+I) = C
               SS(VLEN-K+I) = -S
  240       CONTINUE
C        UPDATE VECTORS
            IPOS = VLEN - K + ISEG
            CALL F06QXF('Right','Variable','Backward',N,K-ISEG+2,1,
     *                  K-ISEG+2,CC(IPOS),SS(IPOS),Z(1,ISEG),IZ)
            CALL F06QXF('Right','Variable','Backward',N,K-ISEG+2,1,
     *                  K-ISEG+2,CC(IPOS),SS(IPOS),W(1,ISEG),IW)
  260    CONTINUE
         E(L) = S*P
         D(L) = C*P
         IF (ABS(E(L)).GT.B) GO TO 100
  280    D(L) = D(L) + F
  300 CONTINUE
C     ORDER EIGENVALUES AND EIGENVECTORS
      DO 380 I = 1, N
         K = I
         P = D(I)
         I1 = I + 1
         IF (I1.GT.N) GO TO 340
         DO 320 J = I1, N
            IF (D(J).GE.P) GO TO 320
            K = J
            P = D(J)
  320    CONTINUE
  340    IF (K.EQ.I) GO TO 380
         D(K) = D(I)
         D(I) = P
         DO 360 J = 1, N
            P = Z(J,I)
            Z(J,I) = Z(J,K)
            Z(J,K) = P
            P = W(J,I)
            W(J,I) = W(J,K)
            W(J,K) = P
  360    CONTINUE
  380 CONTINUE
      IFAIL = 0
      RETURN
  400 IFAIL = P01ABF(ISAVE,1,SRNAME,0,P01REC)
      RETURN
      END
C
      SUBROUTINE DSCAL( N, ALPHA, X, INCX )
C     .. Scalar Arguments ..
      REAL               ALPHA
      INTEGER            INCX, N
C     .. Array Arguments ..
      DIMENSION          X( * )
C     ..
C
C  F06EDF performs the operation
C
C     x := alpha*x
C
C     .. Parameters ..
      REAL               ONE         , ZERO
      PARAMETER        ( ONE = 1.0E+0, ZERO = 0.0E+0 )
C     .. Local Scalars ..
      INTEGER            IX
C     ..
C     .. Executable Statements ..
      IF( N.GT.0 )THEN
         IF( ALPHA.EQ.ZERO )THEN
            DO 10, IX = 1, 1 + ( N - 1 )*INCX, INCX
               X( IX ) = ZERO
   10       CONTINUE
         ELSE IF( ALPHA.EQ.( -ONE ) )THEN
            DO 20, IX = 1, 1 + ( N - 1 )*INCX, INCX
               X( IX ) = -X( IX )
   20       CONTINUE
         ELSE IF( ALPHA.NE.ONE )THEN
            DO 30, IX = 1, 1 + ( N - 1 )*INCX, INCX
               X( IX ) = ALPHA*X( IX )
   30       CONTINUE
         END IF
      END IF
C
      RETURN
      END
