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