hqrvec Subroutine

public subroutine hqrvec(n, low, high, h, wr, wi, eivec, rc)

Computes the eigenvectors for the eigenvalues found in hqr2

Input parameters

n : int n; ( n > 0 ) : Dimension of mat and eivec, number of eigenvalues.

low : int low;

high : int high; see balance

h : n x n upper Hessenberg matrix

wr : vector of size n; : Real parts of the n eigenvalues.

wi : vector of size n; Imaginary parts of the n eigenvalues.

Output parameter:

eivec : n x n matrix, whose columns are the eigenvectors w x | | h[i][en] | | -r | | | | = | | y z | | h[i+1][en] | | -s | w+iq x | | h[i][na] + ih[i][en] | | -ra+isa | | | | = | | y z+iq| | h[i+1][na]+ih[i+1][en]| | -r+is |

Arguments

Type IntentOptional AttributesName
integer, intent(in) :: n
integer, intent(in) :: low
integer, intent(in) :: high
real(kind=double) :: h(0:n-1,0:n-1)
real(kind=double), intent(in) :: wr(0:n-1)
real(kind=double), intent(in) :: wi(0:n-1)
real(kind=double), intent(out) :: eivec(0:n-1,0:n-1)
integer :: rc

Called by

proc~~hqrvec~~CalledByGraph proc~hqrvec hqrvec proc~hqr2 hqr2 proc~hqr2->proc~hqrvec proc~eigen eigen proc~eigen->proc~hqr2 proc~colrad colrad proc~colrad->proc~eigen proc~fida_weights_los fida_weights_los proc~fida_weights_los->proc~colrad proc~attenuate attenuate proc~attenuate->proc~colrad proc~npa_weights npa_weights proc~npa_weights->proc~attenuate proc~npa_f npa_f proc~npa_f->proc~attenuate program~fidasim fidasim program~fidasim->proc~fida_weights_los program~fidasim->proc~npa_weights program~fidasim->proc~npa_f proc~npa_mc npa_mc program~fidasim->proc~npa_mc proc~npa_mc->proc~attenuate

Contents

Source Code


