      SUBROUTINE A02AAF(XXR,XXI,YR,YI)
C     MARK 2A RELEASE.  NAG COPYRIGHT 1973
C     MARK 4.5 REVISED
C     MARK 5C REVISED
C     MARK 11C REVISED. IER-467 (MAR 1985)
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     COMPUTES THE SQUARE ROOT OF A COMPLEX NUMBER
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  XXI, XXR, YI, YR
C     .. Local Scalars ..
      DOUBLE PRECISION  H, HALF, ONE, XI, XR, ZERO
C     .. External Functions ..
      DOUBLE PRECISION  A02ABF
      EXTERNAL          A02ABF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, SQRT
C     .. Data statements ..
      DATA              ZERO/0.0D0/, HALF/0.5D0/, ONE/1.0D0/
C     .. Executable Statements ..
C
      XR = ABS(XXR)
      XI = XXI
      IF (XR.GT.ONE) H = SQRT(XR*HALF+A02ABF(XR*HALF,XI*HALF))
      IF (XR.LE.ONE) H = SQRT(XR+A02ABF(XR,XI))*SQRT(HALF)
      IF (XI.NE.ZERO) XI = XI/(H+H)
      IF (XXR.LT.ZERO) GO TO 20
      YR = H
      YI = XI
      RETURN
   20 IF (XI.LT.ZERO) GO TO 40
      YR = XI
      YI = H
      RETURN
   40 YR = -XI
      YI = -H
      RETURN
      END
      SUBROUTINE C02AGU(A0,B0,C0,ZSM,ZLG)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C     BASED ON THE ROUTINE  CQDRTC, WRITTEN BY BRIAN T. SMITH
C
C     THIS SUBROUTINE DETERMINES THE ROOTS OF THE QUADRATIC EQUATION
C
C         A0*X**2 + B0*X + C0
C
C     WHERE A0, B0, AND C0 ARE COMPLEX COEFFICIENTS, AND ZSM AND ZLG
C     THE SMALLEST AND LARGEST ROOT IN MAGNITUDE, RESPECTIVELY.
C
C     THE ROOTS ARE COMPUTED TO WITHIN A RELATIVE ERROR OF A FEW
C     UNITS IN THE LAST PLACE (DEPENDING ON THE ACCURACY OF THE
C     BASIC ARITHMETIC OPERATIONS) EXCEPT WHEN UNDERFLOW OR OVERFLOW
C     OCCURS IN WHICH CASE THE TRUE ROOTS ARE WITHIN A FEW UNITS IN
C     THE LAST PLACE OF THE UNDERFLOW OR OVERFLOW THRESHOLD.
C
C     IF THE LEADING COEFFICIENT IS ZERO, THE LARGER ROOT IS
C     SET TO THE LARGEST MACHINE REPRESENTABLE NUMBER AND THE
C     OVERFLOW FLAG OVFLOW IS SET TRUE.  IF ALL THREE COEFFICIENTS ARE
C     ZERO, THE OVERFLOW FLAG IS SET, BOTH ROOTS ARE SET TO THE LARGEST
C     REPRESENTABLE NUMBER, BUT NO DIVIDE CHECK IS CREATED.
C
C     THIS PROGRAM IS DESIGNED TO TAKE ADVANTAGE OF SYSTEMS THAT REPORT
C     OVERFLOW AND UNDERFLOW CONDITIONS IN AN EFFICIENT WAY.  THAT IS,
C     IF, WHENEVER AN OVERFLOW OR UNDERFLOW OCCURS, CERTAIN FLAGS ARE
C     SET (THAT IS, THE LOGICAL VARIABLES OVFLOW AND UNFLOW IN THE
C     COMMON BLOCK AC02AG), C02AGU CAN USE THESE INDICATORS TO INDICATE
C     THAT THE ROOTS OVERFLOW OR UNDERFLOW AND CANNOT BE REPRESENTED.
C
C     HOWEVER, AS IMPLEMENTED IN THE NAG LIBRARY, THE ROUTINE SIMPLY
C     ASSUMES THAT THE MACHINE TERMINATES ON OVERFLOW AND IGNORES
C     UNDERFLOW.
C
C     C02AGX -- DETERMINE THE EXPONENT OF A NUMBER IN TERMS OF THE
C               MODEL.
C     C02AGR -- FORM A NUMBER WITH A GIVEN MANTISSA AND EXPONENT
C               PRECISION.
C
C     .. Parameters ..
      DOUBLE PRECISION  HALF, ONE, ZERO, TWO
      PARAMETER         (HALF=0.5D0,ONE=1.0D0,ZERO=0.0D0,TWO=2.0D0)
C     .. Array Arguments ..
      DOUBLE PRECISION  A0(2), B0(2), C0(2), ZLG(2), ZSM(2)
C     .. Scalars in Common ..
      DOUBLE PRECISION  DEPS, FINITY, SQRTFY, SQRTTY, TINY
      INTEGER           EMAXM1, EMINM1, EXPDEP, LRGEXP
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C nag routines
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      LOGICAL           OVFLOW, UNFLOW
C     .. Local Scalars ..
      DOUBLE PRECISION  SC, XC1, XC2
      INTEGER           EXPBSQ, SCLEXP
C     .. Local Arrays ..
      DOUBLE PRECISION  A(2), B(2), C(2), CT(2), D(2)
C     .. External Functions ..
      DOUBLE PRECISION  C02AGR
      INTEGER           C02AGX
      EXTERNAL          C02AGR, C02AGX
C     .. External Subroutines ..
      EXTERNAL          A02AAF, A02ACF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, MIN
C     .. Common blocks ..
      COMMON            /AC02AG/OVFLOW, UNFLOW
      COMMON            /BC02AG/FINITY, SQRTFY, SQRTTY, TINY, DEPS,
     *                  EMINM1, EMAXM1, EXPDEP, LRGEXP
C     .. Statement Functions ..
      DOUBLE PRECISION  APPABS
C     .. Save statement ..
      SAVE              /AC02AG/, /BC02AG/
C     .. Statement Function definitions ..
      APPABS(XC1,XC2) = MAX(ABS(XC1),ABS(XC2))
C     .. Executable Statements ..
C
C     INITIALIZE LOCAL VARIABLES WITH THE INPUT COEFFICIENTS.
C
      A(1) = A0(1)
      A(2) = A0(2)
      B(1) = -B0(1)
      B(2) = -B0(2)
      C(1) = C0(1)
      C(2) = C0(2)
C
C     CHECK FOR  A = CMPLX(ZERO, ZERO)  OR  C = CMPLX(ZERO, ZERO).
C
      IF (APPABS(A(1),A(2)).NE.ZERO) THEN
         IF (APPABS(C(1),C(2)).NE.ZERO) THEN
C
C           AT THIS POINT, A AND C ARE NON-ZERO.
C
C           SCALE THE COEFFICIENTS SO THAT THE PRODUCT A * C IS NEAR
C           1.0D0 IN MAGNITUDE.  THIS AVOIDS SPURIOUS UNDERFLOW/OVERFLOW
C           CONDITIONS WHEN THE TRUE RESULTS ARE WITHIN RANGE.
C
C           THE SCALE FACTOR IS A POWER OF THE BASE NEAR TO
C           SQRT(ABS(A*C)).  THIS CHOICE AVOIDS UNNECESSARY ROUNDING
C           ERRORS BUT IS EXPENSIVE TO COMPUTE WHEN FLOATING POINT
C           MANIPULATIVE FUNCTIONS ARE NOT AVAILABLE IN MACHINE CODE.
C
            SCLEXP = (C02AGX(APPABS(A(1),A(2)))+C02AGX(APPABS(C(1),C(2))
     *               ))/2
C
C           THE SCALE FACTOR IS  BASE ** SCLEXP.  IF A AND C ARE SCALED
C           USING THIS SCALE FACTOR AS A DIVIDEND, THEN THE
C           THE SCALED PRODUCT A'*C' IS BETWEEN BASE**(-2) AND
C           BASE IN MAGNITUDE, WHERE BASE IS THE BASE FOR MODEL NUMBERS
C           OF THE TYPE OF A.
C
C           BUT BEFORE PERFORMING THE SCALING, CHECK TO SEE IF IT IS
C           NECESSARY -- THAT IS, IF B IS SO LARGE IN MAGNITUDE THAT
C           B**2 EXCEEDS ABS(4*A*C) BY MORE THAN THE RELATIVE MACHINE
C           PRECISION FOR THE DOUBLE PRECISION DATA TYPE,
C           THE DISCRIMINANT IS IN EFFECT B AND NO SCALING IS REQUIRED.
C           HOWEVER, IF B IS SO SMALL IN MAGNITUDE THAT ABS(4*A*C)
C           EXCEEDS B**2 IN MAGNITUDE BY MORE THAN THIS SAME RELATIVE
C           MACHINE PRECISION, B IS IN EFFECT ZERO, BUT A AND C ARE
C           STILL SCALED TO AVOID SPURIOUS UNDERFLOWS/OVERFLOWS.
C
C           COMPUTE THE EXPONENT OF THE SQUARE OF THE SCALED B.
C
            IF (APPABS(B(1),B(2)).NE.ZERO) THEN
               EXPBSQ = 2*(C02AGX(APPABS(B(1),B(2)))-SCLEXP)
            ELSE
               EXPBSQ = -2*EXPDEP
            END IF
C
C           CHECK IF B**2 IS TOO BIG.
C
            IF (EXPBSQ.LE.EXPDEP) THEN
C
C              B**2 IS NOT TOO BIG.  SCALING WILL BE PERFORMED.
C
C              A AND C SHOULD BE SCALED USING THE USUAL SCALE
C              MANIPULATION FUNCTION BUT FOR EFFICIENCY, THE
C              SCALING IS PERFORMED BY DIVISION.
C
               SCLEXP = MIN(SCLEXP+1,EMAXM1)
               SCLEXP = MAX(SCLEXP,EMINM1)
               SC = C02AGR(ONE,SCLEXP)
C
C              CHECK IF IT IS TOO SMALL.
C
               IF (EXPBSQ.LT.-EXPDEP) THEN
C
C                 B IS TOO SMALL.  SET IT TO ZERO.
C
                  B(1) = ZERO
                  B(2) = ZERO
               ELSE
C
C                 B IS NEITHER TOO LARGE NOR TOO SMALL.  SCALE IT.
C
                  B(1) = (B(1)/SC)*HALF
                  B(2) = (B(2)/SC)*HALF
               END IF
               A(1) = A(1)/SC
               A(2) = A(2)/SC
               C(1) = C(1)/SC
               C(2) = C(2)/SC
C
C              THE MAGNITUDE OF THE DISCRIMINANT WILL NOT UNDERFLOW
C              OR OVERFLOW -- HOWEVER, A COMPONENT OF IT MAY UNDERFLOW.
C
               CT(1) = B(1)*B(1) - B(2)*B(2) - A(1)*C(1) + A(2)*C(2)
               CT(2) = TWO*B(2)*B(1) - A(2)*C(1) - A(1)*C(2)
               CALL A02AAF(CT(1),CT(2),D(1),D(2))
C
C              IN ORDER TO ENSURE THAT THE LARGER ROOT IS ASSIGNED TO
C              ZLG, SELECT THE SIGN OF D SO THAT B+D IS LARGER IN
C              MAGNITUDE THAN B-D.  (THIS CONDITION REDUCES TO THE
C              CONDITION THAT REAL(B)*REAL(D)+AIMAG(B)*AIMAG(D)>0.)
C
               IF (D(1)*B(1)+D(2)*B(2).LE.ZERO) THEN
                  D(1) = -D(1)
                  D(2) = -D(2)
               END IF
               B(1) = B(1) + D(1)
               B(2) = B(2) + D(2)
            END IF
C
C           AT THIS POINT, B IS EITHER VERY LARGE OR MODERATE; IN CASE
C           IT IS MODERATE, THE COEFFICIENTS HAVE BEEN SCALED AND B
C           REPRESENTS THE SUM OF THE SCALED INPUT COEFFICENT B0 AND THE
C           DISCRIMINANT; IN CASE B IS VERY LARGE, THE COEFFICIENTS NEED
C           NOT BE SCALED (THE DISCRIMINANT IS ESSENTIALLY B), AND THE
C           ROOTS CAN BE COMPUTED WITH THE SAME QUOTIENTS AS IN THE
C           FORMER CASE.
C
            CALL A02ACF(B(1),B(2),A(1),A(2),ZLG(1),ZLG(2))
            CALL A02ACF(C(1),C(2),B(1),B(2),ZSM(1),ZSM(2))
         ELSE
C
C           C IS ZERO, BUT A IS NOT.
C
            ZSM(1) = ZERO
            ZSM(2) = ZERO
            CALL A02ACF(B(1),B(2),A(1),A(2),ZLG(1),ZLG(2))
         END IF
      ELSE
C
C        A IS ZERO.  INDICATE THAT AT LEAST ONE ROOT HAS OVERFLOWED.
C
         OVFLOW = .TRUE.
         ZLG(1) = FINITY
         ZLG(2) = FINITY
         IF (APPABS(B(1),B(2)).EQ.ZERO .AND. APPABS(C(1),C(2)).NE.ZERO)
     *       THEN
C
C           A AND B ARE ZERO, BUT C IS NOT.  SET THE ROOTS TO INFINITY
C           BUT OF OPPOSITE SIGN TO INDICATE THIS.
C
            ZSM(1) = -ZLG(1)
            ZSM(2) = -ZLG(2)
         ELSE
            IF (APPABS(B(1),B(2)).EQ.ZERO) THEN
C
C              ALL COEFFICIENTS ARE ZERO.  SET BOTH ROOTS TO + INFINITY.
C
               ZSM(1) = ZLG(1)
               ZSM(2) = ZLG(2)
            ELSE
C
C              A IS ZERO, BUT B IS NOT.  COMPUTE THE SMALLER ROOT.
C
               CALL A02ACF(C(1),C(2),B(1),B(2),ZSM(1),ZSM(2))
            END IF
         END IF
      END IF
      RETURN
      END
      SUBROUTINE C02AGV(DX,NDEG,A,P,PPRIME,PDPRIM,ERROR,DEFLAT)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C     BASED ON THE ROUTINE  RPOLYR, WRITTEN BY BRIAN T. SMITH
