comqr Subroutine

private subroutine comqr(nm, n, low, igh, hr, hi, wr, wi, ierr)

this subroutine finds the eigenvalues of a complex upper hessenberg matrix by the qr method.

Notes

  • calls cdiv for complex division.
  • calls csroot for complex square root.
  • calls pythag for sqrt(aa + bb) .

Reference

  • this subroutine is a translation of a unitary analogue of the algol procedure comlr, num. math. 12, 369-376(1968) by martin and wilkinson. handbook for auto. comp., vol.ii-linear algebra, 396-403(1971). the unitary analogue substitutes the qr algorithm of francis (comp. jour. 4, 332-345(1962)) for the lr algorithm.

History

  • From EISPACK. this version dated august 1983. questions and comments should be directed to burton s. garbow, mathematics and computer science div, argonne national laboratory
  • Jacob Williams, 9/14/2022 : modernized this code

Arguments

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

row dimension of two-dimensional array parameters
integer, intent(in) :: n

the order of the matrix
integer, intent(in) :: low

integer determined by the balancing subroutine cbal. if cbal has not been used, set low=1
integer, intent(in) :: igh

integer determined by the balancing subroutine cbal. if cbal has not been used, igh=n.
real(kind=wp), intent(inout) :: hr(nm,n)

On input: hr and hi contain the real and imaginary parts, respectively, of the complex upper hessenberg matrix. their lower triangles below the subdiagonal contain information about the unitary transformations used in the reduction by corth, if performed.

On Output: the upper hessenberg portions of hr and hi have been destroyed. therefore, they must be saved before calling comqr if subsequent calculation of eigenvectors is to be performed.
real(kind=wp), intent(inout) :: hi(nm,n)

See hr description
real(kind=wp), intent(out) :: wr(n)

the real parts of the eigenvalues. if an error exit is made, the eigenvalues should be correct for indices ierr+1,...,n.
real(kind=wp), intent(out) :: wi(n)

the imaginary parts of the eigenvalues. if an error exit is made, the eigenvalues should be correct for indices ierr+1,...,n.
integer, intent(out) :: ierr

is set to:

  • 0 -- for normal return
  • j -- if the limit of 30*n iterations is exhausted while the j-th eigenvalue is being sought.

Calls

proc~~comqr~~CallsGraph proc~comqr polyroots_module::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

Called by

proc~~comqr~~CalledByGraph proc~comqr polyroots_module::comqr proc~cpqr79 polyroots_module::cpqr79 proc~cpqr79->proc~comqr

Source Code

