Computes the eigenvectors for the eigenvalues found in hqr2
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.
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 |
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| 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 |
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