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