C
C     THIS SUBROUTINE EVALUATES A POLYNOMIAL OF DEGREE NDEG WITH REAL
C     COEFFICIENTS AT A REAL POINT, ITS FIRST DERIVATIVE AND HALF ITS
C     SECOND DERIVATIVE AT THAT POINT AND AN ESTIMATE OF THE ERROR IN
C     THE POLYNOMIAL VALUE.  THE ESTIMATE IS A GUARANTEED UPPER BOUND,
C     ASSUMING THE COEFFICIENTS ARE EXACT, AND THE FLOATING POINT
C     ARITHMETIC IS 'WELL-BEHAVED'.
C     THAT IS,
C
C           FL(C+D) = (C+D)(1+E),  ABS(E) <= DEPS * (SMALL INTEGER)
C           FL(C*D) = (C*D)(1+F),  ABS(F) <= DEPS * (SMALL INTEGER)
C
C     WHERE FL(.) IS THE MACHINE VALUE FOR THE ENCLOSED OPERATION
C           AND  DEPS  IS A SMALL CONSTANT, USUALLY EQUAL TO THE
C           RELATIVE MACHINE PRECISION.
C
C     THE POLYNOMIAL IS ASSUMED TO BE OF THE FORM
C
C        P(X) = SUM(A(NDEG-I)*X**I, I=0,1,...,NDEG )
C
C     WHERE  A(I) AND  X  ARE DOUBLE PRECISION.
C     NOTE: A(0) IS THE COEFFICIENT OF X**NDEG AND  A(NDEG) IS THE
C           CONSTANT COEFFICIENT.
C
C     USAGE: THIS PROGRAM CAN BE USED TO EVALUATE A REAL POLYNOMIAL
C            AT A REAL POINT AND TO CHECK WHETHER THE POLYNOMIAL VALUE
C            IS ESSENTIALLY ZERO.  THAT IS, IF  ABS(P) <= ERROR, THE
C            POINT X IS INDISTINGUISHABLE FROM A ZERO OF P, BECAUSE OF
C            ROUNDING ERROR INTRODUCED BY THE FLOATING POINT ARITHMETIC.
C
C            AS A SECOND USAGE, THE CIRCLE CENTERED AT X AND OF RADIUS
C            R = NDEG*(ABS(P)+ERROR)/ABS(PPRIME)  IS KNOWN TO CONTAIN A
C            ZERO OF P.  IF IT IS KNOWN THAT THERE IS NO COMPLEX ZERO
C            IN THIS REGION, THEN THERE IS A REAL ZERO IN THE CLOSED
C            INTERVAL  [X-R,X+R].
C
C     HISTORY: THIS CODE WAS TAKEN FROM THE PROGRAM  ZERPOL, A ZERO-
C              FINDING ALGORITHM FOR POLYNOMIALS, UNIV. OF TORONTO,
C              COMPUTER SCIENCE DEPARTMENT, 1965.  THE CODE HAS BEEN
C              MODIFIED BY B. T. SMITH TO BE PORTABLE TO FORTRAN 77
C              SYSTEMS.
C
C     THIS PROGRAM IS DESIGNED TO TAKE ADVANTAGE OF SYSTEMS THAT REPORT
C     OVERFLOW AND UNDERFLOW CONDITIONS IN AN EFFICIENT WAY.  THAT IS,
C     IF, WHENEVER AN OVERFLOW OR UNDERFLOW OCCURS, CERTAIN FLAGS ARE
C     SET (THAT IS, THE LOGICAL VARIABLES OVFLOW AND UNFLOW IN THE
C     COMMON BLOCK AC02AG), C02AGV CAN USE THESE INDICATORS TO INDICATE
C     THAT THE COEFFICIENTS OF THE POLYNOMIAL NEED TO BE SCALED.  SUCH
C     SCALING PERMITS THE DETERMINATION OF THE ROOTS OF THE POLYNOMIAL
C     WITHOUT INTERMEDIATE UNDERFLOW/OVERFLOW CONTAMINATING THE
C     COMPUTED ROOTS.
C
C     HOWEVER, AS IMPLEMENTED IN THE NAG LIBRARY, THE ROUTINE SIMPLY
C     ASSUMES THAT THE MACHINE TERMINATES ON OVERFLOW AND IGNORES
C     UNDERFLOW.
C
C     .. Parameters ..
      DOUBLE PRECISION  SXTEEN, THREE, TWO, TO3RDS, ZERO
      PARAMETER         (SXTEEN=16.0D0,THREE=3.0D0,TWO=2.0D0,
     *                  TO3RDS=TWO/THREE,ZERO=0.0D0)
C     .. Scalar Arguments ..
      DOUBLE PRECISION  DX, ERROR, P, PDPRIM, PPRIME
      INTEGER           NDEG
C     .. Array Arguments ..
      DOUBLE PRECISION  A(0:NDEG), DEFLAT(0:NDEG)
C     .. Scalars in Common ..
      DOUBLE PRECISION  DEPS, FINITY, SQRTFY, SQRTTY, TINY
      INTEGER           EMAXM1, EMINM1, EXPDEP, LRGEXP
      LOGICAL           OVFLOW, UNFLOW
C     .. Local Scalars ..
      DOUBLE PRECISION  ABSP, ABSX, DV, SX, W
      INTEGER           I
C     .. Intrinsic Functions ..
      INTRINSIC         ABS
C     .. Common blocks ..
      COMMON            /AC02AG/OVFLOW, UNFLOW
      COMMON            /BC02AG/FINITY, SQRTFY, SQRTTY, TINY, DEPS,
     *                  EMINM1, EMAXM1, EXPDEP, LRGEXP
C     .. Save statement ..
      SAVE              /AC02AG/, /BC02AG/
C     .. Executable Statements ..
      IF (NDEG.GE.1) THEN
C
C        ASSUME THE POLYNOMIAL IS SCALED SO THAT NO OVER/UNDERFLOW
C        OCCURS.  IF THE LARGEST COEFFICIENT IS LARGER THAN THE
C        SQRT(FINITY), UNDERFLOW CAN PROBABLY BE IGNORED.  HOWEVER,
C        OVERFLOW CANNOT BE IGNORED IN THESE CALCULATIONS IN GENERAL.
C
         SX = DX
         ABSX = ABS(SX)
         W = ZERO
         DV = A(0)
         DEFLAT(0) = DV
         DEFLAT(1) = A(1) + DX*DEFLAT(0)
         DO 20 I = 2, NDEG
            W = DV + SX*W
            DV = DEFLAT(I-1) + DX*DV
            DEFLAT(I) = A(I) + DX*DEFLAT(I-1)
   20    CONTINUE
         P = DEFLAT(NDEG)
C
C        IF OVFLOW, THE COEFFICIENTS MUST BE SCALED DOWN.
C        THIS RESCALING IS TO BE PERFORMED BY THE CALLING PROGRAM.
C        THE NEED TO RESCALE IS INDICATED BY THE FLAG  OVFLOW.
C
C        FOR A MACHINE THAT USES IEEE ARITHMETIC, AN EFFECTIVE BUT
C        INCOMPLETE TEST FOR OVERFLOW IS TO CHECK WHETHER THE RESULT
C        P IS EQUAL TO INFINITY.  THIS MAY WORK ON OTHER MACHINES IF
C        AN OVERFLOW DOES NOT TERMINATE THE COMPUTATION AND IS REPLACED
C        WITH A VALUE THAT PERSISTS THROUGH SUBSEQUENT COMPUTATIONS.
C
         ABSP = ABS(P)
         OVFLOW = ABSP .GE. FINITY
         IF ( .NOT. OVFLOW) THEN
            ERROR = TO3RDS*ABS(A(0))
            DO 40 I = 1, NDEG - 1
               ERROR = ABS(DEFLAT(I)) + ABSX*ERROR
   40       CONTINUE
            ERROR = SXTEEN*DEPS*(ABS(DEFLAT(NDEG))+THREE*ABSX*ERROR)
C
C           CHECK TO SEE IF THE ERROR BOUND COMPUTATION HAS OVERFLOWED.
C           IF SO, MULTIPLY IT BY NDEG TO INSURE IT IS LARGE AND IGNORE
C           THE OVERFLOW CONDITION.
C
C           FOR A MACHINE THAT USES IEEE ARITHMETIC, AN EFFECTIVE BUT
C           INCOMPLETE TEST FOR OVERFLOW IS TO CHECK WHETHER THE RESULT
C           ERROR IS EQUAL TO INFINITY.  THIS MAY WORK ON OTHER MACHINES
C           IF AN OVERFLOW DOES NOT TERMINATE THE COMPUTATION AND IS
C           REPLACED WITH A VALUE THAT PERSISTS THROUGH SUBSEQUENT
C           COMPUTATIONS.
C
            OVFLOW = ERROR .GE. FINITY
            IF (OVFLOW) THEN
               P = FINITY
               PPRIME = FINITY
               PDPRIM = FINITY
               ERROR = TINY
               OVFLOW = .FALSE.
            ELSE
C
C              COMPUTE THE FIRST AND SECOND DERIVATIVES.
C
               PPRIME = DV
               PDPRIM = W
            END IF
         END IF
      ELSE IF (NDEG.EQ.0) THEN
         P = A(0)
         ERROR = ZERO
         PPRIME = ZERO
         PDPRIM = ZERO
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION A02ABF(XXR,XXI)
C     NAG COPYRIGHT 1975
C     MARK 4.5 REVISED
C     MARK 5C REVISED
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     RETURNS THE ABSOLUTE VALUE OF A COMPLEX NUMBER VIA ROUTINE
C     NAME
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION                 XXI, XXR
C     .. Local Scalars ..
      DOUBLE PRECISION                 H, ONE, XI, XR, ZERO
C     .. Intrinsic Functions ..
      INTRINSIC                        ABS, SQRT
C     .. Data statements ..
      DATA                             ZERO/0.0D0/, ONE/1.0D0/
C     .. Executable Statements ..
C
      XR = ABS(XXR)
      XI = ABS(XXI)
      IF (XI.LE.XR) GO TO 20
      H = XR
      XR = XI
      XI = H
   20 IF (XI.NE.ZERO) GO TO 40
      A02ABF = XR
      RETURN
   40 H = XR*SQRT(ONE+(XI/XR)**2)
      A02ABF = H
      RETURN
      END
      SUBROUTINE C02AGW(XR,XI,NDEG,A,P,PPRIME,PDPRIM,ERROR,DEFLAT)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C     BASED ON THE ROUTINE  CPOLYR, WRITTEN BY BRIAN T. SMITH
C
C     THIS SUBROUTINE EVALUATES A POLYNOMIAL OF DEGREE NDEG WITH REAL
C     COEFFICIENTS AT A COMPLEX POINT (XR, XI), ITS FIRST DERIVATIVE
C     AND HALF ITS SECOND DERIVATIVE AT THAT POINT AND AN ESTIMATE OF
C     THE ERROR IN THE POLYNOMIAL VALUE.  THE ESTIMATE IS A GUARANTEED
C     UPPER BOUND, ASSUMING THE COEFFICIENTS ARE EXACT, AND THE FLOATING
C     POINT ARITHMETIC IS 'WELL-BEHAVED'.
C     THAT IS,
C
C           FL(C+D) = (C+D)(1+E),  ABS(E) <= DEPS * (SMALL INTEGER)
C           FL(C*D) = (C*D)(1+F),  ABS(F) <= DEPS * (SMALL INTEGER)
C
C     WHERE FL(.) IS THE MACHINE VALUE FOR THE ENCLOSED OPERATION
C           AND  DEPS  IS A SMALL CONSTANT, USUALLY EQUAL TO THE
C           RELATIVE MACHINE PRECISION.
C
C     THE POLYNOMIAL IS ASSUMED TO BE OF THE FORM
C
C        P(X) = SUM(A(NDEG-I)*X**I, I=0,1,...,NDEG )
C
C     WHERE  A(I)  AND  X  (THE INDETERMINATE) ARE OF TYPE COMPLEX.
C     NOTE: A(0) IS THE COEFFICIENT OF X**NDEG AND  A(NDEG) IS THE
C           CONSTANT COEFFICIENT.
C
C     USAGE: THIS PROGRAM CAN BE USED TO EVALUATE A REAL POLYNOMIAL
C            AT A COMPLEX POINT AND TO CHECK WHETHER THE POLYNOMIAL
C            VALUE IS ESSENTIALLY ZERO.  THAT IS, IF  ABS(P) <= ERROR,
C            THE POINT X IS INDISTINGUISHABLE FROM A ZERO OF P, BECAUSE
C            OF ROUNDING ERROR INTRODUCED BY THE FLOATING POINT
C            ARITHMETIC.
C
C            AS A SECOND USAGE, THE CIRCLE CENTERED AT X AND OF RADIUS
C            R = NDEG*(ABS(P)+ERROR)/ABS(PPRIME)  IS KNOWN TO CONTAIN A
C            ZERO OF P.  IF IT IS KNOWN THAT THERE IS NO COMPLEX ZERO
C            IN THIS REGION, THEN THERE IS A REAL ZERO IN THE CLOSED
C            INTERVAL  [X-R,X+R].
C
C     HISTORY: THIS CODE WAS TAKEN FROM THE PROGRAM  ZERPOL, A ZERO-
C              FINDING ALGORITHM FOR POLYNOMIALS, UNIV. OF TORONTO,
C              COMPUTER SCIENCE DEPARTMENT, 1965.  THE CODE HAS BEEN
C              MODIFIED BY B. T. SMITH TO BE PORTABLE TO FORTRAN 77
C              SYSTEMS.
C
C     THIS PROGRAM IS DESIGNED TO TAKE ADVANTAGE OF SYSTEMS THAT REPORT
C     OVERFLOW AND UNDERFLOW CONDITIONS IN AN EFFICIENT WAY.  THAT IS,
C     IF, WHENEVER AN OVERFLOW OR UNDERFLOW OCCURS, CERTAIN FLAGS ARE
C     SET (THAT IS, THE LOGICAL VARIABLES OVFLOW AND UNFLOW IN THE
C     COMMON BLOCK AC02AG), C02AGW CAN USE THESE INDICATORS TO INDICATE
C     THAT THE COEFFICIENTS OF THE POLYNOMIAL NEED TO BE SCALED. SUCH
C     SCALING PERMITS THE DETERMINATION OF THE ROOTS OF THE POLYNOMIAL
C     WITHOUT INTERMEDIATE UNDERFLOW/OVERFLOW CONTAMINATING THE
C     COMPUTED ROOTS.
C
C     HOWEVER, AS IMPLEMENTED IN THE NAG LIBRARY, THE ROUTINE SIMPLY
C     ASSUMES THAT THE MACHINE TERMINATES ON OVERFLOW AND IGNORES
C     UNDERFLOW.
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO, TWO, NINE, SXTEEN, THREE
      PARAMETER         (ZERO=0.0D0,TWO=2.0D0,NINE=9.0D0,SXTEEN=16.0D0,
     *                  THREE=3.0D0)
C     .. Scalar Arguments ..
      DOUBLE PRECISION  ERROR, XI, XR
      INTEGER           NDEG
C     .. Array Arguments ..
      DOUBLE PRECISION  A(0:NDEG), DEFLAT(0:NDEG), P(2), PDPRIM(2),
     *                  PPRIME(2)
C     .. Scalars in Common ..
      DOUBLE PRECISION  DEPS, FINITY, SQRTFY, SQRTTY, TINY
      INTEGER           EMAXM1, EMINM1, EXPDEP, LRGEXP
      LOGICAL           OVFLOW, UNFLOW
C     .. Local Scalars ..
      DOUBLE PRECISION  ABSP, ABSX, DR, DSC, DT, DT1, DV, DX, DX2, DY,
     *                  S, S1, SC, V, XC1, XC2
      INTEGER           I
C     .. Local Arrays ..
      DOUBLE PRECISION  CX(2)
C     .. External Functions ..
      DOUBLE PRECISION  A02ABF
      EXTERNAL          A02ABF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS
C     .. Common blocks ..
      COMMON            /AC02AG/OVFLOW, UNFLOW
      COMMON            /BC02AG/FINITY, SQRTFY, SQRTTY, TINY, DEPS,
     *                  EMINM1, EMAXM1, EXPDEP, LRGEXP
C     .. Statement Functions ..
      DOUBLE PRECISION  APXABS
C     .. Save statement ..
      SAVE              /AC02AG/, /BC02AG/
C     .. Statement Function definitions ..
      APXABS(XC1,XC2) = ABS(XC1) + ABS(XC2)
C     .. Executable Statements ..
      IF (NDEG.GE.1) THEN
C
C        ASSUME THE POLYNOMIAL IS SCALED SO THAT NO OVER/UNDERFLOW
C        OCCURS.  IF THE LARGEST COEFFICIENT IS LARGER THAN THE
C        SQRT(FINITY), UNDERFLOW CAN PROBABLY BE IGNORED.  HOWEVER,
C        OVERFLOW CANNOT BE IGNORED IN THESE CALCULATIONS IN GENERAL.
C
         S = ZERO
         S1 = ZERO
         DT1 = ZERO
         DT = A(0)
         DEFLAT(0) = A(0)
         DX = XR
         DY = XI
C
C        SC  IS COMPUTED TO CHECK IF SCALING IS NEEDED TO AVOID
C        SPURIOUS OVERFLOW OR UNDERFLOW.  THE SCALE FACTOR NEEDS TO BE
C        AN UPPER BOUND FOR THE MODULUS OF THE ROOT TO AVOID SUCH RANGE
C        EXCEPTIONS.  HOWEVER, THE EXACT MODULUS OF THE ROOT IS NEEDED
C        TO COMPUTE AN ESTIMATE OF THE ERROR IN THE POLYNOMIAL VALUE.
C
         CX(1) = DX
         CX(2) = DY
         ABSX = A02ABF(CX(1),CX(2))
         SC = APXABS(CX(1),CX(2))