subroutine comqr(nm, n, low, igh, hr, hi, wr, wi, ierr)
    implicit none

    integer, intent(in) :: nm !! row dimension of two-dimensional array parameters
    integer, intent(in) :: n !! the order of the matrix
    integer, intent(in) :: low !! integer determined by the balancing
                               !! subroutine  cbal.  if  cbal  has not been used,
                               !! set low=1
    integer, intent(in) :: igh !! integer determined by the balancing
                               !! subroutine  cbal.  if  cbal  has not been used,
                               !! igh=n.
    real(wp), intent(inout) :: hr(nm, n) !! On input: hr and hi contain the real and imaginary parts,
                                         !! respectively, of the complex upper hessenberg matrix.
                                         !! their lower triangles below the subdiagonal contain
                                         !! information about the unitary transformations used in
                                         !! the reduction by  corth, if performed.
                                         !!
                                         !! On Output: the upper hessenberg portions of hr and hi have been
                                         !! destroyed.  therefore, they must be saved before
                                         !! calling  comqr  if subsequent calculation of
                                         !! eigenvectors is to be performed.
    real(wp), intent(inout) :: hi(nm, n) !! See `hr` description
    real(wp), intent(out) :: wr(n) !! the real parts of the eigenvalues.  if an error
                                   !! exit is made, the eigenvalues should be correct
                                   !! for indices `ierr+1,...,n`.
    real(wp), intent(out) :: wi(n) !! the imaginary parts of the eigenvalues.  if an error
                                   !! exit is made, the eigenvalues should be correct
                                   !! for indices `ierr+1,...,n`.
    integer, intent(out) :: ierr !! is set to:
                                 !!
                                 !!  * 0 -- for normal return
                                 !!  * j -- if the limit of 30*n iterations is exhausted
                                 !!    while the j-th eigenvalue is being sought.

    integer :: i, j, l, en, ll, itn, its, lp1, enm1
    real(wp) :: si, sr, ti, tr, xi, xr, yi, yr, zzi, &
                zzr, norm, tst1, tst2, xr2, xi2

    ierr = 0
    if (low /= igh) then
        ! create real subdiagonal elements
        l = low + 1
        do i = l, igh
            ll = min(i + 1, igh)
            if (hi(i, i - 1) /= 0.0_wp) then
                norm = pythag(hr(i, i - 1), hi(i, i - 1))
                yr = hr(i, i - 1)/norm
                yi = hi(i, i - 1)/norm
                hr(i, i - 1) = norm
                hi(i, i - 1) = 0.0_wp
                do j = i, igh
                    si = yr*hi(i, j) - yi*hr(i, j)
                    hr(i, j) = yr*hr(i, j) + yi*hi(i, j)
                    hi(i, j) = si
                end do
                do j = low, ll
                    si = yr*hi(j, i) + yi*hr(j, i)
                    hr(j, i) = yr*hr(j, i) - yi*hi(j, i)
                    hi(j, i) = si
                end do
            end if
        end do
    end if
    ! store roots isolated by cbal
    do i = 1, n
        if (i < low .or. i > igh) then
            wr(i) = hr(i, i)
            wi(i) = hi(i, i)
        end if
    end do

    en = igh
    tr = 0.0_wp
    ti = 0.0_wp
    itn = 30*n

    main: do

        if (en < low) return

        ! search for next eigenvalue
        its = 0
        enm1 = en - 1

        do

            ! look for single small sub-diagonal element
            ! for l=en step -1 until low d0 --
            do ll = low, en
                l = en + low - ll
                if (l == low) exit
                tst1 = abs(hr(l - 1, l - 1)) + abs(hi(l - 1, l - 1)) + abs(hr(l, l)) + abs(hi(l, l))
                tst2 = tst1 + abs(hr(l, l - 1))
                if (tst2 == tst1) exit
            end do

            ! form shift
            if (l == en) then
                ! a root found
                wr(en) = hr(en, en) + tr
                wi(en) = hi(en, en) + ti
                en = enm1
                cycle main
            elseif (itn == 0) then
                ! set error -- all eigenvalues have not converged after 30*n iterations
                ierr = en
                return
            else
                if (its == 10 .or. its == 20) then
                    ! form exceptional shift
                    sr = abs(hr(en, enm1)) + abs(hr(enm1, en - 2))
                    si = 0.0_wp
                else
                    sr = hr(en, en)
                    si = hi(en, en)
                    xr = hr(enm1, en)*hr(en, enm1)
                    xi = hi(enm1, en)*hr(en, enm1)
                    if (xr /= 0.0_wp .or. xi /= 0.0_wp) then
                        yr = (hr(enm1, enm1) - sr)/2.0_wp
                        yi = (hi(enm1, enm1) - si)/2.0_wp
                        call csroot(yr**2 - yi**2 + xr, 2.0_wp*yr*yi + xi, zzr, zzi)
                        if (yr*zzr + yi*zzi < 0.0_wp) then
                            zzr = -zzr
                            zzi = -zzi
                        end if
                        call cdiv(xr, xi, yr + zzr, yi + zzi, xr2, xi2)
                        xr = xr2
                        xi = xi2
                        sr = sr - xr
                        si = si - xi
                    end if
                end if

                do i = low, en
                    hr(i, i) = hr(i, i) - sr
                    hi(i, i) = hi(i, i) - si
                end do

                tr = tr + sr
                ti = ti + si
                its = its + 1
                itn = itn - 1
                ! reduce to triangle (rows)
                lp1 = l + 1

                do i = lp1, en
                    sr = hr(i, i - 1)
                    hr(i, i - 1) = 0.0_wp
                    norm = pythag(pythag(hr(i - 1, i - 1), hi(i - 1, i - 1)), sr)
                    xr = hr(i - 1, i - 1)/norm
                    wr(i - 1) = xr
                    xi = hi(i - 1, i - 1)/norm
                    wi(i - 1) = xi
                    hr(i - 1, i - 1) = norm
                    hi(i - 1, i - 1) = 0.0_wp
                    hi(i, i - 1) = sr/norm

                    do j = i, en
                        yr = hr(i - 1, j)
                        yi = hi(i - 1, j)
                        zzr = hr(i, j)
                        zzi = hi(i, j)
                        hr(i - 1, j) = xr*yr + xi*yi + hi(i, i - 1)*zzr
                        hi(i - 1, j) = xr*yi - xi*yr + hi(i, i - 1)*zzi
                        hr(i, j) = xr*zzr - xi*zzi - hi(i, i - 1)*yr
                        hi(i, j) = xr*zzi + xi*zzr - hi(i, i - 1)*yi
                    end do

                end do

                si = hi(en, en)
                if (si /= 0.0_wp) then
                    norm = pythag(hr(en, en), si)
                    sr = hr(en, en)/norm
                    si = si/norm
                    hr(en, en) = norm
                    hi(en, en) = 0.0_wp
                end if
                ! inverse operation (columns)
                do j = lp1, en
                    xr = wr(j - 1)
                    xi = wi(j - 1)

                    do i = l, j
                        yr = hr(i, j - 1)
                        yi = 0.0_wp
                        zzr = hr(i, j)
                        zzi = hi(i, j)
                        if (i /= j) then
                            yi = hi(i, j - 1)
                            hi(i, j - 1) = xr*yi + xi*yr + hi(j, j - 1)*zzi
                        end if
                        hr(i, j - 1) = xr*yr - xi*yi + hi(j, j - 1)*zzr
                        hr(i, j) = xr*zzr + xi*zzi - hi(j, j - 1)*yr
                        hi(i, j) = xr*zzi - xi*zzr - hi(j, j - 1)*yi
                    end do

                end do

                if (si /= 0.0_wp) then
                    do i = l, en
                        yr = hr(i, en)
                        yi = hi(i, en)
                        hr(i, en) = sr*yr - si*yi
                        hi(i, en) = sr*yi + si*yr
                    end do
                end if
            end if

        end do

    end do main

end subroutine comqr