*     C02AFF Example Program Text
*     Mark 20 Revised. NAG Copyright 2001.
*     .. Parameters ..
      INTEGER          NOUT
      PARAMETER        (NOUT=6)
*     .. External Subroutines ..
      EXTERNAL         EX1, EX2
*     .. Executable Statements ..
      WRITE (NOUT,*) 'C02AFF Example Program Results'
      CALL EX1
      CALL EX2
      END
*
      SUBROUTINE EX1
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      INTEGER          MAXDEG
      PARAMETER        (MAXDEG=100)
      LOGICAL          SCAL
      PARAMETER        (SCAL=.TRUE.)
*     .. Local Scalars ..
      INTEGER          I, IFAIL, N
*     .. Local Arrays ..
      DOUBLE PRECISION A(2,0:MAXDEG), W(4*MAXDEG+4), Z(2,MAXDEG)
*     .. External Subroutines ..
      EXTERNAL         C02AFF
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Example 1'
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*)
      READ (NIN,*)
      READ (NIN,*) N
      IF (N.GT.0 .AND. N.LE.MAXDEG) THEN
         READ (NIN,*) (A(1,I),A(2,I),I=0,N)
         IFAIL = 1
*
         CALL C02AFF(A,N,SCAL,Z,W,IFAIL)
*
         IF (IFAIL.EQ.0) THEN
            WRITE (NOUT,*)
            WRITE (NOUT,99999) 'Degree of polynomial = ', N
            WRITE (NOUT,*)
            WRITE (NOUT,*) 'Computed roots of polynomial'
            WRITE (NOUT,*)
            DO 20 I = 1, N
               WRITE (NOUT,99998) 'z = ', Z(1,I), Z(2,I), '*i'
   20       CONTINUE
         ELSE
            WRITE (NOUT,*)
            WRITE (NOUT,99997) ' ** C02AFF returned with IFAIL = ',
     +        IFAIL
         END IF
      ELSE
         WRITE (NOUT,*) 'N is out of range'
      END IF
*
99999 FORMAT (1X,A,I4)
99998 FORMAT (1X,A,1P,E12.4,SP,E12.4,A)
99997 FORMAT (1X,A,I5)
      END
*
      SUBROUTINE EX2
*     .. Parameters ..
      INTEGER          NIN, NOUT
      PARAMETER        (NIN=5,NOUT=6)
      DOUBLE PRECISION ZERO, ONE, THREE
      PARAMETER        (ZERO=0.0D0,ONE=1.0D0,THREE=3.0D0)
      INTEGER          MAXDEG
      PARAMETER        (MAXDEG=100)
      LOGICAL          SCAL
      PARAMETER        (SCAL=.TRUE.)
*     .. Local Scalars ..
      DOUBLE PRECISION DELTAC, DELTAI, DI, EPS, EPSBAR, F, R1, R2, R3,
     +                 RMAX
      INTEGER          I, IFAIL, J, JMIN, N
*     .. Local Arrays ..
      DOUBLE PRECISION A(2,0:MAXDEG), ABAR(2,0:MAXDEG), R(MAXDEG),
     +                 W(4*MAXDEG+4), Z(2,MAXDEG), ZBAR(2,MAXDEG)
      INTEGER          M(MAXDEG)
*     .. External Functions ..
      DOUBLE PRECISION A02ABF, X02AJF, X02ALF
      EXTERNAL         A02ABF, X02AJF, X02ALF
*     .. External Subroutines ..
      EXTERNAL         C02AFF
*     .. Intrinsic Functions ..
      INTRINSIC        ABS, MAX, MIN
*     .. Executable Statements ..
      WRITE (NOUT,*)
      WRITE (NOUT,*)
      WRITE (NOUT,*) 'Example 2'
*     Skip heading in data file
      READ (NIN,*)
      READ (NIN,*)
      READ (NIN,*) N
      IF (N.GT.0 .AND. N.LE.MAXDEG) THEN