C
C        IF  ABS(CMPLX(X,Y)).LE.SQRT(SMALLEST NO.), SCALE UP   X AND Y.
C        IF  ABS(CMPLX(X,Y)).GE.SQRT(LARGEST  NO.), SCALE DOWN X AND Y.
C
         IF (SC.LT.SQRTFY .AND. SC.GE.SQRTTY) THEN
C
C           SCALING OF  DX2  AND  DR  IS UNNECESSARY.
C
            DX2 = DX + DX
            DR = DX*DX + DY*DY
            IF (NDEG.GE.3) THEN
               DEFLAT(1) = A(1) + DX2*A(0)
               DEFLAT(2) = A(2) + (DX2*DEFLAT(1)-DR*A(0))
               DO 20 I = 3, NDEG - 1
                  V = S1*DR
                  S1 = S
                  S = DT1 + (DX2*S-V)
                  DV = DT1*DR
                  DT1 = DT
                  DT = (DX2*DT-DV) + DEFLAT(I-2)
                  DEFLAT(I) = A(I) + (DX2*DEFLAT(I-1)-DR*DEFLAT(I-2))
   20          CONTINUE
               V = S1*DR
               S1 = S
               S = DT1 + (DX*S-V)
               DV = DT1*DR
               DT1 = DT
               DT = (DX*DT-DV) + DEFLAT(NDEG-2)
               DEFLAT(NDEG) = A(NDEG) + (DX*DEFLAT(NDEG-1)
     *                        -DR*DEFLAT(NDEG-2))
            ELSE IF (NDEG.EQ.2) THEN
               DEFLAT(1) = A(1) + DX2*A(0)
               DEFLAT(2) = A(2) + (DX*DEFLAT(1)-DR*A(0))
            ELSE
C
C              NDEG IS 1.
C
               DT = ZERO
               DEFLAT(1) = A(1) + DX*A(0)
            END IF
         ELSE
C
C           SCALE  DX  AND  DY  LEST  DR  OVERFLOWS OR UNDERFLOWS.
C
            DSC = SC
            DX = DX/DSC
            DY = DY/DSC
C
C           DR  CANNOT OVERFLOW, FORTUNATELY.
C
            DR = (DX*DX+DY*DY)*DSC
            DX2 = DX + DX
            IF (NDEG.GE.3) THEN
               DEFLAT(1) = A(1) + (DX2*A(0))*DSC
               DEFLAT(2) = A(2) + (DX2*DEFLAT(1)-DR*A(0))*DSC
               DO 40 I = 3, NDEG - 1
                  V = S1*DR
                  S1 = S
                  S = DT1 + (DX2*S-V)*DSC
                  DV = DT1*DR
                  DT1 = DT
                  DV = DX2*DT - DV
                  DT = DEFLAT(I-2) + DV*DSC
                  DEFLAT(I) = A(I) + (DX2*DEFLAT(I-1)-DR*DEFLAT(I-2))
     *                        *DSC
   40          CONTINUE
               V = S1*DR
               S1 = S
               S = DT1 + (DX*S-V)*DSC
               DV = DT1*DR
               DT1 = DT
               DV = DX*DT - DV
               DT = DEFLAT(NDEG-2) + DV*DSC
               DEFLAT(NDEG) = A(NDEG) + (DX*DEFLAT(NDEG-1)
     *                        -DR*DEFLAT(NDEG-2))*DSC
            ELSE IF (NDEG.EQ.2) THEN
               DEFLAT(1) = A(1) + (DX2*A(0))*DSC
               DEFLAT(2) = A(2) + (DX*DEFLAT(1)-DR*A(0))*DSC
            ELSE
C
C              NDEG IS 1.
C
               DT = ZERO
               DEFLAT(1) = A(1) + (DX*A(0))*DSC
            END IF
            DY = DY*DSC
         END IF
C
C        COMPUTE THE POLYNOMIAL VALUE.
C
         P(1) = DEFLAT(NDEG)
         P(2) = DY*DEFLAT(NDEG-1)
C
C        IF OVFLOW, THE COEFFICIENTS MUST BE SCALED DOWN.
C        THIS RESCALING IS TO BE PERFORMED BY THE CALLING PROGRAM.
C        THE NEED TO RESCALE IS INDICATED BY THE FLAG  OVFLOW.
C
C        FOR A MACHINE THAT USES IEEE ARITHMETIC, AN EFFECTIVE BUT
C        INCOMPLETE TEST FOR OVERFLOW IS TO CHECK WHETHER THE RESULT
C        P IS EQUAL TO INFINITY.  THIS MAY WORK ON OTHER MACHINES IF
C        AN OVERFLOW DOES NOT TERMINATE THE COMPUTATION AND IS REPLACED
C        WITH A VALUE THAT PERSISTS THROUGH SUBSEQUENT COMPUTATIONS.
C
         ABSP = A02ABF(P(1),P(2))
         OVFLOW = ABSP .GE. FINITY
         IF ( .NOT. OVFLOW) THEN
            ERROR = ABS(A(0))
            DO 60 I = 2, NDEG - 1
               ERROR = ABS(DEFLAT(I-1)) + ABSX*ERROR
   60       CONTINUE
            ERROR = SXTEEN*DEPS*((NINE*ERROR*ABSX+THREE*ABS(DEFLAT(NDEG-
     *              1)))*ABSX+ABS(DEFLAT(NDEG)))
C
C           CHECK TO SEE IF THE ERROR BOUND COMPUTATION HAS OVERFLOWED.
C           IF SO, MULTIPLY IT BY NDEG TO INSURE IT IS LARGE AND IGNORE
C           THE OVERFLOW CONDITION.
C
C           FOR A MACHINE THAT USES IEEE ARITHMETIC, AN EFFECTIVE BUT
C           INCOMPLETE TEST FOR OVERFLOW IS TO CHECK WHETHER THE RESULT
C           ERROR IS EQUAL TO INFINITY.  THIS MAY WORK ON OTHER MACHINES
C           IF AN OVERFLOW DOES NOT TERMINATE THE COMPUTATION AND IS
C           REPLACED WITH A VALUE THAT PERSISTS THROUGH SUBSEQUENT
C           COMPUTATIONS.
C
            OVFLOW = ERROR .GE. FINITY
            IF (OVFLOW) THEN
               ERROR = NDEG*ERROR
               OVFLOW = .FALSE.
            END IF
C
C           COMPUTE THE FIRST AND SECOND DERIVATIVES.  IF OVERFLOW
C           OCCURS IN THIS COMPUTATION, THE FLAG OVFLOW WILL BE SET.
C
            DV = TWO*DY
            PPRIME(1) = DEFLAT(NDEG-1) - DV*DT1*DY
            PPRIME(2) = DV*DT
            PDPRIM(1) = DT - DV*(DV*S)
            PDPRIM(2) = DY*(THREE*DT1-DV*(DV*S1))
         END IF
      ELSE IF (NDEG.EQ.0) THEN
         ERROR = ZERO
         P(1) = A(0)
         P(2) = ZERO
         PPRIME(1) = ZERO
         PPRIME(2) = ZERO
         PDPRIM(1) = ZERO
         PDPRIM(2) = ZERO
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION C02AGR(X,EXP)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C     BASED ON THE ROUTINE  SETEXP, WRITTEN BY BRIAN T. SMITH
C
C     THIS FUNCTION COMPUTES THE VALUE  FRACTION(X) * B ** EXP
C     WHERE  B  IS THE BASE OF ENTITIES OF TYPE DOUBLE PRECISION.
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION                 X
      INTEGER                          EXP
C     .. External Functions ..
      DOUBLE PRECISION                 C02AGY
      INTEGER                          C02AGX
      EXTERNAL                         C02AGY, C02AGX
C     .. Executable Statements ..
      C02AGR = C02AGY(X,EXP-C02AGX(X))
      RETURN
      END
      SUBROUTINE A02ACF(XXR,XXI,YYR,YYI,ZR,ZI)
C     MARK 2A RELEASE.  NAG COPYRIGHT 1973
C     MARK 4.5 REVISED
C     MARK 5C REVISED
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C
C     DIVIDES ONE COMPLEX NUMBER BY A SECOND
C
C     .. Scalar Arguments ..
      DOUBLE PRECISION  XXI, XXR, YYI, YYR, ZI, ZR
C     .. Local Scalars ..
      DOUBLE PRECISION  A, H, ONE
C     .. Intrinsic Functions ..
      INTRINSIC         ABS
C     .. Data statements ..
      DATA              ONE/1.0D0/
C     .. Executable Statements ..
C
      IF (ABS(YYR).LE.ABS(YYI)) GO TO 20
      H = YYI/YYR
      A = ONE/(H*YYI+YYR)
      ZR = (XXR+H*XXI)*A
      ZI = (XXI-H*XXR)*A
      RETURN
   20 H = YYR/YYI
      A = ONE/(H*YYR+YYI)
      ZR = (H*XXR+XXI)*A
      ZI = (H*XXI-XXR)*A
      RETURN
      END
      SUBROUTINE C02AGT(A0,B0,C0,ZSM,ZLG)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C     BASED ON THE ROUTINE  QDRTC, WRITTEN BY BRIAN T. SMITH
C
C     THIS SUBROUTINE DETERMINES THE ROOTS OF THE QUADRATIC EQUATION
C
C         A0*X**2 + B0*X + C0
C
C     WHERE A0, B0, AND C0 ARE REAL COEFFICIENTS, AND ZSM AND ZLG
C     THE SMALLEST AND LARGEST ROOT IN MAGNITUDE, RESPECTIVELY.
C
C     THE ROOTS ARE COMPUTED TO WITHIN A RELATIVE ERROR OF A FEW
C     UNITS IN THE LAST PLACE (DEPENDING ON THE ACCURACY OF THE
C     BASIC ARITHMETIC OPERATIONS) EXCEPT WHEN UNDERFLOW OR OVERFLOW
C     OCCURS IN WHICH CASE THE TRUE ROOTS ARE WITHIN A FEW UNITS IN
C     THE LAST PLACE OF THE UNDERFLOW OR OVERFLOW THRESHOLD.
C
C     IF THE LEADING COEFFICIENT IS ZERO, THE LARGER ROOT IS
C     SET TO THE LARGEST MACHINE REPRESENTABLE NUMBER AND THE
C     OVERFLOW FLAG OVFLOW IS SET TRUE.  IF ALL THREE COEFFICIENTS ARE
C     ZERO, THE OVERFLOW FLAG IS SET, BOTH ROOTS ARE SET TO THE LARGEST
C     REPRESENTABLE NUMBER, BUT NO DIVIDE CHECK IS CREATED.
C
C     THIS PROGRAM IS DESIGNED TO TAKE ADVANTAGE OF SYSTEMS THAT REPORT
C     OVERFLOW AND UNDERFLOW CONDITIONS IN AN EFFICIENT WAY.  THAT IS,
C     IF, WHENEVER AN OVERFLOW OR UNDERFLOW OCCURS, CERTAIN FLAGS ARE
C     SET (THAT IS, THE LOGICAL VARIABLES OVFLOW AND UNFLOW IN THE
C     COMMON BLOCK AC02AG), C02AGT CAN USE THESE INDICATORS TO INDICATE
C     THAT THE ROOTS OVERFLOW OR UNDERFLOW AND CANNOT BE REPRESENTED.
C
C     HOWEVER, AS IMPLEMENTED IN THE NAG LIBRARY, THE ROUTINE SIMPLY
C     ASSUMES THAT THE MACHINE TERMINATES ON OVERFLOW AND IGNORES
C     UNDERFLOW.
C
C     C02AGX -- DETERMINE THE EXPONENT OF A NUMBER IN TERMS OF THE
C               MODEL.
C     C02AGR -- FORM A NUMBER WITH A GIVEN MANTISSA AND EXPONENT
C               PRECISION.
C
C     .. Parameters ..
      DOUBLE PRECISION  HALF, ONE, ZERO
      PARAMETER         (HALF=0.5D0,ONE=1.0D0,ZERO=0.0D0)
C     .. Scalar Arguments ..
      DOUBLE PRECISION  A0, B0, C0
C     .. Array Arguments ..
      DOUBLE PRECISION  ZLG(2), ZSM(2)
C     .. Scalars in Common ..
      DOUBLE PRECISION  DEPS, FINITY, SQRTFY, SQRTTY, TINY
      INTEGER           EMAXM1, EMINM1, EXPDEP, LRGEXP
      LOGICAL           OVFLOW, UNFLOW
C     .. Local Scalars ..
      DOUBLE PRECISION  A, B, C, D, SC, SQRTD
      INTEGER           EXPBSQ, SCLEXP
C     .. External Functions ..
      DOUBLE PRECISION  C02AGR
      INTEGER           C02AGX
      EXTERNAL          C02AGR, C02AGX
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MAX, MIN, SIGN, SQRT
C     .. Common blocks ..
      COMMON            /AC02AG/OVFLOW, UNFLOW
      COMMON            /BC02AG/FINITY, SQRTFY, SQRTTY, TINY, DEPS,
     *                  EMINM1, EMAXM1, EXPDEP, LRGEXP
C     .. Save statement ..
      SAVE              /AC02AG/, /BC02AG/
C     .. Executable Statements ..
C
C     INITIALIZE LOCAL VARIABLES WITH THE INPUT COEFFICIENTS.
C
      A = A0
      B = -B0
      C = C0
C
C     CHECK FOR  A = ZERO  OR  C = ZERO.
C
      IF (A.NE.ZERO) THEN
         IF (C.NE.ZERO) THEN
C
C           AT THIS POINT, A AND C ARE NON-ZERO.
C
C           SCALE THE COEFFICIENTS SO THAT THE PRODUCT A * C IS NEAR
C           1.0D0 IN MAGNITUDE.  THIS AVOIDS SPURIOUS UNDERFLOW/OVERFLOW
C           CONDITIONS WHEN THE TRUE RESULTS ARE WITHIN RANGE.
C
C           THE SCALE FACTOR IS A POWER OF THE BASE NEAR TO
C           SQRT(ABS(A*C)).  THIS CHOICE AVOIDS UNNECESSARY ROUNDING
C           ERRORS BUT IS EXPENSIVE TO COMPUTE WHEN FLOATING POINT
C           MANIPULATIVE FUNCTIONS ARE NOT AVAILABLE IN MACHINE CODE.
C
            SCLEXP = (C02AGX(A)+C02AGX(C))/2
C
C           THE SCALE FACTOR IS  BASE ** SCLEXP.  IF A AND C ARE SCALED
C           USING THIS SCALE FACTOR AS A DIVIDEND, THEN THE
C           THE SCALED PRODUCT A'*C' IS BETWEEN BASE**(-2) AND
C           BASE IN MAGNITUDE, WHERE BASE IS THE BASE FOR MODEL NUMBERS
C           OF THE TYPE OF A.
C
C           BUT BEFORE PERFORMING THE SCALING, CHECK TO SEE IF IT IS
C           NECESSARY -- THAT IS, IF B IS SO LARGE IN MAGNITUDE THAT
C           B**2 EXCEEDS ABS(4*A*C) BY MORE THAN THE RELATIVE MACHINE
C           PRECISION FOR THE DOUBLE PRECISION DATA TYPE,
C           THE DISCRIMINANT IS IN EFFECT B AND NO SCALING IS REQUIRED.
C           HOWEVER, IF B IS SO SMALL IN MAGNITUDE THAT ABS(4*A*C)
C           EXCEEDS B**2 IN MAGNITUDE BY MORE THAN THIS SAME RELATIVE
C           MACHINE PRECISION, B IS IN EFFECT ZERO, BUT A AND C ARE
C           STILL SCALED TO AVOID SPURIOUS UNDERFLOWS/OVERFLOWS.
C
C           COMPUTE THE EXPONENT OF THE SQUARE OF THE SCALED B.
C
            IF (ABS(B).NE.ZERO) THEN
               EXPBSQ = 2*(C02AGX(B)-SCLEXP)
            ELSE
               EXPBSQ = -2*EXPDEP
            END IF
C
C           CHECK IF B**2 IS TOO BIG.
C
            IF (EXPBSQ.LE.EXPDEP) THEN
