start Subroutine

private subroutine start(n, a, y, radius, nz, small, big)

compute the starting approximations of the roots

this routines selects starting approximations along circles center at 0 and having suitable radii. the computation of the number of circles and of the corresponding radii is performed by computing the upper convex hull of the set (i,log(a(i))), i=1,...,n+1.

Arguments

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

number of the coefficients of the polynomial
real(kind=wp), intent(inout) :: a(n+1)

moduli of the coefficients of the polynomial
complex(kind=wp), intent(out) :: y(n)

starting approximations
real(kind=wp), intent(out) :: radius(n)

if a component is -1 then the corresponding root has a too big or too small modulus in order to be represented as double float with no overflow/underflow
integer, intent(out) :: nz

number of roots which cannot be represented without overflow/underflow
real(kind=wp), intent(in) :: small

the min positive real(wp), small=2**(-1074) for the ieee.
real(kind=wp), intent(in) :: big

the max real(wp), big=2**1023 for the ieee standard.

Calls

proc~~start~~CallsGraph proc~start polyroots_module::start proc~cnvex polyroots_module::cnvex proc~start->proc~cnvex proc~cmerge polyroots_module::cmerge proc~cnvex->proc~cmerge proc~ctest polyroots_module::ctest proc~cmerge->proc~ctest proc~left polyroots_module::left proc~cmerge->proc~left proc~right polyroots_module::right proc~cmerge->proc~right

Called by

proc~~start~~CalledByGraph proc~start polyroots_module::start proc~polzeros polyroots_module::polzeros proc~polzeros->proc~start

Source Code

    subroutine start(n, a, y, radius, nz, small, big)

        !! compute the starting approximations of the roots
        !!
        !! this routines selects starting approximations along circles center at
        !! 0 and having suitable radii. the computation of the number of circles
        !! and of the corresponding radii is performed by computing the upper
        !! convex hull of the set (i,log(a(i))), i=1,...,n+1.

        implicit none

        integer,intent(in) :: n !! number of the coefficients of the polynomial
        real(wp),intent(inout) :: a(n + 1) !! moduli of the coefficients of the polynomial
        complex(wp),intent(out) :: y(n) !! starting approximations
        real(wp),intent(out) :: radius(n) !! if a component is -1 then the corresponding root has a
                                          !! too big or too small modulus in order to be represented
                                          !! as double float with no overflow/underflow
        integer,intent(out) :: nz !! number of roots which cannot be represented without
                                  !! overflow/underflow
        real(wp),intent(in) :: small !! the min positive real(wp), small=2**(-1074) for the ieee.
        real(wp),intent(in) :: big !! the max real(wp), big=2**1023 for the ieee standard.

        logical :: h(n + 1) !! auxiliary variable: needed for the computation of the convex hull

        integer :: i, iold, nzeros, j, jj
        real(wp) :: r, th, ang, temp
        real(wp) :: xsmall, xbig

        real(wp),parameter :: pi2 = 2.0_wp * pi
        real(wp),parameter :: sigma = 0.7_wp

        xsmall = log(small)
        xbig = log(big)
        nz = 0
        ! compute the logarithm a(i) of the moduli of the coefficients of
        ! the polynomial and then the upper covex hull of the set (a(i),i)
        do i = 1, n + 1
            if (a(i) /= 0.0_wp) then
                a(i) = log(a(i))
            else
                a(i) = -1.0e30_wp ! maybe replace with -huge(1.0_wp) ?? -JW
            end if
        end do
        call cnvex(n + 1, a, h)
        ! given the upper convex hull of the set (a(i),i) compute the moduli
        ! of the starting approximations by means of rouche's theorem
        iold = 1
        th = pi2/n
        do i = 2, n + 1
            if (h(i)) then
                nzeros = i - iold
                temp = (a(iold) - a(i))/nzeros
                ! check if the modulus is too small
                if ((temp < -xbig) .and. (temp >= xsmall)) then
                    write (*, *) 'warning:', nzeros, ' zero(s) are too small to'
                    write (*, *) 'represent their inverses as complex(wp), they'
                    write (*, *) 'are replaced by small numbers, the corresponding'
                    write (*, *) 'radii are set to -1'
                    nz = nz + nzeros
                    r = 1.0_wp/big
                end if
                if (temp < xsmall) then
                    nz = nz + nzeros
                    write (*, *) 'warning: ', nzeros, ' zero(s) are too small to be'
                    write (*, *) 'represented as complex(wp), they are set to 0'
                    write (*, *) 'the corresponding radii are set to -1'
                end if
                ! check if the modulus is too big
                if (temp > xbig) then
                    r = big
                    nz = nz + nzeros
                    write (*, *) 'warning: ', nzeros, ' zeros(s) are too big to be'
                    write (*, *) 'represented as complex(wp),'
                    write (*, *) 'the corresponding radii are set to -1'
                end if
                if ((temp <= xbig) .and. (temp > max(-xbig, xsmall))) r = exp(temp)
                ! compute nzeros approximations equally distributed in the disk of
                ! radius r
                ang = pi2/nzeros
                do j = iold, i - 1
                    jj = j - iold + 1
                    if ((r <= (1.0_wp/big)) .or. (r == big)) radius(j) = -1.0_wp
                    y(j) = r*(cos(ang*jj + th*i + sigma) + cmplx(0.0_wp, 1.0_wp, wp)*sin(ang*jj + th*i + sigma))
                end do
                iold = i
            end if
        end do

    end subroutine start