cpevl Subroutine

private subroutine cpevl(n, m, a, z, c, b, kbd)

Evaluate a complex polynomial and its derivatives. Optionally compute error bounds for these values.

REVISION HISTORY (YYMMDD)

  • 810223 DATE WRITTEN
  • 890531 Changed all specific intrinsics to generic. (WRB)
  • 890831 Modified array declarations. (WRB)
  • 891214 Prologue converted to Version 4.0 format. (BAB)
  • 900402 Added TYPE section. (WRB)
  • Jacob Williams, 9/13/2022 : modernized this routine

Arguments

Type IntentOptional Attributes Name
integer, intent(in) :: n

Degree of the polynomial
integer, intent(in) :: m

Number of derivatives to be calculated:

  • M=0 evaluates only the function
  • M=1 evaluates the function and first derivative, etc.

if M > N+1 function and all N derivatives will be calculated.
complex(kind=wp), intent(in) :: a(*)

vector containing the N+1 coefficients of polynomial. A(I) = coefficient of Z**(N+1-I)
complex(kind=wp), intent(in) :: z

point at which the evaluation is to take place
complex(kind=wp), intent(out) :: c(*)

Array of 2(M+1) words: C(I+1) contains the complex value of the I-th derivative at Z, I=0,...,M
complex(kind=wp), intent(out) :: b(*)

Array of 2(M+1) words: B(I) contains the bounds on the real and imaginary parts of C(I) if they were requested. only needed if bounds are to be calculated. It is not used otherwise.
logical, intent(in) :: kbd

A logical variable, e.g. .TRUE. or .FALSE. which is to be set .TRUE. if bounds are to be computed.

Called by

proc~~cpevl~~CalledByGraph proc~cpevl polyroots_module::cpevl proc~cpzero polyroots_module::cpzero proc~cpzero->proc~cpevl proc~rpzero polyroots_module::rpzero proc~rpzero->proc~cpzero

Source Code

subroutine cpevl(n, m, a, z, c, b, kbd)

    implicit none

    integer, intent(in) :: n !! Degree of the polynomial
    integer, intent(in) :: m !! Number of derivatives to be calculated:
                             !!
                             !!  * `M=0` evaluates only the function
                             !!  * `M=1` evaluates the function and first derivative, etc.
                             !!
                             !! if `M > N+1` function and all `N` derivatives will be calculated.
    complex(wp), intent(in) :: a(*) !! vector containing the `N+1` coefficients of polynomial.
                                    !! `A(I)` = coefficient of `Z**(N+1-I)`
    complex(wp), intent(in) :: z !! point at which the evaluation is to take place
    complex(wp), intent(out) :: c(*) !! Array of `2(M+1)` words: `C(I+1)` contains the complex value of the I-th
                                     !! derivative at `Z, I=0,...,M`
    complex(wp), intent(out) :: b(*) !! Array of `2(M+1)` words: `B(I)` contains the bounds on the real and imaginary parts
                                     !! of `C(I)` if they were requested. only needed if bounds are to be calculated.
                                     !! It is not used otherwise.
    logical, intent(in) :: kbd !! A logical variable, e.g. .TRUE. or .FALSE. which is
                               !! to be set .TRUE. if bounds are to be computed.

    real(wp) :: r, s
    integer :: i, j, mini, np1
    complex(wp) :: ci, cim1, bi, bim1, t, za, q

    za(q) = cmplx(abs(real(q, wp)), abs(aimag(q)), wp)
    np1 = n + 1
    do j = 1, np1
        ci = 0.0_wp
        cim1 = a(j)
        bi = 0.0_wp
        bim1 = 0.0_wp
        mini = min(m + 1, n + 2 - j)
        do i = 1, mini
            if (j /= 1) ci = c(i)
            if (i /= 1) cim1 = c(i - 1)
            c(i) = cim1 + z*ci
            if (kbd) then
                if (j /= 1) bi = b(i)
                if (i /= 1) bim1 = b(i - 1)
                t = bi + (3.0_wp*eps + 4.0_wp*eps*eps)*za(ci)
                r = real(za(z)*cmplx(real(t, wp), -aimag(t), wp), wp)
                s = aimag(za(z)*t)
                b(i) = (1.0_wp + 8.0_wp*eps)*(bim1 + eps*za(cim1) + cmplx(r, s, wp))
                if (j == 1) b(i) = 0.0_wp
            end if
        end do
    end do

end subroutine cpevl