1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304
|
C/MEMBR ADD NAME=QVALZ,SSI=0
subroutine qvalz(nm,n,a,b,epsb,alfr,alfi,beta,matq,q,matz,z)
c
integer i,j,n,en,na,nm,nn,isw
double precision a(nm,n),b(nm,n),alfr(n),alfi(n),beta(n)
double precision z(nm,n),q(nm,n)
double precision c,d,e,r,s,t,an,a1,a2,bn,cq,cz,di,dr,ei,ti,tr
double precision u1,u2,v1,v2,a1i,a11,a12,a2i,a21,a22
double precision b11,b12,b22,sqi,sqr
double precision ssi,ssr,szi,szr,a11i,a11r,a12i,a12r,a22i,a22r
double precision epsb
logical matz,matq
c
c! purpose
c this subroutine accepts a pair of real matrices, one of them
c in quasi-triangular form and the other in upper triangular form.
c it reduces the quasi-triangular matrix further, so that any
c remaining 2-by-2 blocks correspond to pairs of complex
c eigenvalues, and returns quantities whose ratios give the
c generalized eigenvalues. it is usually preceded by qzhes
c and qzit and may be followed by qzvec.
c
c MODIFIED FROM EISPACK ROUTINE ``QZVAL'' TO ALSO RETURN THE Q
c MATRIX. IN ADDITION, THE TOLERANCE epsb IS DIRECTLY PASSED IN
c THE CALLING LIST INSTEAD OF VIA b(n,1)
c
c! calling sequence
c
c subroutine qvalz(nm,n,a,b,epsb,alfr,alfi,beta,matq,q,matz,z)
c on input:
c
c nm must be set to the row dimension of two-dimensional
c array parameters as declared in the calling program
c dimension statement;
c
c n is the order of the matrices;
c
c a contains a real upper quasi-triangular matrix;
c
c b contains a real upper triangular matrix.
c
c epsb: tolerance computed and saved in qitz (qzit)
c
c matz (resp matq) should be set to .true. if the right
c (resp left) hand transformations are to be accumulated
c for later use in computing eigenvectors, and to .false.
c otherwise;
c
c z (resp q) contains, if matz (resp matq) has been set
c to .true., the transformation matrix produced in the
c reductions by qzhes and qzit, if performed, or else the
c identity matrix. if matz has been set to .false., z is not
c referenced.
c
c on output:
c
c a has been reduced further to a quasi-triangular matrix
c in which all nonzero subdiagonal elements correspond to
c pairs of complex eigenvalues;
c
c b is still in upper triangular form, although its elements
c have been altered. b(n,1) is unaltered;
c
c alfr and alfi contain the real and imaginary parts of the
c diagonal elements of the triangular matrix that would be
c obtained if a were reduced completely to triangular form
c by unitary transformations. non-zero values of alfi occur
c in pairs, the first member positive and the second negative;
c
c beta contains the diagonal elements of the corresponding b,
c normalized to be real and non-negative. the generalized
c eigenvalues are then the ratios ((alfr+i*alfi)/beta);
c
c z (resp q) contains the product of the right resp left hand
c (for all three steps) if matz (resp, matq) has been set
c to .true.
c
c! originator
c
c this subroutine is the third step of the qz algorithm
c for solving generalized matrix eigenvalue problems,
c siam j. numer. anal. 10, 241-256(1973) by moler and stewart.
c modification de la routine qzval de eispack pour avoir la matrice
c q en option
c!
c questions and comments should be directed to b. s. garbow,
c applied mathematics division, argonne national laboratory
c
c ------------------------------------------------------------------
c
isw = 1
c :::::::::: find eigenvalues of quasi-triangular matrices.
c for en=n step -1 until 1 do -- ::::::::::
do 510 nn = 1, n
en = n + 1 - nn
na = en - 1
if (isw .eq. 2) go to 505
if (en .eq. 1) go to 410
if (a(en,na) .ne. 0.0d+0) go to 420
c :::::::::: 1-by-1 block, one real root ::::::::::
410 alfr(en) = a(en,en)
if (b(en,en) .lt. 0.0d+0) alfr(en) = -alfr(en)
beta(en) = abs(b(en,en))
alfi(en) = 0.0d+0
go to 510
c :::::::::: 2-by-2 block ::::::::::
420 if (abs(b(na,na)) .le. epsb) go to 455
if (abs(b(en,en)) .gt. epsb) go to 430
a1 = a(en,en)
a2 = a(en,na)
bn = 0.0d+0
go to 435
430 an = abs(a(na,na)) + abs(a(na,en)) + abs(a(en,na))
& + abs(a(en,en))
bn = abs(b(na,na)) + abs(b(na,en)) + abs(b(en,en))
a11 = a(na,na) / an
a12 = a(na,en) / an
a21 = a(en,na) / an
a22 = a(en,en) / an
b11 = b(na,na) / bn
b12 = b(na,en) / bn
b22 = b(en,en) / bn
e = a11 / b11
ei = a22 / b22
s = a21 / (b11 * b22)
t = (a22 - e * b22) / b22
if (abs(e) .le. abs(ei)) go to 431
e = ei
t = (a11 - e * b11) / b11
431 c = 0.50d+0 * (t - s * b12)
d = c * c + s * (a12 - e * b12)
if (d .lt. 0.0d+0) go to 480
c :::::::::: two real roots.
c zero both a(en,na) and b(en,na) ::::::::::
e = e + (c + sign(sqrt(d),c))
a11 = a11 - e * b11
a12 = a12 - e * b12
a22 = a22 - e * b22
if (abs(a11) + abs(a12) .lt.
x abs(a21) + abs(a22)) go to 432
a1 = a12
a2 = a11
go to 435
432 a1 = a22
a2 = a21
c :::::::::: choose and apply real z ::::::::::
435 s = abs(a1) + abs(a2)
u1 = a1 / s
u2 = a2 / s
r = sign(sqrt(u1*u1+u2*u2),u1)
v1 = -(u1 + r) / r
v2 = -u2 / r
u2 = v2 / v1
c
do 440 i = 1, en
t = a(i,en) + u2 * a(i,na)
a(i,en) = a(i,en) + t * v1
a(i,na) = a(i,na) + t * v2
t = b(i,en) + u2 * b(i,na)
b(i,en) = b(i,en) + t * v1
b(i,na) = b(i,na) + t * v2
440 continue
c
if (.not. matz) go to 450
c
do 445 i = 1, n
t = z(i,en) + u2 * z(i,na)
z(i,en) = z(i,en) + t * v1
z(i,na) = z(i,na) + t * v2
445 continue
c
450 if (bn .eq. 0.0d+0) go to 475
if (an .lt. abs(e) * bn) go to 455
a1 = b(na,na)
a2 = b(en,na)
go to 460
455 a1 = a(na,na)
a2 = a(en,na)
c :::::::::: choose and apply real q ::::::::::
460 s = abs(a1) + abs(a2)
if (s .eq. 0.0d+0) go to 475
u1 = a1 / s
u2 = a2 / s
r = sign(sqrt(u1*u1+u2*u2),u1)
v1 = -(u1 + r) / r
v2 = -u2 / r
u2 = v2 / v1
c
do 470 j = na, n
t = a(na,j) + u2 * a(en,j)
a(na,j) = a(na,j) + t * v1
a(en,j) = a(en,j) + t * v2
t = b(na,j) + u2 * b(en,j)
b(na,j) = b(na,j) + t * v1
b(en,j) = b(en,j) + t * v2
470 continue
ccccccccccccccccccccccccccccccccccccccccc
c MODIFIED TO ACCUMULATE Q AS WELL
ccccccccccccccccccccccccccccccccccccccc
if(.not.matq) goto 471
do 472 j=1,n
t=q(na,j)+u2*q(en,j)
q(na,j)=q(na,j)+t*v1
q(en,j)=q(en,j)+t*v2
472 continue
cccccccccccccccccccccccccccccccccccccccc
471 continue
c
475 a(en,na) = 0.0d+0
b(en,na) = 0.0d+0
alfr(na) = a(na,na)
alfr(en) = a(en,en)
if (b(na,na) .lt. 0.0d+0) alfr(na) = -alfr(na)
if (b(en,en) .lt. 0.0d+0) alfr(en) = -alfr(en)
beta(na) = abs(b(na,na))
beta(en) = abs(b(en,en))
alfi(en) = 0.0d+0
alfi(na) = 0.0d+0
go to 505
c :::::::::: two complex roots ::::::::::
480 e = e + c
ei = sqrt(-d)
a11r = a11 - e * b11
a11i = ei * b11
a12r = a12 - e * b12
a12i = ei * b12
a22r = a22 - e * b22
a22i = ei * b22
if (abs(a11r) + abs(a11i) + abs(a12r) + abs(a12i) .lt.
x abs(a21) + abs(a22r) + abs(a22i)) go to 482
a1 = a12r
a1i = a12i
a2 = -a11r
a2i = -a11i
go to 485
482 a1 = a22r
a1i = a22i
a2 = -a21
a2i = 0.0d+0
c :::::::::: choose complex z ::::::::::
485 cz = sqrt(a1*a1+a1i*a1i)
if (cz .eq. 0.0d+0) go to 487
szr = (a1 * a2 + a1i * a2i) / cz
szi = (a1 * a2i - a1i * a2) / cz
r = sqrt(cz*cz+szr*szr+szi*szi)
cz = cz / r
szr = szr / r
szi = szi / r
go to 490
487 szr = 1.0d+0
szi = 0.0d+0
490 if (an .lt. (abs(e) + ei) * bn) go to 492
a1 = cz * b11 + szr * b12
a1i = szi * b12
a2 = szr * b22
a2i = szi * b22
go to 495
492 a1 = cz * a11 + szr * a12
a1i = szi * a12
a2 = cz * a21 + szr * a22
a2i = szi * a22
c :::::::::: choose complex q ::::::::::
495 cq = sqrt(a1*a1+a1i*a1i)
if (cq .eq. 0.0d+0) go to 497
sqr = (a1 * a2 + a1i * a2i) / cq
sqi = (a1 * a2i - a1i * a2) / cq
r = sqrt(cq*cq+sqr*sqr+sqi*sqi)
cq = cq / r
sqr = sqr / r
sqi = sqi / r
go to 500
497 sqr = 1.0d+0
sqi = 0.0d+0
c :::::::::: compute diagonal elements that would result
c if transformations were applied ::::::::::
500 ssr = sqr * szr + sqi * szi
ssi = sqr * szi - sqi * szr
i = 1
tr = cq * cz * a11 + cq * szr * a12 + sqr * cz * a21
x + ssr * a22
ti = cq * szi * a12 - sqi * cz * a21 + ssi * a22
dr = cq * cz * b11 + cq * szr * b12 + ssr * b22
di = cq * szi * b12 + ssi * b22
go to 503
502 i = 2
tr = ssr * a11 - sqr * cz * a12 - cq * szr * a21
x + cq * cz * a22
ti = -ssi * a11 - sqi * cz * a12 + cq * szi * a21
dr = ssr * b11 - sqr * cz * b12 + cq * cz * b22
di = -ssi * b11 - sqi * cz * b12
503 t = ti * dr - tr * di
j = na
if (t .lt. 0.0d+0) j = en
r = sqrt(dr*dr+di*di)
beta(j) = bn * r
alfr(j) = an * (tr * dr + ti * di) / r
alfi(j) = an * t / r
if (i .eq. 1) go to 502
505 isw = 3 - isw
510 continue
c
return
c :::::::::: last card of qzval ::::::::::
end
|