rpqr79 Subroutine

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

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

REVISION HISTORY (YYMMDD)

  • 800601 DATE WRITTEN. Vandevender, W. H., (SNLA)
  • 890505 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)
  • 911010 Code reworked and simplified. (RWC and WRB)
  • Jacob Williams, 9/13/2022 : modernized this routine

Arguments

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

degree of polynomial
real(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~~rpqr79~~CallsGraph proc~rpqr79 polyroots_module::rpqr79 proc~hqr polyroots_module::hqr proc~rpqr79->proc~hqr

Source Code

subroutine rpqr79(ndeg, coeff, root, ierr)

    implicit none

    integer, intent(in) :: ndeg !! degree of polynomial
    real(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)

    real(wp) :: scale
    integer :: k, kh, kwr, kwi, kcol, km1, kwend
    real(wp), dimension(:), allocatable :: work !! work array of dimension `NDEG*(NDEG+2)`

    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) = cmplx(-coeff(2)/coeff(1), 0.0_wp, wp)
        return
    end if

    allocate (work(ndeg*(ndeg + 2))) ! work array

    scale = 1.0_wp/coeff(1)
    kh = 1
    kwr = kh + ndeg*ndeg
    kwi = kwr + ndeg
    kwend = kwi + ndeg - 1

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

    do k = 1, ndeg
        kcol = (k - 1)*ndeg + 1
        work(kcol) = -coeff(k + 1)*scale
        if (k /= ndeg) work(kcol + k) = 1.0_wp
    end do

    call hqr(ndeg, ndeg, 1, ndeg, work(kh), work(kwr), work(kwi), ierr)

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

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

end subroutine rpqr79