cpqr79 Subroutine

public subroutine cpqr79(ndeg, coeff, root, ierr)

This routine computes all zeros of a polynomial of degree NDEG with complex coefficients by computing the eigenvalues of the companion matrix.

REVISION HISTORY (YYMMDD)

  • 791201 DATE WRITTEN. Vandevender, W. H., (SNLA)
  • 890531 Changed all specific intrinsics to generic. (WRB)
  • 890531 REVISION DATE from Version 3.2
  • 891214 Prologue converted to Version 4.0 format. (BAB)
  • 900315 CALLs to XERROR changed to CALLs to XERMSG. (THJ)
  • 900326 Removed duplicate information from DESCRIPTION section. (WRB)
  • 911010 Code reworked and simplified. (RWC and WRB)
  • Jacob Williams, 9/14/2022 : modernized this code

Arguments

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

degree of polynomial
complex(kind=wp), intent(in), dimension(ndeg+1) :: coeff

(NDEG+1) coefficients in descending order. i.e., P(Z)= COEFF(1)*(Z**NDEG) + COEFF(NDEG)*Z + COEFF(NDEG+1)
complex(kind=wp), intent(out), dimension(ndeg) :: root

(NDEG) vector of roots
integer, intent(out) :: ierr

Output Error Code.

Normal Code:

  • 0 -- means the roots were computed.

Abnormal Codes:

  • 1 -- more than 30 QR iterations on some eigenvalue of the companion matrix
  • 2 -- COEFF(1)=0.0
  • 3 -- NDEG is invalid (less than or equal to 0)

Calls

proc~~cpqr79~~CallsGraph proc~cpqr79 polyroots_module::cpqr79 proc~comqr polyroots_module::comqr proc~cpqr79->proc~comqr proc~cdiv polyroots_module::cdiv proc~comqr->proc~cdiv proc~csroot polyroots_module::csroot proc~comqr->proc~csroot proc~pythag polyroots_module::pythag proc~comqr->proc~pythag proc~csroot->proc~pythag

Source Code

subroutine cpqr79(ndeg, coeff, root, ierr)
    implicit none

    integer, intent(in) :: ndeg !! degree of polynomial
    complex(wp), dimension(ndeg+1), intent(in) :: coeff !! `(NDEG+1)` coefficients in descending order.  i.e.,
                                                        !! `P(Z)= COEFF(1)*(Z**NDEG) + COEFF(NDEG)*Z + COEFF(NDEG+1)`
    complex(wp), dimension(ndeg), intent(out) :: root !! `(NDEG)` vector of roots
    integer, intent(out) :: ierr !! Output Error Code.
                                 !!
                                 !!### Normal Code:
                                 !!
                                 !!  * 0 -- means the roots were computed.
                                 !!
                                 !!### Abnormal Codes:
                                 !!
                                 !!  * 1 --  more than 30 QR iterations on some eigenvalue of the companion matrix
                                 !!  * 2 --  COEFF(1)=0.0
                                 !!  * 3 --  NDEG is invalid (less than or equal to 0)

    complex(wp) :: scale, c
    integer :: k, khr, khi, kwr, kwi, kad, kj, km1
    real(wp), dimension(:), allocatable :: work !! work array of dimension `2*NDEG*(NDEG+1)`

    ierr = 0
    if (abs(coeff(1)) == 0.0_wp) then
        ierr = 2
        write (*, *) 'leading coefficient is zero.'
        return
    end if

    if (ndeg <= 0) then
        ierr = 3
        write (*, *) 'degree invalid.'
        return
    end if

    if (ndeg == 1) then
        root(1) = -coeff(2)/coeff(1)
        return
    end if

    ! allocate work array:
    allocate (work(2*NDEG*(NDEG + 1)))

    scale = 1.0_wp/coeff(1)
    khr = 1
    khi = khr + ndeg*ndeg
    kwr = khi + khi - khr
    kwi = kwr + ndeg

    do k = 1, kwr
        work(k) = 0.0_wp
    end do

    do k = 1, ndeg
        kad = (k - 1)*ndeg + 1
        c = scale*coeff(k + 1)
        work(kad) = -real(c, wp)
        kj = khi + kad - 1
        work(kj) = -aimag(c)
        if (k /= ndeg) work(kad + k) = 1.0_wp
    end do

    call comqr(ndeg, ndeg, 1, ndeg, work(khr), work(khi), work(kwr), work(kwi), ierr)

    if (ierr /= 0) then
        write (*, *) 'no convergence in 30 qr iterations. ierr = ', ierr
        ierr = 1
        return
    end if

    do k = 1, ndeg
        km1 = k - 1
        root(k) = cmplx(work(kwr + km1), work(kwi + km1), wp)
    end do

end subroutine cpqr79