C
C              B**2 IS NOT TOO BIG.  SCALING WILL BE PERFORMED.
C
C              A AND C SHOULD BE SCALED USING THE USUAL SCALE
C              MANIPULATION FUNCTION BUT FOR EFFICIENCY, THE
C              SCALING IS PERFORMED BY DIVISION.
C
               SCLEXP = MIN(SCLEXP+1,EMAXM1)
               SCLEXP = MAX(SCLEXP,EMINM1)
               SC = C02AGR(ONE,SCLEXP)
C
C              CHECK IF IT IS TOO SMALL.
C
               IF (EXPBSQ.LT.-EXPDEP) THEN
C
C                 B IS TOO SMALL.  SET IT TO ZERO.
C
                  B = ZERO
               ELSE
C
C                 B IS NEITHER TOO LARGE NOR TOO SMALL.  SCALE IT.
C
                  B = (B/SC)*HALF
               END IF
               A = A/SC
               C = C/SC
               D = B*B - A*C
               SQRTD = SQRT(ABS(D))
               IF (D.LE.ZERO) THEN
C
C                 THE ROOTS ARE COMPLEX.
C
                  ZLG(1) = B/A
                  ZLG(2) = ABS(SQRTD/A)
                  ZSM(1) = ZLG(1)
                  ZSM(2) = -ZLG(2)
               ELSE
C
C                 THE ROOTS ARE REAL AND SQRTD IS NOT ZERO.
C
                  B = SIGN(SQRTD,B) + B
                  ZSM(1) = C/B
                  ZSM(2) = ZERO
                  ZLG(1) = B/A
                  ZLG(2) = ZERO
C
C                 BECAUSE OF ROUNDING ERRORS IN THE SQUARE ROOT AND
C                 DIVISIONS ABOVE (PARTICULARLY ON MACHINES THAT
C                 TRUNCATE AND ONLY WHEN B IS SMALL), THE REAL ROOTS MAY
C                 BE IMPROPERLY ORDERED -- SET THEM SO THAT THE SMALLER
C                 ONE IS OPPOSITE IN SIGN TO THE LARGER ONE.
C
                  IF (ABS(ZLG(1)).LT.ABS(ZSM(1))) THEN
                     ZSM(1) = -ZLG(1)
                     ZSM(2) = -ZLG(2)
                  END IF
               END IF
            ELSE
C
C              AT THIS POINT, B IS VERY LARGE; IN THIS CASE, THE
C              COEFFICIENTS NEED NOT BE SCALED AS THE DISCRIMINANT
C              IS ESSENTIALLY B.
C
               ZSM(1) = C/B
               ZSM(2) = ZERO
               ZLG(1) = B/A
               ZLG(2) = ZERO
C
C              BECAUSE OF ROUNDING ERRORS IN THE SQUARE ROOT AND
C              DIVISIONS ABOVE (PARTICULARLY ON MACHINES THAT TRUNCATE
C              AND ONLY WHEN B IS SMALL), THE REAL ROOTS MAY BE
C              IMPROPERLY ORDERED -- SET THEM SO THAT THE SMALLER ONE
C              IS OPPOSITE IN SIGN TO THE LARGER ONE.
C
               IF (ABS(ZLG(1)).LT.ABS(ZSM(1))) THEN
                  ZSM(1) = -ZLG(1)
                  ZSM(2) = -ZLG(2)
               END IF
            END IF
         ELSE
C
C           C IS ZERO, BUT A IS NOT.
C
            ZSM(1) = ZERO
            ZSM(2) = ZERO
            ZLG(1) = B/A
            ZLG(2) = ZERO
         END IF
      ELSE
C
C        A IS ZERO.  INDICATE THAT AT LEAST ONE ROOT HAS OVERFLOWED.
C
         OVFLOW = .TRUE.
         ZLG(1) = FINITY
         ZLG(2) = ZERO
         IF (B.EQ.ZERO .AND. C.NE.ZERO) THEN
C
C           A AND B ARE ZERO, BUT C IS NOT.  SET THE ROOTS TO INFINITY
C           BUT OF OPPOSITE SIGN TO INDICATE THIS.
C
            ZSM(1) = -ZLG(1)
            ZSM(2) = -ZLG(2)
         ELSE
            IF (B.EQ.ZERO) THEN
C
C              ALL COEFFICIENTS ARE ZERO.  SET BOTH ROOTS TO + INFINITY.
C
               ZSM(1) = ZLG(1)
               ZSM(2) = ZLG(2)
            ELSE
C
C              A IS ZERO, BUT B IS NOT.  COMPUTE THE SMALLER ROOT.
C
               ZSM(1) = C/B
               ZSM(2) = ZERO
            END IF
         END IF
      END IF
      RETURN
      END
      DOUBLE PRECISION FUNCTION F06BLF( A, B, FAIL )
C     MARK 12 RELEASE. NAG COPYRIGHT 1986.
C     .. Scalar Arguments ..
      DOUBLE PRECISION                  A, B
      LOGICAL                           FAIL
