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

    Module c02agfe_mod

!     C02AGF 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 c02agfe_mod
    Program c02agfe

!     C02AGF Example Main Program

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

      Call ex1

      Call ex2

    Contains
      Subroutine ex1

!       .. Use Statements ..
        Use c02agfe_mod, Only: nin, scal
        Use nag_library, Only: c02agf, nag_wp
!       .. Local Scalars ..
        Real (Kind=nag_wp)             :: zi, zr
        Integer                        :: i, ifail, n, nroot
!       .. Local Arrays ..
        Real (Kind=nag_wp), Allocatable :: a(:), w(:), z(:,:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs
!       .. 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(0:n),w(2*(n+1)),z(2,n))

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

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

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

        Write (nout,99998) 'Computed roots of polynomial'

        nroot = 1

        Do While (nroot<=n)

          zr = z(1,nroot)
          zi = z(2,nroot)
          If (zi==0.0E0_nag_wp) Then
            Write (nout,99997) 'z = ', zr
            nroot = nroot + 1
          Else
            Write (nout,99997) 'z = ', zr, ' +/- ', abs(zi), '*i'
            nroot = nroot + 2
          End If

        End Do

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

!       .. Use Statements ..
        Use c02agfe_mod, Only: nin, scal
        Use nag_library, Only: a02abf, c02agf, 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(0:n),abar(0:n),r(n),w(2*(n+1)),z(2,n),zbar(2,n),m(n))

!       Read in the coefficients of the original polynomial.

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

!       Compute the roots of the original polynomial.

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

!       Form the coefficients of the perturbed polynomial.

        eps = x02ajf()
        epsbar = 3.0_nag_wp*eps

        Do i = 0, n

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

        End Do

!       Compute the roots of the perturbed polynomial.

        ifail = 0
        Call c02agf(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.0_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 c02agfe