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**T)*T = w*(y**T)
49: *>
50: *> where y**T denotes the 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: *> \date November 2011
204: *
205: *> \ingroup doubleOTHERcomputational
206: *
207: *> \par Further Details:
208: * =====================
209: *>
210: *> \verbatim
211: *>
212: *> The algorithm used in this program is basically backward (forward)
213: *> substitution, with scaling to make the the code robust against
214: *> possible overflow.
215: *>
216: *> Each eigenvector is normalized so that the element of largest
217: *> magnitude has magnitude 1; here the magnitude of a complex number
218: *> (x,y) is taken to be |x| + |y|.
219: *> \endverbatim
220: *>
221: * =====================================================================
222: SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
223: $ LDVR, MM, M, WORK, INFO )
224: *
225: * -- LAPACK computational routine (version 3.4.0) --
226: * -- LAPACK is a software package provided by Univ. of Tennessee, --
227: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
228: * November 2011
229: *
230: * .. Scalar Arguments ..
231: CHARACTER HOWMNY, SIDE
232: INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
233: * ..
234: * .. Array Arguments ..
235: LOGICAL SELECT( * )
236: DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
237: $ WORK( * )
238: * ..
239: *
240: * =====================================================================
241: *
242: * .. Parameters ..
243: DOUBLE PRECISION ZERO, ONE
244: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
245: * ..
246: * .. Local Scalars ..
247: LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
248: INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
249: DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
250: $ SMIN, SMLNUM, ULP, UNFL, VCRIT, VMAX, WI, WR,
251: $ XNORM
252: * ..
253: * .. External Functions ..
254: LOGICAL LSAME
255: INTEGER IDAMAX
256: DOUBLE PRECISION DDOT, DLAMCH
257: EXTERNAL LSAME, IDAMAX, DDOT, DLAMCH
258: * ..
259: * .. External Subroutines ..
260: EXTERNAL DAXPY, DCOPY, DGEMV, DLALN2, DSCAL, XERBLA
261: * ..
262: * .. Intrinsic Functions ..
263: INTRINSIC ABS, MAX, SQRT
264: * ..
265: * .. Local Arrays ..
266: DOUBLE PRECISION X( 2, 2 )
267: * ..
268: * .. Executable Statements ..
269: *
270: * Decode and test the input parameters
271: *
272: BOTHV = LSAME( SIDE, 'B' )
273: RIGHTV = LSAME( SIDE, 'R' ) .OR. BOTHV
274: LEFTV = LSAME( SIDE, 'L' ) .OR. BOTHV
275: *
276: ALLV = LSAME( HOWMNY, 'A' )
277: OVER = LSAME( HOWMNY, 'B' )
278: SOMEV = LSAME( HOWMNY, 'S' )
279: *
280: INFO = 0
281: IF( .NOT.RIGHTV .AND. .NOT.LEFTV ) THEN
282: INFO = -1
283: ELSE IF( .NOT.ALLV .AND. .NOT.OVER .AND. .NOT.SOMEV ) THEN
284: INFO = -2
285: ELSE IF( N.LT.0 ) THEN
286: INFO = -4
287: ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
288: INFO = -6
289: ELSE IF( LDVL.LT.1 .OR. ( LEFTV .AND. LDVL.LT.N ) ) THEN
290: INFO = -8
291: ELSE IF( LDVR.LT.1 .OR. ( RIGHTV .AND. LDVR.LT.N ) ) THEN
292: INFO = -10
293: ELSE
294: *
295: * Set M to the number of columns required to store the selected
296: * eigenvectors, standardize the array SELECT if necessary, and
297: * test MM.
298: *
299: IF( SOMEV ) THEN
300: M = 0
301: PAIR = .FALSE.
302: DO 10 J = 1, N
303: IF( PAIR ) THEN
304: PAIR = .FALSE.
305: SELECT( J ) = .FALSE.
306: ELSE
307: IF( J.LT.N ) THEN
308: IF( T( J+1, J ).EQ.ZERO ) THEN
309: IF( SELECT( J ) )
310: $ M = M + 1
311: ELSE
312: PAIR = .TRUE.
313: IF( SELECT( J ) .OR. SELECT( J+1 ) ) THEN
314: SELECT( J ) = .TRUE.
315: M = M + 2
316: END IF
317: END IF
318: ELSE
319: IF( SELECT( N ) )
320: $ M = M + 1
321: END IF
322: END IF
323: 10 CONTINUE
324: ELSE
325: M = N
326: END IF
327: *
328: IF( MM.LT.M ) THEN
329: INFO = -11
330: END IF
331: END IF
332: IF( INFO.NE.0 ) THEN
333: CALL XERBLA( 'DTREVC', -INFO )
334: RETURN
335: END IF
336: *
337: * Quick return if possible.
338: *
339: IF( N.EQ.0 )
340: $ RETURN
341: *
342: * Set the constants to control overflow.
343: *
344: UNFL = DLAMCH( 'Safe minimum' )
345: OVFL = ONE / UNFL
346: CALL DLABAD( UNFL, OVFL )
347: ULP = DLAMCH( 'Precision' )
348: SMLNUM = UNFL*( N / ULP )
349: BIGNUM = ( ONE-ULP ) / SMLNUM
350: *
351: * Compute 1-norm of each column of strictly upper triangular
352: * part of T to control overflow in triangular solver.
353: *
354: WORK( 1 ) = ZERO
355: DO 30 J = 2, N
356: WORK( J ) = ZERO
357: DO 20 I = 1, J - 1
358: WORK( J ) = WORK( J ) + ABS( T( I, J ) )
359: 20 CONTINUE
360: 30 CONTINUE
361: *
362: * Index IP is used to specify the real or complex eigenvalue:
363: * IP = 0, real eigenvalue,
364: * 1, first of conjugate complex pair: (wr,wi)
365: * -1, second of conjugate complex pair: (wr,wi)
366: *
367: N2 = 2*N
368: *
369: IF( RIGHTV ) THEN
370: *
371: * Compute right eigenvectors.
372: *
373: IP = 0
374: IS = M
375: DO 140 KI = N, 1, -1
376: *
377: IF( IP.EQ.1 )
378: $ GO TO 130
379: IF( KI.EQ.1 )
380: $ GO TO 40
381: IF( T( KI, KI-1 ).EQ.ZERO )
382: $ GO TO 40
383: IP = -1
384: *
385: 40 CONTINUE
386: IF( SOMEV ) THEN
387: IF( IP.EQ.0 ) THEN
388: IF( .NOT.SELECT( KI ) )
389: $ GO TO 130
390: ELSE
391: IF( .NOT.SELECT( KI-1 ) )
392: $ GO TO 130
393: END IF
394: END IF
395: *
396: * Compute the KI-th eigenvalue (WR,WI).
397: *
398: WR = T( KI, KI )
399: WI = ZERO
400: IF( IP.NE.0 )
401: $ WI = SQRT( ABS( T( KI, KI-1 ) ) )*
402: $ SQRT( ABS( T( KI-1, KI ) ) )
403: SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
404: *
405: IF( IP.EQ.0 ) THEN
406: *
407: * Real right eigenvector
408: *
409: WORK( KI+N ) = ONE
410: *
411: * Form right-hand side
412: *
413: DO 50 K = 1, KI - 1
414: WORK( K+N ) = -T( K, KI )
415: 50 CONTINUE
416: *
417: * Solve the upper quasi-triangular system:
418: * (T(1:KI-1,1:KI-1) - WR)*X = SCALE*WORK.
419: *
420: JNXT = KI - 1
421: DO 60 J = KI - 1, 1, -1
422: IF( J.GT.JNXT )
423: $ GO TO 60
424: J1 = J
425: J2 = J
426: JNXT = J - 1
427: IF( J.GT.1 ) THEN
428: IF( T( J, J-1 ).NE.ZERO ) THEN
429: J1 = J - 1
430: JNXT = J - 2
431: END IF
432: END IF
433: *
434: IF( J1.EQ.J2 ) THEN
435: *
436: * 1-by-1 diagonal block
437: *
438: CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
439: $ LDT, ONE, ONE, WORK( J+N ), N, WR,
440: $ ZERO, X, 2, SCALE, XNORM, IERR )
441: *
442: * Scale X(1,1) to avoid overflow when updating
443: * the right-hand side.
444: *
445: IF( XNORM.GT.ONE ) THEN
446: IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
447: X( 1, 1 ) = X( 1, 1 ) / XNORM
448: SCALE = SCALE / XNORM
449: END IF
450: END IF
451: *
452: * Scale if necessary
453: *
454: IF( SCALE.NE.ONE )
455: $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
456: WORK( J+N ) = X( 1, 1 )
457: *
458: * Update right-hand side
459: *
460: CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
461: $ WORK( 1+N ), 1 )
462: *
463: ELSE
464: *
465: * 2-by-2 diagonal block
466: *
467: CALL DLALN2( .FALSE., 2, 1, SMIN, ONE,
468: $ T( J-1, J-1 ), LDT, ONE, ONE,
469: $ WORK( J-1+N ), N, WR, ZERO, X, 2,
470: $ SCALE, XNORM, IERR )
471: *
472: * Scale X(1,1) and X(2,1) to avoid overflow when
473: * updating the right-hand side.
474: *
475: IF( XNORM.GT.ONE ) THEN
476: BETA = MAX( WORK( J-1 ), WORK( J ) )
477: IF( BETA.GT.BIGNUM / XNORM ) THEN
478: X( 1, 1 ) = X( 1, 1 ) / XNORM
479: X( 2, 1 ) = X( 2, 1 ) / XNORM
480: SCALE = SCALE / XNORM
481: END IF
482: END IF
483: *
484: * Scale if necessary
485: *
486: IF( SCALE.NE.ONE )
487: $ CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
488: WORK( J-1+N ) = X( 1, 1 )
489: WORK( J+N ) = X( 2, 1 )
490: *
491: * Update right-hand side
492: *
493: CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
494: $ WORK( 1+N ), 1 )
495: CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
496: $ WORK( 1+N ), 1 )
497: END IF
498: 60 CONTINUE
499: *
500: * Copy the vector x or Q*x to VR and normalize.
501: *
502: IF( .NOT.OVER ) THEN
503: CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS ), 1 )
504: *
505: II = IDAMAX( KI, VR( 1, IS ), 1 )
506: REMAX = ONE / ABS( VR( II, IS ) )
507: CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
508: *
509: DO 70 K = KI + 1, N
510: VR( K, IS ) = ZERO
511: 70 CONTINUE
512: ELSE
513: IF( KI.GT.1 )
514: $ CALL DGEMV( 'N', N, KI-1, ONE, VR, LDVR,
515: $ WORK( 1+N ), 1, WORK( KI+N ),
516: $ VR( 1, KI ), 1 )
517: *
518: II = IDAMAX( N, VR( 1, KI ), 1 )
519: REMAX = ONE / ABS( VR( II, KI ) )
520: CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
521: END IF
522: *
523: ELSE
524: *
525: * Complex right eigenvector.
526: *
527: * Initial solve
528: * [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
529: * [ (T(KI,KI-1) T(KI,KI) ) ]
530: *
531: IF( ABS( T( KI-1, KI ) ).GE.ABS( T( KI, KI-1 ) ) ) THEN
532: WORK( KI-1+N ) = ONE
533: WORK( KI+N2 ) = WI / T( KI-1, KI )
534: ELSE
535: WORK( KI-1+N ) = -WI / T( KI, KI-1 )
536: WORK( KI+N2 ) = ONE
537: END IF
538: WORK( KI+N ) = ZERO
539: WORK( KI-1+N2 ) = ZERO
540: *
541: * Form right-hand side
542: *
543: DO 80 K = 1, KI - 2
544: WORK( K+N ) = -WORK( KI-1+N )*T( K, KI-1 )
545: WORK( K+N2 ) = -WORK( KI+N2 )*T( K, KI )
546: 80 CONTINUE
547: *
548: * Solve upper quasi-triangular system:
549: * (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
550: *
551: JNXT = KI - 2
552: DO 90 J = KI - 2, 1, -1
553: IF( J.GT.JNXT )
554: $ GO TO 90
555: J1 = J
556: J2 = J
557: JNXT = J - 1
558: IF( J.GT.1 ) THEN
559: IF( T( J, J-1 ).NE.ZERO ) THEN
560: J1 = J - 1
561: JNXT = J - 2
562: END IF
563: END IF
564: *
565: IF( J1.EQ.J2 ) THEN
566: *
567: * 1-by-1 diagonal block
568: *
569: CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
570: $ LDT, ONE, ONE, WORK( J+N ), N, WR, WI,
571: $ X, 2, SCALE, XNORM, IERR )
572: *
573: * Scale X(1,1) and X(1,2) to avoid overflow when
574: * updating the right-hand side.
575: *
576: IF( XNORM.GT.ONE ) THEN
577: IF( WORK( J ).GT.BIGNUM / XNORM ) THEN
578: X( 1, 1 ) = X( 1, 1 ) / XNORM
579: X( 1, 2 ) = X( 1, 2 ) / XNORM
580: SCALE = SCALE / XNORM
581: END IF
582: END IF
583: *
584: * Scale if necessary
585: *
586: IF( SCALE.NE.ONE ) THEN
587: CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
588: CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
589: END IF
590: WORK( J+N ) = X( 1, 1 )
591: WORK( J+N2 ) = X( 1, 2 )
592: *
593: * Update the right-hand side
594: *
595: CALL DAXPY( J-1, -X( 1, 1 ), T( 1, J ), 1,
596: $ WORK( 1+N ), 1 )
597: CALL DAXPY( J-1, -X( 1, 2 ), T( 1, J ), 1,
598: $ WORK( 1+N2 ), 1 )
599: *
600: ELSE
601: *
602: * 2-by-2 diagonal block
603: *
604: CALL DLALN2( .FALSE., 2, 2, SMIN, ONE,
605: $ T( J-1, J-1 ), LDT, ONE, ONE,
606: $ WORK( J-1+N ), N, WR, WI, X, 2, SCALE,
607: $ XNORM, IERR )
608: *
609: * Scale X to avoid overflow when updating
610: * the right-hand side.
611: *
612: IF( XNORM.GT.ONE ) THEN
613: BETA = MAX( WORK( J-1 ), WORK( J ) )
614: IF( BETA.GT.BIGNUM / XNORM ) THEN
615: REC = ONE / XNORM
616: X( 1, 1 ) = X( 1, 1 )*REC
617: X( 1, 2 ) = X( 1, 2 )*REC
618: X( 2, 1 ) = X( 2, 1 )*REC
619: X( 2, 2 ) = X( 2, 2 )*REC
620: SCALE = SCALE*REC
621: END IF
622: END IF
623: *
624: * Scale if necessary
625: *
626: IF( SCALE.NE.ONE ) THEN
627: CALL DSCAL( KI, SCALE, WORK( 1+N ), 1 )
628: CALL DSCAL( KI, SCALE, WORK( 1+N2 ), 1 )
629: END IF
630: WORK( J-1+N ) = X( 1, 1 )
631: WORK( J+N ) = X( 2, 1 )
632: WORK( J-1+N2 ) = X( 1, 2 )
633: WORK( J+N2 ) = X( 2, 2 )
634: *
635: * Update the right-hand side
636: *
637: CALL DAXPY( J-2, -X( 1, 1 ), T( 1, J-1 ), 1,
638: $ WORK( 1+N ), 1 )
639: CALL DAXPY( J-2, -X( 2, 1 ), T( 1, J ), 1,
640: $ WORK( 1+N ), 1 )
641: CALL DAXPY( J-2, -X( 1, 2 ), T( 1, J-1 ), 1,
642: $ WORK( 1+N2 ), 1 )
643: CALL DAXPY( J-2, -X( 2, 2 ), T( 1, J ), 1,
644: $ WORK( 1+N2 ), 1 )
645: END IF
646: 90 CONTINUE
647: *
648: * Copy the vector x or Q*x to VR and normalize.
649: *
650: IF( .NOT.OVER ) THEN
651: CALL DCOPY( KI, WORK( 1+N ), 1, VR( 1, IS-1 ), 1 )
652: CALL DCOPY( KI, WORK( 1+N2 ), 1, VR( 1, IS ), 1 )
653: *
654: EMAX = ZERO
655: DO 100 K = 1, KI
656: EMAX = MAX( EMAX, ABS( VR( K, IS-1 ) )+
657: $ ABS( VR( K, IS ) ) )
658: 100 CONTINUE
659: *
660: REMAX = ONE / EMAX
661: CALL DSCAL( KI, REMAX, VR( 1, IS-1 ), 1 )
662: CALL DSCAL( KI, REMAX, VR( 1, IS ), 1 )
663: *
664: DO 110 K = KI + 1, N
665: VR( K, IS-1 ) = ZERO
666: VR( K, IS ) = ZERO
667: 110 CONTINUE
668: *
669: ELSE
670: *
671: IF( KI.GT.2 ) THEN
672: CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
673: $ WORK( 1+N ), 1, WORK( KI-1+N ),
674: $ VR( 1, KI-1 ), 1 )
675: CALL DGEMV( 'N', N, KI-2, ONE, VR, LDVR,
676: $ WORK( 1+N2 ), 1, WORK( KI+N2 ),
677: $ VR( 1, KI ), 1 )
678: ELSE
679: CALL DSCAL( N, WORK( KI-1+N ), VR( 1, KI-1 ), 1 )
680: CALL DSCAL( N, WORK( KI+N2 ), VR( 1, KI ), 1 )
681: END IF
682: *
683: EMAX = ZERO
684: DO 120 K = 1, N
685: EMAX = MAX( EMAX, ABS( VR( K, KI-1 ) )+
686: $ ABS( VR( K, KI ) ) )
687: 120 CONTINUE
688: REMAX = ONE / EMAX
689: CALL DSCAL( N, REMAX, VR( 1, KI-1 ), 1 )
690: CALL DSCAL( N, REMAX, VR( 1, KI ), 1 )
691: END IF
692: END IF
693: *
694: IS = IS - 1
695: IF( IP.NE.0 )
696: $ IS = IS - 1
697: 130 CONTINUE
698: IF( IP.EQ.1 )
699: $ IP = 0
700: IF( IP.EQ.-1 )
701: $ IP = 1
702: 140 CONTINUE
703: END IF
704: *
705: IF( LEFTV ) THEN
706: *
707: * Compute left eigenvectors.
708: *
709: IP = 0
710: IS = 1
711: DO 260 KI = 1, N
712: *
713: IF( IP.EQ.-1 )
714: $ GO TO 250
715: IF( KI.EQ.N )
716: $ GO TO 150
717: IF( T( KI+1, KI ).EQ.ZERO )
718: $ GO TO 150
719: IP = 1
720: *
721: 150 CONTINUE
722: IF( SOMEV ) THEN
723: IF( .NOT.SELECT( KI ) )
724: $ GO TO 250
725: END IF
726: *
727: * Compute the KI-th eigenvalue (WR,WI).
728: *
729: WR = T( KI, KI )
730: WI = ZERO
731: IF( IP.NE.0 )
732: $ WI = SQRT( ABS( T( KI, KI+1 ) ) )*
733: $ SQRT( ABS( T( KI+1, KI ) ) )
734: SMIN = MAX( ULP*( ABS( WR )+ABS( WI ) ), SMLNUM )
735: *
736: IF( IP.EQ.0 ) THEN
737: *
738: * Real left eigenvector.
739: *
740: WORK( KI+N ) = ONE
741: *
742: * Form right-hand side
743: *
744: DO 160 K = KI + 1, N
745: WORK( K+N ) = -T( KI, K )
746: 160 CONTINUE
747: *
748: * Solve the quasi-triangular system:
749: * (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
750: *
751: VMAX = ONE
752: VCRIT = BIGNUM
753: *
754: JNXT = KI + 1
755: DO 170 J = KI + 1, N
756: IF( J.LT.JNXT )
757: $ GO TO 170
758: J1 = J
759: J2 = J
760: JNXT = J + 1
761: IF( J.LT.N ) THEN
762: IF( T( J+1, J ).NE.ZERO ) THEN
763: J2 = J + 1
764: JNXT = J + 2
765: END IF
766: END IF
767: *
768: IF( J1.EQ.J2 ) THEN
769: *
770: * 1-by-1 diagonal block
771: *
772: * Scale if necessary to avoid overflow when forming
773: * the right-hand side.
774: *
775: IF( WORK( J ).GT.VCRIT ) THEN
776: REC = ONE / VMAX
777: CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
778: VMAX = ONE
779: VCRIT = BIGNUM
780: END IF
781: *
782: WORK( J+N ) = WORK( J+N ) -
783: $ DDOT( J-KI-1, T( KI+1, J ), 1,
784: $ WORK( KI+1+N ), 1 )
785: *
786: * Solve (T(J,J)-WR)**T*X = WORK
787: *
788: CALL DLALN2( .FALSE., 1, 1, SMIN, ONE, T( J, J ),
789: $ LDT, ONE, ONE, WORK( J+N ), N, WR,
790: $ ZERO, X, 2, SCALE, XNORM, IERR )
791: *
792: * Scale if necessary
793: *
794: IF( SCALE.NE.ONE )
795: $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
796: WORK( J+N ) = X( 1, 1 )
797: VMAX = MAX( ABS( WORK( J+N ) ), VMAX )
798: VCRIT = BIGNUM / VMAX
799: *
800: ELSE
801: *
802: * 2-by-2 diagonal block
803: *
804: * Scale if necessary to avoid overflow when forming
805: * the right-hand side.
806: *
807: BETA = MAX( WORK( J ), WORK( J+1 ) )
808: IF( BETA.GT.VCRIT ) THEN
809: REC = ONE / VMAX
810: CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
811: VMAX = ONE
812: VCRIT = BIGNUM
813: END IF
814: *
815: WORK( J+N ) = WORK( J+N ) -
816: $ DDOT( J-KI-1, T( KI+1, J ), 1,
817: $ WORK( KI+1+N ), 1 )
818: *
819: WORK( J+1+N ) = WORK( J+1+N ) -
820: $ DDOT( J-KI-1, T( KI+1, J+1 ), 1,
821: $ WORK( KI+1+N ), 1 )
822: *
823: * Solve
824: * [T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
825: * [T(J+1,J) T(J+1,J+1)-WR] ( WORK2 )
826: *
827: CALL DLALN2( .TRUE., 2, 1, SMIN, ONE, T( J, J ),
828: $ LDT, ONE, ONE, WORK( J+N ), N, WR,
829: $ ZERO, X, 2, SCALE, XNORM, IERR )
830: *
831: * Scale if necessary
832: *
833: IF( SCALE.NE.ONE )
834: $ CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
835: WORK( J+N ) = X( 1, 1 )
836: WORK( J+1+N ) = X( 2, 1 )
837: *
838: VMAX = MAX( ABS( WORK( J+N ) ),
839: $ ABS( WORK( J+1+N ) ), VMAX )
840: VCRIT = BIGNUM / VMAX
841: *
842: END IF
843: 170 CONTINUE
844: *
845: * Copy the vector x or Q*x to VL and normalize.
846: *
847: IF( .NOT.OVER ) THEN
848: CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
849: *
850: II = IDAMAX( N-KI+1, VL( KI, IS ), 1 ) + KI - 1
851: REMAX = ONE / ABS( VL( II, IS ) )
852: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
853: *
854: DO 180 K = 1, KI - 1
855: VL( K, IS ) = ZERO
856: 180 CONTINUE
857: *
858: ELSE
859: *
860: IF( KI.LT.N )
861: $ CALL DGEMV( 'N', N, N-KI, ONE, VL( 1, KI+1 ), LDVL,
862: $ WORK( KI+1+N ), 1, WORK( KI+N ),
863: $ VL( 1, KI ), 1 )
864: *
865: II = IDAMAX( N, VL( 1, KI ), 1 )
866: REMAX = ONE / ABS( VL( II, KI ) )
867: CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
868: *
869: END IF
870: *
871: ELSE
872: *
873: * Complex left eigenvector.
874: *
875: * Initial solve:
876: * ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
877: * ((T(KI+1,KI) T(KI+1,KI+1)) )
878: *
879: IF( ABS( T( KI, KI+1 ) ).GE.ABS( T( KI+1, KI ) ) ) THEN
880: WORK( KI+N ) = WI / T( KI, KI+1 )
881: WORK( KI+1+N2 ) = ONE
882: ELSE
883: WORK( KI+N ) = ONE
884: WORK( KI+1+N2 ) = -WI / T( KI+1, KI )
885: END IF
886: WORK( KI+1+N ) = ZERO
887: WORK( KI+N2 ) = ZERO
888: *
889: * Form right-hand side
890: *
891: DO 190 K = KI + 2, N
892: WORK( K+N ) = -WORK( KI+N )*T( KI, K )
893: WORK( K+N2 ) = -WORK( KI+1+N2 )*T( KI+1, K )
894: 190 CONTINUE
895: *
896: * Solve complex quasi-triangular system:
897: * ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
898: *
899: VMAX = ONE
900: VCRIT = BIGNUM
901: *
902: JNXT = KI + 2
903: DO 200 J = KI + 2, N
904: IF( J.LT.JNXT )
905: $ GO TO 200
906: J1 = J
907: J2 = J
908: JNXT = J + 1
909: IF( J.LT.N ) THEN
910: IF( T( J+1, J ).NE.ZERO ) THEN
911: J2 = J + 1
912: JNXT = J + 2
913: END IF
914: END IF
915: *
916: IF( J1.EQ.J2 ) THEN
917: *
918: * 1-by-1 diagonal block
919: *
920: * Scale if necessary to avoid overflow when
921: * forming the right-hand side elements.
922: *
923: IF( WORK( J ).GT.VCRIT ) THEN
924: REC = ONE / VMAX
925: CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
926: CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
927: VMAX = ONE
928: VCRIT = BIGNUM
929: END IF
930: *
931: WORK( J+N ) = WORK( J+N ) -
932: $ DDOT( J-KI-2, T( KI+2, J ), 1,
933: $ WORK( KI+2+N ), 1 )
934: WORK( J+N2 ) = WORK( J+N2 ) -
935: $ DDOT( J-KI-2, T( KI+2, J ), 1,
936: $ WORK( KI+2+N2 ), 1 )
937: *
938: * Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
939: *
940: CALL DLALN2( .FALSE., 1, 2, SMIN, ONE, T( J, J ),
941: $ LDT, ONE, ONE, WORK( J+N ), N, WR,
942: $ -WI, X, 2, SCALE, XNORM, IERR )
943: *
944: * Scale if necessary
945: *
946: IF( SCALE.NE.ONE ) THEN
947: CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
948: CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
949: END IF
950: WORK( J+N ) = X( 1, 1 )
951: WORK( J+N2 ) = X( 1, 2 )
952: VMAX = MAX( ABS( WORK( J+N ) ),
953: $ ABS( WORK( J+N2 ) ), VMAX )
954: VCRIT = BIGNUM / VMAX
955: *
956: ELSE
957: *
958: * 2-by-2 diagonal block
959: *
960: * Scale if necessary to avoid overflow when forming
961: * the right-hand side elements.
962: *
963: BETA = MAX( WORK( J ), WORK( J+1 ) )
964: IF( BETA.GT.VCRIT ) THEN
965: REC = ONE / VMAX
966: CALL DSCAL( N-KI+1, REC, WORK( KI+N ), 1 )
967: CALL DSCAL( N-KI+1, REC, WORK( KI+N2 ), 1 )
968: VMAX = ONE
969: VCRIT = BIGNUM
970: END IF
971: *
972: WORK( J+N ) = WORK( J+N ) -
973: $ DDOT( J-KI-2, T( KI+2, J ), 1,
974: $ WORK( KI+2+N ), 1 )
975: *
976: WORK( J+N2 ) = WORK( J+N2 ) -
977: $ DDOT( J-KI-2, T( KI+2, J ), 1,
978: $ WORK( KI+2+N2 ), 1 )
979: *
980: WORK( J+1+N ) = WORK( J+1+N ) -
981: $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
982: $ WORK( KI+2+N ), 1 )
983: *
984: WORK( J+1+N2 ) = WORK( J+1+N2 ) -
985: $ DDOT( J-KI-2, T( KI+2, J+1 ), 1,
986: $ WORK( KI+2+N2 ), 1 )
987: *
988: * Solve 2-by-2 complex linear equation
989: * ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B
990: * ([T(j+1,j) T(j+1,j+1)] )
991: *
992: CALL DLALN2( .TRUE., 2, 2, SMIN, ONE, T( J, J ),
993: $ LDT, ONE, ONE, WORK( J+N ), N, WR,
994: $ -WI, X, 2, SCALE, XNORM, IERR )
995: *
996: * Scale if necessary
997: *
998: IF( SCALE.NE.ONE ) THEN
999: CALL DSCAL( N-KI+1, SCALE, WORK( KI+N ), 1 )
1000: CALL DSCAL( N-KI+1, SCALE, WORK( KI+N2 ), 1 )
1001: END IF
1002: WORK( J+N ) = X( 1, 1 )
1003: WORK( J+N2 ) = X( 1, 2 )
1004: WORK( J+1+N ) = X( 2, 1 )
1005: WORK( J+1+N2 ) = X( 2, 2 )
1006: VMAX = MAX( ABS( X( 1, 1 ) ), ABS( X( 1, 2 ) ),
1007: $ ABS( X( 2, 1 ) ), ABS( X( 2, 2 ) ), VMAX )
1008: VCRIT = BIGNUM / VMAX
1009: *
1010: END IF
1011: 200 CONTINUE
1012: *
1013: * Copy the vector x or Q*x to VL and normalize.
1014: *
1015: IF( .NOT.OVER ) THEN
1016: CALL DCOPY( N-KI+1, WORK( KI+N ), 1, VL( KI, IS ), 1 )
1017: CALL DCOPY( N-KI+1, WORK( KI+N2 ), 1, VL( KI, IS+1 ),
1018: $ 1 )
1019: *
1020: EMAX = ZERO
1021: DO 220 K = KI, N
1022: EMAX = MAX( EMAX, ABS( VL( K, IS ) )+
1023: $ ABS( VL( K, IS+1 ) ) )
1024: 220 CONTINUE
1025: REMAX = ONE / EMAX
1026: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS ), 1 )
1027: CALL DSCAL( N-KI+1, REMAX, VL( KI, IS+1 ), 1 )
1028: *
1029: DO 230 K = 1, KI - 1
1030: VL( K, IS ) = ZERO
1031: VL( K, IS+1 ) = ZERO
1032: 230 CONTINUE
1033: ELSE
1034: IF( KI.LT.N-1 ) THEN
1035: CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
1036: $ LDVL, WORK( KI+2+N ), 1, WORK( KI+N ),
1037: $ VL( 1, KI ), 1 )
1038: CALL DGEMV( 'N', N, N-KI-1, ONE, VL( 1, KI+2 ),
1039: $ LDVL, WORK( KI+2+N2 ), 1,
1040: $ WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
1041: ELSE
1042: CALL DSCAL( N, WORK( KI+N ), VL( 1, KI ), 1 )
1043: CALL DSCAL( N, WORK( KI+1+N2 ), VL( 1, KI+1 ), 1 )
1044: END IF
1045: *
1046: EMAX = ZERO
1047: DO 240 K = 1, N
1048: EMAX = MAX( EMAX, ABS( VL( K, KI ) )+
1049: $ ABS( VL( K, KI+1 ) ) )
1050: 240 CONTINUE
1051: REMAX = ONE / EMAX
1052: CALL DSCAL( N, REMAX, VL( 1, KI ), 1 )
1053: CALL DSCAL( N, REMAX, VL( 1, KI+1 ), 1 )
1054: *
1055: END IF
1056: *
1057: END IF
1058: *
1059: IS = IS + 1
1060: IF( IP.NE.0 )
1061: $ IS = IS + 1
1062: 250 CONTINUE
1063: IF( IP.EQ.-1 )
1064: $ IP = 0
1065: IF( IP.EQ.1 )
1066: $ IP = -1
1067: *
1068: 260 CONTINUE
1069: *
1070: END IF
1071: *
1072: RETURN
1073: *
1074: * End of DTREVC
1075: *
1076: 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