cpolz Subroutine

public subroutine cpolz(a, ndeg, z, ierr)

In the discussion below, the notation A([*,],k} should be interpreted as the complex number A(k) if A is declared complex, and should be interpreted as the complex number A(1,k) + i * A(2,k) if A is not declared to be of type complex. Similar statements apply for Z(k).

Given complex coefficients A([,[1),...,A([,]NDEG+1) this subr computes the NDEG roots of the polynomial A([,]1)XNDEG + ... + A([,]NDEG+1) storing the roots as complex numbers in the array Z( ). Require NDEG >= 1 and A([,]1) /= 0.

Reference

License

Copyright (c) 1996 California Institute of Technology, Pasadena, CA. ALL RIGHTS RESERVED. Based on Government Sponsored Research NAS7-03001.

History

  • C.L.Lawson & S.Y.Chan, JPL, June 3,1986.
  • 1987-02-25 CPOLZ Lawson Initial code.
  • 1989-10-20 CLL Delcared all variables.
  • 1992-05-11 CLL IERR was not being set when N = 0 or 1. Fixed this.
  • 1995-01-18 CPOLZ Krogh More M77CON for conversion to C.
  • 1995-11-17 CPOLZ Krogh Added M77CON statements for conversion to C
  • 1995-11-17 CPOLZ Krogh Converted SFTRAN to Fortran 77.
  • 1996-03-30 CPOLZ Krogh Added external statement.
  • 1996-04-27 CPOLZ Krogh Changes to use .C. and C%%.
  • 2001-05-25 CPOLZ Krogh Minor change for making .f90 version.
  • 2022-10-06, Jacob Williams modernized this routine

Arguments

Type IntentOptional Attributes Name
complex(kind=wp), intent(in) :: a(ndeg+1)

contains the complex coefficients of a polynomial high order coefficient first, with a([,]1)/=0. the real and imaginary parts of the jth coefficient must be provided in a([],j). the contents of this array will not be modified by the subroutine.
integer, intent(in) :: ndeg

degree of the polynomial
complex(kind=wp), intent(out) :: z(ndeg)

contains the polynomial roots stored as complex numbers. the real and imaginary parts of the jth root will be stored in z([*,]j).
integer, intent(out) :: ierr

error flag. set by the subroutine to 0 on normal termination. set to -1 if a([*,]1)=0. set to -2 if ndeg <= 0. set to j > 0 if the iteration count limit has been exceeded and roots 1 through j have not been determined.

Calls

proc~~cpolz~~CallsGraph proc~cpolz polyroots_module::cpolz proc~scomqr polyroots_module::scomqr proc~cpolz->proc~scomqr

Source Code

    subroutine cpolz(a,ndeg,z,ierr)

    integer,intent(in) :: ndeg !! degree of the polynomial
    complex(wp),intent(in) :: a(ndeg+1) !! contains the complex coefficients of a polynomial
                                        !! high order coefficient first, with a([*,]1)/=0. the
                                        !! real and imaginary parts of the jth coefficient must
                                        !! be provided in a([*],j). the contents of this array will
                                        !! not be modified by the subroutine.
    complex(wp),intent(out) :: z(ndeg) !! contains the polynomial roots stored as complex
                                       !! numbers.  the real and imaginary parts of the jth root
                                       !! will be stored in z([*,]j).
    integer,intent(out) :: ierr !! error flag. set by the subroutine to 0 on normal
                                !! termination. set to -1 if a([*,]1)=0. set to -2 if ndeg
                                !! <= 0. set to  j > 0 if the iteration count limit
                                !! has been exceeded and roots 1 through j have not been
                                !! determined.

    complex(wp) :: temp
    integer :: i, j, n
    real(wp) :: c, f, g, r, s
    logical :: more, first
    real(wp) :: h(ndeg,ndeg,2) !! array of work space

    real(wp),parameter :: zero = 0.0_wp
    real(wp),parameter :: one = 1.0_wp
    real(wp),parameter :: c95 = 0.95_wp
    integer,parameter :: base = radix(1.0_wp)   !! i1mach(10)
    integer,parameter :: b2 = base * base

    if (ndeg <= 0) then
        ierr = -2
        write(*,*) 'ndeg <= 0'
        return
    end if

    if (a(1) == cmplx(zero, zero, wp)) then
        ierr = -1
        write(*,*) 'a(*,1) == zero'
        return
    end if

    n = ndeg
    ierr = 0

    ! build first row of companion matrix.

    do i = 2,n+1
        temp = -(a(i)/a(1))
        h(1,i-1,1) = real(temp,wp)
        h(1,i-1,2) = aimag(temp)
    end do

    ! extract any exact zero roots and set n = degree of
    ! remaining polynomial.

    do j = ndeg,1,-1
        if (h(1,j,1)/=zero .or. h(1,j,2)/=zero) exit
        z(j) = zero
        n = n - 1
    end do

    ! special for n = 0 or 1.

    if (n == 0) return
    if (n == 1) then
        z(1) = cmplx(h(1,1,1),h(1,1,2),wp)
        return
    end if

    ! build rows 2 thru n of the companion matrix.

    do i = 2,n
        do j = 1,n
            if (j == i-1) then
                h(i,j,1) = one
                h(i,j,2) = zero
            else
                h(i,j,1) = zero
                h(i,j,2) = zero
            end if
        end do
    end do

    ! ***************** balance the matrix ***********************
    !
    !     this is an adaption of the eispack subroutine balanc to
    !     the special case of a complex companion matrix. the eispack
    !     balance is a translation of the algol procedure balance,
    !     num. math. 13, 293-304(1969) by parlett and reinsch.
    !     handbook for auto. comp., vol.ii-linear algebra, 315-326(1971).

    ! ********** iterative loop for norm reduction **********
    do
        more = .false.
        do i = 1, n
          ! compute r = sum of magnitudes in row i skipping diagonal.
          !         c = sum of magnitudes in col i skipping diagonal.
          if (i == 1) then
            r = abs(h(1,2,1)) + abs(h(1,2,2))
            do j = 3,n
              r = r + abs(h(1,j,1)) + abs(h(1,j,2))
            end do
            c = abs(h(2,1,1)) + abs(h(2,1,2))
          else
            r = abs(h(i,i-1,1)) + abs(h(i,i-1,2))
            c = abs(h(1,i,1)) + abs(h(1,i,2))
            if (i /= n) then
              c = c + abs(h(i+1,i,1)) + abs(h(i+1,i,2))
            end if
          end if

          ! determine column scale factor, f.

          g = r / base
          f = one
          s = c + r

          do
            if (c >= g) exit
            f = f * base
            c = c * b2
          end do
          g = r * base
          do
            if (c < g) exit
            f = f / base
            c = c / b2
          end do
          ! will the factor f have a significant effect ?

          if ((c + r) / f < c95 * s) then

            ! yes, so do the scaling.

            g = one / f
            more = .true.

            ! scale row i

            if (i == 1) then
              do j = 1,n
                h(1,j,1) = h(1,j,1)*g
                h(1,j,2) = h(1,j,2)*g
              end do
            else
              h(i,i-1,1) = h(i,i-1,1)*g
              h(i,i-1,2) = h(i,i-1,2)*g
            end if

            ! scale column i

            h(1,i,1) = h(1,i,1) * f
            h(1,i,2) = h(1,i,2) * f
            if (i /= n) then
              h(i+1,i,1) = h(i+1,i,1) * f
              h(i+1,i,2) = h(i+1,i,2) * f
            end if

          end if
        end do
        if (.not. more) exit
    end do

    call scomqr(ndeg,n,1,n,h(1,1,1),h(1,1,2),z,ierr)

    if (ierr /= 0) write(*,*) 'Convergence failure in cpolz'

    end subroutine cpolz