LAPACK 3.12.1
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*> Download DTREVC + dependencies
9*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc.f">
10*> [TGZ]</a>
11*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc.f">
12*> [ZIP]</a>
13*> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc.f">
14*> [TXT]</a>
15*
16* Definition:
17* ===========
18*
19* SUBROUTINE DTREVC( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR,
20* LDVR, MM, M, WORK, INFO )
21*
22* .. Scalar Arguments ..
23* CHARACTER HOWMNY, SIDE
24* INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
25* ..
26* .. Array Arguments ..
27* LOGICAL SELECT( * )
28* DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
29* $ WORK( * )
30* ..
31*
32*
33*> \par Purpose:
34* =============
35*>
36*> \verbatim
37*>
38*> DTREVC computes some or all of the right and/or left eigenvectors of
39*> a real upper quasi-triangular matrix T.
40*> Matrices of this type are produced by the Schur factorization of
41*> a real general matrix: A = Q*T*Q**T, as computed by DHSEQR.
42*>
43*> The right eigenvector x and the left eigenvector y of T corresponding
44*> to an eigenvalue w are defined by:
45*>
46*> T*x = w*x, (y**H)*T = w*(y**H)
47*>
48*> where y**H denotes the conjugate transpose of y.
49*> The eigenvalues are not input to this routine, but are read directly
50*> from the diagonal blocks of T.
51*>
52*> This routine returns the matrices X and/or Y of right and left
53*> eigenvectors of T, or the products Q*X and/or Q*Y, where Q is an
54*> input matrix. If Q is the orthogonal factor that reduces a matrix
55*> A to Schur form T, then Q*X and Q*Y are the matrices of right and
56*> left eigenvectors of A.
57*> \endverbatim
58*
59* Arguments:
60* ==========
61*
62*> \param[in] SIDE
63*> \verbatim
64*> SIDE is CHARACTER*1
65*> = 'R': compute right eigenvectors only;
66*> = 'L': compute left eigenvectors only;
67*> = 'B': compute both right and left eigenvectors.
68*> \endverbatim
69*>
70*> \param[in] HOWMNY
71*> \verbatim
72*> HOWMNY is CHARACTER*1
73*> = 'A': compute all right and/or left eigenvectors;
74*> = 'B': compute all right and/or left eigenvectors,
75*> backtransformed by the matrices in VR and/or VL;
76*> = 'S': compute selected right and/or left eigenvectors,
77*> as indicated by the logical array SELECT.
78*> \endverbatim
79*>
80*> \param[in,out] SELECT
81*> \verbatim
82*> SELECT is LOGICAL array, dimension (N)
83*> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
84*> computed.
85*> If w(j) is a real eigenvalue, the corresponding real
86*> eigenvector is computed if SELECT(j) is .TRUE..
87*> If w(j) and w(j+1) are the real and imaginary parts of a
88*> complex eigenvalue, the corresponding complex eigenvector is
89*> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
90*> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
91*> .FALSE..
92*> Not referenced if HOWMNY = 'A' or 'B'.
93*> \endverbatim
94*>
95*> \param[in] N
96*> \verbatim
97*> N is INTEGER
98*> The order of the matrix T. N >= 0.
99*> \endverbatim
100*>
101*> \param[in] T
102*> \verbatim
103*> T is DOUBLE PRECISION array, dimension (LDT,N)
104*> The upper quasi-triangular matrix T in Schur canonical form.
105*> \endverbatim
106*>
107*> \param[in] LDT
108*> \verbatim
109*> LDT is INTEGER
110*> The leading dimension of the array T. LDT >= max(1,N).
111*> \endverbatim
112*>
113*> \param[in,out] VL
114*> \verbatim
115*> VL is DOUBLE PRECISION array, dimension (LDVL,MM)
116*> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
117*> contain an N-by-N matrix Q (usually the orthogonal matrix Q
118*> of Schur vectors returned by DHSEQR).
119*> On exit, if SIDE = 'L' or 'B', VL contains:
120*> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
121*> if HOWMNY = 'B', the matrix Q*Y;
122*> if HOWMNY = 'S', the left eigenvectors of T specified by
123*> SELECT, stored consecutively in the columns
124*> of VL, in the same order as their
125*> eigenvalues.
126*> A complex eigenvector corresponding to a complex eigenvalue
127*> is stored in two consecutive columns, the first holding the
128*> real part, and the second the imaginary part.
129*> Not referenced if SIDE = 'R'.
130*> \endverbatim
131*>
132*> \param[in] LDVL
133*> \verbatim
134*> LDVL is INTEGER
135*> The leading dimension of the array VL. LDVL >= 1, and if
136*> SIDE = 'L' or 'B', LDVL >= N.
137*> \endverbatim
138*>
139*> \param[in,out] VR
140*> \verbatim
141*> VR is DOUBLE PRECISION array, dimension (LDVR,MM)
142*> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
143*> contain an N-by-N matrix Q (usually the orthogonal matrix Q
144*> of Schur vectors returned by DHSEQR).
145*> On exit, if SIDE = 'R' or 'B', VR contains:
146*> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
147*> if HOWMNY = 'B', the matrix Q*X;
148*> if HOWMNY = 'S', the right eigenvectors of T specified by
149*> SELECT, stored consecutively in the columns
150*> of VR, in the same order as their
151*> eigenvalues.
152*> A complex eigenvector corresponding to a complex eigenvalue
153*> is stored in two consecutive columns, the first holding the
154*> real part and the second the imaginary part.
155*> Not referenced if SIDE = 'L'.
156*> \endverbatim
157*>
158*> \param[in] LDVR
159*> \verbatim
160*> LDVR is INTEGER
161*> The leading dimension of the array VR. LDVR >= 1, and if
162*> SIDE = 'R' or 'B', LDVR >= N.
163*> \endverbatim
164*>
165*> \param[in] MM
166*> \verbatim
167*> MM is INTEGER
168*> The number of columns in the arrays VL and/or VR. MM >= M.
169*> \endverbatim
170*>
171*> \param[out] M
172*> \verbatim
173*> M is INTEGER
174*> The number of columns in the arrays VL and/or VR actually
175*> used to store the eigenvectors.
176*> If HOWMNY = 'A' or 'B', M is set to N.
177*> Each selected real eigenvector occupies one column and each
178*> selected complex eigenvector occupies two columns.
179*> \endverbatim
180*>
181*> \param[out] WORK
182*> \verbatim
183*> WORK is DOUBLE PRECISION array, dimension (3*N)
184*> \endverbatim
185*>
186*> \param[out] INFO
187*> \verbatim
188*> INFO is INTEGER
189*> = 0: successful exit
190*> < 0: if INFO = -i, the i-th argument had an illegal value
191*> \endverbatim
192*
193* Authors:
194* ========
195*
196*> \author Univ. of Tennessee
197*> \author Univ. of California Berkeley
198*> \author Univ. of Colorado Denver
199*> \author NAG Ltd.
200*
201*> \ingroup trevc
202*
203*> \par Further Details:
204* =====================
205*>
206*> \verbatim
207*>
208*> The algorithm used in this program is basically backward (forward)
209*> substitution, with scaling to make the the code robust against
210*> possible overflow.
211*>
212*> Each eigenvector is normalized so that the element of largest
213*> magnitude has magnitude 1; here the magnitude of a complex number
214*> (x,y) is taken to be |x| + |y|.
215*> \endverbatim
216*>
217* =====================================================================
218 SUBROUTINE dtrevc( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
219 $ VR,
220 $ LDVR, MM, M, WORK, INFO )
221*
222* -- LAPACK computational routine --
223* -- LAPACK is a software package provided by Univ. of Tennessee, --
224* -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
225*
226* .. Scalar Arguments ..
227 CHARACTER HOWMNY, SIDE
228 INTEGER INFO, LDT, LDVL, LDVR, M, MM, N
229* ..
230* .. Array Arguments ..
231 LOGICAL SELECT( * )
232 DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ),
233 $ work( * )
234* ..
235*
236* =====================================================================
237*
238* .. Parameters ..
239 DOUBLE PRECISION ZERO, ONE
240 PARAMETER ( ZERO = 0.0d+0, one = 1.0d+0 )
241* ..
242* .. Local Scalars ..
243 LOGICAL ALLV, BOTHV, LEFTV, OVER, PAIR, RIGHTV, SOMEV
244 INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI, N2
245 DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
246 $ smin, smlnum, ulp, unfl, vcrit, vmax, wi, wr,
247 $ xnorm
248* ..
249* .. External Functions ..
250 LOGICAL LSAME
251 INTEGER IDAMAX
252 DOUBLE PRECISION DDOT, DLAMCH
253 EXTERNAL lsame, idamax, ddot, dlamch
254* ..
255* .. External Subroutines ..
256 EXTERNAL daxpy, dcopy, dgemv, dlaln2, dscal,
257 $ 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,
435 $ j ),
436 $ ldt, one, one, work( j+n ), n, wr,
437 $ zero, x, 2, scale, xnorm, ierr )
438*
439* Scale X(1,1) to avoid overflow when updating
440* the right-hand side.
441*
442 IF( xnorm.GT.one ) THEN
443 IF( work( j ).GT.bignum / xnorm ) THEN
444 x( 1, 1 ) = x( 1, 1 ) / xnorm
445 scale = scale / xnorm
446 END IF
447 END IF
448*
449* Scale if necessary
450*
451 IF( scale.NE.one )
452 $ CALL dscal( ki, scale, work( 1+n ), 1 )
453 work( j+n ) = x( 1, 1 )
454*
455* Update right-hand side
456*
457 CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
458 $ work( 1+n ), 1 )
459*
460 ELSE
461*
462* 2-by-2 diagonal block
463*
464 CALL dlaln2( .false., 2, 1, smin, one,
465 $ t( j-1, j-1 ), ldt, one, one,
466 $ work( j-1+n ), n, wr, zero, x, 2,
467 $ scale, xnorm, ierr )
468*
469* Scale X(1,1) and X(2,1) to avoid overflow when
470* updating the right-hand side.
471*
472 IF( xnorm.GT.one ) THEN
473 beta = max( work( j-1 ), work( j ) )
474 IF( beta.GT.bignum / xnorm ) THEN
475 x( 1, 1 ) = x( 1, 1 ) / xnorm
476 x( 2, 1 ) = x( 2, 1 ) / xnorm
477 scale = scale / xnorm
478 END IF
479 END IF
480*
481* Scale if necessary
482*
483 IF( scale.NE.one )
484 $ CALL dscal( ki, scale, work( 1+n ), 1 )
485 work( j-1+n ) = x( 1, 1 )
486 work( j+n ) = x( 2, 1 )
487*
488* Update right-hand side
489*
490 CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
491 $ work( 1+n ), 1 )
492 CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
493 $ work( 1+n ), 1 )
494 END IF
495 60 CONTINUE
496*
497* Copy the vector x or Q*x to VR and normalize.
498*
499 IF( .NOT.over ) THEN
500 CALL dcopy( ki, work( 1+n ), 1, vr( 1, is ), 1 )
501*
502 ii = idamax( ki, vr( 1, is ), 1 )
503 remax = one / abs( vr( ii, is ) )
504 CALL dscal( ki, remax, vr( 1, is ), 1 )
505*
506 DO 70 k = ki + 1, n
507 vr( k, is ) = zero
508 70 CONTINUE
509 ELSE
510 IF( ki.GT.1 )
511 $ CALL dgemv( 'N', n, ki-1, one, vr, ldvr,
512 $ work( 1+n ), 1, work( ki+n ),
513 $ vr( 1, ki ), 1 )
514*
515 ii = idamax( n, vr( 1, ki ), 1 )
516 remax = one / abs( vr( ii, ki ) )
517 CALL dscal( n, remax, vr( 1, ki ), 1 )
518 END IF
519*
520 ELSE
521*
522* Complex right eigenvector.
523*
524* Initial solve
525* [ (T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I* WI)]*X = 0.
526* [ (T(KI,KI-1) T(KI,KI) ) ]
527*
528 IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) ) THEN
529 work( ki-1+n ) = one
530 work( ki+n2 ) = wi / t( ki-1, ki )
531 ELSE
532 work( ki-1+n ) = -wi / t( ki, ki-1 )
533 work( ki+n2 ) = one
534 END IF
535 work( ki+n ) = zero
536 work( ki-1+n2 ) = zero
537*
538* Form right-hand side
539*
540 DO 80 k = 1, ki - 2
541 work( k+n ) = -work( ki-1+n )*t( k, ki-1 )
542 work( k+n2 ) = -work( ki+n2 )*t( k, ki )
543 80 CONTINUE
544*
545* Solve upper quasi-triangular system:
546* (T(1:KI-2,1:KI-2) - (WR+i*WI))*X = SCALE*(WORK+i*WORK2)
547*
548 jnxt = ki - 2
549 DO 90 j = ki - 2, 1, -1
550 IF( j.GT.jnxt )
551 $ GO TO 90
552 j1 = j
553 j2 = j
554 jnxt = j - 1
555 IF( j.GT.1 ) THEN
556 IF( t( j, j-1 ).NE.zero ) THEN
557 j1 = j - 1
558 jnxt = j - 2
559 END IF
560 END IF
561*
562 IF( j1.EQ.j2 ) THEN
563*
564* 1-by-1 diagonal block
565*
566 CALL dlaln2( .false., 1, 2, smin, one, t( j,
567 $ 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 ),
678 $ 1 )
679 CALL dscal( n, work( ki+n2 ), vr( 1, ki ), 1 )
680 END IF
681*
682 emax = zero
683 DO 120 k = 1, n
684 emax = max( emax, abs( vr( k, ki-1 ) )+
685 $ abs( vr( k, ki ) ) )
686 120 CONTINUE
687 remax = one / emax
688 CALL dscal( n, remax, vr( 1, ki-1 ), 1 )
689 CALL dscal( n, remax, vr( 1, ki ), 1 )
690 END IF
691 END IF
692*
693 is = is - 1
694 IF( ip.NE.0 )
695 $ is = is - 1
696 130 CONTINUE
697 IF( ip.EQ.1 )
698 $ ip = 0
699 IF( ip.EQ.-1 )
700 $ ip = 1
701 140 CONTINUE
702 END IF
703*
704 IF( leftv ) THEN
705*
706* Compute left eigenvectors.
707*
708 ip = 0
709 is = 1
710 DO 260 ki = 1, n
711*
712 IF( ip.EQ.-1 )
713 $ GO TO 250
714 IF( ki.EQ.n )
715 $ GO TO 150
716 IF( t( ki+1, ki ).EQ.zero )
717 $ GO TO 150
718 ip = 1
719*
720 150 CONTINUE
721 IF( somev ) THEN
722 IF( .NOT.SELECT( ki ) )
723 $ GO TO 250
724 END IF
725*
726* Compute the KI-th eigenvalue (WR,WI).
727*
728 wr = t( ki, ki )
729 wi = zero
730 IF( ip.NE.0 )
731 $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
732 $ sqrt( abs( t( ki+1, ki ) ) )
733 smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
734*
735 IF( ip.EQ.0 ) THEN
736*
737* Real left eigenvector.
738*
739 work( ki+n ) = one
740*
741* Form right-hand side
742*
743 DO 160 k = ki + 1, n
744 work( k+n ) = -t( ki, k )
745 160 CONTINUE
746*
747* Solve the quasi-triangular system:
748* (T(KI+1:N,KI+1:N) - WR)**T*X = SCALE*WORK
749*
750 vmax = one
751 vcrit = bignum
752*
753 jnxt = ki + 1
754 DO 170 j = ki + 1, n
755 IF( j.LT.jnxt )
756 $ GO TO 170
757 j1 = j
758 j2 = j
759 jnxt = j + 1
760 IF( j.LT.n ) THEN
761 IF( t( j+1, j ).NE.zero ) THEN
762 j2 = j + 1
763 jnxt = j + 2
764 END IF
765 END IF
766*
767 IF( j1.EQ.j2 ) THEN
768*
769* 1-by-1 diagonal block
770*
771* Scale if necessary to avoid overflow when forming
772* the right-hand side.
773*
774 IF( work( j ).GT.vcrit ) THEN
775 rec = one / vmax
776 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
777 vmax = one
778 vcrit = bignum
779 END IF
780*
781 work( j+n ) = work( j+n ) -
782 $ ddot( j-ki-1, t( ki+1, j ), 1,
783 $ work( ki+1+n ), 1 )
784*
785* Solve (T(J,J)-WR)**T*X = WORK
786*
787 CALL dlaln2( .false., 1, 1, smin, one, t( j,
788 $ 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 ),
849 $ 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 ),
863 $ ldvl,
864 $ work( ki+1+n ), 1, work( ki+n ),
865 $ vl( 1, ki ), 1 )
866*
867 ii = idamax( n, vl( 1, ki ), 1 )
868 remax = one / abs( vl( ii, ki ) )
869 CALL dscal( n, remax, vl( 1, ki ), 1 )
870*
871 END IF
872*
873 ELSE
874*
875* Complex left eigenvector.
876*
877* Initial solve:
878* ((T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI))*X = 0.
879* ((T(KI+1,KI) T(KI+1,KI+1)) )
880*
881 IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) ) THEN
882 work( ki+n ) = wi / t( ki, ki+1 )
883 work( ki+1+n2 ) = one
884 ELSE
885 work( ki+n ) = one
886 work( ki+1+n2 ) = -wi / t( ki+1, ki )
887 END IF
888 work( ki+1+n ) = zero
889 work( ki+n2 ) = zero
890*
891* Form right-hand side
892*
893 DO 190 k = ki + 2, n
894 work( k+n ) = -work( ki+n )*t( ki, k )
895 work( k+n2 ) = -work( ki+1+n2 )*t( ki+1, k )
896 190 CONTINUE
897*
898* Solve complex quasi-triangular system:
899* ( T(KI+2,N:KI+2,N) - (WR-i*WI) )*X = WORK1+i*WORK2
900*
901 vmax = one
902 vcrit = bignum
903*
904 jnxt = ki + 2
905 DO 200 j = ki + 2, n
906 IF( j.LT.jnxt )
907 $ GO TO 200
908 j1 = j
909 j2 = j
910 jnxt = j + 1
911 IF( j.LT.n ) THEN
912 IF( t( j+1, j ).NE.zero ) THEN
913 j2 = j + 1
914 jnxt = j + 2
915 END IF
916 END IF
917*
918 IF( j1.EQ.j2 ) THEN
919*
920* 1-by-1 diagonal block
921*
922* Scale if necessary to avoid overflow when
923* forming the right-hand side elements.
924*
925 IF( work( j ).GT.vcrit ) THEN
926 rec = one / vmax
927 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
928 CALL dscal( n-ki+1, rec, work( ki+n2 ), 1 )
929 vmax = one
930 vcrit = bignum
931 END IF
932*
933 work( j+n ) = work( j+n ) -
934 $ ddot( j-ki-2, t( ki+2, j ), 1,
935 $ work( ki+2+n ), 1 )
936 work( j+n2 ) = work( j+n2 ) -
937 $ ddot( j-ki-2, t( ki+2, j ), 1,
938 $ work( ki+2+n2 ), 1 )
939*
940* Solve (T(J,J)-(WR-i*WI))*(X11+i*X12)= WK+I*WK2
941*
942 CALL dlaln2( .false., 1, 2, smin, one, t( j,
943 $ j ),
944 $ ldt, one, one, work( j+n ), n, wr,
945 $ -wi, x, 2, scale, xnorm, ierr )
946*
947* Scale if necessary
948*
949 IF( scale.NE.one ) THEN
950 CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
951 CALL dscal( n-ki+1, scale, work( ki+n2 ), 1 )
952 END IF
953 work( j+n ) = x( 1, 1 )
954 work( j+n2 ) = x( 1, 2 )
955 vmax = max( abs( work( j+n ) ),
956 $ abs( work( j+n2 ) ), vmax )
957 vcrit = bignum / vmax
958*
959 ELSE
960*
961* 2-by-2 diagonal block
962*
963* Scale if necessary to avoid overflow when forming
964* the right-hand side elements.
965*
966 beta = max( work( j ), work( j+1 ) )
967 IF( beta.GT.vcrit ) THEN
968 rec = one / vmax
969 CALL dscal( n-ki+1, rec, work( ki+n ), 1 )
970 CALL dscal( n-ki+1, rec, work( ki+n2 ), 1 )
971 vmax = one
972 vcrit = bignum
973 END IF
974*
975 work( j+n ) = work( j+n ) -
976 $ ddot( j-ki-2, t( ki+2, j ), 1,
977 $ work( ki+2+n ), 1 )
978*
979 work( j+n2 ) = work( j+n2 ) -
980 $ ddot( j-ki-2, t( ki+2, j ), 1,
981 $ work( ki+2+n2 ), 1 )
982*
983 work( j+1+n ) = work( j+1+n ) -
984 $ ddot( j-ki-2, t( ki+2, j+1 ), 1,
985 $ work( ki+2+n ), 1 )
986*
987 work( j+1+n2 ) = work( j+1+n2 ) -
988 $ ddot( j-ki-2, t( ki+2, j+1 ),
989 $ 1,
990 $ work( ki+2+n2 ), 1 )
991*
992* Solve 2-by-2 complex linear equation
993* ([T(j,j) T(j,j+1) ]**T-(wr-i*wi)*I)*X = SCALE*B
994* ([T(j+1,j) T(j+1,j+1)] )
995*
996 CALL dlaln2( .true., 2, 2, smin, one, t( j, j ),
997 $ ldt, one, one, work( j+n ), n, wr,
998 $ -wi, x, 2, scale, xnorm, ierr )
999*
1000* Scale if necessary
1001*
1002 IF( scale.NE.one ) THEN
1003 CALL dscal( n-ki+1, scale, work( ki+n ), 1 )
1004 CALL dscal( n-ki+1, scale, work( ki+n2 ), 1 )
1005 END IF
1006 work( j+n ) = x( 1, 1 )
1007 work( j+n2 ) = x( 1, 2 )
1008 work( j+1+n ) = x( 2, 1 )
1009 work( j+1+n2 ) = x( 2, 2 )
1010 vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1011 $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ), vmax )
1012 vcrit = bignum / vmax
1013*
1014 END IF
1015 200 CONTINUE
1016*
1017* Copy the vector x or Q*x to VL and normalize.
1018*
1019 IF( .NOT.over ) THEN
1020 CALL dcopy( n-ki+1, work( ki+n ), 1, vl( ki, is ),
1021 $ 1 )
1022 CALL dcopy( n-ki+1, work( ki+n2 ), 1, vl( ki,
1023 $ is+1 ),
1024 $ 1 )
1025*
1026 emax = zero
1027 DO 220 k = ki, n
1028 emax = max( emax, abs( vl( k, is ) )+
1029 $ abs( vl( k, is+1 ) ) )
1030 220 CONTINUE
1031 remax = one / emax
1032 CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
1033 CALL dscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1034*
1035 DO 230 k = 1, ki - 1
1036 vl( k, is ) = zero
1037 vl( k, is+1 ) = zero
1038 230 CONTINUE
1039 ELSE
1040 IF( ki.LT.n-1 ) THEN
1041 CALL dgemv( 'N', n, n-ki-1, one, vl( 1, ki+2 ),
1042 $ ldvl, work( ki+2+n ), 1, work( ki+n ),
1043 $ vl( 1, ki ), 1 )
1044 CALL dgemv( 'N', n, n-ki-1, one, vl( 1, ki+2 ),
1045 $ ldvl, work( ki+2+n2 ), 1,
1046 $ work( ki+1+n2 ), vl( 1, ki+1 ), 1 )
1047 ELSE
1048 CALL dscal( n, work( ki+n ), vl( 1, ki ), 1 )
1049 CALL dscal( n, work( ki+1+n2 ), vl( 1, ki+1 ),
1050 $ 1 )
1051 END IF
1052*
1053 emax = zero
1054 DO 240 k = 1, n
1055 emax = max( emax, abs( vl( k, ki ) )+
1056 $ abs( vl( k, ki+1 ) ) )
1057 240 CONTINUE
1058 remax = one / emax
1059 CALL dscal( n, remax, vl( 1, ki ), 1 )
1060 CALL dscal( n, remax, vl( 1, ki+1 ), 1 )
1061*
1062 END IF
1063*
1064 END IF
1065*
1066 is = is + 1
1067 IF( ip.NE.0 )
1068 $ is = is + 1
1069 250 CONTINUE
1070 IF( ip.EQ.-1 )
1071 $ ip = 0
1072 IF( ip.EQ.1 )
1073 $ ip = -1
1074*
1075 260 CONTINUE
1076*
1077 END IF
1078*
1079 RETURN
1080*
1081* End of DTREVC
1082*
1083 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:216
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:221