C     ..
C
C  F06BLF returns the value div given by
C
C     div = ( a/b                 if a/b does not overflow,
C           (
C           ( 0.0                 if a .eq. 0.0,
C           (
C           ( sign( a/b )*flmax   if a .ne. 0.0  and a/b would overflow,
C
C  where  flmax  is a large value, via the function name. In addition if
C  a/b would overflow then  fail is returned as true, otherwise  fail is
C  returned as false.
C
C  Note that when  a and b  are both zero, fail is returned as true, but
C  div  is returned as  0.0. In all other cases of overflow  div is such
C  that  abs( div ) = flmax.
C
C  When  b = 0  then  sign( a/b )  is taken as  sign( a ).
C
C  Nag Fortran 77 O( 1 ) basic linear algebra routine.
C
C  -- Written on 26-October-1982.
C     Sven Hammarling, Nag Central Office.
C
C
C     .. Parameters ..
      DOUBLE PRECISION      ONE         , ZERO
      PARAMETER           ( ONE = 1.0D+0, ZERO = 0.0D+0 )
C     .. Local Scalars ..
      DOUBLE PRECISION      ABSB, DIV, FLMAX, FLMIN
      LOGICAL               FIRST
C     .. External Functions ..
      DOUBLE PRECISION      X02AMF
      EXTERNAL              X02AMF
C     .. Intrinsic Functions ..
      INTRINSIC             ABS, SIGN
C     .. Save statement ..
      SAVE                  FIRST, FLMIN, FLMAX
C     .. Data statements ..
      DATA                  FIRST/ .TRUE. /
C     ..
C     .. Executable Statements ..
      IF( A.EQ.ZERO )THEN
         DIV = ZERO
         IF( B.EQ.ZERO )THEN
            FAIL = .TRUE.
         ELSE
            FAIL = .FALSE.
         END IF
      ELSE
C
         IF( FIRST )THEN
            FIRST = .FALSE.
            FLMIN =  X02AMF( )
            FLMAX =  1/FLMIN
         END IF
C
         IF( B.EQ.ZERO )THEN
            DIV  =  SIGN( FLMAX, A )
            FAIL = .TRUE.
         ELSE
            ABSB = ABS( B )
            IF( ABSB.GE.ONE )THEN
               FAIL = .FALSE.
               IF( ABS( A ).GE.ABSB*FLMIN )THEN
                  DIV = A/B
               ELSE
                  DIV = ZERO
               END IF
            ELSE
               IF( ABS( A ).LE.ABSB*FLMAX )THEN
                  FAIL = .FALSE.
                  DIV  =  A/B
               ELSE
                  FAIL = .TRUE.
                  DIV  = FLMAX
                  IF( ( ( A.LT.ZERO ).AND.( B.GT.ZERO ) ).OR.
     $                ( ( A.GT.ZERO ).AND.( B.LT.ZERO ) )     )
     $               DIV = -DIV
               END IF
            END IF
         END IF
      END IF
C
      F06BLF = DIV
      RETURN
C
C     End of F06BLF. ( SDIV )
C
      END
      LOGICAL FUNCTION C02AGS(X)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C     BASED ON THE ROUTINE  ISUNRM, WRITTEN BY BRIAN T. SMITH
C
C     THIS FUNCTION RETURNS TRUE IF X IS UNNORMALIZED OR ZERO.
C
C     .. Parameters ..
      DOUBLE PRECISION        ZERO
      PARAMETER               (ZERO=0.0D0)
C     .. Scalar Arguments ..
      DOUBLE PRECISION        X
C     .. Executable Statements ..
      C02AGS = (X+ZERO) .EQ. ZERO
      RETURN
      END
      DOUBLE PRECISION FUNCTION C02AGY(DX,EXP)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C     MARK 15 REVISED. IER-893 (APR 1991).
C     BASED ON THE ROUTINE  DSCALE, WRITTEN BY BRIAN T. SMITH
C
C     THIS FUNCTION COMPUTES THE SCALED PRODUCT  DX * B ** EXP
C     WHERE  B  IS THE BASE OF ENTITIES OF TYPE  DOUBLE PRECISION.
C
C     .. Parameters ..
      DOUBLE PRECISION                 ONE
      PARAMETER                        (ONE=1.0D0)
C     .. Scalar Arguments ..
      DOUBLE PRECISION                 DX
      INTEGER                          EXP
C     .. Scalars in Common ..
      DOUBLE PRECISION                 DEPS, DPNEWL, DPNEWU, FACT, TEMP
      INTEGER                          DBASE, MNEXP, MXEXP, NEWL, NEWU
C     .. Local Scalars ..
      DOUBLE PRECISION                 DPE, DSC, POWER
      INTEGER                          E
C     .. Intrinsic Functions ..
      INTRINSIC                        MOD
C     .. Common blocks ..
      COMMON                           /CC02AG/DPNEWL, DPNEWU, DEPS,
     *                                 TEMP, FACT, DBASE, MNEXP, MXEXP,
     *                                 NEWL, NEWU
C     .. Save statement ..
      SAVE                             /CC02AG/
C     .. Executable Statements ..
      E = EXP
      DSC = DX
C
C     IF THE EXPONENT SCALING IS OUT OF RANGE FOR THE PRECOMPUTED
C     POWERS OF THE BASE, SCALE REPETITIVELY BY THE LARGEST
C     PRECOMPUTED POWER UNTIL E IS WITHIN RANGE.
C
C     CHECK FOR E TOO LARGE.
C
   20 IF (E.GT.NEWU) THEN
         DSC = DSC*FACT
         E = E - NEWU
         GO TO 20
      END IF
C
C     CHECK FOR E TOO SMALL.
C
   40 IF (E.LT.NEWL) THEN
         DSC = DSC*DPNEWL
         E = E - NEWL
         GO TO 40
      END IF
C
C     SCALE BY THE REMAINING SCALING FACTOR.
C     SET DPE = DBASE**E.
C
      IF (E.EQ.0) THEN
         DPE = ONE
      ELSE
         IF (E.LT.0) THEN
            E = -E
            POWER = ONE/DBASE
         ELSE
            POWER = DBASE
         END IF
         DPE = ONE
   60    IF (MOD(E,2).EQ.1) DPE = DPE*POWER
         E = E/2
         IF (E.GT.0) THEN
            POWER = POWER*POWER
            GO TO 60
         END IF
      END IF
      C02AGY = DSC*DPE
      RETURN
      END
      INTEGER FUNCTION X02BHF()
C     MARK 12 RELEASE. NAG COPYRIGHT 1986.
C
C     RETURNS THE MODEL PARAMETER, B.
C
C     .. Executable Statements ..
      X02BHF =     2
      RETURN
      END
      INTEGER FUNCTION C02AGX(DX)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C     MARK 15 REVISED. IER-892 (APR 1991).
C     BASED ON THE ROUTINE  DEXPNT, WRITTEN BY BRIAN T. SMITH
C
C     THIS FUNCTION COMPUTES THE EXPONENT E WHERE X IS REPRESENTED
C     AS  0.0 OR S * F * B**E  WHERE  S  IS A SIGN (PLUS OR MINUS ONE),
C     F  IS A FRACTION, EITHER ZERO OR SATISFIES  1/DBASE <= F < 1,
C     AND  E  SATISFIES  X02BKF( X ) <= E <= X02BLF( X ).
C
C     .. Parameters ..
      DOUBLE PRECISION        ONE, ZERO
      PARAMETER               (ONE=1.0D0,ZERO=0.0D0)
C     .. Scalar Arguments ..
      DOUBLE PRECISION        DX
C     .. Scalars in Common ..
      DOUBLE PRECISION        DEPS, DPNEWL, DPNEWU, FACT, TEMP
      INTEGER                 DBASE, MNEXP, MXEXP, NEWL, NEWU
C     .. Local Scalars ..
      DOUBLE PRECISION        A, ABSX, P
      INTEGER                 E, LP, UP
      LOGICAL                 FIRST
C     .. External Functions ..
      DOUBLE PRECISION        C02AGY, X02AJF, X02AKF, X02ALF
      INTEGER                 X02BHF, X02BKF, X02BLF
      EXTERNAL                C02AGY, X02AJF, X02AKF, X02ALF, X02BHF,
     *                        X02BKF, X02BLF
C     .. Intrinsic Functions ..
      INTRINSIC               ABS, INT, LOG, DBLE
C     .. Common blocks ..
      COMMON                  /CC02AG/DPNEWL, DPNEWU, DEPS, TEMP, FACT,
     *                        DBASE, MNEXP, MXEXP, NEWL, NEWU
C     .. Save statement ..
      SAVE                    /CC02AG/, FIRST
C     .. Data statements ..
      DATA                    FIRST/.TRUE./
C     .. Executable Statements ..
      IF (FIRST) THEN
         FIRST = .FALSE.
         DBASE = X02BHF()
         DEPS = X02AJF()
         MNEXP = X02BKF()
         MXEXP = X02BLF()
         DPNEWL = X02AKF()
         DPNEWU = X02ALF()
         NEWU = MXEXP - 1
         NEWL = MNEXP - 1
         TEMP = DBLE(DBASE)*(ONE-DEPS)
         FACT = DPNEWU/TEMP
      END IF
      IF (DX.NE.ZERO) THEN
C
C        DX IS IN THE RANGE
C
C        DBASE**(MNEXP-1) <= ABS(DX) < DBASE**MXEXP.
C
         ABSX = ABS(DX)
         IF (ABSX.LT.DPNEWU) THEN
            P = (LOG(ABSX)/LOG(DBLE(DBASE))) + ONE
            LP = INT(P)
            UP = LP - 1
            IF (LP.GE.1) THEN
               A = ABSX/C02AGY(ONE,UP)
               IF (A.GE.ONE) THEN
                  E = LP
               ELSE
                  E = UP
               END IF
            ELSE
               A = ABSX/C02AGY(ONE,LP)
               IF (ONE/DBLE(DBASE).GT.A) THEN
                  E = UP
               ELSE
                  E = LP
               END IF
            END IF
         ELSE
            E = MXEXP
         END IF
      ELSE
         E = 0
      END IF
      C02AGX = E
      RETURN
      END
      INTEGER FUNCTION X02BLF()
C     MARK 12 RELEASE. NAG COPYRIGHT 1986.
C
C     RETURNS THE MODEL PARAMETER, EMAX.
C
C     .. Executable Statements ..
      X02BLF =  1024
      RETURN
      END
      INTEGER FUNCTION X02BKF()
C     MARK 12 RELEASE. NAG COPYRIGHT 1986.
C
C     RETURNS THE MODEL PARAMETER, EMIN.
C
C     .. Executable Statements ..
      X02BKF =  -1021
      RETURN
      END
      DOUBLE PRECISION FUNCTION X02ALF()
C     MARK 12 RELEASE. NAG COPYRIGHT 1986.
C
C     RETURNS  (1 - B**(-P)) * B**EMAX  (THE LARGEST POSITIVE MODEL
C     NUMBER)
C
      DOUBLE PRECISION X02CON
      DATA X02CON /1.79769313486231D+308 /
C     .. Executable Statements ..
      X02ALF = X02CON
      RETURN
      END
      INTEGER FUNCTION X02BJF()
C     MARK 12 RELEASE. NAG COPYRIGHT 1986.
C
C     RETURNS THE MODEL PARAMETER, p.
C
C     .. Executable Statements ..
      X02BJF =    53
      RETURN
      END
      DOUBLE PRECISION FUNCTION X02AKF()
C     MARK 12 RELEASE. NAG COPYRIGHT 1986.
C
C     RETURNS  B**(EMIN-1)  (THE SMALLEST POSITIVE MODEL NUMBER)
C
      DOUBLE PRECISION X02CON
      DATA X02CON /2.22507385850721D-308 /
C     .. Executable Statements ..
      X02AKF = X02CON
      RETURN
      END
      SUBROUTINE P01ABZ
C     MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986.
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
      SUBROUTINE X04BAF(NOUT,REC)
C     MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986.
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
      SUBROUTINE X04AAF(I,NERR)
C     MARK 7 RELEASE. NAG COPYRIGHT 1978
C     MARK 7C REVISED IER-190 (MAY 1979)
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     MARK 14 REVISED. IER-829 (DEC 1989).
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/0/
C     .. Executable Statements ..
      IF (I.EQ.0) NERR = NERR1
      IF (I.EQ.1) NERR1 = NERR
      RETURN
      END
      INTEGER FUNCTION P01ABF(IFAIL,IERROR,SRNAME,NREC,REC)
C     MARK 11.5(F77) RELEASE. NAG COPYRIGHT 1986.
C     MARK 13 REVISED. IER-621 (APR 1988).
C     MARK 13B REVISED. IER-668 (AUG 1988).
C
C     P01ABF is the error-handling routine for the NAG Library.
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     For compatibility with certain routines included before Mark 12
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     S.Hammarling, M.P.Hooper and J.J.du Croz, NAG Central Office.
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 NAG Library routine ',A,': IFAIL',
     *  ' =',I6)
      END
      SUBROUTINE C02AGZ(A,NDEG,SCALE,Z,DU,DEFLAT,IER)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C     MARK 14 REVISED. IER-709 (DEC 1989).
C     BASED ON THE ROUTINE  ZERPOL, WRITTEN BY BRIAN T. SMITH
C
C     THIS SUBROUTINE COMPUTES THE N ZEROS OF THE REAL POLYNOMIAL
C
C            A(0)*Z**N + A(1)*Z**(N-1) + ... + A(N-1)*Z + A(N) = 0
C
C     GIVEN BY CMPLX(Z(1,J),Z(2,J)), WHERE J = 1,2,...,N.
C
C     GAMA, THETA AND PHI ARE ARBITRARY PARAMETERS WHICH FOR THIS
C     IMPLEMENTATION HAVE BEEN SET TO 1.0, 2.0 AND 0.2 RESPECTIVELY.
C     THERE IS NO INHERENT LIMITATION ON THE DEGREE OTHER THAN AS THE
C     DEGREE OF A POLYNOMIAL INCREASES, ITS ROOTS BECOME ILL-CONDITIONED
C     AS A FUNCTION OF THE COEFFICIENTS.
C
C     THIS PROGRAM IS DESIGNED TO TAKE ADVANTAGE OF SYSTEMS THAT REPORT
C     OVERFLOW AND UNDERFLOW CONDITIONS IN AN EFFICIENT WAY.  THAT IS,
C     IF, WHENEVER AN OVERFLOW OR UNDERFLOW OCCURS, CERTAIN FLAGS ARE
C     SET (THAT IS, THE LOGICAL VARIABLES OVFLOW AND UNFLOW IN THE
C     COMMON BLOCK AC02AG), C02AGZ CAN USE THESE INDICATORS TO OPTIMALLY
C     SCALE THE COEFFICIENTS OF THE POLYNOMIAL.  THE OPTIMAL SCALING
C     PERMITS THE DETERMINATION OF THE ROOTS OF THE POLYNOMIAL WITHOUT
C     INTERMEDIATE UNDERFLOW/OVERFLOW CONTAMINATING THE COMPUTED ROOTS.
C
C     HOWEVER, AS IMPLEMENTED IN THE NAG LIBRARY, THE ROUTINE SIMPLY
C     ASSUMES THAT THE MACHINE TERMINATES ON OVERFLOW AND IGNORES
C     UNDERFLOW.
C
C     X02BJF -- NUMBER OF DIGITS IN THE MANTISSA OF MODEL NUMBERS OF
C               TYPE DOUBLE PRECISION.
C     X02AJF -- RELATIVE MACHINE PRECISION FOR ENTITIES OF TYPE
C               DOUBLE PRECISION.
C     X02ALF -- LARGEST POSITIVE MACHINE REPRESENTABLE NUMBER OF TYPE
C               DOUBLE PRECISION.
C     X02BLF -- MAXIMUM EXPONENT OF ENTITIES OF TYPE DOUBLE PRECISION.
C     X02AKF -- SMALLEST POSITIVE MACHINE REPRESENTABLE NUMBER OF TYPE
C               DOUBLE PRECISION.
C     C02AGY -- SCALE THE FIRST ARGUMENT BY A VALUE WITH AN EXPONENT
C               EQUAL TO THE SECOND ARGUMENT.
C     C02AGX -- DETERMINE THE EXPONENT OF A NUMBER IN TERMS OF THE
C               MODEL.
C     C02AGW -- COMPUTE THE POLYNOMIAL VALUE, FIRST DERIVATIVE, AND
C               SECOND DERIVATIVE AT A COMPLEX POINT.
C     C02AGU -- DETERMINE THE ROOTS OF A QUADRATIC EQUATION WITH COMPLEX
C               COEFFICIENTS.
C     C02AGT -- DETERMINE THE ROOTS OF A QUADRATIC EQUATION WITH REAL
C               COEFFICIENTS.
C     C02AGV -- COMPUTE THE POLYNOMIAL VALUE, FIRST DERIVATIVE, AND
C               SECOND DERIVATIVE AT A REAL POINT.
C     C02AGS -- RETURNS TRUE IF THE ARGUMENT IS UNNORMALIZED OR ZERO.
C
C     .. Parameters ..
      DOUBLE PRECISION  GAMA, THETA, PHI
      PARAMETER         (GAMA=1.0D0,THETA=2.0D0,PHI=0.2D0)
      DOUBLE PRECISION  HALF, ONE, SMALL, BIGONE, SMLONE, RCONST,
     *                  ONEPQT, ZERO, TWO
      PARAMETER         (HALF=0.5D0,ONE=1.0D0,SMALL=1.0D-3,
     *                  BIGONE=1.0001D0,SMLONE=0.99999D0,RCONST=1.445D0,
     *                  ONEPQT=1.25D0,ZERO=0.0D0,TWO=2.0D0)
C     .. Scalar Arguments ..
      INTEGER           IER, NDEG
      LOGICAL           SCALE
C     .. Array Arguments ..
      DOUBLE PRECISION  A(0:NDEG), DEFLAT(0:NDEG), DU(0:NDEG), Z(2,NDEG)
C     .. Scalars in Common ..
      DOUBLE PRECISION  DEPS, FINITY, SQRTFY, SQRTTY, TINY
      INTEGER           EMAXM1, EMINM1, EXPDEP, LRGEXP
      LOGICAL           OVFLOW, UNFLOW
C     .. Local Scalars ..
      DOUBLE PRECISION  ABDIR, ABDIRO, ABSCL, DX, DZ0I, DZ0R, DZNI,
     *                  DZNR, E, F, F0, FEJER, FN, G, LOWERB, MXCOEF, R,
     *                  RATIO, RTN, S, T, UPPERB, V, W, X2N, X2N1, XC1,
     *                  XC2, XN, XN1, XN2, XN2N
      INTEGER           I, IERS, IHALF, ISPIR, ITER, K, MXCFEX, N, NERR,
     *                  SCBYEX
      LOGICAL           CAUCHY, CONTIN, OVF, SAVO, SAVU, SPIRAL, STARTD,
     *                  UNF
C     .. Local Arrays ..
      DOUBLE PRECISION  C(2), CDIR(2), CDIRO(2), CF(2), CF1(2), CF2(2),
     *                  CL(2), CSPIR(2), CTEMP(2)
      CHARACTER*80      REC(2)
C     .. External Functions ..
      DOUBLE PRECISION  A02ABF, C02AGY, F06BLF, X02AJF, X02AKF, X02ALF
      INTEGER           C02AGX, X02BJF, X02BKF, X02BLF
      LOGICAL           C02AGS
      EXTERNAL          A02ABF, C02AGY, F06BLF, X02AJF, X02AKF, X02ALF,
     *                  C02AGX, X02BJF, X02BKF, X02BLF, C02AGS
C     .. External Subroutines ..
      EXTERNAL          A02ACF, C02AGT, C02AGU, C02AGV, C02AGW, X04AAF,
     *                  X04BAF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, DBLE, EXP, LOG, MAX, MIN, SQRT
C     .. Common blocks ..
      COMMON            /AC02AG/OVFLOW, UNFLOW
      COMMON            /BC02AG/FINITY, SQRTFY, SQRTTY, TINY, DEPS,
     *                  EMINM1, EMAXM1, EXPDEP, LRGEXP
C     .. Statement Functions ..
      DOUBLE PRECISION  APXABS
C     .. Save statement ..
      SAVE              /AC02AG/, /BC02AG/
C     .. Statement Function definitions ..
      APXABS(XC1,XC2) = ABS(XC1) + ABS(XC2)
C     .. Executable Statements ..
      TINY = X02AKF()
      SQRTTY = SQRT(TINY)
      FINITY = X02ALF()
      SQRTFY = SQRT(FINITY)
      EXPDEP = X02BJF() + 1
      EMINM1 = X02BKF() - 1
      EMAXM1 = X02BLF() - 1
      LRGEXP = EMAXM1 + 1 - EXPDEP
      DEPS = X02AJF()
      IERS = IER
      IER = 0
      ITER = 0
      IHALF = 0
      ISPIR = 0
      N = NDEG
C
C     SAVE OVERFLOW/UNDERFLOW INDICATORS OF THE CALLING PROGRAM.
C
      SAVO = OVFLOW
      SAVU = UNFLOW
      OVF = .FALSE.
      UNF = .FALSE.
C
C     MOVE THE COEFFICIENTS A(I) TO DU(I) AND DETERMINE THE LARGEST
C     COEFFICIENT IN MAGNITUDE.
C
      MXCOEF = ZERO
      DO 20 I = 0, N
         DU(I) = A(I)
         MXCOEF = MAX(MXCOEF,ABS(A(I)))
   20 CONTINUE
      IF (MXCOEF.EQ.ZERO) THEN
         DO 40 I = 1, N
            Z(1,I) = FINITY
            Z(2,I) = ZERO
   40    CONTINUE
         N = 0
         OVF = .TRUE.
      ELSE
C
C        DETERMINE A SCALING FOR THE COEFFICIENTS SO THAT THE LARGEST
C        COEFFICIENT IN MAGNITUDE IS LARGE IN MAGNITUDE -- THAT IS, NEAR
C        DEPS * BASE ** X02BLF(), UNLESS THE LARGEST COEFFICIENT IN
C        MAGNITUDE IS LARGER THAN THIS QUANTITY.  IN THIS CASE, SET
C        SCALE TO .FALSE. AND DO NOT SCALE THE COEFFICIENTS.
C
         MXCFEX = C02AGX(MXCOEF)
         IF (MXCFEX.GT.LRGEXP) THEN
            SCBYEX = 0
            SCALE = .FALSE.
         ELSE
            SCBYEX = LRGEXP - MXCFEX
         END IF
      END IF
C
C     INDICATE THAT THE CAUCHY REGION CONTAINING THE SMALLEST ZEROS
C     OF THE CURRENT POLYNOMIAL HAS NOT BEEN COMPUTED.
C
      CAUCHY = .FALSE.
C
C     DO LOOP WHILE N>2.
C
   60 IF (N.GT.2) THEN
C
C        IF SCALE = .TRUE., SCALE THE COEFFICIENTS SO THAT THE LARGEST
C        COEFFICIENT IN MAGNITUDE IS LARGE.
C
         IF (SCALE) THEN
            IF (SCBYEX.NE.0) THEN
               DO 80 I = 0, N
                  DU(I) = C02AGY(DU(I),SCBYEX)
   80          CONTINUE
               SCBYEX = 0
            END IF
         END IF
         UNF = UNFLOW .OR. UNF
C
C        FIND THE NUMBER I OF CONSECUTIVE LEADING COEFFICIENTS EQUAL TO
C        ZERO.
C
         DO 100 I = 0, N - 1
            IF (C02AGS(DU(I))) THEN
C
C              EACH VANISHED LEADING COEFFICIENT YIELDS AN INFINITE
C              ZERO.
C
               Z(1,N-I) = FINITY
               Z(2,N-I) = ZERO
            ELSE
C
C              EXIT THE LOOP ON I -- THE FIRST NON-ZERO COEFFICIENT IS
C              THE I-TH COEFFICIENT.
C
               GO TO 120
            END IF
  100    CONTINUE
  120    IF (I.NE.0) THEN
C
C           SLIDE BACK THE COEFFICIENTS AND DECLARE OVERFLOW.
C
            DO 140 K = I, N
               DU(K-I) = DU(K)
  140       CONTINUE
            N = N - I
C
C           GIVE AN ERROR MESSAGE IF THE VANISHING LEADING COEFFICIENTS
C           HAVE OCCURRED BECAUSE THE POLYNOMIAL HAS BEEN SCALED DOWN.
C           THIS CAN ONLY HAPPEN IF OVERFLOW WAS DETECTED DURING
C           THE COMPUTATION WITH THE ORIGINAL INPUT COEFFICIENTS.
C
            IF (SCBYEX.EQ.-EXPDEP) THEN
               IER = 3
               IF (IERS.NE.1) THEN
                  CALL X04AAF(0,NERR)
                  WRITE (REC,FMT=99999)
                  CALL X04BAF(NERR,REC(1))
                  CALL X04BAF(NERR,REC(2))
               END IF
               GO TO 300
            END IF
C
C           SIGNAL OVERFLOW AND CYCLE THE DO LOOP ON N.
C
            OVF = .TRUE.
            GO TO 60
         END IF
C
C        FIND THE NUMBER I OF CONSECUTIVE TRAILING COEFFICIENTS EQUAL TO
C        ZERO.
C
         DO 160 I = N, 1, -1
            IF (C02AGS(DU(I))) THEN
C
C              EXTRACT ROOTS (IF ANY) FROM THE ORIGIN = (0., 0.)
C
               Z(1,I) = ZERO
               Z(2,I) = ZERO
            ELSE
C
C              EXIT THE LOOP ON I -- THE FIRST NON-ZERO COEFFICIENT IS
C              THE I-TH COEFFICIENT.
C
               GO TO 180
            END IF
  160    CONTINUE
  180    IF (I.NE.N) THEN
C
C           REDUCE THE DEGREE BY THE NUMBER OF ZERO ROOTS AND
C           THEN CYCLE THE DO LOOP ON N.
C
            N = I
            GO TO 60
         END IF
C
C        INITIALIZE LOGICAL UNDERFLOW/OVERFLOW CONDITION STATUS
C        VARIABLES.
C
         OVFLOW = .FALSE.
         UNFLOW = .FALSE.
C
C        HENCEFORTH  N .GT. 2,  DU0 .NE. 0. , AND  DU(N) .NE. 0.
C        CHECK TO SEE WHETHER THE CAUCHY BOUNDS NEED TO BE COMPUTED.
C
         IF ( .NOT. CAUCHY) THEN
C
C           INITIALIZE SOME USEFUL CONSTANTS.
C
            XN = DBLE(N)
            XN1 = DBLE(N-1)
            XN2 = DBLE(N-2)
            X2N = TWO/XN
            X2N1 = X2N/XN1
            XN2N = XN2/XN
            RTN = SQRT(XN)
C
C           CALCULATE  G, AN UPPER BOUND FOR THE SMALLEST ZERO.
C           START WITH  G = ABS( GEOMETRIC MEAN OF THE ZEROS).
C
            G = EXP((LOG(ABS(DU(N)))-LOG(ABS(DU(0))))/XN+SMALL)
C
C           CALCULATE LAGUERRE-STEP  CDIR  AND  FEJER-BOUND, WHICH IS
C           AN UPPER BOUND FOR THE SMALLEST ZERO.
C           CALCULATION OF THE LAGUERRE STEP INVOLVES THE SQUARE OF
C           RECIPROCAL OF NEWTON'S STEP.  SINCE IT CAN EASILY OVERFLOW,
C           THE FEJER BOUND IS CALCULATED WITH NO SUCH OVERFLOWS AND THE
C           LAGUERRE STEP IS CALCULATED FROM IT.
C
            OVFLOW = .FALSE.
            R = F06BLF(DU(N-1),DU(N),OVFLOW)
C
C           IF OVFLOW, A ROOT OF POLYNOMIAL IS WITHIN
C           N * BASE ** (X02BKF()-1) OF  ZERO.
C
            IF (OVFLOW) THEN
C
C              THUS, ASSUME A ROOT IS ZERO, BY ASSUMING DU(N) IS ZERO.
C
               Z(1,N) = ZERO
               Z(2,N) = ZERO
               N = N - 1
C
C              CYCLE THE DO LOOP ON N.
C
               GO TO 60
            END IF
C
C           THE LAGUERRE STEP AND FEJER BOUNDS ARE COMPUTED FROM THE
C           SMALLER ROOT OF A QUADRATIC POLYNOMIAL.
C
            CALL C02AGT(X2N1*DU(N-2),X2N*DU(N-1),DU(N),C,CF1)
            R = XN2N*R
            CTEMP(1) = C(1)*R + XN1
            CTEMP(2) = C(2)*R
            CALL A02ACF(C(1),C(2),CTEMP(1),CTEMP(2),CDIRO(1),CDIRO(2))
            ABDIRO = APXABS(CDIRO(1),CDIRO(2))
            G = MIN(G,BIGONE*MIN(APXABS(C(1),C(2)),RTN*ABDIRO))
C
C           CALCULATE THE CAUCHY-LOWER BOUND  R  FOR THE SMALLEST ZERO
C           BY SOLVING FOR THE ROOT  R  OF THE POLYNOMIAL EQUATION
C             ABS(DU(N)) = SUM( ABS(DU(I))*R**(N-I), I = 0, N-1 )
C           USING  NEWTON'S METHOD.
C
            R = G
            S = BIGONE*G
            UNFLOW = .FALSE.
C
C           NEWTON ITERATION LOOP FOR THE CAUCHY LOWER BOUND R.
C
  200       IF (R.LT.S) THEN
               T = ABS(DU(0))
               S = ZERO
               OVFLOW = .FALSE.
               DO 220 I = 1, N - 1
                  S = R*S + T
                  T = R*T + ABS(DU(I))
  220          CONTINUE
               S = R*S + T
C
C              IT CAN BE PROVED THAT S CANNOT UNDERFLOW.
C
               T = (R*T-ABS(DU(N)))/S
               S = R
               R = R - T
               GO TO 200
            END IF
C
            IF (OVFLOW) THEN
C
C              THE COEFFICIENTS ARE TOO LARGE;  SCALE THEM DOWN AND
C              THEN CYCLE THE DO LOOP ON N.
C
               SCBYEX = -EXPDEP
               GO TO 60
            END IF
C
C           ABS( SMALLEST ROOT ) < R/(2**(1/N) - 1 ) <  1.445*N*R.
C           THUS, 1.445*N*R IS ANOTHER UPPER BOUND AND THE CAUCHY BOUND
C           HAS BEEN COMPUTED, SO SET
C
            CAUCHY = .TRUE.
            UPPERB = MIN(RCONST*XN*R,G)
            LOWERB = SMLONE*S
            UNF = UNFLOW .OR. UNF
         END IF
C
C        NOW   LOWERB < ABS( SMALLEST ZERO ) < UPPERB
C        INITIALIZE THE ITERATION TO BEGIN AT THE ORIGIN.
C        (IN THE CODE BELOW, F0 IS INITIALIZED BUT ITS VALUE NEVER
C        USEFULLY REFERENCED -- IT AVOIDS REFERENCE TO AN UNINITIALIZED
C        VARIABLE IN THE TEST TO ACCEPT THE NEXT ITERATE WHEN THE
C        ITERATION IS NOT STARTED.)
C
         FEJER = UPPERB
         G = UPPERB
         CDIR(1) = CDIRO(1)
         CDIR(2) = CDIRO(2)
         ABDIR = ABDIRO
         RATIO = ABDIR/G
         DZNR = ZERO
         DZNI = ZERO
         FN = ABS(DU(N))
         F0 = FN
         SPIRAL = .FALSE.
         STARTD = .FALSE.
         CONTIN = .TRUE.
C
C        DO WHILE (CONTIN) LOOP, SEARCHING FOR A REAL ROOT,
C        OR PAIR OF COMPLEX ROOTS.  THE NEXT ITERATE IS
C             ZN=CMPLX(DZNR , DZNI).
C
  240    IF (CONTIN) THEN
            ITER = ITER + 1
C
C           RE-ENTRY POINT TO ACCEPT, MODIFY, OR REJECT THE
C           LAGUERRE STEP.
C
C           REJECT  CDIR  IF  ABS(CDIR) > THETA*G .
C
            IF (RATIO.GT.THETA) THEN
C
C              CURRENT LAGUERRE STEP IS NOT ACCEPTABLE.
C              IF STARTD, REDUCE PREVIOUS LAGUERRE STEP BY HALF.
C
               IF (STARTD) THEN
                  IHALF = IHALF + 1
                  ABSCL = HALF*ABSCL
                  CL(1) = HALF*CL(1)
                  CL(2) = HALF*CL(2)
C
C                 HAS THE STEP BECOME NEGLIGIBLE ?
C
                  DX = ABS(DZNR) + ABS(DZNI)
                  IF (DX+ABSCL.NE.DX) THEN
                     DZNR = DZ0R + CL(1)
                     DZNI = DZ0I + CL(2)
                  ELSE
C
C                    OTHERWISE, C02AGF HAS HUNG-UP.
C
                     IF (FN.GE.E*XN**2) THEN
                        IER = 2
                        IF (IERS.NE.1) THEN
                           CALL X04AAF(0,NERR)
                           WRITE (REC,FMT=99997)
                           CALL X04BAF(NERR,REC(1))
                           CALL X04BAF(NERR,REC(2))
                        END IF
                        GO TO 300
                     END IF
C
C                    EXIT THE ITERATION LOOP  DO WHILE(CONTIN).
C
                     CONTIN = .FALSE.
                     GO TO 240
                  END IF
               ELSE
C
C                 IF .NOT. STARTD, HAS ZN BEEN ON THE INNER CAUCHY
C                 RADIUS ?
C
                  ISPIR = ISPIR + 1
                  IF (SPIRAL) THEN
                     C(1) = CSPIR(1)*DZNR - CSPIR(2)*DZNI
                     C(2) = CSPIR(2)*DZNR + CSPIR(1)*DZNI
                  ELSE
C
C                    SET SPIRAL TO  .TRUE..  PUT  ZN  ON THE INNER
C                    CIRCLE OF THE ANNULUS CONTAINING THE SMALLEST
C                    ZERO IN THE DIRECTION OF THE LAGUERRE STEP.
C
                     SPIRAL = .TRUE.
                     CSPIR(1) = -ONEPQT/XN
                     CSPIR(2) = ONE
                     ABSCL = LOWERB/XN**2
                     CTEMP(1) = CDIR(1)/ABDIR
                     CTEMP(2) = CDIR(2)/ABDIR
                     C(1) = CTEMP(1)*LOWERB
                     C(2) = CTEMP(2)*LOWERB
                  END IF
C
C                 SET  ZN  TO THE NEXT POINT ON THE SPIRAL.
C
                  DZNR = C(1)
                  DZNI = C(2)
               END IF
            ELSE
C
C              CDIR  AT THE ORIGIN IS IN THE DIRECTION OF DECREASING
C              FUNCTION VALUE, SO
C
               STARTD = .TRUE.
C
C              ACCEPT  CDIR  IF  ABS(CDIR) <= GAMA*G.
C
               IF (RATIO.GT.GAMA .AND. (STARTD .OR. SPIRAL .OR.
     *             LOWERB.LE.GAMA*G)) THEN
                  RATIO = GAMA/RATIO
                  CDIR(1) = CDIR(1)*RATIO
                  CDIR(2) = CDIR(2)*RATIO
                  ABDIR = ABDIR*RATIO
               END IF
C
C              ACCEPT THE PREVIOUS ITERATE.
C              SAVE THE DATA ASSOCIATED WITH THE CURRENT ITERATE.
C
               G = FEJER
               CL(1) = CDIR(1)
               CL(2) = CDIR(2)
               ABSCL = ABDIR
               F0 = FN
               DZ0R = DZNR
               DZ0I = DZNI
               DZNR = DZ0R + CL(1)
               DZNI = DZ0I + CL(2)
            END IF
C
C           IS  ZN  CLOSE TO THE REAL AXIS RELATIVE TO STEP SIZE ?
C
            IF (ABS(DZNI).GT.PHI*ABSCL) THEN
C
C              ZN  IS COMPLEX.  THUS, ATTEMPT TO DIVIDE THE POLYNOMIAL
C              BY THE QUADRATIC FACTOR  (Z**2-X2*Z+R); THAT IS,
C
C              SUM(DU(I)*Z**(N-I)) =
C                  (Z**2-X2*Z+R) * SUM( D(I)*Z**(N-I-2), I=0,N-2) +
C                                     D(N-1)*(Z-X) + D(N)   FOR ALL Z,
C
C              WHERE (X,Y) AND (X,-Y) ARE ZEROS OF  Z**2-X2*Z+R.
C              IN THE CODE BELOW,
C              E  IS ERROR BOUND FOR THE VALUE OF POLYNOMIAL,
C              CF IS THE  VALUE OF THE POLYNOMIAL AT (X,Y),
C              CF1 IS THE FIRST  DERIVITIVE OF POLYNOMIAL AT (X,Y),
C              2*CF2 IS THE SECOND DERIVATIVE OF POLYNOMIAL AT (X,Y),
C              AND D(I) ARE THE COEFFICIENTS OF QUOTIENT POLYNOMIAL
C              DEFLAT(I).
C
C              BE SURE THAT THE OVERFLOW INDICATOR IS TURNED OFF.
C
               UNFLOW = .FALSE.
               OVFLOW = .FALSE.
               CALL C02AGW(DZNR,DZNI,N,DU,CF,CF1,CF2,E,DEFLAT)
               FN = A02ABF(CF(1),CF(2))
C
C              CHECK FOR OVERFLOW.
C
               IF (OVFLOW) THEN
C
C                 INDICATE THAT THE POLYNOMIAL NEEDS TO BE SCALED DOWN
C                 AND CYCLE THE DO LOOP ON N.  NOTE: CAUCHY IS NOT RESET
C                 AS THE CAUCHY BOUNDS NEED NOT BE RECOMPUTED.
C
                  SCBYEX = -EXPDEP
                  GO TO 60
               END IF
C
C              CHECK TO SEE IF  ZN  IS A ZERO OR IF UNDERFLOW HAS
C              OCCURRED.
C
               IF (FN.LE.E .OR. UNFLOW) THEN
                  IF (UNFLOW) THEN
                     IER = 3
                     IF (IERS.NE.1) THEN
                        CALL X04AAF(0,NERR)
                        WRITE (REC,FMT=99998)
                        CALL X04BAF(NERR,REC(1))
                        CALL X04BAF(NERR,REC(2))
                     END IF
                     GO TO 300
                  END IF
C
C                 A ROOT HAS BEEN FOUND -- EXIT THE ITERATION
C                 LOOP  DO WHILE(CONTIN).
C
                  CONTIN = .FALSE.
                  GO TO 240
               END IF
               UNF = UNFLOW .OR. UNF
C
C              HAS THE FUNCTION VALUE DECREASED ?
C
               IF (FN.GE.F0 .AND. STARTD) THEN
C
C                 NO, IT HAS NOT.  INDICATE THAT THE LAGUERRE STEP IS
C                 UNACCEPTABLE.  (A RATIO LARGER THAN THETA INDICATES
C                 THAT THE LAGUERRE STEP SHOULD BE SHORTENED.)
C
                  RATIO = BIGONE*THETA
C
C                 CYCLE ITERATION LOOP  DO WHILE(CONTIN).
C
                  GO TO 240
               END IF
C
C              FIND THE LAGUERRE STEP AT  ZN.
C
               OVFLOW = .FALSE.
               CALL A02ACF(CF1(1),CF1(2),CF(1),CF(2),C(1),C(2))
C
C              IF OVFLOW, A ROOT OF POLYNOMIAL IS WITHIN  N*2**X02BKF()
C              OF  ZERO.
C
               IF (OVFLOW) THEN
                  UNF = .TRUE.
C
C                 A ROOT HAS BEEN FOUND -- EXIT THE ITERATION LOOP.
C
                  CONTIN = .FALSE.
                  GO TO 240
               END IF
C
C              COMPUTE THE LAGUERRE STEP  CDIR  AND THE BOUND  FEJER
C              AT  ZN.  THE LAGUERRE STEP AND FEJER BOUNDS ARE COMPUTED
C              FROM THE SMALLER ROOT OF A QUADRATIC POLYNOMIAL.
C
               CF2(1) = CF2(1)*X2N1
               CF2(2) = CF2(2)*X2N1
               CTEMP(1) = CF1(1)*X2N
               CTEMP(2) = CF1(2)*X2N
               CALL C02AGU(CF2,CTEMP,CF,CDIR,CF1)
               FEJER = APXABS(CDIR(1),CDIR(2))
               CTEMP(1) = C(1)*XN2N
               CTEMP(2) = C(2)*XN2N
               C(1) = CTEMP(1)*CDIR(1) - CTEMP(2)*CDIR(2)
               C(2) = CTEMP(2)*CDIR(1) + CTEMP(1)*CDIR(2)
               C(1) = C(1) + XN1
               CTEMP(1) = CDIR(1)
               CTEMP(2) = CDIR(2)
               CALL A02ACF(CTEMP(1),CTEMP(2),C(1),C(2),CDIR(1),CDIR(2))
               ABDIR = APXABS(CDIR(1),CDIR(2))
               RATIO = ABDIR/G
               FEJER = MIN(RTN*ABDIR,FEJER)
C
C              IS THE STEP SIZE NEGLIGIBLE ?  (THIS TEST MAY BE PROVABLY
C              REDUNDANT ON SOME WELL-BEHAVED ARITHMETICS.)
C
               DX = ABS(DZNR) + ABS(DZNI)
               IF (DX+ABDIR.EQ.DX) THEN
C
C                 THE STEP IS NEGLIGIBLE.  ASSUME  ZN=(DZNR,DZNI) IS A
C                 ROOT. EXIT THE ITERATION LOOP  DO WHILE(CONTIN).
C
                  CONTIN = .FALSE.
                  GO TO 240
               END IF
C
C              NOW DETERMINE WHETHER  CDIR  IS ACCEPTABLE.
C
            ELSE
C
C              ZN  IS REAL.  THUS, ATTEMPT TO DIVIDE THE POLYNOMIAL
C              BY THE LINEAR FACTOR (Z-X); THAT IS,
C
C              SUM(DU(I)*Z**(N-I)) =
C                  (Z-X) * SUM(D(I)*Z**(N-I-1), I=0,N-1) + D(N)
C                               FOR ALL Z,
C
C              WHERE X IS A ZERO OF  Z-X.  IN THE CODE BELOW,
C              E    IS ERROR BOUND FOR THE VALUE OF POLYNOMIAL AT X,
C              F    IS THE VALUE OF THE POLYNOMIAL AT X,
C              V    IS THE FIRST DERIVITIVE OF POLYNOMIAL AT X,
C              2*W  IS THE SECOND DERIVATIVE OF POLYNOMIAL AT X, AND
C              D(I) ARE THE COEFFICIENTS OF QUOTIENT POLYNOMIAL
C              DEFLAT(I).
C
C              BE SURE THAT THE OVERFLOW INDICATOR IS TURNED OFF.
C
               OVFLOW = .FALSE.
               UNFLOW = .FALSE.
C
C              THE ITERATE IS TAKEN AS REAL. SET THE IMAGINARY PART TO
C              ZERO.
C
               DZNI = ZERO
               CALL C02AGV(DZNR,N,DU,F,V,W,E,DEFLAT)
               FN = ABS(F)
C
C              CHECK FOR OVERFLOW.
C
               IF (OVFLOW) THEN
C
C                 INDICATE THAT THE POLYNOMIAL NEEDS TO BE SCALED DOWN
C                 AND CYCLE THE DO LOOP ON N.  NOTE: CAUCHY IS NOT RESET
C                 AS THE CAUCHY BOUNDS NEED NOT BE RECOMPUTED.
C
                  SCBYEX = -EXPDEP
                  GO TO 60
               END IF
C
C              CHECK TO SEE IF  ZN  IS A ZERO OR IF UNDERFLOW HAS
C              OCCURRED.
C
               IF (FN.LE.E .OR. UNFLOW) THEN
                  IF (UNFLOW) THEN
                     IER = 3
                     IF (IERS.NE.1) THEN
                        CALL X04AAF(0,NERR)
                        WRITE (REC,FMT=99998)
                        CALL X04BAF(NERR,REC(1))
                        CALL X04BAF(NERR,REC(2))
                     END IF
                     GO TO 300
                  END IF
C
C                 A ROOT HAS BEEN FOUND -- EXIT THE ITERATION
C                 LOOP  DO WHILE(CONTIN).
C
                  CONTIN = .FALSE.
                  GO TO 240
               END IF
               UNF = UNFLOW .OR. UNF
C
C              HAS THE FUNCTION VALUE DECREASED ?
C
               IF (FN.GE.F0 .AND. STARTD) THEN
C
C                 NO, IT HAS NOT.  INDICATE THAT THE LAGUERRE STEP IS
C                 UNACCEPTABLE.  (A RATIO LARGER THAN THETA INDICATES
C                 THAT THE LAGUERRE STEP SHOULD BE SHORTENED.)
C
                  RATIO = BIGONE*THETA
C
C                 CYCLE ITERATION LOOP  DO WHILE(CONTIN).
C
                  GO TO 240
               END IF
               OVFLOW = .FALSE.
C
C              FIND THE LAGUERRE STEP AT  DZNR.
C
               R = F06BLF(V,F,OVFLOW)
C
C              IF OVFLOW,  A ROOT OF POLYNOMIAL IS WITHIN
C              4 * N * BASE ** (X02BKF()-1)  OF  ZN.
C
               IF (OVFLOW) THEN
                  UNF = .TRUE.
C
C                 A ROOT HAS BEEN FOUND -- EXIT THE ITERATION
C                 LOOP  DO WHILE(CONTIN).
C
                  CONTIN = .FALSE.
                  GO TO 240
               END IF
C
C              COMPUTE THE LAGUERRE STEP  CDIR  AND THE BOUND  FEJER
C              AT  ZN .  THE LAGUERRE STEP AND FEJER BOUNDS ARE COMPUTED
C              FROM THE SMALLER ROOT OF A QUADRATIC POLYNOMIAL.
C
               CALL C02AGT(X2N1*W,X2N*V,F,C,CF1)
               FEJER = APXABS(C(1),C(2))
               R = XN2N*R
               CTEMP(1) = C(1)*R + XN1
               CTEMP(2) = C(2)*R
               CALL A02ACF(C(1),C(2),CTEMP(1),CTEMP(2),CDIR(1),CDIR(2))
               ABDIR = APXABS(CDIR(1),CDIR(2))
               RATIO = ABDIR/G
               FEJER = MIN(RTN*ABDIR,FEJER)
C
C              IS THE STEP SIZE NEGLIGIBLE ?
C
               DX = ABS(DZNR)
               IF (DX+ABDIR.EQ.DX) THEN
C
C                 THE STEP IS NEGLIGIBLE.  ASSUME  ZN=(DZNR,DZNI) IS
C                 A ROOT. EXIT THE ITERATION LOOP  DO WHILE(CONTIN).
C
                  CONTIN = .FALSE.
                  GO TO 240
               END IF
C
C              NOW DETERMINE WHETHER  CDIR  IS ACCEPTABLE.
C
            END IF
C
C           REPEAT THE ITERATION LOOP DO WHILE(CONTIN).
C
            GO TO 240
         END IF
C
C        A ROOT HAS BEEN COMPUTED.  DEFLATE THE POLYNOMIAL.
C
         IF (DZNI.NE.ZERO) THEN
C
C           ACCEPT ZN AS A COMPLEX ROOT AND DEFLATE FOR A COMPLEX ROOT.
C           PUT COEFFICIENTS OF THE QUOTIENT POLYNOMIAL IN THE DU ARRAY.
C           DU(0) IS UNCHANGED FOR THE DEFLATED POLYNOMIAL.
C
            DO 260 I = 1, N - 2
               DU(I) = DEFLAT(I)
  260       CONTINUE
            Z(1,N) = DZNR
            Z(2,N) = DZNI
            Z(1,N-1) = Z(1,N)
            Z(2,N-1) = -Z(2,N)
            N = N - 2
         ELSE
C
C           ACCEPT ZN AS A REAL ROOT AND DEFLATE FOR A REAL ROOT.
C           PUT COEFFICIENTS OF THE QUOTIENT POLYNOMIAL IN THE DU ARRAY.
C           DU(0) IS UNCHANGED FOR THE DEFLATED POLYNOMIAL.
C
            DO 280 I = 1, N - 1
               DU(I) = DEFLAT(I)
  280       CONTINUE
            Z(1,N) = DZNR
            Z(2,N) = ZERO
            N = N - 1
         END IF
C
C        INDICATE THAT THE CAUCHY REGION CONTAINING THE SMALLEST ZEROS
C        OF THE CURRENT POLYNOMIAL HAS NOT BEEN COMPUTED.
C
         CAUCHY = .FALSE.
C
C        REPEAT THE LOOP WHILE N>2 FOR DECREASING N.
C
         GO TO 60
      END IF
C
C     THE POLYNOMIAL IS NOW OF DEGREE 2 OR LESS.  DETERMINE THE
C     REMAINING ROOTS DIRECTLY RATHER THAN ITERATIVELY.
C
      OVFLOW = .FALSE.
      UNFLOW = .FALSE.
      IF (N.EQ.2) THEN
         CALL C02AGT(DU(0),DU(1),DU(2),CTEMP,C)
         Z(1,1) = C(1)
         Z(2,1) = C(2)
         Z(1,2) = CTEMP(1)
         Z(2,2) = CTEMP(2)
      ELSE IF (N.EQ.1) THEN
         Z(1,1) = -DU(1)/DU(0)
         Z(2,1) = ZERO
      ELSE
         OVF = OVF .OR. OVFLOW
         UNF = UNF .OR. UNFLOW
C
C        RESTORE OVERFLOW AND UNDERFLOW INDICATORS AND ENABLE MESSAGE.
C
         OVFLOW = SAVO
         UNFLOW = SAVU
C
C        PROVIDE ONLY THE RELEVANT OVER/UNDERFLOW MESSAGES.
C
         IF (OVF) R = FINITY*FINITY
         IF (UNF) R = TINY*TINY
      END IF
  300 RETURN
C
99999 FORMAT (' ** C02AGF cannot evaluate p(z) near some of its zeros ',
     *       'without overflow.',/' ** If this message occurs please c',
     *       'ontact NAG Central Office.')
99998 FORMAT (' ** C02AGF cannot evaluate p(z) near some of its zeros ',
     *       'without underflow.',/' ** If this message occurs please ',
     *       'contact NAG Central Office.')
99997 FORMAT (' ** The method has failed. This error is very unlikely ',
     *       'to occur.',/' ** Please contact NAG Central Office.')
      END
      DOUBLE PRECISION FUNCTION X02AMF()
C     MARK 12 RELEASE. NAG COPYRIGHT 1986.
C
C     RETURNS THE 'SAFE RANGE' PARAMETER
C     I.E. THE SMALLEST POSITIVE MODEL NUMBER Z SUCH THAT
C     FOR ANY X WHICH SATISFIES X.GE.Z AND X.LE.1/Z
C     THE FOLLOWING CAN BE COMPUTED WITHOUT OVERFLOW, UNDERFLOW OR OTHER
C     ERROR
C
C        -X
C        1.0/X
C        SQRT(X)
C        LOG(X)
C        EXP(LOG(X))
C        Y**(LOG(X)/LOG(Y)) FOR ANY Y
C
      DOUBLE PRECISION X02CON
      DATA X02CON /2.22507385850721D-308 /
C     .. Executable Statements ..
      X02AMF = X02CON
      RETURN
      END
      INTEGER FUNCTION X02BBF(X)
C     NAG COPYRIGHT 1975
C     MARK 4.5 RELEASE
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     * MAXINT *
C     RETURNS THE LARGEST INTEGER REPRESENTABLE ON THE COMPUTER
C     THE X PARAMETER IS NOT USED
C     .. Scalar Arguments ..
      DOUBLE PRECISION        X
C     .. Executable Statements ..
C     FOR ICL 1900
C      X02BBF = 8388607
      X02BBF =           2147483647
      RETURN
      END
      DOUBLE PRECISION FUNCTION X02AJF()
C     MARK 12 RELEASE. NAG COPYRIGHT 1986.
C
C     RETURNS  (1/2)*B**(1-P)  IF ROUNDS IS .TRUE.
C     RETURNS  B**(1-P)  OTHERWISE
C
      DOUBLE PRECISION X02CON
      DATA X02CON /1.11022302462516D-16 /
C     .. Executable Statements ..
      X02AJF = X02CON
      RETURN
      END
      SUBROUTINE C02AGF(A,N,SCALE,Z,WORK,IFAIL)
C     MARK 13 RELEASE. NAG COPYRIGHT 1988.
C
C     C02AGF ATTEMPTS TO FIND ALL THE ROOTS OF THE NTH ORDER REAL
C     POLYNOMIAL EQUATION
C
C        A(0)*Z**N + A(1)*Z**(N-1) + ... + A(N-1)*Z + A(N) = 0.
C
C     THE ZEROS OF POLYNOMIALS OF DEGREE 1 AND 2 ARE CALCULATED USING
C     THE "STANDARD" CLOSED FORMULAS
C           Z = -A(1)/A(0) AND
C           Z = (-A(1) +/- SQRT(DISC))/(2*A(0)) RESPECTIVELY, WHERE
C        DISC = A(1)**2 - 4*A(0)*A(2).
C     FOR N >= 3, THE ROOTS ARE LOCATED ITERATIVELY USING A VARIANT OF
C     LAGUERRE'S METHOD, WHICH IS CUBICALLY CONVERGENT FOR ISOLATED
C     ZEROS (REAL OR COMPLEX) AND LINEARLY CONVERGENT FOR MULTIPLE
C     ZEROS.
C
C     C02AGF ITSELF IS ESSENTIALLY A DUMMY ROUTINE WHOSE FUNCTION IS TO
C     PARTITION THE WORK ARRAY WORK FOR USE BY C02AGZ.
C     WORK IS PARTITIONED INTO 2 ARRAYS EACH OF SIZE (N + 1).
C
C     .. Parameters ..
      DOUBLE PRECISION  ZERO
      PARAMETER         (ZERO=0.0D0)
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='C02AGF')
C     .. Scalar Arguments ..
      INTEGER           IFAIL, N
      LOGICAL           SCALE
C     .. Array Arguments ..
      DOUBLE PRECISION  A(0:N), WORK(2*(N+1)), Z(2,N)
C     .. Local Scalars ..
      INTEGER           IER, NDEG, NREC
      LOGICAL           SC
C     .. Local Arrays ..
      CHARACTER*80      REC(2)
C     .. External Functions ..
      INTEGER           P01ABF
      EXTERNAL          P01ABF
C     .. External Subroutines ..
      EXTERNAL          C02AGZ
C     .. Executable Statements ..
      IER = IFAIL
      SC = SCALE
      NREC = 0
      NDEG = N
      IF (N.LT.1 .OR. A(0).EQ.ZERO) THEN
         IER = 1
         WRITE (REC,FMT=99999) N, A(0)
         NREC = 2
         GO TO 20
      END IF
      CALL C02AGZ(A,NDEG,SC,Z,WORK(1),WORK(N+2),IER)
      IF (IER.EQ.0) THEN
         IFAIL = 0
         GO TO 40
      END IF
   20 IFAIL = P01ABF(IFAIL,IER,SRNAME,NREC,REC)
   40 RETURN
C
99999 FORMAT (' ** On entry, N.lt.1 or A(0).eq.0:',/'    N = ',I8,'  A',
     *       '(0) = ',1P,D13.5)
      END
      SUBROUTINE D01FCF(NDIM,A,B,MINPTS,MAXPTS,FUNCTN,EPS,ACC,LENWRK,
     *                  WRKSTR,FINVAL,IFAIL)
C     MARK 8 RELEASE. NAG COPYRIGHT 1979.
C     MARK 10B REVISED. IER-399 (JAN 1983).
C     MARK 11C REVISED. IER-464 (MAR 1985).
C     MARK 11.5(F77) REVISED. (SEPT 1985.)
C     MARK 12A REVISED. IER-498 (AUG 1986).
C     MARK 13 REVISED. USE OF MARK 12 X02 FUNCTIONS (APR 1988).
C
C     ADAPTIVE MULTIDIMENSIONAL INTEGRATION SUBROUTINE
C
C     *********  PARAMETERS FOR D01FCF  ****************************
C
C      INPUT PARAMETERS
C
C     NDIM    INTEGER NUMBER OF VARIABLES, MUST EXCEED 1 BUT
C             NOT EXCEED 15.
C
C     A       REAL ARRAY OF LOWER LIMITS, WITH DIMENSION NDIM
C
C     B       REAL ARRAY OF UPPER LIMITS, WITH DIMENSION NDIM
C
C     MINPTS  INTEGER MINIMUM NUMBER OF INTEGRAND VALUES TO BE
C             ALLOWED, WHICH MUST NOT EXCEED MAXPTS.
C
C     MAXPTS  INTEGER MAXIMUM NUMBER OF INTEGRAND VALUES TO BE
C             ALLOWED, WHICH MUST BE AT LEAST
C             2**NDIM+2*NDIM**2+2*NDIM+1.
C
C     FUNCTN  EXTERNALLY DECLARED USER DEFINED REAL FUNCTION
C             INTEGRAND. IT MUST HAVE PARAMETERS (NDIM,Z),
C             WHERE Z IS A REAL ARRAY OF DIMENSION NDIM.
C
C     EPS     REAL REQUIRED RELATIVE ACCURACY, MUST BE GREATER
C             THAN ZERO
C
C     LENWRK  INTEGER LENGTH OF ARRAY WRKSTR, MUST BE AT LEAST
C             2*NDIM+4.
C
C     IFAIL   INTEGER NAG FAILURE PARAMETER
C             IFAIL=0 FOR HARD FAIL
C             IFAIL=1 FOR SOFT FAIL
C
C      OUTPUT PARAMETERS
C
C     MINPTS  INTEGER NUMBER OF INTEGRAND VALUES USED BY THE
C             ROUTINE
C
C     WRKSTR  REAL ARRAY OF WORKING STORAGE OF DIMENSION (LENWRK).
C
C     ACC     REAL ESTIMATED RELATIVE ACCURACY OF FINVAL
C
C     FINVAL  REAL ESTIMATED VALUE OF INTEGRAL
C
C     IFAIL   IFAIL=0 FOR NORMAL EXIT, WHEN ESTIMATED RELATIVE
C                  LESS INTEGACCURACY RAND VALUES USED.
C
C          IFAIL=1 IF NDIM.LT.2, NDIM.GT.15, MINPTS.GT.MAXPTS,
C                  MAXPTS.LT.2**NDIM+2*NDIM*(NDIM+1)+1, EPS.LE.0
C                  OR LENWRK.LT.2*NDIM+4.
C
C          IFAIL=2 IF MAXPTS WAS TOO SMALL FOR D01FCF TO OBTAIN THE
C                  REQUIRED RELATIVE ACCURACY EPS.  IN THIS
C                  CASE D01FCF RETURNS A VALUE OF FINVAL
C                  WITH ESTIMATED RELATIVE ACCURACY ACC.
C
C          IFAIL=3 IF LENWRK TOO SMALL FOR MAXPTS INTEGRAND
C                  VALUES.  IN THIS CASE D01FCF RETURNS A
C                  VALUE OF FINVAL WITH ESTIMATED ACCURACY
C                  ACC USING THE WORKING STORAGE
C                  AVAILABLE, BUT ACC WILL BE GREATER
C                  THAN EPS.
C
C     **************************************************************
C
C     .. Parameters ..
      CHARACTER*6       SRNAME
      PARAMETER         (SRNAME='D01FCF')
C     .. Scalar Arguments ..
      DOUBLE PRECISION  ACC, EPS, FINVAL
      INTEGER           IFAIL, LENWRK, MAXPTS, MINPTS, NDIM
C     .. Array Arguments ..
      DOUBLE PRECISION  A(NDIM), B(NDIM), WRKSTR(LENWRK)
C     .. Function Arguments ..
      DOUBLE PRECISION  FUNCTN
      EXTERNAL          FUNCTN
C     .. Local Scalars ..
      DOUBLE PRECISION  ABSERR, DF1, DF2, DIFMAX, F1, F2, F3, F4, HALF,
     *                  LAMDA2, LAMDA4, LAMDA5, ONE, RATIO, RGNCMP,
     *                  RGNERR, RGNERT, RGNVAL, RGNVLT, RGNVOL, RLNDIM,
     *                  SUM1, SUM2, SUM3, SUM4, SUM5, TWO, TWONDM,
     *                  WEIT1, WEIT2, WEIT3, WEIT4, WEIT5, WEITP1,
     *                  WEITP2, WEITP3, WEITP4, ZERO
      INTEGER           DVAXES, DVAXIS, DVFLAG, FUNCLS, IERROR, J, K,
     *                  MAXAXS, MXRGNS, POINTR, RGNCLS, RULCLS, SBRGNS,
     *                  SUBRGN, SUBTMP, TPONTP, TPONTR
C     .. Local Arrays ..
      DOUBLE PRECISION  CENTER(15), DIF(15), OLDCNT(15), WIDTH(15),
     *                  Z(15)
      INTEGER           DVCNTL(15), DVCNTR(15)
      CHARACTER*1       P01REC(1)
C     .. External Functions ..
      DOUBLE PRECISION  X02AJF, X02AMF
      INTEGER           P01ABF, X02BBF
      EXTERNAL          X02AJF, X02AMF, P01ABF, X02BBF
C     .. Intrinsic Functions ..
      INTRINSIC         ABS, MIN, DBLE, SQRT, INT
C     .. Data statements ..
      DATA              ZERO, ONE, TWO, HALF/0.0D0, 1.0D0, 2.0D0, 0.5D0/
C     .. Executable Statements ..
C
C       SUBROUTINE INITIALISATION AND PARAMETER CHECKING
C
      FUNCLS = 0
      IF (NDIM.LT.2 .OR. NDIM.GT.15) GO TO 580
      IF (MINPTS.GT.MAXPTS) GO TO 580
      IF (EPS.LE.ZERO) GO TO 580
      IF (LENWRK.LT.2*NDIM+4) GO TO 580
      TWONDM = TWO**NDIM
      RGNVOL = TWONDM
      DVFLAG = 1
      MAXAXS = INT(MIN(DBLE(X02BBF(ONE)),1.0D0/X02AJF()))
      MAXAXS = (MAXAXS-NDIM)/(NDIM+1)
      MXRGNS = LENWRK/(2*NDIM+4)
      SBRGNS = 0
      RGNVLT = ZERO
      RGNERT = ZERO
      DO 20 J = 1, NDIM
         CENTER(J) = (A(J)+B(J))*HALF
         DIF(J) = ZERO
         WIDTH(J) = (B(J)-A(J))*HALF
         DVCNTL(J) = 1
         DVCNTR(J) = 1
         OLDCNT(J) = CENTER(J)
         RGNVOL = RGNVOL*WIDTH(J)
   20 CONTINUE
C
C       END SUBROUTINE INITIALISATION
C       BASIC RULE INITIALISATION
C
      RULCLS = 2**NDIM + 2*NDIM*NDIM + 2*NDIM + 1
      IF (MAXPTS.LT.RULCLS) GO TO 580
      FUNCLS = RULCLS
      RLNDIM = NDIM
      LAMDA2 = SQRT(9.0D0/70.0D0)
      LAMDA4 = SQRT(9.0D0/10.0D0)
      LAMDA5 = SQRT(9.0D0/19.0D0)
      WEIT1 = (12824.0D0-9120.0D0*RLNDIM+400.0D0*RLNDIM*RLNDIM)
     *        /19683.0D0
      WEIT2 = 980.0D0/6561.0D0
      WEIT3 = (1820.0D0-400.0D0*RLNDIM)/19683.0D0
      WEIT4 = 200.0D0/19683.0D0
      WEIT5 = 6859.0D0/19683.0D0/TWONDM
      WEITP1 = (729.0D0-950.0D0*RLNDIM+50.0D0*RLNDIM**2)/729.0D0
      WEITP2 = 245.0D0/486.0D0
      WEITP3 = (265.0D0-100.0D0*RLNDIM)/1458.0D0
      WEITP4 = 25.0D0/729.0D0
      RATIO = (LAMDA2/LAMDA4)**2
C
C       END BASIC RULE INITIALISATION
      GO TO 100
C       DIVIDE SUBREGION WITH LARGEST ERROR AND PREPARE TO USE
C       BASIC RULE ON EACH PORTION
C
   40 SUBRGN = 1
      POINTR = WRKSTR(1)
      RGNCLS = RULCLS
      RGNVOL = TWONDM
      TPONTR = POINTR + 2
      DO 60 J = 1, NDIM
         TPONTR = TPONTR + 2
         CENTER(J) = WRKSTR(TPONTR-1)
         WIDTH(J) = WRKSTR(TPONTR)
         DVCNTR(J) = 1
         DVCNTL(J) = 1
         OLDCNT(J) = CENTER(J)
         RGNVOL = RGNVOL*WIDTH(J)
   60 CONTINUE
      DVAXES = WRKSTR(POINTR+2)
      IF (DVAXES.LT.0) GO TO 620
   80 DVAXIS = DVAXES
      DVAXES = DVAXIS/(NDIM+1)
      DVAXIS = DVAXIS - (NDIM+1)*DVAXES
      DVCNTL(DVAXIS) = 2*DVCNTL(DVAXIS)
      RGNCLS = RGNCLS*2
      IF (DVAXES.GT.0) GO TO 80
      IF (FUNCLS+RGNCLS.GT.MAXPTS) GO TO 600
      IF (RGNCLS/RULCLS+SBRGNS-1.GT.MXRGNS) DVFLAG = 2
      FUNCLS = FUNCLS + RGNCLS
C
C       BEGIN BASIC RULE
  100 DO 120 J = 1, NDIM
         Z(J) = CENTER(J)
  120 CONTINUE
      SUM1 = FUNCTN(NDIM,Z)
      SUM2 = ZERO
      SUM3 = ZERO
      DO 140 J = 1, NDIM
         Z(J) = CENTER(J) - LAMDA2*WIDTH(J)
         F1 = FUNCTN(NDIM,Z)
         Z(J) = CENTER(J) + LAMDA2*WIDTH(J)
         F2 = FUNCTN(NDIM,Z)
         Z(J) = CENTER(J) - LAMDA4*WIDTH(J)
         F3 = FUNCTN(NDIM,Z)
         Z(J) = CENTER(J) + LAMDA4*WIDTH(J)
         F4 = FUNCTN(NDIM,Z)
         SUM2 = SUM2 + F1 + F2
         SUM3 = SUM3 + F3 + F4
         DF1 = F1 + F2 - TWO*SUM1
         DF2 = F3 + F4 - TWO*SUM1
         DIF(J) = DIF(J) + ABS(DF1-RATIO*DF2)
         Z(J) = CENTER(J)
  140 CONTINUE
      SUM4 = ZERO
      DO 200 J = 2, NDIM
         Z(J-1) = CENTER(J-1) - LAMDA4*WIDTH(J-1)
         DO 160 K = J, NDIM
            Z(K) = CENTER(K) - LAMDA4*WIDTH(K)
            SUM4 = SUM4 + FUNCTN(NDIM,Z)
            Z(K) = CENTER(K) + LAMDA4*WIDTH(K)
            SUM4 = SUM4 + FUNCTN(NDIM,Z)
            Z(K) = CENTER(K)
  160    CONTINUE
         Z(J-1) = CENTER(J-1) + LAMDA4*WIDTH(J-1)
         DO 180 K = J, NDIM
            Z(K) = CENTER(K) - LAMDA4*WIDTH(K)
            SUM4 = SUM4 + FUNCTN(NDIM,Z)
            Z(K) = CENTER(K) + LAMDA4*WIDTH(K)
            SUM4 = SUM4 + FUNCTN(NDIM,Z)
            Z(K) = CENTER(K)
  180    CONTINUE
         Z(J-1) = CENTER(J-1)
  200 CONTINUE
      SUM5 = ZERO
      DO 220 J = 1, NDIM
         Z(J) = CENTER(J) - LAMDA5*WIDTH(J)
  220 CONTINUE
  240 DO 260 J = 2, NDIM
         IF (Z(J-1).LT.CENTER(J-1)+WIDTH(J-1)) GO TO 280
         Z(J-1) = CENTER(J-1) - LAMDA5*WIDTH(J-1)
         Z(J) = Z(J) + TWO*LAMDA5*WIDTH(J)
  260 CONTINUE
      IF (Z(NDIM).GT.CENTER(NDIM)+WIDTH(NDIM)) GO TO 300
  280 SUM5 = SUM5 + FUNCTN(NDIM,Z)
      Z(1) = Z(1) + TWO*LAMDA5*WIDTH(1)
      GO TO 240
  300 RGNVAL = RGNVOL*(WEIT1*SUM1+WEIT2*SUM2+WEIT3*SUM3+WEIT4*SUM4+
     *         WEIT5*SUM5)
      RGNCMP = RGNVOL*(WEITP1*SUM1+WEITP2*SUM2+WEITP3*SUM3+WEITP4*SUM4)
      RGNERR = ABS(RGNVAL-RGNCMP)
C
C       END BASIC RULE
C       STORE RESULTS OF BASIC RULE APPLICATION
C
      RGNVLT = RGNVLT + RGNVAL
      RGNERT = RGNERT + RGNERR
      IF (DVFLAG.EQ.0) GO TO 340
      IF (DVFLAG.EQ.2) GO TO 500
      POINTR = MXRGNS + SBRGNS*(2*NDIM+3) + 1
      SBRGNS = SBRGNS + 1
      WRKSTR(SBRGNS) = POINTR
      SUBRGN = SBRGNS
      TPONTR = POINTR + 2
      DO 320 J = 1, NDIM
         TPONTR = TPONTR + 2
         WRKSTR(TPONTR-1) = CENTER(J)
         WRKSTR(TPONTR) = WIDTH(J)
  320 CONTINUE
  340 WRKSTR(POINTR) = RGNERT
      WRKSTR(POINTR+1) = RGNVLT
C       DETERMINE AXIS ALONG WHICH FOURTH DIFFERENCE IS LARGEST
      DIFMAX = ZERO
      DO 380 J = 1, NDIM
         IF (DIFMAX.GT.DIF(J)) GO TO 360
         DIFMAX = DIF(J)
         DVAXIS = J
  360    DIF(J) = ZERO
  380 CONTINUE
      TPONTR = POINTR + 2*(DVAXIS+1)
      WRKSTR(TPONTR) = WIDTH(DVAXIS)*HALF
      WRKSTR(TPONTR-1) = CENTER(DVAXIS) - WRKSTR(TPONTR)
      IF (DVFLAG.NE.2) GO TO 400
      DVAXES = WRKSTR(POINTR+2)
      IF (DVAXES.GT.MAXAXS) DVAXES = -1
      DVAXIS = DVAXIS + (NDIM+1)*DVAXES
  400 WRKSTR(POINTR+2) = DVAXIS
      IF (DVFLAG.EQ.1) GO TO 460
C       DETERMINE THE POSITION IN THE PARIALLY ORDERED LIST OF
C       THE SUBREGION WHICH REPLACES MOST RECENTLY DIVIDED SUBREGION
  420 SUBTMP = 2*SUBRGN
      IF (SUBTMP.GT.SBRGNS) GO TO 480
      TPONTR = WRKSTR(SUBTMP)
      IF (SUBTMP.EQ.SBRGNS) GO TO 440
      TPONTP = WRKSTR(SUBTMP+1)
      IF (WRKSTR(TPONTR).GE.WRKSTR(TPONTP)) GO TO 440
      SUBTMP = SUBTMP + 1
      TPONTR = TPONTP
  440 IF (RGNERT.GE.WRKSTR(TPONTR)) GO TO 480
      WRKSTR(SUBTMP) = POINTR
      WRKSTR(SUBRGN) = TPONTR
      SUBRGN = SUBTMP
      GO TO 420
C       WHEN WORKING STORAGE IS NOT USED UP, DETERMINE THE
C       POSITION IN THE PARTIALLY ORDERED LIST FOR THE DESCRIPTION
C       OF OTHER PORTION(S) OF MOST RECENTLY DIVIDED SUBREGION
  460 SUBTMP = SUBRGN/2
      IF (SUBTMP.LT.1) GO TO 480
      TPONTR = WRKSTR(SUBTMP)
      IF (RGNERT.LE.WRKSTR(TPONTR)) GO TO 480
      WRKSTR(SUBTMP) = POINTR
      WRKSTR(SUBRGN) = TPONTR
      SUBRGN = SUBTMP
      GO TO 460
  480 RGNVLT = ZERO
      RGNERT = ZERO
      IF (DVFLAG.EQ.2) GO TO 540
      DVFLAG = 1 - DVFLAG
C       COUNT TO DETERMINE THE NEXT PART OF THE RECENTLY DIVIDED
C       SUBREGION FOR APPLICATION OF THE BASIC RULE
  500 CENTER(1) = CENTER(1) + TWO*WIDTH(1)
      DVCNTR(1) = DVCNTR(1) + 1
      DO 520 J = 2, NDIM
         IF (DVCNTR(J-1).LE.DVCNTL(J-1)) GO TO 100
         DVCNTR(J-1) = 1
         CENTER(J-1) = OLDCNT(J-1)
         DVCNTR(J) = DVCNTR(J) + 1
         CENTER(J) = CENTER(J) + TWO*WIDTH(J)
  520 CONTINUE
      IF (DVCNTR(NDIM).LE.DVCNTL(NDIM)) GO TO 100
      CENTER(NDIM) = OLDCNT(NDIM)
      IF (DVFLAG.EQ.2) GO TO 340
C
C       END ORDERING OF BASIC RULE RESULTS
C       MAKE CHECKS FOR POSSIBLE TERMINATION OF ROUTINE
C
  540 ACC = ONE
      ABSERR = ZERO
      FINVAL = ZERO
      POINTR = MXRGNS + 1
      TPONTR = 2*NDIM + 3
      DO 560 J = 1, SBRGNS
         ABSERR = ABSERR + WRKSTR(POINTR)
         FINVAL = FINVAL + WRKSTR(POINTR+1)
         POINTR = POINTR + TPONTR
  560 CONTINUE
      IF (ABS(FINVAL).GT.X02AMF()) ACC = ABSERR/ABS(FINVAL)
      IF (ABSERR.EQ.ZERO) ACC = ZERO
      IF (ACC.GT.EPS .OR. FUNCLS.LT.MINPTS) GO TO 40
C
C       LOOP BACK TO APPLY BASIC RULE
C
C       TERMINATION POINT, SET IFAIL AND RETURN
C
      IERROR = 0
      GO TO 640
  580 IERROR = 1
      GO TO 640
  600 IERROR = 2
      GO TO 640
  620 IERROR = 3
  640 MINPTS = FUNCLS
      IFAIL = P01ABF(IFAIL,IERROR,SRNAME,0,P01REC)
      RETURN
      END

