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