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