*
*        Read in the coefficients of the original polynomial.
*
         READ (NIN,*) (A(1,I),A(2,I),I=0,N)
*
*        Compute the roots of the original polynomial.
*
         IFAIL = 1
*
         CALL C02AFF(A,N,SCAL,Z,W,IFAIL)
*
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,*)
            WRITE (NOUT,99997) ' ** C02AFF returned with IFAIL = ',
     +        IFAIL
            GO TO 120
         END IF
*
*        Form the coefficients of the perturbed polynomial.
*
         EPS = X02AJF()
         EPSBAR = THREE*EPS
         DO 20 I = 0, N
            IF (A(1,I).NE.ZERO) THEN
               F = ONE + EPSBAR
               EPSBAR = -EPSBAR
               ABAR(1,I) = F*A(1,I)
               IF (A(2,I).NE.ZERO) THEN
                  ABAR(2,I) = F*A(2,I)
               ELSE
                  ABAR(2,I) = ZERO
               END IF
            ELSE
               ABAR(1,I) = ZERO
               IF (A(2,I).NE.ZERO) THEN
                  F = ONE + EPSBAR
                  EPSBAR = -EPSBAR
                  ABAR(2,I) = F*A(2,I)
               ELSE
                  ABAR(2,I) = ZERO
               END IF
            END IF
   20    CONTINUE
*
*        Compute the roots of the perturbed polynomial.
*
         IFAIL = 1
*
         CALL C02AFF(ABAR,N,SCAL,ZBAR,W,IFAIL)
*
         IF (IFAIL.NE.0) THEN
            WRITE (NOUT,*)
            WRITE (NOUT,99997) ' ** C02AFF returned with IFAIL = ',
     +        IFAIL
            GO TO 120
         END IF
*
*        Perform error analysis.
*
         DO 40 I = 1, N
*           Initialize markers to 0 (unmarked).
            M(I) = 0
   40    CONTINUE
         RMAX = X02ALF()
*        Loop over all unperturbed roots (stored in Z).
         DO 80 I = 1, N
            DELTAI = RMAX
            R1 = A02ABF(Z(1,I),Z(2,I))
*           Loop over all perturbed roots (stored in ZBAR).
            DO 60 J = 1, N
*              Compare the current unperturbed root to all unmarked
*              perturbed roots.
               IF (M(J).EQ.0) THEN
                  R2 = A02ABF(ZBAR(1,J),ZBAR(2,J))
                  DELTAC = ABS(R1-R2)
                  IF (DELTAC.LT.DELTAI) THEN
                     DELTAI = DELTAC
                     JMIN = J
                  END IF
               END IF
   60       CONTINUE
*           Mark the selected perturbed root.
            M(JMIN) = 1
*           Compute the relative error.
            IF (R1.NE.ZERO) THEN
               R3 = A02ABF(ZBAR(1,JMIN),ZBAR(2,JMIN))
               DI = MIN(R1,R3)
               R(I) = MAX(DELTAI/MAX(DI,DELTAI/RMAX),EPS)
            ELSE
               R(I) = ZERO
            END IF
   80    CONTINUE
*
         WRITE (NOUT,*)
         WRITE (NOUT,99999) 'Degree of polynomial = ', N
         WRITE (NOUT,*)
         WRITE (NOUT,*) 'Computed roots of polynomial   ',
     +     '  Error estimates'
         WRITE (NOUT,*) '                            ',
     +     '   (machine-dependent)'
         WRITE (NOUT,*)
         DO 100 I = 1, N
            WRITE (NOUT,99998) 'z = ', Z(1,I), Z(2,I), '*i', R(I)
  100    CONTINUE
      ELSE
         WRITE (NOUT,*) 'N is out of range'
      END IF
  120 CONTINUE
*
99999 FORMAT (1X,A,I4)
99998 FORMAT (1X,A,1P,E12.4,SP,E12.4,A,5X,SS,E9.1)
99997 FORMAT (1X,A,I5)
      END