!   C02AFF Example Program Text
!   Mark 26.1 Release. NAG Copyright 2016.

    Module c02affe_mod

!     C02AFF Example Program Module:
!            Parameters

!     .. Implicit None Statement ..
      Implicit None
!     .. Accessibility Statements ..
      Private
!     .. Parameters ..
      Integer, Parameter, Public       :: nin = 5, nout = 6
      Logical, Parameter, Public       :: scal = .True.
    End Module c02affe_mod
    Program c02affe

!     C02AFF Example Main Program

!     .. Use Statements ..
      Use c02affe_mod, Only: nout
!     .. Implicit None Statement ..
      Implicit None
!     .. Executable Statements ..
      Write (nout,*) 'C02AFF Example Program Results'

      Call ex1

      Call ex2

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use c02affe_mod, Only: nin, scal
        Use nag_library, Only: c02aff, nag_wp
!       .. Local Scalars ..
        Integer                        :: i, ifail, n
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: a(:,:), w(:), z(:,:)
!       .. Executable Statements ..
        Write (nout,*)
        Write (nout,*)
        Write (nout,*) 'Example 1'

!       Skip heading in data file
        Read (nin,*)
        Read (nin,*)
        Read (nin,*)

        Read (nin,*) n
        Allocate (a(2,0:n),w(4*(n+1)),z(2,n))

        Read (nin,*)(a(1,i),a(2,i),i=0,n)

        ifail = 0
        Call c02aff(a,n,scal,z,w,ifail)

        Write (nout,*)
        Write (nout,99999) 'Degree of polynomial = ', n
        Write (nout,*)
        Write (nout,*) 'Computed roots of polynomial'
        Write (nout,*)

        Do i = 1, n
          Write (nout,99998) 'z = ', z(1,i), z(2,i), '*i'
        End Do

99999   Format (1X,A,I4)
99998   Format (1X,A,1P,E12.4,Sp,E12.4,A)
      End Subroutine ex1
      Subroutine ex2

!       .. Use Statements ..
        Use c02affe_mod, Only: nin, scal
        Use nag_library, Only: a02abf, c02aff, nag_wp, x02ajf, x02alf
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: deltac, deltai, di, eps, epsbar, f,  &
                                          r1, r2, r3, rmax
        Integer                        :: i, ifail, j, jmin, n
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: a(:,:), abar(:,:), r(:), w(:),      &
                                          z(:,:), zbar(:,:)
        Integer, Allocatable           :: m(:)
!       .. Intrinsic Procedures ..
        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
        Allocate (a(2,0:n),abar(2,0:n),r(n),w(4*(n+1)),z(2,n),zbar(2,n),m(n))

!       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 = 0
        Call c02aff(a,n,scal,z,w,ifail)

!       Form the coefficients of the perturbed polynomial.

        eps = x02ajf()
        epsbar = 3.0E0_nag_wp*eps

        Do i = 0, n

          If (a(1,i)/=0.0E0_nag_wp) Then
            f = 1.0E0_nag_wp + epsbar
            epsbar = -epsbar
            abar(1,i) = f*a(1,i)

            If (a(2,i)/=0.0E0_nag_wp) Then
              abar(2,i) = f*a(2,i)
            Else
              abar(2,i) = 0.0E0_nag_wp
            End If

          Else
            abar(1,i) = 0.0E0_nag_wp

            If (a(2,i)/=0.0E0_nag_wp) Then
              f = 1.0E0_nag_wp + epsbar
              epsbar = -epsbar
              abar(2,i) = f*a(2,i)
            Else
              abar(2,i) = 0.0E0_nag_wp
            End If
          End If

        End Do

!       Compute the roots of the perturbed polynomial.

        ifail = 0
        Call c02aff(abar,n,scal,zbar,w,ifail)

!       Perform error analysis.

!       Initialize markers to 0 (unmarked).

        m(1:n) = 0

        rmax = x02alf()

!       Loop over all unperturbed roots (stored in Z).

        Do i = 1, n
          deltai = rmax
          r1 = a02abf(z(1,i),z(2,i))

!         Loop over all perturbed roots (stored in ZBAR).

          Do j = 1, n

!           Compare the current unperturbed root to all unmarked
!           perturbed roots.

            If (m(j)==0) Then
              r2 = a02abf(zbar(1,j),zbar(2,j))
              deltac = abs(r1-r2)

              If (deltac<deltai) Then
                deltai = deltac
                jmin = j
              End If

            End If

          End Do

!         Mark the selected perturbed root.

          m(jmin) = 1

!         Compute the relative error.

          If (r1/=0.0E0_nag_wp) 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) = 0.0E0_nag_wp
          End If

        End Do

        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 i = 1, n
          Write (nout,99998) 'z = ', z(1,i), z(2,i), '*i', r(i)
        End Do

99999   Format (1X,A,I4)
99998   Format (1X,A,1P,E12.4,Sp,E12.4,A,5X,Ss,E9.1)
      End Subroutine ex2
    End Program c02affe