cpzero Subroutine

subroutine cpzero(in, a, r, iflg, s)

    implicit none

    integer, intent(in) :: in !! `N`, the degree of `P(Z)`
    complex(wp), dimension(in+1), intent(in) :: a !! complex vector containing coefficients of `P(Z)`,
                                                  !! `A(I)` = coefficient of `Z**(N+1-i)`
    complex(wp), dimension(in), intent(inout) :: r !! `N` word complex vector. On input: containing initial
                                                   !! estimates for zeros if these are known. On output: Ith zero
    integer, intent(inout) :: iflg !!### On Input:
                                   !!
                                   !! flag to indicate if initial estimates of zeros are input:
                                   !!
                                   !!  * If `IFLG == 0`, no estimates are input.
                                   !!  * If `IFLG /= 0`, the vector `R` contains estimates of the zeros
                                   !!
                                   !! ** WARNING ****** If estimates are input, they must
                                   !!                   be separated, that is, distinct or
                                   !!                   not repeated.
                                   !!### On Output:
                                   !!
                                   !! Error Diagnostics:
                                   !!
                                   !! * If `IFLG == 0` on return, all is well
                                   !! * If `IFLG == 1` on return, `A(1)=0.0` or `N=0` on input
                                   !! * If `IFLG == 2` on return, the program failed to converge
                                   !!   after `25*N` iterations.  Best current estimates of the
                                   !!   zeros are in `R(I)`.  Error bounds are not calculated.
    real(wp), intent(out) :: s(in) !! an `N` word array. `S(I)` = bound for `R(I)`

    integer :: i, imax, j, n, n1, nit, nmax, nr
    real(wp) :: u, v, x
    complex(wp) :: pn, temp
    complex(wp) :: ctmp(1), btmp(1)
    complex(wp), dimension(:), allocatable :: t !! `4(N+1)` word array used for temporary storage

    if (in <= 0 .or. abs(a(1)) == 0.0_wp) then
        iflg = 1
    else
        ! work array:
        allocate(t(4*(in+1)))
        ! check for easily obtained zeros
        n = in
        n1 = n + 1
        if (iflg == 0) then
            do
                n1 = n + 1
                if (n <= 1) then
                    r(1) = -a(2)/a(1)
                    s(1) = 0.0_wp
                    return
                elseif (abs(a(n1)) /= 0.0_wp) then
                    ! if initial estimates for zeros not given, find some
                    temp = -a(2)/(a(1)*n)
                    call cpevl(n, n, a, temp, t, t, .false.)
                    imax = n + 2
                    t(n1) = abs(t(n1))
                    do i = 2, n1
                        t(n + i) = -abs(t(n + 2 - i))
                        if (real(t(n + i), wp) < real(t(imax), wp)) imax = n + i
                    end do
                    x = (-real(t(imax), wp)/real(t(n1), wp))**(1.0_wp/(imax - n1))
                    do
                        x = 2.0_wp*x
                        call cpevl(n, 0, t(n1), cmplx(x, 0.0_wp, wp), ctmp, btmp, .false.)
                        pn = ctmp(1)
                        if (real(pn, wp) >= 0.0_wp) exit
                    end do
                    u = 0.5_wp*x
                    v = x
                    do
                        x = 0.5_wp*(u + v)
                        call cpevl(n, 0, t(n1), cmplx(x, 0.0_wp, wp), ctmp, btmp, .false.)
                        pn = ctmp(1)
                        if (real(pn, wp) > 0.0_wp) v = x
                        if (real(pn, wp) <= 0.0_wp) u = x
                        if ((v - u) <= 0.001_wp*(1.0_wp + v)) exit
                    end do
                    do i = 1, n
                        u = (pi/n)*(2*i - 1.5_wp)
                        r(i) = max(x, 0.001_wp*abs(temp))*cmplx(cos(u), sin(u), wp) + temp
                    end do
                    exit
                else
                    r(n) = 0.0_wp
                    s(n) = 0.0_wp
                    n = n - 1
                end if
            end do
        end if

        ! main iteration loop starts here
        nr = 0
        nmax = 25*n
        do nit = 1, nmax
            do i = 1, n
                if (nit == 1 .or. abs(t(i)) /= 0.0_wp) then
                    call cpevl(n, 0, a, r(i), ctmp, btmp, .true.)
                    pn = ctmp(1)
                    temp = btmp(1)
                    if (abs(real(pn, wp)) + abs(aimag(pn)) > real(temp, wp) + aimag(temp)) then
                        temp = a(1)
                        do j = 1, n
                            if (j /= i) temp = temp*(r(i) - r(j))
                        end do
                        t(i) = pn/temp
                    else
                        t(i) = 0.0_wp
                        nr = nr + 1
                    end if
                end if
            end do
            do i = 1, n
                r(i) = r(i) - t(i)
            end do
            if (nr == n) then
                ! calculate error bounds for zeros
                do nr = 1, n
                    call cpevl(n, n, a, r(nr), t, t(n + 2), .true.)
                    x = abs(cmplx(abs(real(t(1), wp)), abs(aimag(t(1))), wp) + t(n + 2))
                    s(nr) = 0.0_wp
                    do i = 1, n
                        x = x*real(n1 - i, wp)/i
                        temp = cmplx(max(abs(real(t(i + 1), wp)) - real(t(n1 + i), wp), 0.0_wp), &
                                     max(abs(aimag(t(i + 1))) - aimag(t(n1 + i)), 0.0_wp), wp)
                        s(nr) = max(s(nr), (abs(temp)/x)**(1.0_wp/i))
                    end do
                    s(nr) = 1.0_wp/s(nr)
                end do
                return
            end if
        end do
        iflg = 2  ! error exit
    end if

end subroutine cpzero