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