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