Source Code

  subroutine  hqrvec(n,     & !Dimension of matrix .......
       low,   & !first nonzero row .........
       high,  & !last nonzero row ..........
       h,     & !upper Hessenberg matrix ...
       wr,    & !Real parts of evalues .....
       wi,    & !Imaginary parts of evalues
       eivec, & !Eigenvectors ..............
       rc  )   !return code ...............
    !+Computes the eigenvectors for the eigenvalues found in hqr2
    !+
    !+###Input parameters
    !+   n     :   int n;  ( n > 0 )
    !+         :   Dimension of  mat and eivec, number of eigenvalues.
    !+
    !+   low   :   int low;
    !+
    !+   high  :   int high; see  balance
    !+
    !+   h     :   n x n upper Hessenberg matrix
    !+
    !+   wr    :   vector of size n;
    !+         :   Real parts of the n eigenvalues.
    !+
    !+   wi    :   vector of size n; Imaginary parts of the n eigenvalues.
    !+
    !+###Output parameter:
    !+   eivec :  n x n matrix, whose columns are the eigenvectors
    integer,intent(in)         :: n
    integer,intent(in)         :: high, low
    real(double), intent(in)   :: wr(0:n-1),wi(0:n-1)
    real(double), intent(out)  :: eivec(0:n-1,0:n-1)
    real(double)  :: h(0:n-1,0:n-1)
    integer :: rc
    integer :: i, j, m, k, na, l
    integer :: code, en
    real(double)  :: p, q, r, s, t, w, x, y, z, ra, sa, vr, vi, norm, temp

    r=ZERO; s=ZERO; z=ZERO; norm=ZERO
    do i = 0, n-1                               !find norm of h
       do j = i, n-1
          norm = norm + DABS(h(i,j))
       enddo
    enddo
    if (norm == ZERO) then
       rc = 1                                    !zero matrix
       return
    endif
    do en = n-1, 0, -1                          !transform back
       p = wr(en)
       q = wi(en)
       na = en - 1
       if (q == ZERO) then
          m = en
          h(en,en) = ONE
          do i = na, 0, -1
             w = h(i,i) - p
             r = h(i,en)
             do j = m, na
                r = r + h(i,j) * h(j,en)
             enddo
             if (wi(i) < ZERO) then
                z = w
                s = r
             else
                m = i
                if (wi(i) == ZERO) then
                   if (w.ne.ZERO) then
                      temp = w
                   else
                      temp=XMACH_EPS * norm
                   endif
                   h(i,en) = -r/temp
                else
                   !Solve the linear system:
                   !| w   x |  | h[i][en]   |   | -r |
                   !|       |  |            | = |    |
                   !| y   z |  | h[i+1][en] |   | -s |
                   x = h(i,i+1)
                   y = h(i+1,i)
                   q = (wr(i) - p)**2 + wi(i)**2
                   h(i,en) = (x * s - z * r) / q
                   t = h(i,en)
                   if (DABS(x) > DABS(z)) then
                      temp = (-r -w * t) / x
                   else
                      temp = (-s -y * t) / z
                   endif
                   h(i+1,en) = temp
                endif
             endif !wi[i] < 0
          enddo !i loop
       else if (q < ZERO) then
          m = na
          if (DABS(h(en,na)) > DABS(h(na,en))) then
             h(na,na) = - (h(en,en) - p) / h(en,na)
             h(na,en) = - q / h(en,na)
          else
             call Comdiv(-h(na,en),0.d0, h(na,na)-p, q, h(na,na), h(na,en),code)
          endif
          h(en,na) = ONE
          h(en,en) = ZERO
          do i = na - 1, 0, -1
             w = h(i,i) - p
             ra = h(i,en)
             sa = ZERO
             do j = m, na
                ra = ra + h(i,j) * h(j,na)
                sa = sa + h(i,j) * h(j,en)
             enddo
             if (wi(i) < ZERO) then
                z = w
                r = ra
                s = sa
             else
                m = i
                if (wi(i) == ZERO) then
                   call Comdiv(-ra, -sa, w, q, h(i,na), h(i,en),code)
                else
            !  solve complex linear system:
                   !| w+i*q     x | | h[i][na] + i*h[i][en]  |   | -ra+i*sa |
                   !|             | |                        | = |          |
            !|   y    z+i*q| | h[i+1][na]+i*h[i+1][en]|   | -r+i*s   |
                   x = h(i,i+1)
                   y = h(i+1,i)
                   vr = (wr(i) - p)**2 + wi(i)**2 - q*q
                   vi = TWO * q * (wr(i) - p)
                   if (vr == ZERO.AND.vi == ZERO) then
                      vr = XMACH_EPS * norm * (DABS(w) + DABS(q)  &
                           + DABS(x) + DABS(y) + DABS(z))
                   endif

                   call Comdiv (x*r-z*ra+q*sa,x*s-z*sa-q*ra &
                        ,vr,vi,h(i,na),h(i,en),code)
                   if (DABS(x) > DABS(z) + DABS(q)) then
                      h(i+1,na) = (-ra - w * h(i,na) + q * h(i,en)) / x
                      h(i+1,en) = (-sa - w * h(i,en) - q * h(i,na)) / x
                   else
                      call Comdiv (-r - y * h(i,na), -s - y * h(i,en) &
                           , z, q, h(i+1,na), h(i+1,en),code)
                   endif
                endif !wi[i] = 0
             endif !wi[i] < 0
          enddo !i loop
       endif !else if q < 0
    enddo !en loop
    do i = 0, n-1                        !Eigenvectors for the evalues for
       if (i < low.or.i > high) then      !rows < low  and rows  > high
          do k = i + 1, n-1
             eivec(i,k) = h(i,k)
          enddo
       endif
    enddo
    j = n-1
    do while (j>=low)
       if(j<=high)then
          m =j
       else
          j = high
       endif
       if(j < 0) exit
       if (wi(j) < ZERO) then
          l=j-1
          do i = low, high
             y=ZERO; z=ZERO
             do k = low, m
                y = y + eivec(i,k) * h(k,l)
                z = z + eivec(i,k) * h(k,j)
             enddo
             eivec(i,l) = y
             eivec(i,j) = z
          enddo
       else
          if (wi(j) == ZERO) then
             do i = low, high
                z = ZERO
                do k = low, m
                   z = z + eivec(i,k) * h(k,j)
                enddo
                eivec(i,j) = z
             enddo
          endif
       endif
       j = j - 1
    enddo !j loop
    rc = 0
  end subroutine hqrvec