Computes the eigenvalues and (if vec = True) the eigenvectors of an n * n upper Hessenberg matrix.
n : integer; ( n > 0 ) Dimension of h and eivec, length of the real parts vector wr and of the imaginary parts vector wi of the eigenvalues.
low : integer;
high : integer; see balance
h : n x n matrix; upper Hessenberg matrix as output of Elmhes (destroyed in the process).
eivec : n x n matrix; (only if vec = 1) Matrix, which for vec = 1 contains the eigenvectors as follows: For real eigebvalues the corresponding column contains the corresponding eigenvactor, while for complex eigenvalues the corresponding column contains the real part of the eigenvactor with its imaginary part is stored in the subsequent column of eivec. The eigenvactor for the complex conjugate eigenvactor is given by the complex conjugate eigenvactor.
wr : vector of size n; Real part of the n eigenvalues.
wi : vector of size n; Imaginary parts of the eigenvalues
cnt : Integer vector of size n; vector of iterations used for each eigenvalue. For a complex conjugate eigenvalue pair the second entry is negative.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n | |||
| integer, | intent(in) | :: | low | |||
| integer, | intent(in) | :: | high | |||
| real(kind=double), | intent(out) | :: | h(0:n-1,0:n-1) | |||
| real(kind=double), | intent(out) | :: | wr(0:n-1) | |||
| real(kind=double), | intent(out) | :: | wi(0:n-1) | |||
| real(kind=double), | intent(out) | :: | eivec(0:n-1,0:n-1) | |||
| integer, | intent(out) | :: | cnt(0:n-1) | |||
| integer, | intent(out) | :: | rc |
subroutine hqr2(n, & !Dimension of matrix .........
low, & !first nonzero row ...........
high, & !last nonzero row ............
h, & !Hessenberg matrix ...........
wr, & !Real parts of eigenvalues ...
wi, & !Imaginary parts of evalues ..
eivec,& !Matrix of eigenvectors ......
cnt, & !Iteration counter ...........
rc ) !return code .................
!+Computes the eigenvalues and (if vec = True) the eigenvectors
!+of an n * n upper Hessenberg matrix.
!+
!+###Input parameters
!+ n : integer; ( n > 0 )
!+ Dimension of h and eivec,
!+ length of the real parts vector wr and of the
!+ imaginary parts vector wi of the eigenvalues.
!+
!+ low : integer;
!+
!+ high : integer; see balance
!+
!+ h : n x n matrix;
!+ upper Hessenberg matrix as output of Elmhes
!+ (destroyed in the process).
!+###Output parameters
!+ eivec : n x n matrix; (only if vec = 1)
!+ Matrix, which for vec = 1 contains the
!+ eigenvectors as follows:
!+ For real eigebvalues the corresponding column
!+ contains the corresponding eigenvactor, while for
!+ complex eigenvalues the corresponding column contains
!+ the real part of the eigenvactor with its imaginary
!+ part is stored in the subsequent column of eivec.
!+ The eigenvactor for the complex conjugate eigenvactor
!+ is given by the complex conjugate eigenvactor.
!+
!+ wr : vector of size n;
!+ Real part of the n eigenvalues.
!+
!+ wi : vector of size n;
!+ Imaginary parts of the eigenvalues
!+
!+ cnt : Integer vector of size n;
!+ vector of iterations used for each eigenvalue.
!+ For a complex conjugate eigenvalue pair the second
!+ entry is negative.
integer,intent(in) :: n
integer,intent(in) :: high, low
real(double) ,intent(out) :: h(0:n-1,0:n-1)
real(double), intent(out) :: wr(0:n-1),wi(0:n-1)
real(double), intent(out) :: eivec(0:n-1,0:n-1)
integer,intent(out) :: rc
integer,intent(out) :: cnt(0:n-1)
integer :: en
integer :: i, j, na, iter, l, ll, m, k
real(double) :: p, q, r, s, t, w, x, y, z
real(double) :: r_p,r_r,r_s,r_x,r_z
p=ZERO; q=ZERO; r=ZERO
do i = 0, n-1
if (i < low.or.i > high) then
wr(i) = h(i,i)
wi(i) = ZERO
cnt(i) = 0
endif
enddo
en = high
t = ZERO
do while (en >= low)
iter = 0
na = en - 1
do while(1<2)
ll=999
do l = en, low+1, -1 !search for small
!subdiagonal element
if(DABS(h(l,l-1))<=XMACH_EPS*(DABS(h(l-1,l-1))+DABS(h(l,l))))then
ll=l; !save current index
goto 10 !exit l loop
endif
enddo
10 if(ll.ne.999)then
l=ll
else
l=0 !restore l
endif
x = h(en,en)
if (l == en) then !found one evalue
wr(en) = x + t
h(en,en) = x + t
wi(en) = ZERO
cnt(en) = iter
en = en - 1
goto 15 !exit from loop while(True)
endif
y = h(na,na)
w = h(en,na) * h(na,en)
if (l == na) then !found two evalues
p = (y - x) * 0.5d0
q = p * p + w
z = DSQRT(DABS(q))
x = x + t
h(en,en) = x + t
h(na,na) = y + t
cnt(en) = -iter
cnt(na) = iter
if (q >= ZERO) then !real eigenvalues
if (p<ZERO) then
z=p-z
else
z=p+z
endif
r_z = 1.0d0/z
wr(na) = x + z
wr(en) = x - w * r_z
s = w - w * r_z
wi(na) = ZERO
wi(en) = ZERO
x = h(en,na)
r = DSQRT (x * x + z * z)
r_r = 1.0d0/r
p = x * r_r
q = z * r_r
do j = na, n-1
z = h(na,j)
h(na,j) = q * z + p * h(en,j)
h(en,j) = q * h(en,j) - p * z
enddo
do i = 0, en
z = h(i,na)
h(i,na) = q * z + p * h(i,en)
h(i,en) = q * h(i,en) - p * z
enddo
do i = low, high
z = eivec(i,na)
eivec(i,na) = q * z + p * eivec(i,en)
eivec(i,en) = q * eivec(i,en) - p * z
enddo
else !pair of complex
wr(na) = x + p
wr(en) = x + p
wi(na) = z
wi(en) = - z
endif !if q>=ZERO
en = en - 2
goto 15 !exit while(1<2)
endif !if l = na
if (iter >= MAXIT) then
cnt(en) = MAXIT + 1
rc = en
write(*,*) ' stop at iter >= MAXIT.'
return
endif
if (iter.ne.0.and.MOD(iter,10) == 0) then
t = t + x
do i = low, en
h(i,i) = h(i,i) - x
enddo
s = DABS(h(en,na)) + DABS(h(na,en-2))
x = 0.75d0 * s; y = x
w = -0.4375d0 * s * s
endif
iter = iter + 1
do m = en - 2, l, -1
z = h(m,m)
r = x - z
s = y - z
p = ( r * s - w ) / h(m+1,m) + h(m,m+1)
q = h(m + 1,m + 1) - z - r - s
r = h(m + 2,m + 1)
s = DABS(p) + DABS(q) + DABS (r)
r_s = 1.0d0/s
p = p * r_s
q = q * r_s
r = r * r_s
if (m == l) goto 12
if (DABS(h(m,m-1)) * (DABS(q) + DABS(r)) <= XMACH_EPS * DABS(p) &
* (DABS(h(m-1,m-1)) + DABS(z) + DABS(h(m+1,m+1)))) then
goto 12 !exit m loop
endif
enddo
12 do i = m + 2, en
h(i,i-2) = ZERO
enddo
do i = m + 3, en
h(i,i-3) = ZERO
enddo
do k = m, na
if(k.ne.m)then!double QR step, for rows l to en and columns m to en
p = h(k,k-1)
q = h(k+1,k-1)
if (k.ne.na) then
r = h(k+2,k-1)
else
r = ZERO
endif
x = DABS(p) + DABS(q) + DABS(r)
if (x == ZERO) goto 30 !next k
r_x = 1.0d0/x
p = p * r_x
q = q * r_x
r = r * r_x
endif
s = DSQRT(p * p + q * q + r * r)
if (p < ZERO) s = -s
r_s = 1.0d0/s
if (k.ne.m) then
h(k,k-1) = -s * x
else if (l.ne.m) then
h(k,k-1) = -h(k,k-1)
endif
p = p + s
r_p = 1.0d0/p
x = p * r_s
y = q * r_s
z = r * r_s
q = q * r_p
r = r * r_p
do j = k, n-1 !modify rows
p = h(k,j) + q * h(k+1,j)
if (k.ne.na) then
p = p + r * h(k+2,j)
h(k+2,j) = h(k+2,j) - p * z
endif
h(k+1,j) = h(k+1,j) - p * y
h(k,j) = h(k,j) - p * x
enddo
if (k+3 < en) then
j=k+3
else
j=en
endif
do i = 0, j !modify columns
p = x * h(i,k) + y * h(i,k+1)
if (k.ne.na) then
p = p + z * h(i,k+2)
h(i,k+2) = h(i,k+2) - p * r
endif
h(i,k+1) = h(i,k+1) - p * q
h(i,k) = h(i,k) - p
enddo
do i = low, high
p = x * eivec(i,k) + y * eivec(i,k+1)
if (k.ne.na) then
p = p + z * eivec(i,k+2)
eivec(i,k+2) = eivec(i,k+2) - p * r
endif
eivec(i,k+1) = eivec(i,k+1) - p * q
eivec(i,k) = eivec(i,k) - p
enddo
30 continue
enddo !k loop
enddo !while(1<2)
15 continue
enddo !while en >= low All evalues found
!transform evectors back
call hqrvec (n, low, high, h, wr, wi, eivec,rc)
end subroutine hqr2