scomqr – polyroots-fortran

private subroutine scomqr(nm, n, low, igh, hr, hi, z, ierr)

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

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.

Reference

Original code from JPL MATH77 Library

License

History

1987-11-12 SCOMQR Lawson Initial code.
1992-03-13 SCOMQR FTK Removed implicit statements.
1995-01-03 SCOMQR WVS Added EXTERNAL CQUO, CSQRT so VAX won't gripe
1996-01-18 SCOMQR Krogh Added M77CON statements for conv. to C.
1996-03-30 SCOMQR Krogh Added external statement.
1996-04-27 SCOMQR Krogh Changes to use .C. and C%%.
2001-01-24 SCOMQR Krogh CSQRT -> CSQRTX to avoid C lib. conflicts.
2022-10-06, Jacob Williams modernized this routine

Arguments

Type	Intent		Name
integer,	intent(in)	::	nm	the row dimension of two-dimensional array parameters as declared in the calling program dimension statement
integer,	intent(in)	::	n	the order of the matrix
integer,	intent(in)	::	low	low and igh are integers determined by the balancing subroutine cbal. if cbal has not been used, set low=1, igh=n
integer,	intent(in)	::	igh	low and igh are integers determined by the balancing subroutine cbal. if cbal has not been used, set low=1, igh=n
real(kind=wp),	intent(inout)	::	hr(nm,n)	see `hi` description
real(kind=wp),	intent(inout)	::	hi(nm,n)	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. 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,
complex(kind=wp),	intent(out)	::	z(n)	the real and imaginary parts, respectively, 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: zero -- for normal return, j -- if the j-th eigenvalue has not been determined after 30 iterations.

Called by

Help

Source Code

    subroutine scomqr(nm,n,low,igh,hr,hi,z,ierr)

    integer,intent(in) :: nm !! the row dimension of two-dimensional array
                             !! parameters as declared in the calling program
                             !! dimension statement
    integer,intent(in) :: n !! the order of the matrix
    integer,intent(in) :: low !! low and igh are integers determined by the balancing
                              !! subroutine  cbal.  if  cbal  has not been used,
                              !! set low=1, igh=n
    integer,intent(in) :: igh !! low and igh are integers determined by the balancing
                              !! subroutine  cbal.  if  cbal  has not been used,
                              !! set low=1, igh=n
    real(wp),intent(inout) :: hi(nm,n) !! * 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.
                                       !!
                                       !! * 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) :: hr(nm,n) !! see `hi` description
    complex(wp),intent(out) :: z(n) !! the real and imaginary parts,
                                    !! respectively, 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:
                                !!
                                !!  * zero -- for normal return,
                                !!  * j -- if the j-th eigenvalue has not been
                                !!    determined after 30 iterations.

    integer :: en,enm1,i,its,j,l,ll,lp1
    real(wp) :: norm,si,sr,ti,tr,xi,xr,yi,yr,zzi,zzr
    complex(wp) :: z3

    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) cycle
            norm = abs(cmplx(hr(i,i-1),hi(i,i-1),wp))
            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 do
    end if

      ! store roots isolated by cbal
      do i = 1, n
         if (i >= low .and. i <= igh) cycle
         z(i) = cmplx(hr(i,i),hi(i,i),wp)
      end do

      en = igh
      tr = 0.0_wp
      ti = 0.0_wp

    main : do
          ! search for next eigenvalue
          if (en < low) return
          its = 0
          enm1 = en - 1

          do
              ! look for single small sub-diagonal element
              ! for l=en step -1 until low
              do ll = low, en
                 l = en + low - ll
                 if (l == low) exit
                 if (abs(hr(l,l-1)) <= &
                          eps * (abs(hr(l-1,l-1)) + abs(hi(l-1,l-1)) &
                          + abs(hr(l,l)) +abs(hi(l,l)))) exit
              end do
              ! form shift
              if (l == en) then
              ! a root found
                z(en) = cmplx(hr(en,en)+tr,hi(en,en)+ti,wp)
                en = enm1
                cycle main
              end if
              if (its == 30) exit main
              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
                    z3 = sqrt(cmplx(yr**2-yi**2+xr,2.0_wp*yr*yi+xi,wp))
                    zzr = real(z3,wp)
                    zzi = aimag(z3)
                    if (yr * zzr + yi * zzi < 0.0_wp) then
                        zzr = -zzr
                        zzi = -zzi
                    end if
                    z3 = cmplx(xr,xi,wp) / cmplx(yr+zzr,yi+zzi,wp)
                    sr = sr - real(z3,wp)
                    si = si - aimag(z3)
                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
              ! reduce to triangle (rows)
              lp1 = l + 1

              do i = lp1, en
                 sr = hr(i,i-1)
                 hr(i,i-1) = 0.0_wp
                 norm = sqrt(hr(i-1,i-1)*hr(i-1,i-1)+hi(i-1,i-1)*hi(i-1,i-1)+sr*sr)
                 xr = hr(i-1,i-1) / norm
                 xi = hi(i-1,i-1) / norm
                 z(i-1) = cmplx(xr,xi,wp)
                 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 = abs(cmplx(hr(en,en),si,wp))
                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 = real(z(j-1),wp)
                 xi = aimag(z(j-1))
                 do i = l, j
                    yr = hr(i,j-1)
                    yi = 0.0
                    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 do

    end do main

    ! set error -- no convergence to an
    ! eigenvalue after 30 iterations
    ierr = en

    end subroutine scomqr

scomqr Subroutine

Contents

Source Code

private subroutine scomqr(nm, n, low, igh, hr, hi, z, ierr)

Reference

License

History

Arguments

Called by

Source Code