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