|
2329
|
1 SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, |
|
|
2 $ LDVR, MM, M, WORK, INFO ) |
|
|
3 * |
|
3333
|
4 * -- LAPACK routine (version 3.0) -- |
|
2329
|
5 * Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd., |
|
|
6 * Courant Institute, Argonne National Lab, and Rice University |
|
3333
|
7 * June 30, 1999 |
|
2329
|
8 * |
|
|
9 * .. Scalar Arguments .. |
|
|
10 CHARACTER HOWMNY, SIDE |
|
|
11 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N |
|
|
12 * .. |
|
|
13 * .. Array Arguments .. |
|
|
14 LOGICAL SELECT( * ) |
|
|
15 DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), |
|
|
16 $ WORK( * ) |
|
|
17 * .. |
|
|
18 * |
|
|
19 * Purpose |
|
|
20 * ======= |
|
|
21 * |
|
|
22 * DTREVC computes some or all of the right and/or left eigenvectors of |
|
|
23 * a real upper quasi-triangular matrix T. |
|
|
24 * |
|
|
25 * The right eigenvector x and the left eigenvector y of T corresponding |
|
|
26 * to an eigenvalue w are defined by: |
|
|
27 * |
|
|
28 * T*x = w*x, y'*T = w*y' |
|
|
29 * |
|
|
30 * where y' denotes the conjugate transpose of the vector y. |
|
|
31 * |
|
|
32 * If all eigenvectors are requested, the routine may either return the |
|
|
33 * matrices X and/or Y of right or left eigenvectors of T, or the |
|
|
34 * products Q*X and/or Q*Y, where Q is an input orthogonal |
|
|
35 * matrix. If T was obtained from the real-Schur factorization of an |
|
|
36 * original matrix A = Q*T*Q', then Q*X and Q*Y are the matrices of |
|
|
37 * right or left eigenvectors of A. |
|
|
38 * |
|
|
39 * T must be in Schur canonical form (as returned by DHSEQR), that is, |
|
|
40 * block upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each |
|
|
41 * 2-by-2 diagonal block has its diagonal elements equal and its |
|
|
42 * off-diagonal elements of opposite sign. Corresponding to each 2-by-2 |
|
|
43 * diagonal block is a complex conjugate pair of eigenvalues and |
|
|
44 * eigenvectors; only one eigenvector of the pair is computed, namely |
|
|
45 * the one corresponding to the eigenvalue with positive imaginary part. |
|
|
46 * |
|
|
47 * Arguments |
|
|
48 * ========= |
|
|
49 * |
|
|
50 * SIDE (input) CHARACTER*1 |
|
|
51 * = 'R': compute right eigenvectors only; |
|
|
52 * = 'L': compute left eigenvectors only; |
|
|
53 * = 'B': compute both right and left eigenvectors. |
|
|
54 * |
|
|
55 * HOWMNY (input) CHARACTER*1 |
|
|
56 * = 'A': compute all right and/or left eigenvectors; |
|
|
57 * = 'B': compute all right and/or left eigenvectors, |
|
|
58 * and backtransform them using the input matrices |
|
|
59 * supplied in VR and/or VL; |
|
|
60 * = 'S': compute selected right and/or left eigenvectors, |
|
|
61 * specified by the logical array SELECT. |
|
|
62 * |
|
|
63 * SELECT (input/output) LOGICAL array, dimension (N) |
|
|
64 * If HOWMNY = 'S', SELECT specifies the eigenvectors to be |
|
|
65 * computed. |
|
|
66 * If HOWMNY = 'A' or 'B', SELECT is not referenced. |
|
|
67 * To select the real eigenvector corresponding to a real |
|
|
68 * eigenvalue w(j), SELECT(j) must be set to .TRUE.. To select |
|
|
69 * the complex eigenvector corresponding to a complex conjugate |
|
|
70 * pair w(j) and w(j+1), either SELECT(j) or SELECT(j+1) must be |
|
|
71 * set to .TRUE.; then on exit SELECT(j) is .TRUE. and |
|
|
72 * SELECT(j+1) is .FALSE.. |
|
|
73 * |
|
|
74 * N (input) INTEGER |
|
|
75 * The order of the matrix T. N >= 0. |
|
|
76 * |
|
|
77 * T (input) DOUBLE PRECISION array, dimension (LDT,N) |
|
|
78 * The upper quasi-triangular matrix T in Schur canonical form. |
|
|
79 * |
|
|
80 * LDT (input) INTEGER |
|
|
81 * The leading dimension of the array T. LDT >= max(1,N). |
|
|
82 * |
|
|
83 * VL (input/output) DOUBLE PRECISION array, dimension (LDVL,MM) |
|
|
84 * On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must |
|
|
85 * contain an N-by-N matrix Q (usually the orthogonal matrix Q |
|
|
86 * of Schur vectors returned by DHSEQR). |
|
|
87 * On exit, if SIDE = 'L' or 'B', VL contains: |
|
|
88 * if HOWMNY = 'A', the matrix Y of left eigenvectors of T; |
|
3333
|
89 * VL has the same quasi-lower triangular form |
|
|
90 * as T'. If T(i,i) is a real eigenvalue, then |
|
|
91 * the i-th column VL(i) of VL is its |
|
|
92 * corresponding eigenvector. If T(i:i+1,i:i+1) |
|
|
93 * is a 2-by-2 block whose eigenvalues are |
|
|
94 * complex-conjugate eigenvalues of T, then |
|
|
95 * VL(i)+sqrt(-1)*VL(i+1) is the complex |
|
|
96 * eigenvector corresponding to the eigenvalue |
|
|
97 * with positive real part. |
|
2329
|
98 * if HOWMNY = 'B', the matrix Q*Y; |
|
|
99 * if HOWMNY = 'S', the left eigenvectors of T specified by |
|
|
100 * SELECT, stored consecutively in the columns |
|
|
101 * of VL, in the same order as their |
|
|
102 * eigenvalues. |
|
|
103 * A complex eigenvector corresponding to a complex eigenvalue |
|
|
104 * is stored in two consecutive columns, the first holding the |
|
|
105 * real part, and the second the imaginary part. |
|
|
106 * If SIDE = 'R', VL is not referenced. |
|
|
107 * |
|
|
108 * LDVL (input) INTEGER |
|
|
109 * The leading dimension of the array VL. LDVL >= max(1,N) if |
|
|
110 * SIDE = 'L' or 'B'; LDVL >= 1 otherwise. |
|
|
111 * |
|
|
112 * VR (input/output) DOUBLE PRECISION array, dimension (LDVR,MM) |
|
|
113 * On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must |
|
|
114 * contain an N-by-N matrix Q (usually the orthogonal matrix Q |
|
|
115 * of Schur vectors returned by DHSEQR). |
|
|
116 * On exit, if SIDE = 'R' or 'B', VR contains: |
|
|
117 * if HOWMNY = 'A', the matrix X of right eigenvectors of T; |
|
3333
|
118 * VR has the same quasi-upper triangular form |
|
|
119 * as T. If T(i,i) is a real eigenvalue, then |
|
|
120 * the i-th column VR(i) of VR is its |
|
|
121 * corresponding eigenvector. If T(i:i+1,i:i+1) |
|
|
122 * is a 2-by-2 block whose eigenvalues are |
|
|
123 * complex-conjugate eigenvalues of T, then |
|
|
124 * VR(i)+sqrt(-1)*VR(i+1) is the complex |
|
|
125 * eigenvector corresponding to the eigenvalue |
|
|
126 * with positive real part. |
|
2329
|
127 * if HOWMNY = 'B', the matrix Q*X; |
|
|
128 * if HOWMNY = 'S', the right eigenvectors of T specified by |
|
|
129 * SELECT, stored consecutively in the columns |
|
|
130 * of VR, in the same order as their |
|
|
131 * eigenvalues. |
|
|
132 * A complex eigenvector corresponding to a complex eigenvalue |
|
|
133 * is stored in two consecutive columns, the first holding the |
|
|
134 * real part and the second the imaginary part. |
|
|
135 * If SIDE = 'L', VR is not referenced. |
|
|
136 * |
|
|
137 * LDVR (input) INTEGER |
|
|
138 * The leading dimension of the array VR. LDVR >= max(1,N) if |
|
|
139 * SIDE = 'R' or 'B'; LDVR >= 1 otherwise. |
|
|
140 * |
|
|
141 * MM (input) INTEGER |
|
|
142 * The number of columns in the arrays VL and/or VR. MM >= M. |
|
|
143 * |
|
|
144 * M (output) INTEGER |
|
|
145 * The number of columns in the arrays VL and/or VR actually |
|
|
146 * used to store the eigenvectors. |
|
|
147 * If HOWMNY = 'A' or 'B', M is set to N. |
|
|
148 * Each selected real eigenvector occupies one column and each |
|
|
149 * selected complex eigenvector occupies two columns. |
|
|
150 * |
|
|
151 * WORK (workspace) DOUBLE PRECISION array, dimension (3*N) |
|
|
152 * |
|
|
153 * INFO (output) INTEGER |
|
|
154 * = 0: successful exit |
|
|
155 * < 0: if INFO = -i, the i-th argument had an illegal value |
|
|
156 * |
|
|
157 * Further Details |
|
|
158 * =============== |
|
|
159 * |
|
|
160 * The algorithm used in this program is basically backward (forward) |
|
|
161 * substitution, with scaling to make the the code robust against |
|
|
162 * possible overflow. |
|
|
163 * |
|
|
164 * Each eigenvector is normalized so that the element of largest |
|
|
165 * magnitude has magnitude 1; here the magnitude of a complex number |
|
|
166 * (x,y) is taken to be |x| + |y|. |
|
|
167 * |
|
|
168 * ===================================================================== |
|
|
169 * |
|
|
170 * .. Parameters .. |
|
|
171 DOUBLE PRECISION ZERO, ONE |
|
|
172 PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 ) |
|
|
173 * .. |
|
|
174 * .. Local Scalars .. |
|
|
175 LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV |
|
|
176 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2 |
|
|
177 DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE, |
|
|
178 $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR, |
|
|
179 $ XNORM |
|
|
180 * .. |
|
|
181 * .. External Functions .. |
|
|
182 LOGICAL LSAME |
|
|
183 INTEGER IDAMAX |
|
|
184 DOUBLE PRECISION DDOT, DLAMCH |
|
|
185 EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH |
|
|
186 * .. |
|
|
187 * .. External Subroutines .. |
|
3333
|
188 EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA |
|
2329
|
189 * .. |
|
|
190 * .. Intrinsic Functions .. |
|
|
191 INTRINSIC ABS, MAX, SQRT |
|
|
192 * .. |
|
|
193 * .. Local Arrays .. |
|
|
194 DOUBLE PRECISION X( 2, 2 ) |
|
|
195 * .. |
|
|
196 * .. Executable Statements .. |
|
|
197 * |
|
|
198 * Decode and test the input parameters |
|
|
199 * |
|
|
200 BOTHV = LSAME( SIDE, 'B' ) |
|
|
201 RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV |
|
|
202 LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV |
|
|
203 * |
|
|
204 ALLV = LSAME( HOWMNY, 'A' ) |
|
3333
|
205 OVER = LSAME( HOWMNY, 'B' ) |
|
2329
|
206 SOMEV = LSAME( HOWMNY, 'S' ) |
|
|
207 * |
|
|
208 INFO = 0 |
|
|
209 IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN |
|
|
210 INFO = -1 |
|
|
211 ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN |
|
|
212 INFO = -2 |
|
|
213 ELSE IF( N.LT.0 ) THEN |
|
|
214 INFO = -4 |
|
|
215 ELSE IF( LDT.LT.MAX( 1, N ) ) THEN |
|
|
216 INFO = -6 |
|
|
217 ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN |
|
|
218 INFO = -8 |
|
|
219 ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN |
|
|
220 INFO = -10 |
|
|
221 ELSE |
|
|
222 * |
|
|
223 * Set M to the number of columns required to store the selected |
|
|
224 * eigenvectors, standardize the array SELECT if necessary, and |
|
|
225 * test MM. |
|
|
226 * |
|
|
227 IF( SOMEV ) THEN |
|
|
228 M = 0 |
|
|
229 PAIR = .FALSE. |
|
|
230 DO 10 J = 1, N |
|
|
231 IF( PAIR ) THEN |
|
|
232 PAIR = .FALSE. |
|
|
233 SELECT( J ) = .FALSE. |
|
|
234 ELSE |
|
|
235 IF( J.LT.N ) THEN |
|
|
236 IF( T( J+1, J ).EQ.ZERO ) THEN |
|
|
237 IF( SELECT( J ) ) |
|
|
238 $ M = M + 1 |
|
|
239 ELSE |
|
|
240 PAIR = .TRUE. |
|
|
241 IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN |
|
|
242 SELECT( J ) = .TRUE. |
|
|
243 M = M + 2 |
|
|
244 END IF |
|
|
245 END IF |
|
|
246 ELSE |
|
|
247 IF( SELECT( N ) ) |
|
|
248 $ M = M + 1 |
|
|
249 END IF |
|
|
250 END IF |
|
|
251 10 CONTINUE |
|
|
252 ELSE |
|
|
253 M = N |
|
|
254 END IF |
|
|
255 * |
|
|
256 IF( MM.LT.M ) THEN |
|
|
257 INFO = -11 |
|
|
258 END IF |
|
|
259 END IF |
|
|
260 IF( INFO.NE.0 ) THEN |
|
|
261 CALL XERBLA( 'DTREVC', -INFO ) |
|
|
262 RETURN |
|
|
263 END IF |
|
|
264 * |
|
|
265 * Quick return if possible. |
|
|
266 * |
|
|
267 IF( N.EQ.0 ) |
|
|
268 $ RETURN |
|
|
269 * |
|
|
270 * Set the constants to control overflow. |
|
|
271 * |
|
|
272 UNFL = DLAMCH( 'Safe minimum' ) |
|
|
273 OVFL = ONE / UNFL |
|
|
274 CALL DLABAD( UNFL, OVFL ) |
|
|
275 ULP = DLAMCH( 'Precision' ) |
|
|
276 SMLNUM = UNFL*( N / ULP ) |
|
|
277 BIGNUM = ( ONE-ULP ) / SMLNUM |
|
|
278 * |
|
|
279 * Compute 1-norm of each column of strictly upper triangular |
|
|
280 * part of T to control overflow in triangular solver. |
|
|
281 * |
|
|
282 WORK( 1 ) = ZERO |
|
|
283 DO 30 J = 2, N |
|
|
284 WORK( J ) = ZERO |
|
|
285 DO 20 I = 1, J - 1 |
|
|
286 WORK( J ) = WORK( J ) + ABS( T( I, J ) ) |
|
|
287 20 CONTINUE |
|
|
288 30 CONTINUE |
|
|
289 * |
|
|
290 * Index IP is used to specify the real or complex eigenvalue: |
|
|
291 * IP = 0, real eigenvalue, |
|
|
292 * 1, first of conjugate complex pair: (wr,wi) |
|
|
293 * -1, second of conjugate complex pair: (wr,wi) |
|
|
294 * |
|
|
295 N2 = 2*N |
|
|
296 * |
|
|
297 IF( RIGHTV ) THEN |
|
|
298 * |
|
|
299 * Compute right eigenvectors. |
|
|
300 * |
|
|
301 IP = 0 |
|
|
302 IS = M |
|
|
303 DO 140 KI = N, 1, -1 |
|
|
304 * |
|
|
305 IF( IP.EQ.1 ) |
|
|
306 $ GO TO 130 |
|
|
307 IF( KI.EQ.1 ) |
|
|
308 $ GO TO 40 |
|
|
309 IF( T( KI, KI-1 ).EQ.ZERO ) |
|
|
310 $ GO TO 40 |
|
|
311 IP = -1 |
|
|
312 * |
|
|
313 40 CONTINUE |
|
|
314 IF( SOMEV ) THEN |
|
|
315 IF( IP.EQ.0 ) THEN |
|
|
316 IF( .NOT.SELECT( KI ) ) |
|
|
317 $ GO TO 130 |
|
|
318 ELSE |
|
|
319 IF( .NOT.SELECT( KI-1 ) ) |
|
|
320 $ GO TO 130 |
|
|
321 END IF |
|
|
322 END IF |
|
|
323 * |
|
|
324 * Compute the KI-th eigenvalue (WR,WI). |
|
|
325 * |
|
|
326 WR = T( KI, KI ) |
|
|
327 WI = ZERO |
|
|
328 IF( IP.NE.0 ) |
|
|
329 $ WI = SQRT( ABS( T( KI, KI-1 ) ) )* |
|
|
330 $ SQRT( ABS( T( KI-1, KI ) ) ) |
|
|
331 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) |
|
|
332 * |
|
|
333 IF( IP.EQ.0 ) THEN |
|
|
334 * |
|
|
335 * Real right eigenvector |
|
|
336 * |
|
|
337 WORK( KI+N ) = ONE |
|
|
338 * |
|
|
339 * Form right-hand side |
|
|
340 * |
|
|
341 DO 50 K = 1, KI - 1 |
|
|
342 WORK( K+N ) = -T( K, KI ) |
|
|
343 50 CONTINUE |
|
|
344 * |
|
|
345 * Solve the upper quasi-triangular system: |
|
|
346 * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK. |
|
|
347 * |
|
|
348 JNXT = KI - 1 |
|
|
349 DO 60 J = KI - 1, 1, -1 |
|
|
350 IF( J.GT.JNXT ) |
|
|
351 $ GO TO 60 |
|
|
352 J1 = J |
|
|
353 J2 = J |
|
|
354 JNXT = J - 1 |
|
|
355 IF( J.GT.1 ) THEN |
|
|
356 IF( T( J, J-1 ).NE.ZERO ) THEN |
|
|
357 J1 = J - 1 |
|
|
358 JNXT = J - 2 |
|
|
359 END IF |
|
|
360 END IF |
|
|
361 * |
|
|
362 IF( J1.EQ.J2 ) THEN |
|
|
363 * |
|
|
364 * 1-by-1 diagonal block |
|
|
365 * |
|
|
366 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), |
|
|
367 $ LDT, ONE, ONE, WORK( J+N ), N, WR, |
|
|
368 $ ZERO, X, 2, SCALE, XNORM, IERR ) |
|
|
369 * |
|
|
370 * Scale X(1,1) to avoid overflow when updating |
|
|
371 * the right-hand side. |
|
|
372 * |
|
|
373 IF( XNORM.GT.ONE ) THEN |
|
|
374 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN |
|
|
375 X( 1, 1 ) = X( 1, 1 ) / XNORM |
|
|
376 SCALE = SCALE / XNORM |
|
|
377 END IF |
|
|
378 END IF |
|
|
379 * |
|
|
380 * Scale if necessary |
|
|
381 * |
|
|
382 IF( SCALE.NE.ONE ) |
|
|
383 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) |
|
|
384 WORK( J+N ) = X( 1, 1 ) |
|
|
385 * |
|
|
386 * Update right-hand side |
|
|
387 * |
|
|
388 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, |
|
|
389 $ WORK( 1+N ), 1 ) |
|
|
390 * |
|
|
391 ELSE |
|
|
392 * |
|
|
393 * 2-by-2 diagonal block |
|
|
394 * |
|
|
395 CALL DLALN2( .FALSE., 2, 1, SMIN, ONE, |
|
|
396 $ T( J-1, J-1 ), LDT, ONE, ONE, |
|
|
397 $ WORK( J-1+N ), N, WR, ZERO, X, 2, |
|
|
398 $ SCALE, XNORM, IERR ) |
|
|
399 * |
|
|
400 * Scale X(1,1) and X(2,1) to avoid overflow when |
|
|
401 * updating the right-hand side. |
|
|
402 * |
|
|
403 IF( XNORM.GT.ONE ) THEN |
|
|
404 BETA = MAX( WORK( J-1 ), WORK( J ) ) |
|
|
405 IF( BETA.GT.BIGNUM / XNORM ) THEN |
|
|
406 X( 1, 1 ) = X( 1, 1 ) / XNORM |
|
|
407 X( 2, 1 ) = X( 2, 1 ) / XNORM |
|
|
408 SCALE = SCALE / XNORM |
|
|
409 END IF |
|
|
410 END IF |
|
|
411 * |
|
|
412 * Scale if necessary |
|
|
413 * |
|
|
414 IF( SCALE.NE.ONE ) |
|
|
415 $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) |
|
|
416 WORK( J-1+N ) = X( 1, 1 ) |
|
|
417 WORK( J+N ) = X( 2, 1 ) |
|
|
418 * |
|
|
419 * Update right-hand side |
|
|
420 * |
|
|
421 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, |
|
|
422 $ WORK( 1+N ), 1 ) |
|
|
423 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, |
|
|
424 $ WORK( 1+N ), 1 ) |
|
|
425 END IF |
|
|
426 60 CONTINUE |
|
|
427 * |
|
|
428 * Copy the vector x or Q*x to VR and normalize. |
|
|
429 * |
|
|
430 IF( .NOT.OVER ) THEN |
|
|
431 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 ) |
|
|
432 * |
|
|
433 II = IDAMAX( KI, VR( 1, IS ), 1 ) |
|
|
434 REMAX = ONE / ABS( VR( II, IS ) ) |
|
|
435 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 ) |
|
|
436 * |
|
|
437 DO 70 K = KI + 1, N |
|
|
438 VR( K, IS ) = ZERO |
|
|
439 70 CONTINUE |
|
|
440 ELSE |
|
|
441 IF( KI.GT.1 ) |
|
|
442 $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR, |
|
|
443 $ WORK( 1+N ), 1, WORK( KI+N ), |
|
|
444 $ VR( 1, KI ), 1 ) |
|
|
445 * |
|
|
446 II = IDAMAX( N, VR( 1, KI ), 1 ) |
|
|
447 REMAX = ONE / ABS( VR( II, KI ) ) |
|
|
448 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 ) |
|
|
449 END IF |
|
|
450 * |
|
|
451 ELSE |
|
|
452 * |
|
|
453 * Complex right eigenvector. |
|
|
454 * |
|
|
455 * Initial solve |
|
|
456 * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0. |
|
|
457 * [ (T(KI,KI-1) T(KI,KI) ) ] |
|
|
458 * |
|
|
459 IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN |
|
|
460 WORK( KI-1+N ) = ONE |
|
|
461 WORK( KI+N2 ) = WI / T( KI-1, KI ) |
|
|
462 ELSE |
|
|
463 WORK( KI-1+N ) = -WI / T( KI, KI-1 ) |
|
|
464 WORK( KI+N2 ) = ONE |
|
|
465 END IF |
|
|
466 WORK( KI+N ) = ZERO |
|
|
467 WORK( KI-1+N2 ) = ZERO |
|
|
468 * |
|
|
469 * Form right-hand side |
|
|
470 * |
|
|
471 DO 80 K = 1, KI - 2 |
|
|
472 WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 ) |
|
|
473 WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI ) |
|
|
474 80 CONTINUE |
|
|
475 * |
|
|
476 * Solve upper quasi-triangular system: |
|
|
477 * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2) |
|
|
478 * |
|
|
479 JNXT = KI - 2 |
|
|
480 DO 90 J = KI - 2, 1, -1 |
|
|
481 IF( J.GT.JNXT ) |
|
|
482 $ GO TO 90 |
|
|
483 J1 = J |
|
|
484 J2 = J |
|
|
485 JNXT = J - 1 |
|
|
486 IF( J.GT.1 ) THEN |
|
|
487 IF( T( J, J-1 ).NE.ZERO ) THEN |
|
|
488 J1 = J - 1 |
|
|
489 JNXT = J - 2 |
|
|
490 END IF |
|
|
491 END IF |
|
|
492 * |
|
|
493 IF( J1.EQ.J2 ) THEN |
|
|
494 * |
|
|
495 * 1-by-1 diagonal block |
|
|
496 * |
|
|
497 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), |
|
|
498 $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI, |
|
|
499 $ X, 2, SCALE, XNORM, IERR ) |
|
|
500 * |
|
|
501 * Scale X(1,1) and X(1,2) to avoid overflow when |
|
|
502 * updating the right-hand side. |
|
|
503 * |
|
|
504 IF( XNORM.GT.ONE ) THEN |
|
|
505 IF( WORK( J ).GT.BIGNUM / XNORM ) THEN |
|
|
506 X( 1, 1 ) = X( 1, 1 ) / XNORM |
|
|
507 X( 1, 2 ) = X( 1, 2 ) / XNORM |
|
|
508 SCALE = SCALE / XNORM |
|
|
509 END IF |
|
|
510 END IF |
|
|
511 * |
|
|
512 * Scale if necessary |
|
|
513 * |
|
|
514 IF( SCALE.NE.ONE ) THEN |
|
|
515 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) |
|
|
516 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) |
|
|
517 END IF |
|
|
518 WORK( J+N ) = X( 1, 1 ) |
|
|
519 WORK( J+N2 ) = X( 1, 2 ) |
|
|
520 * |
|
|
521 * Update the right-hand side |
|
|
522 * |
|
|
523 CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1, |
|
|
524 $ WORK( 1+N ), 1 ) |
|
|
525 CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1, |
|
|
526 $ WORK( 1+N2 ), 1 ) |
|
|
527 * |
|
|
528 ELSE |
|
|
529 * |
|
|
530 * 2-by-2 diagonal block |
|
|
531 * |
|
|
532 CALL DLALN2( .FALSE., 2, 2, SMIN, ONE, |
|
|
533 $ T( J-1, J-1 ), LDT, ONE, ONE, |
|
|
534 $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE, |
|
|
535 $ XNORM, IERR ) |
|
|
536 * |
|
|
537 * Scale X to avoid overflow when updating |
|
|
538 * the right-hand side. |
|
|
539 * |
|
|
540 IF( XNORM.GT.ONE ) THEN |
|
|
541 BETA = MAX( WORK( J-1 ), WORK( J ) ) |
|
|
542 IF( BETA.GT.BIGNUM / XNORM ) THEN |
|
|
543 REC = ONE / XNORM |
|
|
544 X( 1, 1 ) = X( 1, 1 )*REC |
|
|
545 X( 1, 2 ) = X( 1, 2 )*REC |
|
|
546 X( 2, 1 ) = X( 2, 1 )*REC |
|
|
547 X( 2, 2 ) = X( 2, 2 )*REC |
|
|
548 SCALE = SCALE*REC |
|
|
549 END IF |
|
|
550 END IF |
|
|
551 * |
|
|
552 * Scale if necessary |
|
|
553 * |
|
|
554 IF( SCALE.NE.ONE ) THEN |
|
|
555 CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 ) |
|
|
556 CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 ) |
|
|
557 END IF |
|
|
558 WORK( J-1+N ) = X( 1, 1 ) |
|
|
559 WORK( J+N ) = X( 2, 1 ) |
|
|
560 WORK( J-1+N2 ) = X( 1, 2 ) |
|
|
561 WORK( J+N2 ) = X( 2, 2 ) |
|
|
562 * |
|
|
563 * Update the right-hand side |
|
|
564 * |
|
|
565 CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1, |
|
|
566 $ WORK( 1+N ), 1 ) |
|
|
567 CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1, |
|
|
568 $ WORK( 1+N ), 1 ) |
|
|
569 CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1, |
|
|
570 $ WORK( 1+N2 ), 1 ) |
|
|
571 CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1, |
|
|
572 $ WORK( 1+N2 ), 1 ) |
|
|
573 END IF |
|
|
574 90 CONTINUE |
|
|
575 * |
|
|
576 * Copy the vector x or Q*x to VR and normalize. |
|
|
577 * |
|
|
578 IF( .NOT.OVER ) THEN |
|
|
579 CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 ) |
|
|
580 CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 ) |
|
|
581 * |
|
|
582 EMAX = ZERO |
|
|
583 DO 100 K = 1, KI |
|
|
584 EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+ |
|
|
585 $ ABS( VR( K, IS ) ) ) |
|
|
586 100 CONTINUE |
|
|
587 * |
|
|
588 REMAX = ONE / EMAX |
|
|
589 CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 ) |
|
|
590 CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 ) |
|
|
591 * |
|
|
592 DO 110 K = KI + 1, N |
|
|
593 VR( K, IS-1 ) = ZERO |
|
|
594 VR( K, IS ) = ZERO |
|
|
595 110 CONTINUE |
|
|
596 * |
|
|
597 ELSE |
|
|
598 * |
|
|
599 IF( KI.GT.2 ) THEN |
|
|
600 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR, |
|
|
601 $ WORK( 1+N ), 1, WORK( KI-1+N ), |
|
|
602 $ VR( 1, KI-1 ), 1 ) |
|
|
603 CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR, |
|
|
604 $ WORK( 1+N2 ), 1, WORK( KI+N2 ), |
|
|
605 $ VR( 1, KI ), 1 ) |
|
|
606 ELSE |
|
|
607 CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 ) |
|
|
608 CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 ) |
|
|
609 END IF |
|
|
610 * |
|
|
611 EMAX = ZERO |
|
|
612 DO 120 K = 1, N |
|
|
613 EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+ |
|
|
614 $ ABS( VR( K, KI ) ) ) |
|
|
615 120 CONTINUE |
|
|
616 REMAX = ONE / EMAX |
|
|
617 CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 ) |
|
|
618 CALL DSCAL( N, REMAX, VR( 1, KI ), 1 ) |
|
|
619 END IF |
|
|
620 END IF |
|
|
621 * |
|
|
622 IS = IS - 1 |
|
|
623 IF( IP.NE.0 ) |
|
|
624 $ IS = IS - 1 |
|
|
625 130 CONTINUE |
|
|
626 IF( IP.EQ.1 ) |
|
|
627 $ IP = 0 |
|
|
628 IF( IP.EQ.-1 ) |
|
|
629 $ IP = 1 |
|
|
630 140 CONTINUE |
|
|
631 END IF |
|
|
632 * |
|
|
633 IF( LEFTV ) THEN |
|
|
634 * |
|
|
635 * Compute left eigenvectors. |
|
|
636 * |
|
|
637 IP = 0 |
|
|
638 IS = 1 |
|
|
639 DO 260 KI = 1, N |
|
|
640 * |
|
|
641 IF( IP.EQ.-1 ) |
|
|
642 $ GO TO 250 |
|
|
643 IF( KI.EQ.N ) |
|
|
644 $ GO TO 150 |
|
|
645 IF( T( KI+1, KI ).EQ.ZERO ) |
|
|
646 $ GO TO 150 |
|
|
647 IP = 1 |
|
|
648 * |
|
|
649 150 CONTINUE |
|
|
650 IF( SOMEV ) THEN |
|
|
651 IF( .NOT.SELECT( KI ) ) |
|
|
652 $ GO TO 250 |
|
|
653 END IF |
|
|
654 * |
|
|
655 * Compute the KI-th eigenvalue (WR,WI). |
|
|
656 * |
|
|
657 WR = T( KI, KI ) |
|
|
658 WI = ZERO |
|
|
659 IF( IP.NE.0 ) |
|
|
660 $ WI = SQRT( ABS( T( KI, KI+1 ) ) )* |
|
|
661 $ SQRT( ABS( T( KI+1, KI ) ) ) |
|
|
662 SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM ) |
|
|
663 * |
|
|
664 IF( IP.EQ.0 ) THEN |
|
|
665 * |
|
|
666 * Real left eigenvector. |
|
|
667 * |
|
|
668 WORK( KI+N ) = ONE |
|
|
669 * |
|
|
670 * Form right-hand side |
|
|
671 * |
|
|
672 DO 160 K = KI + 1, N |
|
|
673 WORK( K+N ) = -T( KI, K ) |
|
|
674 160 CONTINUE |
|
|
675 * |
|
|
676 * Solve the quasi-triangular system: |
|
|
677 * (T(KI+1:N,KI+1:N) - WR)'*X = SCALE*WORK |
|
|
678 * |
|
|
679 VMAX = ONE |
|
|
680 VCRIT = BIGNUM |
|
|
681 * |
|
|
682 JNXT = KI + 1 |
|
|
683 DO 170 J = KI + 1, N |
|
|
684 IF( J.LT.JNXT ) |
|
|
685 $ GO TO 170 |
|
|
686 J1 = J |
|
|
687 J2 = J |
|
|
688 JNXT = J + 1 |
|
|
689 IF( J.LT.N ) THEN |
|
|
690 IF( T( J+1, J ).NE.ZERO ) THEN |
|
|
691 J2 = J + 1 |
|
|
692 JNXT = J + 2 |
|
|
693 END IF |
|
|
694 END IF |
|
|
695 * |
|
|
696 IF( J1.EQ.J2 ) THEN |
|
|
697 * |
|
|
698 * 1-by-1 diagonal block |
|
|
699 * |
|
|
700 * Scale if necessary to avoid overflow when forming |
|
|
701 * the right-hand side. |
|
|
702 * |
|
|
703 IF( WORK( J ).GT.VCRIT ) THEN |
|
|
704 REC = ONE / VMAX |
|
|
705 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) |
|
|
706 VMAX = ONE |
|
|
707 VCRIT = BIGNUM |
|
|
708 END IF |
|
|
709 * |
|
|
710 WORK( J+N ) = WORK( J+N ) - |
|
|
711 $ DDOT( J-KI-1, T( KI+1, J ), 1, |
|
|
712 $ WORK( KI+1+N ), 1 ) |
|
|
713 * |
|
|
714 * Solve (T(J,J)-WR)'*X = WORK |
|
|
715 * |
|
|
716 CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ), |
|
|
717 $ LDT, ONE, ONE, WORK( J+N ), N, WR, |
|
|
718 $ ZERO, X, 2, SCALE, XNORM, IERR ) |
|
|
719 * |
|
|
720 * Scale if necessary |
|
|
721 * |
|
|
722 IF( SCALE.NE.ONE ) |
|
|
723 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) |
|
|
724 WORK( J+N ) = X( 1, 1 ) |
|
|
725 VMAX = MAX( ABS( WORK( J+N ) ), VMAX ) |
|
|
726 VCRIT = BIGNUM / VMAX |
|
|
727 * |
|
|
728 ELSE |
|
|
729 * |
|
|
730 * 2-by-2 diagonal block |
|
|
731 * |
|
|
732 * Scale if necessary to avoid overflow when forming |
|
|
733 * the right-hand side. |
|
|
734 * |
|
|
735 BETA = MAX( WORK( J ), WORK( J+1 ) ) |
|
|
736 IF( BETA.GT.VCRIT ) THEN |
|
|
737 REC = ONE / VMAX |
|
|
738 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) |
|
|
739 VMAX = ONE |
|
|
740 VCRIT = BIGNUM |
|
|
741 END IF |
|
|
742 * |
|
|
743 WORK( J+N ) = WORK( J+N ) - |
|
|
744 $ DDOT( J-KI-1, T( KI+1, J ), 1, |
|
|
745 $ WORK( KI+1+N ), 1 ) |
|
|
746 * |
|
|
747 WORK( J+1+N ) = WORK( J+1+N ) - |
|
|
748 $ DDOT( J-KI-1, T( KI+1, J+1 ), 1, |
|
|
749 $ WORK( KI+1+N ), 1 ) |
|
|
750 * |
|
|
751 * Solve |
|
|
752 * [T(J,J)-WR T(J,J+1) ]'* X = SCALE*( WORK1 ) |
|
|
753 * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 ) |
|
|
754 * |
|
|
755 CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ), |
|
|
756 $ LDT, ONE, ONE, WORK( J+N ), N, WR, |
|
|
757 $ ZERO, X, 2, SCALE, XNORM, IERR ) |
|
|
758 * |
|
|
759 * Scale if necessary |
|
|
760 * |
|
|
761 IF( SCALE.NE.ONE ) |
|
|
762 $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) |
|
|
763 WORK( J+N ) = X( 1, 1 ) |
|
|
764 WORK( J+1+N ) = X( 2, 1 ) |
|
|
765 * |
|
|
766 VMAX = MAX( ABS( WORK( J+N ) ), |
|
|
767 $ ABS( WORK( J+1+N ) ), VMAX ) |
|
|
768 VCRIT = BIGNUM / VMAX |
|
|
769 * |
|
|
770 END IF |
|
|
771 170 CONTINUE |
|
|
772 * |
|
|
773 * Copy the vector x or Q*x to VL and normalize. |
|
|
774 * |
|
|
775 IF( .NOT.OVER ) THEN |
|
|
776 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) |
|
|
777 * |
|
|
778 II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1 |
|
|
779 REMAX = ONE / ABS( VL( II, IS ) ) |
|
|
780 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) |
|
|
781 * |
|
|
782 DO 180 K = 1, KI - 1 |
|
|
783 VL( K, IS ) = ZERO |
|
|
784 180 CONTINUE |
|
|
785 * |
|
|
786 ELSE |
|
|
787 * |
|
|
788 IF( KI.LT.N ) |
|
|
789 $ CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL, |
|
|
790 $ WORK( KI+1+N ), 1, WORK( KI+N ), |
|
|
791 $ VL( 1, KI ), 1 ) |
|
|
792 * |
|
|
793 II = IDAMAX( N, VL( 1, KI ), 1 ) |
|
|
794 REMAX = ONE / ABS( VL( II, KI ) ) |
|
|
795 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 ) |
|
|
796 * |
|
|
797 END IF |
|
|
798 * |
|
|
799 ELSE |
|
|
800 * |
|
|
801 * Complex left eigenvector. |
|
|
802 * |
|
|
803 * Initial solve: |
|
|
804 * ((T(KI,KI) T(KI,KI+1) )' - (WR - I* WI))*X = 0. |
|
|
805 * ((T(KI+1,KI) T(KI+1,KI+1)) ) |
|
|
806 * |
|
|
807 IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN |
|
|
808 WORK( KI+N ) = WI / T( KI, KI+1 ) |
|
|
809 WORK( KI+1+N2 ) = ONE |
|
|
810 ELSE |
|
|
811 WORK( KI+N ) = ONE |
|
|
812 WORK( KI+1+N2 ) = -WI / T( KI+1, KI ) |
|
|
813 END IF |
|
|
814 WORK( KI+1+N ) = ZERO |
|
|
815 WORK( KI+N2 ) = ZERO |
|
|
816 * |
|
|
817 * Form right-hand side |
|
|
818 * |
|
|
819 DO 190 K = KI + 2, N |
|
|
820 WORK( K+N ) = -WORK( KI+N )*T( KI, K ) |
|
|
821 WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K ) |
|
|
822 190 CONTINUE |
|
|
823 * |
|
|
824 * Solve complex quasi-triangular system: |
|
|
825 * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2 |
|
|
826 * |
|
|
827 VMAX = ONE |
|
|
828 VCRIT = BIGNUM |
|
|
829 * |
|
|
830 JNXT = KI + 2 |
|
|
831 DO 200 J = KI + 2, N |
|
|
832 IF( J.LT.JNXT ) |
|
|
833 $ GO TO 200 |
|
|
834 J1 = J |
|
|
835 J2 = J |
|
|
836 JNXT = J + 1 |
|
|
837 IF( J.LT.N ) THEN |
|
|
838 IF( T( J+1, J ).NE.ZERO ) THEN |
|
|
839 J2 = J + 1 |
|
|
840 JNXT = J + 2 |
|
|
841 END IF |
|
|
842 END IF |
|
|
843 * |
|
|
844 IF( J1.EQ.J2 ) THEN |
|
|
845 * |
|
|
846 * 1-by-1 diagonal block |
|
|
847 * |
|
|
848 * Scale if necessary to avoid overflow when |
|
|
849 * forming the right-hand side elements. |
|
|
850 * |
|
|
851 IF( WORK( J ).GT.VCRIT ) THEN |
|
|
852 REC = ONE / VMAX |
|
|
853 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) |
|
|
854 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) |
|
|
855 VMAX = ONE |
|
|
856 VCRIT = BIGNUM |
|
|
857 END IF |
|
|
858 * |
|
|
859 WORK( J+N ) = WORK( J+N ) - |
|
|
860 $ DDOT( J-KI-2, T( KI+2, J ), 1, |
|
|
861 $ WORK( KI+2+N ), 1 ) |
|
|
862 WORK( J+N2 ) = WORK( J+N2 ) - |
|
|
863 $ DDOT( J-KI-2, T( KI+2, J ), 1, |
|
|
864 $ WORK( KI+2+N2 ), 1 ) |
|
|
865 * |
|
|
866 * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2 |
|
|
867 * |
|
|
868 CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ), |
|
|
869 $ LDT, ONE, ONE, WORK( J+N ), N, WR, |
|
|
870 $ -WI, X, 2, SCALE, XNORM, IERR ) |
|
|
871 * |
|
|
872 * Scale if necessary |
|
|
873 * |
|
|
874 IF( SCALE.NE.ONE ) THEN |
|
|
875 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) |
|
|
876 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) |
|
|
877 END IF |
|
|
878 WORK( J+N ) = X( 1, 1 ) |
|
|
879 WORK( J+N2 ) = X( 1, 2 ) |
|
|
880 VMAX = MAX( ABS( WORK( J+N ) ), |
|
|
881 $ ABS( WORK( J+N2 ) ), VMAX ) |
|
|
882 VCRIT = BIGNUM / VMAX |
|
|
883 * |
|
|
884 ELSE |
|
|
885 * |
|
|
886 * 2-by-2 diagonal block |
|
|
887 * |
|
|
888 * Scale if necessary to avoid overflow when forming |
|
|
889 * the right-hand side elements. |
|
|
890 * |
|
|
891 BETA = MAX( WORK( J ), WORK( J+1 ) ) |
|
|
892 IF( BETA.GT.VCRIT ) THEN |
|
|
893 REC = ONE / VMAX |
|
|
894 CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 ) |
|
|
895 CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 ) |
|
|
896 VMAX = ONE |
|
|
897 VCRIT = BIGNUM |
|
|
898 END IF |
|
|
899 * |
|
|
900 WORK( J+N ) = WORK( J+N ) - |
|
|
901 $ DDOT( J-KI-2, T( KI+2, J ), 1, |
|
|
902 $ WORK( KI+2+N ), 1 ) |
|
|
903 * |
|
|
904 WORK( J+N2 ) = WORK( J+N2 ) - |
|
|
905 $ DDOT( J-KI-2, T( KI+2, J ), 1, |
|
|
906 $ WORK( KI+2+N2 ), 1 ) |
|
|
907 * |
|
|
908 WORK( J+1+N ) = WORK( J+1+N ) - |
|
|
909 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1, |
|
|
910 $ WORK( KI+2+N ), 1 ) |
|
|
911 * |
|
|
912 WORK( J+1+N2 ) = WORK( J+1+N2 ) - |
|
|
913 $ DDOT( J-KI-2, T( KI+2, J+1 ), 1, |
|
|
914 $ WORK( KI+2+N2 ), 1 ) |
|
|
915 * |
|
|
916 * Solve 2-by-2 complex linear equation |
|
|
917 * ([T(j,j) T(j,j+1) ]'-(wr-i*wi)*I)*X = SCALE*B |
|
|
918 * ([T(j+1,j) T(j+1,j+1)] ) |
|
|
919 * |
|
|
920 CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ), |
|
|
921 $ LDT, ONE, ONE, WORK( J+N ), N, WR, |
|
|
922 $ -WI, X, 2, SCALE, XNORM, IERR ) |
|
|
923 * |
|
|
924 * Scale if necessary |
|
|
925 * |
|
|
926 IF( SCALE.NE.ONE ) THEN |
|
|
927 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 ) |
|
|
928 CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 ) |
|
|
929 END IF |
|
|
930 WORK( J+N ) = X( 1, 1 ) |
|
|
931 WORK( J+N2 ) = X( 1, 2 ) |
|
|
932 WORK( J+1+N ) = X( 2, 1 ) |
|
|
933 WORK( J+1+N2 ) = X( 2, 2 ) |
|
|
934 VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ), |
|
|
935 $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX ) |
|
|
936 VCRIT = BIGNUM / VMAX |
|
|
937 * |
|
|
938 END IF |
|
|
939 200 CONTINUE |
|
|
940 * |
|
|
941 * Copy the vector x or Q*x to VL and normalize. |
|
|
942 * |
|
|
943 210 CONTINUE |
|
|
944 IF( .NOT.OVER ) THEN |
|
|
945 CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 ) |
|
|
946 CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ), |
|
|
947 $ 1 ) |
|
|
948 * |
|
|
949 EMAX = ZERO |
|
|
950 DO 220 K = KI, N |
|
|
951 EMAX = MAX( EMAX, ABS( VL( K, IS ) )+ |
|
|
952 $ ABS( VL( K, IS+1 ) ) ) |
|
|
953 220 CONTINUE |
|
|
954 REMAX = ONE / EMAX |
|
|
955 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 ) |
|
|
956 CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 ) |
|
|
957 * |
|
|
958 DO 230 K = 1, KI - 1 |
|
|
959 VL( K, IS ) = ZERO |
|
|
960 VL( K, IS+1 ) = ZERO |
|
|
961 230 CONTINUE |
|
|
962 ELSE |
|
|
963 IF( KI.LT.N-1 ) THEN |
|
|
964 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), |
|
|
965 $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ), |
|
|
966 $ VL( 1, KI ), 1 ) |
|
|
967 CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ), |
|
|
968 $ LDVL, WORK( KI+2+N2 ), 1, |
|
|
969 $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) |
|
|
970 ELSE |
|
|
971 CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 ) |
|
|
972 CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 ) |
|
|
973 END IF |
|
|
974 * |
|
|
975 EMAX = ZERO |
|
|
976 DO 240 K = 1, N |
|
|
977 EMAX = MAX( EMAX, ABS( VL( K, KI ) )+ |
|
|
978 $ ABS( VL( K, KI+1 ) ) ) |
|
|
979 240 CONTINUE |
|
|
980 REMAX = ONE / EMAX |
|
|
981 CALL DSCAL( N, REMAX, VL( 1, KI ), 1 ) |
|
|
982 CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 ) |
|
|
983 * |
|
|
984 END IF |
|
|
985 * |
|
|
986 END IF |
|
|
987 * |
|
|
988 IS = IS + 1 |
|
|
989 IF( IP.NE.0 ) |
|
|
990 $ IS = IS + 1 |
|
|
991 250 CONTINUE |
|
|
992 IF( IP.EQ.-1 ) |
|
|
993 $ IP = 0 |
|
|
994 IF( IP.EQ.1 ) |
|
|
995 $ IP = -1 |
|
|
996 * |
|
|
997 260 CONTINUE |
|
|
998 * |
|
|
999 END IF |
|
|
1000 * |
|
|
1001 RETURN |
|
|
1002 * |
|
|
1003 * End of DTREVC |
|
|
1004 * |
|
|
1005 END |
