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