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