Program c02aafe

!     C02AAF Example Program Text
!     Mark 30.1 Release. NAG Copyright 2024.

!     .. Use Statements ..
      Use nag_library, Only: nag_wp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter               :: nin = 5, nout = 6
      Logical, Parameter               :: polish_example = .False.
!     .. Executable Statements ..

      Write (nout,*) 'C02AAF Example Program Results'

      Call ex1_basic
      If (polish_example) Then
        Call ex2_polishing
      End If

    Contains

      Subroutine ex1_basic

!       .. Use Statements ..
        Use nag_library, Only: c02aaf, x02ajf
!       .. Implicit None Statement ..
        Implicit None
!       .. Local Scalars ..
        Integer                        :: i, ifail, itmax, n, polish
!       .. Local Arrays ..
        Complex (Kind=nag_wp), Allocatable :: a(:), z(:)
        Real (Kind=nag_wp), Allocatable :: berr(:), cond(:)
        Integer, Allocatable           :: conv(:)
!       .. Executable Statements ..

        Write (nout,*) ''
        Write (nout,*) 'Basic Problem'
        Write (nout,*) ''

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

!       Read polynomial degree and allocate
        Read (nin,*) n
        Allocate (a(0:n),berr(n),cond(n),conv(n),z(n))

!       Read polynomial coefficients
        Do i = 0, n
          Read (nin,*) a(i)
        End Do

!       Find roots of the polynomial
        itmax = 30
        polish = 1
        ifail = 0
        Call c02aaf(a,n,itmax,polish,z,berr,cond,conv,ifail)

!       Print output
        Write (nout,*) ' i    z                      conv  berr      cond'
        Write (nout,*) '-----------------------------------------------------'
        Do i = 1, n
          If (berr(i)<x02ajf()) Then
            Write (nout,99999) i, z(i), conv(i), cond(i)
          Else
            Write (nout,99998) i, z(i), conv(i), berr(i), cond(i)
          End If
        End Do

        Deallocate (a,berr,cond,conv,z)

99999   Format (1X,1P,I2,'  (',E10.2,', ',E9.2,')',2X,I3,2X,'< epsilon',1X,    &
          E9.2)
99998   Format (1X,1P,I2,'  (',E10.2,', ',E9.2,')',2X,I3,2X,E9.2,1X,E9.2)

      End Subroutine ex1_basic

      Subroutine ex2_polishing

!       .. Use Statements ..
        Use nag_library, Only: c02aaf, nag_wp, x02ajf, x02alf
!       .. Implicit None Statement ..
        Implicit None
!       .. Parameters ..
        Complex (Kind=nag_wp), Parameter :: one = (1.0_nag_wp,0.0_nag_wp)
        Complex (Kind=nag_wp), Parameter :: zero = (0.0_nag_wp,0.0_nag_wp)
!       .. Local Scalars ..
        Complex (Kind=nag_wp)          :: pz
        Real (Kind=nag_wp)             :: delta, eps, err, fwderr, maxfwderr,  &
                                          maxrelerr, relerr, rmax
        Integer                        :: i, ifail, itmax, j, k, n, polish
!       .. Local Arrays ..
        Complex (Kind=nag_wp), Allocatable :: a(:), z(:), zact(:)
        Real (Kind=nag_wp), Allocatable :: berr(:), cond(:)
        Integer, Allocatable           :: conv(:)
        Logical, Allocatable           :: matched(:)
!       .. Intrinsic Procedures ..
        Intrinsic                      :: abs, cmplx, max
!       .. Executable Statements ..

        Write (nout,*) ''
        Write (nout,*) 'Polishing Processes'
        Write (nout,*) ''

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

!       Read polynomial degree and allocate
        Read (nin,*) n
        Allocate (a(0:n),berr(n),cond(n),conv(n),matched(n),z(n),zact(n))

!       Set known roots
        zact(1:n) = (/(cmplx(i,0,kind=nag_wp),i=1,n)/)

!       Multiply out (z-1)(z-2)...(z-n) for coefficients
        a(0:n-1) = zero
        a(n) = one
        Do i = 1, n
          Do j = 0, n - 1
            a(j) = a(j+1) - a(j)*zact(i)
          End Do
          a(n) = -a(n)*zact(i)
        End Do

        Write (nout,*) ' polish  relerr    fwderr'
        Write (nout,*) '----------------------------'

!       Use different polish modes
        Do polish = 0, 2

          itmax = 30
          eps = x02ajf()
          rmax = x02alf()

!         Find roots
          ifail = 0
          Call c02aaf(a,n,itmax,polish,z,berr,cond,conv,ifail)

!         Calculate the maximum relative errors of the roots, and the maximum
!         forward error evaluating the polynomial at those roots. Errors are
!         capped at machine precision.
          maxrelerr = eps
          maxfwderr = eps
          matched(:) = .False.

          Do i = 1, n

!           Evaluate polynomial at this root
            pz = a(0)
            Do j = 1, n
              pz = z(i)*pz + a(j)
            End Do

!           Match to an expected root
            k = 0
            err = rmax
            Do j = 1, n
              If (.Not. matched(j)) Then
                delta = abs(z(i)-zact(j))
                If (delta<=err) Then
                  err = delta
                  k = j
                End If
              End If
            End Do

!           Mark as matched and update max errors
            matched(k) = .True.
            relerr = err/abs(zact(k))
            fwderr = abs(pz)
            maxrelerr = max(maxrelerr,relerr)
            maxfwderr = max(maxfwderr,fwderr)

          End Do

!         Print output
          Write (nout,99999) polish, maxrelerr, maxfwderr

        End Do

!       Deallocate
        Deallocate (a,berr,cond,conv,matched,z,zact)

99999   Format (1X,1P,I2,5X,E10.2,E10.2)

      End Subroutine ex2_polishing

    End Program c02aafe