LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
dtrevc3.f
1 *> \brief \b DTREVC3
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DTREVC3 + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtrevc3.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtrevc3.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtrevc3.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DTREVC3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
22 * VR, LDVR, MM, M, WORK, LWORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER HOWMNY, SIDE
26 * INTEGER INFO, LDT, LDVL, LDVR, LWORK, 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 *> DTREVC3 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 the vector 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 *>
60 *> This uses a Level 3 BLAS version of the back transformation.
61 *> \endverbatim
62 *
63 * Arguments:
64 * ==========
65 *
66 *> \param[in] SIDE
67 *> \verbatim
68 *> SIDE is CHARACTER*1
69 *> = 'R': compute right eigenvectors only;
70 *> = 'L': compute left eigenvectors only;
71 *> = 'B': compute both right and left eigenvectors.
72 *> \endverbatim
73 *>
74 *> \param[in] HOWMNY
75 *> \verbatim
76 *> HOWMNY is CHARACTER*1
77 *> = 'A': compute all right and/or left eigenvectors;
78 *> = 'B': compute all right and/or left eigenvectors,
79 *> backtransformed by the matrices in VR and/or VL;
80 *> = 'S': compute selected right and/or left eigenvectors,
81 *> as indicated by the logical array SELECT.
82 *> \endverbatim
83 *>
84 *> \param[in,out] SELECT
85 *> \verbatim
86 *> SELECT is LOGICAL array, dimension (N)
87 *> If HOWMNY = 'S', SELECT specifies the eigenvectors to be
88 *> computed.
89 *> If w(j) is a real eigenvalue, the corresponding real
90 *> eigenvector is computed if SELECT(j) is .TRUE..
91 *> If w(j) and w(j+1) are the real and imaginary parts of a
92 *> complex eigenvalue, the corresponding complex eigenvector is
93 *> computed if either SELECT(j) or SELECT(j+1) is .TRUE., and
94 *> on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is set to
95 *> .FALSE..
96 *> Not referenced if HOWMNY = 'A' or 'B'.
97 *> \endverbatim
98 *>
99 *> \param[in] N
100 *> \verbatim
101 *> N is INTEGER
102 *> The order of the matrix T. N >= 0.
103 *> \endverbatim
104 *>
105 *> \param[in] T
106 *> \verbatim
107 *> T is DOUBLE PRECISION array, dimension (LDT,N)
108 *> The upper quasi-triangular matrix T in Schur canonical form.
109 *> \endverbatim
110 *>
111 *> \param[in] LDT
112 *> \verbatim
113 *> LDT is INTEGER
114 *> The leading dimension of the array T. LDT >= max(1,N).
115 *> \endverbatim
116 *>
117 *> \param[in,out] VL
118 *> \verbatim
119 *> VL is DOUBLE PRECISION array, dimension (LDVL,MM)
120 *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
121 *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
122 *> of Schur vectors returned by DHSEQR).
123 *> On exit, if SIDE = 'L' or 'B', VL contains:
124 *> if HOWMNY = 'A', the matrix Y of left eigenvectors of T;
125 *> if HOWMNY = 'B', the matrix Q*Y;
126 *> if HOWMNY = 'S', the left eigenvectors of T specified by
127 *> SELECT, stored consecutively in the columns
128 *> of VL, in the same order as their
129 *> eigenvalues.
130 *> A complex eigenvector corresponding to a complex eigenvalue
131 *> is stored in two consecutive columns, the first holding the
132 *> real part, and the second the imaginary part.
133 *> Not referenced if SIDE = 'R'.
134 *> \endverbatim
135 *>
136 *> \param[in] LDVL
137 *> \verbatim
138 *> LDVL is INTEGER
139 *> The leading dimension of the array VL.
140 *> LDVL >= 1, and if SIDE = 'L' or 'B', LDVL >= N.
141 *> \endverbatim
142 *>
143 *> \param[in,out] VR
144 *> \verbatim
145 *> VR is DOUBLE PRECISION array, dimension (LDVR,MM)
146 *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
147 *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
148 *> of Schur vectors returned by DHSEQR).
149 *> On exit, if SIDE = 'R' or 'B', VR contains:
150 *> if HOWMNY = 'A', the matrix X of right eigenvectors of T;
151 *> if HOWMNY = 'B', the matrix Q*X;
152 *> if HOWMNY = 'S', the right eigenvectors of T specified by
153 *> SELECT, stored consecutively in the columns
154 *> of VR, in the same order as their
155 *> eigenvalues.
156 *> A complex eigenvector corresponding to a complex eigenvalue
157 *> is stored in two consecutive columns, the first holding the
158 *> real part and the second the imaginary part.
159 *> Not referenced if SIDE = 'L'.
160 *> \endverbatim
161 *>
162 *> \param[in] LDVR
163 *> \verbatim
164 *> LDVR is INTEGER
165 *> The leading dimension of the array VR.
166 *> LDVR >= 1, and if SIDE = 'R' or 'B', LDVR >= N.
167 *> \endverbatim
168 *>
169 *> \param[in] MM
170 *> \verbatim
171 *> MM is INTEGER
172 *> The number of columns in the arrays VL and/or VR. MM >= M.
173 *> \endverbatim
174 *>
175 *> \param[out] M
176 *> \verbatim
177 *> M is INTEGER
178 *> The number of columns in the arrays VL and/or VR actually
179 *> used to store the eigenvectors.
180 *> If HOWMNY = 'A' or 'B', M is set to N.
181 *> Each selected real eigenvector occupies one column and each
182 *> selected complex eigenvector occupies two columns.
183 *> \endverbatim
184 *>
185 *> \param[out] WORK
186 *> \verbatim
187 *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
188 *> \endverbatim
189 *>
190 *> \param[in] LWORK
191 *> \verbatim
192 *> LWORK is INTEGER
193 *> The dimension of array WORK. LWORK >= max(1,3*N).
194 *> For optimum performance, LWORK >= N + 2*N*NB, where NB is
195 *> the optimal blocksize.
196 *>
197 *> If LWORK = -1, then a workspace query is assumed; the routine
198 *> only calculates the optimal size of the WORK array, returns
199 *> this value as the first entry of the WORK array, and no error
200 *> message related to LWORK is issued by XERBLA.
201 *> \endverbatim
202 *>
203 *> \param[out] INFO
204 *> \verbatim
205 *> INFO is INTEGER
206 *> = 0: successful exit
207 *> < 0: if INFO = -i, the i-th argument had an illegal value
208 *> \endverbatim
209 *
210 * Authors:
211 * ========
212 *
213 *> \author Univ. of Tennessee
214 *> \author Univ. of California Berkeley
215 *> \author Univ. of Colorado Denver
216 *> \author NAG Ltd.
217 *
218 *> \ingroup trevc3
219 *
220 *> \par Further Details:
221 * =====================
222 *>
223 *> \verbatim
224 *>
225 *> The algorithm used in this program is basically backward (forward)
226 *> substitution, with scaling to make the the code robust against
227 *> possible overflow.
228 *>
229 *> Each eigenvector is normalized so that the element of largest
230 *> magnitude has magnitude 1; here the magnitude of a complex number
231 *> (x,y) is taken to be |x| + |y|.
232 *> \endverbatim
233 *>
234 * =====================================================================
235  SUBROUTINE dtrevc3( SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL,
236  $ VR, LDVR, MM, M, WORK, LWORK, INFO )
237  IMPLICIT NONE
238 *
239 * -- LAPACK computational routine --
240 * -- LAPACK is a software package provided by Univ. of Tennessee, --
241 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
242 *
243 * .. Scalar Arguments ..
244  CHARACTER HOWMNY, SIDE
245  INTEGER INFO, LDT, LDVL, LDVR, LWORK, M, MM, N
246 * ..
247 * .. Array Arguments ..
248  LOGICAL SELECT( * )
249  DOUBLE PRECISION T( ldt, * ), VL( ldvl, * ), VR( ldvr, * ),
250  $ work( * )
251 * ..
252 *
253 * =====================================================================
254 *
255 * .. Parameters ..
256  DOUBLE PRECISION ZERO, ONE
257  parameter( zero = 0.0d+0, one = 1.0d+0 )
258  INTEGER NBMIN, NBMAX
259  parameter( nbmin = 8, nbmax = 128 )
260 * ..
261 * .. Local Scalars ..
262  LOGICAL ALLV, BOTHV, LEFTV, LQUERY, OVER, PAIR,
263  $ rightv, somev
264  INTEGER I, IERR, II, IP, IS, J, J1, J2, JNXT, K, KI,
265  $ iv, maxwrk, nb, ki2
266  DOUBLE PRECISION BETA, BIGNUM, EMAX, OVFL, REC, REMAX, SCALE,
267  $ smin, smlnum, ulp, unfl, vcrit, vmax, wi, wr,
268  $ xnorm
269 * ..
270 * .. External Functions ..
271  LOGICAL LSAME
272  INTEGER IDAMAX, ILAENV
273  DOUBLE PRECISION DDOT, DLAMCH
274  EXTERNAL lsame, idamax, ilaenv, ddot, dlamch
275 * ..
276 * .. External Subroutines ..
277  EXTERNAL daxpy, dcopy, dgemv, dlaln2, dscal, xerbla,
278  $ dgemm, dlaset, dlabad, dlacpy
279 * ..
280 * .. Intrinsic Functions ..
281  INTRINSIC abs, max, sqrt
282 * ..
283 * .. Local Arrays ..
284  DOUBLE PRECISION X( 2, 2 )
285  INTEGER ISCOMPLEX( nbmax )
286 * ..
287 * .. Executable Statements ..
288 *
289 * Decode and test the input parameters
290 *
291  bothv = lsame( side, 'B' )
292  rightv = lsame( side, 'R' ) .OR. bothv
293  leftv = lsame( side, 'L' ) .OR. bothv
294 *
295  allv = lsame( howmny, 'A' )
296  over = lsame( howmny, 'B' )
297  somev = lsame( howmny, 'S' )
298 *
299  info = 0
300  nb = ilaenv( 1, 'DTREVC', side // howmny, n, -1, -1, -1 )
301  maxwrk = n + 2*n*nb
302  work(1) = maxwrk
303  lquery = ( lwork.EQ.-1 )
304  IF( .NOT.rightv .AND. .NOT.leftv ) THEN
305  info = -1
306  ELSE IF( .NOT.allv .AND. .NOT.over .AND. .NOT.somev ) THEN
307  info = -2
308  ELSE IF( n.LT.0 ) THEN
309  info = -4
310  ELSE IF( ldt.LT.max( 1, n ) ) THEN
311  info = -6
312  ELSE IF( ldvl.LT.1 .OR. ( leftv .AND. ldvl.LT.n ) ) THEN
313  info = -8
314  ELSE IF( ldvr.LT.1 .OR. ( rightv .AND. ldvr.LT.n ) ) THEN
315  info = -10
316  ELSE IF( lwork.LT.max( 1, 3*n ) .AND. .NOT.lquery ) THEN
317  info = -14
318  ELSE
319 *
320 * Set M to the number of columns required to store the selected
321 * eigenvectors, standardize the array SELECT if necessary, and
322 * test MM.
323 *
324  IF( somev ) THEN
325  m = 0
326  pair = .false.
327  DO 10 j = 1, n
328  IF( pair ) THEN
329  pair = .false.
330  SELECT( j ) = .false.
331  ELSE
332  IF( j.LT.n ) THEN
333  IF( t( j+1, j ).EQ.zero ) THEN
334  IF( SELECT( j ) )
335  $ m = m + 1
336  ELSE
337  pair = .true.
338  IF( SELECT( j ) .OR. SELECT( j+1 ) ) THEN
339  SELECT( j ) = .true.
340  m = m + 2
341  END IF
342  END IF
343  ELSE
344  IF( SELECT( n ) )
345  $ m = m + 1
346  END IF
347  END IF
348  10 CONTINUE
349  ELSE
350  m = n
351  END IF
352 *
353  IF( mm.LT.m ) THEN
354  info = -11
355  END IF
356  END IF
357  IF( info.NE.0 ) THEN
358  CALL xerbla( 'DTREVC3', -info )
359  RETURN
360  ELSE IF( lquery ) THEN
361  RETURN
362  END IF
363 *
364 * Quick return if possible.
365 *
366  IF( n.EQ.0 )
367  $ RETURN
368 *
369 * Use blocked version of back-transformation if sufficient workspace.
370 * Zero-out the workspace to avoid potential NaN propagation.
371 *
372  IF( over .AND. lwork .GE. n + 2*n*nbmin ) THEN
373  nb = (lwork - n) / (2*n)
374  nb = min( nb, nbmax )
375  CALL dlaset( 'F', n, 1+2*nb, zero, zero, work, n )
376  ELSE
377  nb = 1
378  END IF
379 *
380 * Set the constants to control overflow.
381 *
382  unfl = dlamch( 'Safe minimum' )
383  ovfl = one / unfl
384  CALL dlabad( unfl, ovfl )
385  ulp = dlamch( 'Precision' )
386  smlnum = unfl*( n / ulp )
387  bignum = ( one-ulp ) / smlnum
388 *
389 * Compute 1-norm of each column of strictly upper triangular
390 * part of T to control overflow in triangular solver.
391 *
392  work( 1 ) = zero
393  DO 30 j = 2, n
394  work( j ) = zero
395  DO 20 i = 1, j - 1
396  work( j ) = work( j ) + abs( t( i, j ) )
397  20 CONTINUE
398  30 CONTINUE
399 *
400 * Index IP is used to specify the real or complex eigenvalue:
401 * IP = 0, real eigenvalue,
402 * 1, first of conjugate complex pair: (wr,wi)
403 * -1, second of conjugate complex pair: (wr,wi)
404 * ISCOMPLEX array stores IP for each column in current block.
405 *
406  IF( rightv ) THEN
407 *
408 * ============================================================
409 * Compute right eigenvectors.
410 *
411 * IV is index of column in current block.
412 * For complex right vector, uses IV-1 for real part and IV for complex part.
413 * Non-blocked version always uses IV=2;
414 * blocked version starts with IV=NB, goes down to 1 or 2.
415 * (Note the "0-th" column is used for 1-norms computed above.)
416  iv = 2
417  IF( nb.GT.2 ) THEN
418  iv = nb
419  END IF
420 
421  ip = 0
422  is = m
423  DO 140 ki = n, 1, -1
424  IF( ip.EQ.-1 ) THEN
425 * previous iteration (ki+1) was second of conjugate pair,
426 * so this ki is first of conjugate pair; skip to end of loop
427  ip = 1
428  GO TO 140
429  ELSE IF( ki.EQ.1 ) THEN
430 * last column, so this ki must be real eigenvalue
431  ip = 0
432  ELSE IF( t( ki, ki-1 ).EQ.zero ) THEN
433 * zero on sub-diagonal, so this ki is real eigenvalue
434  ip = 0
435  ELSE
436 * non-zero on sub-diagonal, so this ki is second of conjugate pair
437  ip = -1
438  END IF
439 
440  IF( somev ) THEN
441  IF( ip.EQ.0 ) THEN
442  IF( .NOT.SELECT( ki ) )
443  $ GO TO 140
444  ELSE
445  IF( .NOT.SELECT( ki-1 ) )
446  $ GO TO 140
447  END IF
448  END IF
449 *
450 * Compute the KI-th eigenvalue (WR,WI).
451 *
452  wr = t( ki, ki )
453  wi = zero
454  IF( ip.NE.0 )
455  $ wi = sqrt( abs( t( ki, ki-1 ) ) )*
456  $ sqrt( abs( t( ki-1, ki ) ) )
457  smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
458 *
459  IF( ip.EQ.0 ) THEN
460 *
461 * --------------------------------------------------------
462 * Real right eigenvector
463 *
464  work( ki + iv*n ) = one
465 *
466 * Form right-hand side.
467 *
468  DO 50 k = 1, ki - 1
469  work( k + iv*n ) = -t( k, ki )
470  50 CONTINUE
471 *
472 * Solve upper quasi-triangular system:
473 * [ T(1:KI-1,1:KI-1) - WR ]*X = SCALE*WORK.
474 *
475  jnxt = ki - 1
476  DO 60 j = ki - 1, 1, -1
477  IF( j.GT.jnxt )
478  $ GO TO 60
479  j1 = j
480  j2 = j
481  jnxt = j - 1
482  IF( j.GT.1 ) THEN
483  IF( t( j, j-1 ).NE.zero ) THEN
484  j1 = j - 1
485  jnxt = j - 2
486  END IF
487  END IF
488 *
489  IF( j1.EQ.j2 ) THEN
490 *
491 * 1-by-1 diagonal block
492 *
493  CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
494  $ ldt, one, one, work( j+iv*n ), n, wr,
495  $ zero, x, 2, scale, xnorm, ierr )
496 *
497 * Scale X(1,1) to avoid overflow when updating
498 * the right-hand side.
499 *
500  IF( xnorm.GT.one ) THEN
501  IF( work( j ).GT.bignum / xnorm ) THEN
502  x( 1, 1 ) = x( 1, 1 ) / xnorm
503  scale = scale / xnorm
504  END IF
505  END IF
506 *
507 * Scale if necessary
508 *
509  IF( scale.NE.one )
510  $ CALL dscal( ki, scale, work( 1+iv*n ), 1 )
511  work( j+iv*n ) = x( 1, 1 )
512 *
513 * Update right-hand side
514 *
515  CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
516  $ work( 1+iv*n ), 1 )
517 *
518  ELSE
519 *
520 * 2-by-2 diagonal block
521 *
522  CALL dlaln2( .false., 2, 1, smin, one,
523  $ t( j-1, j-1 ), ldt, one, one,
524  $ work( j-1+iv*n ), n, wr, zero, x, 2,
525  $ scale, xnorm, ierr )
526 *
527 * Scale X(1,1) and X(2,1) to avoid overflow when
528 * updating the right-hand side.
529 *
530  IF( xnorm.GT.one ) THEN
531  beta = max( work( j-1 ), work( j ) )
532  IF( beta.GT.bignum / xnorm ) THEN
533  x( 1, 1 ) = x( 1, 1 ) / xnorm
534  x( 2, 1 ) = x( 2, 1 ) / xnorm
535  scale = scale / xnorm
536  END IF
537  END IF
538 *
539 * Scale if necessary
540 *
541  IF( scale.NE.one )
542  $ CALL dscal( ki, scale, work( 1+iv*n ), 1 )
543  work( j-1+iv*n ) = x( 1, 1 )
544  work( j +iv*n ) = x( 2, 1 )
545 *
546 * Update right-hand side
547 *
548  CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
549  $ work( 1+iv*n ), 1 )
550  CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
551  $ work( 1+iv*n ), 1 )
552  END IF
553  60 CONTINUE
554 *
555 * Copy the vector x or Q*x to VR and normalize.
556 *
557  IF( .NOT.over ) THEN
558 * ------------------------------
559 * no back-transform: copy x to VR and normalize.
560  CALL dcopy( ki, work( 1 + iv*n ), 1, vr( 1, is ), 1 )
561 *
562  ii = idamax( ki, vr( 1, is ), 1 )
563  remax = one / abs( vr( ii, is ) )
564  CALL dscal( ki, remax, vr( 1, is ), 1 )
565 *
566  DO 70 k = ki + 1, n
567  vr( k, is ) = zero
568  70 CONTINUE
569 *
570  ELSE IF( nb.EQ.1 ) THEN
571 * ------------------------------
572 * version 1: back-transform each vector with GEMV, Q*x.
573  IF( ki.GT.1 )
574  $ CALL dgemv( 'N', n, ki-1, one, vr, ldvr,
575  $ work( 1 + iv*n ), 1, work( ki + iv*n ),
576  $ vr( 1, ki ), 1 )
577 *
578  ii = idamax( n, vr( 1, ki ), 1 )
579  remax = one / abs( vr( ii, ki ) )
580  CALL dscal( n, remax, vr( 1, ki ), 1 )
581 *
582  ELSE
583 * ------------------------------
584 * version 2: back-transform block of vectors with GEMM
585 * zero out below vector
586  DO k = ki + 1, n
587  work( k + iv*n ) = zero
588  END DO
589  iscomplex( iv ) = ip
590 * back-transform and normalization is done below
591  END IF
592  ELSE
593 *
594 * --------------------------------------------------------
595 * Complex right eigenvector.
596 *
597 * Initial solve
598 * [ ( T(KI-1,KI-1) T(KI-1,KI) ) - (WR + I*WI) ]*X = 0.
599 * [ ( T(KI, KI-1) T(KI, KI) ) ]
600 *
601  IF( abs( t( ki-1, ki ) ).GE.abs( t( ki, ki-1 ) ) ) THEN
602  work( ki-1 + (iv-1)*n ) = one
603  work( ki + (iv )*n ) = wi / t( ki-1, ki )
604  ELSE
605  work( ki-1 + (iv-1)*n ) = -wi / t( ki, ki-1 )
606  work( ki + (iv )*n ) = one
607  END IF
608  work( ki + (iv-1)*n ) = zero
609  work( ki-1 + (iv )*n ) = zero
610 *
611 * Form right-hand side.
612 *
613  DO 80 k = 1, ki - 2
614  work( k+(iv-1)*n ) = -work( ki-1+(iv-1)*n )*t(k,ki-1)
615  work( k+(iv )*n ) = -work( ki +(iv )*n )*t(k,ki )
616  80 CONTINUE
617 *
618 * Solve upper quasi-triangular system:
619 * [ T(1:KI-2,1:KI-2) - (WR+i*WI) ]*X = SCALE*(WORK+i*WORK2)
620 *
621  jnxt = ki - 2
622  DO 90 j = ki - 2, 1, -1
623  IF( j.GT.jnxt )
624  $ GO TO 90
625  j1 = j
626  j2 = j
627  jnxt = j - 1
628  IF( j.GT.1 ) THEN
629  IF( t( j, j-1 ).NE.zero ) THEN
630  j1 = j - 1
631  jnxt = j - 2
632  END IF
633  END IF
634 *
635  IF( j1.EQ.j2 ) THEN
636 *
637 * 1-by-1 diagonal block
638 *
639  CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
640  $ ldt, one, one, work( j+(iv-1)*n ), n,
641  $ wr, wi, x, 2, scale, xnorm, ierr )
642 *
643 * Scale X(1,1) and X(1,2) to avoid overflow when
644 * updating the right-hand side.
645 *
646  IF( xnorm.GT.one ) THEN
647  IF( work( j ).GT.bignum / xnorm ) THEN
648  x( 1, 1 ) = x( 1, 1 ) / xnorm
649  x( 1, 2 ) = x( 1, 2 ) / xnorm
650  scale = scale / xnorm
651  END IF
652  END IF
653 *
654 * Scale if necessary
655 *
656  IF( scale.NE.one ) THEN
657  CALL dscal( ki, scale, work( 1+(iv-1)*n ), 1 )
658  CALL dscal( ki, scale, work( 1+(iv )*n ), 1 )
659  END IF
660  work( j+(iv-1)*n ) = x( 1, 1 )
661  work( j+(iv )*n ) = x( 1, 2 )
662 *
663 * Update the right-hand side
664 *
665  CALL daxpy( j-1, -x( 1, 1 ), t( 1, j ), 1,
666  $ work( 1+(iv-1)*n ), 1 )
667  CALL daxpy( j-1, -x( 1, 2 ), t( 1, j ), 1,
668  $ work( 1+(iv )*n ), 1 )
669 *
670  ELSE
671 *
672 * 2-by-2 diagonal block
673 *
674  CALL dlaln2( .false., 2, 2, smin, one,
675  $ t( j-1, j-1 ), ldt, one, one,
676  $ work( j-1+(iv-1)*n ), n, wr, wi, x, 2,
677  $ scale, xnorm, ierr )
678 *
679 * Scale X to avoid overflow when updating
680 * the right-hand side.
681 *
682  IF( xnorm.GT.one ) THEN
683  beta = max( work( j-1 ), work( j ) )
684  IF( beta.GT.bignum / xnorm ) THEN
685  rec = one / xnorm
686  x( 1, 1 ) = x( 1, 1 )*rec
687  x( 1, 2 ) = x( 1, 2 )*rec
688  x( 2, 1 ) = x( 2, 1 )*rec
689  x( 2, 2 ) = x( 2, 2 )*rec
690  scale = scale*rec
691  END IF
692  END IF
693 *
694 * Scale if necessary
695 *
696  IF( scale.NE.one ) THEN
697  CALL dscal( ki, scale, work( 1+(iv-1)*n ), 1 )
698  CALL dscal( ki, scale, work( 1+(iv )*n ), 1 )
699  END IF
700  work( j-1+(iv-1)*n ) = x( 1, 1 )
701  work( j +(iv-1)*n ) = x( 2, 1 )
702  work( j-1+(iv )*n ) = x( 1, 2 )
703  work( j +(iv )*n ) = x( 2, 2 )
704 *
705 * Update the right-hand side
706 *
707  CALL daxpy( j-2, -x( 1, 1 ), t( 1, j-1 ), 1,
708  $ work( 1+(iv-1)*n ), 1 )
709  CALL daxpy( j-2, -x( 2, 1 ), t( 1, j ), 1,
710  $ work( 1+(iv-1)*n ), 1 )
711  CALL daxpy( j-2, -x( 1, 2 ), t( 1, j-1 ), 1,
712  $ work( 1+(iv )*n ), 1 )
713  CALL daxpy( j-2, -x( 2, 2 ), t( 1, j ), 1,
714  $ work( 1+(iv )*n ), 1 )
715  END IF
716  90 CONTINUE
717 *
718 * Copy the vector x or Q*x to VR and normalize.
719 *
720  IF( .NOT.over ) THEN
721 * ------------------------------
722 * no back-transform: copy x to VR and normalize.
723  CALL dcopy( ki, work( 1+(iv-1)*n ), 1, vr(1,is-1), 1 )
724  CALL dcopy( ki, work( 1+(iv )*n ), 1, vr(1,is ), 1 )
725 *
726  emax = zero
727  DO 100 k = 1, ki
728  emax = max( emax, abs( vr( k, is-1 ) )+
729  $ abs( vr( k, is ) ) )
730  100 CONTINUE
731  remax = one / emax
732  CALL dscal( ki, remax, vr( 1, is-1 ), 1 )
733  CALL dscal( ki, remax, vr( 1, is ), 1 )
734 *
735  DO 110 k = ki + 1, n
736  vr( k, is-1 ) = zero
737  vr( k, is ) = zero
738  110 CONTINUE
739 *
740  ELSE IF( nb.EQ.1 ) THEN
741 * ------------------------------
742 * version 1: back-transform each vector with GEMV, Q*x.
743  IF( ki.GT.2 ) THEN
744  CALL dgemv( 'N', n, ki-2, one, vr, ldvr,
745  $ work( 1 + (iv-1)*n ), 1,
746  $ work( ki-1 + (iv-1)*n ), vr(1,ki-1), 1)
747  CALL dgemv( 'N', n, ki-2, one, vr, ldvr,
748  $ work( 1 + (iv)*n ), 1,
749  $ work( ki + (iv)*n ), vr( 1, ki ), 1 )
750  ELSE
751  CALL dscal( n, work(ki-1+(iv-1)*n), vr(1,ki-1), 1)
752  CALL dscal( n, work(ki +(iv )*n), vr(1,ki ), 1)
753  END IF
754 *
755  emax = zero
756  DO 120 k = 1, n
757  emax = max( emax, abs( vr( k, ki-1 ) )+
758  $ abs( vr( k, ki ) ) )
759  120 CONTINUE
760  remax = one / emax
761  CALL dscal( n, remax, vr( 1, ki-1 ), 1 )
762  CALL dscal( n, remax, vr( 1, ki ), 1 )
763 *
764  ELSE
765 * ------------------------------
766 * version 2: back-transform block of vectors with GEMM
767 * zero out below vector
768  DO k = ki + 1, n
769  work( k + (iv-1)*n ) = zero
770  work( k + (iv )*n ) = zero
771  END DO
772  iscomplex( iv-1 ) = -ip
773  iscomplex( iv ) = ip
774  iv = iv - 1
775 * back-transform and normalization is done below
776  END IF
777  END IF
778 
779  IF( nb.GT.1 ) THEN
780 * --------------------------------------------------------
781 * Blocked version of back-transform
782 * For complex case, KI2 includes both vectors (KI-1 and KI)
783  IF( ip.EQ.0 ) THEN
784  ki2 = ki
785  ELSE
786  ki2 = ki - 1
787  END IF
788 
789 * Columns IV:NB of work are valid vectors.
790 * When the number of vectors stored reaches NB-1 or NB,
791 * or if this was last vector, do the GEMM
792  IF( (iv.LE.2) .OR. (ki2.EQ.1) ) THEN
793  CALL dgemm( 'N', 'N', n, nb-iv+1, ki2+nb-iv, one,
794  $ vr, ldvr,
795  $ work( 1 + (iv)*n ), n,
796  $ zero,
797  $ work( 1 + (nb+iv)*n ), n )
798 * normalize vectors
799  DO k = iv, nb
800  IF( iscomplex(k).EQ.0 ) THEN
801 * real eigenvector
802  ii = idamax( n, work( 1 + (nb+k)*n ), 1 )
803  remax = one / abs( work( ii + (nb+k)*n ) )
804  ELSE IF( iscomplex(k).EQ.1 ) THEN
805 * first eigenvector of conjugate pair
806  emax = zero
807  DO ii = 1, n
808  emax = max( emax,
809  $ abs( work( ii + (nb+k )*n ) )+
810  $ abs( work( ii + (nb+k+1)*n ) ) )
811  END DO
812  remax = one / emax
813 * else if ISCOMPLEX(K).EQ.-1
814 * second eigenvector of conjugate pair
815 * reuse same REMAX as previous K
816  END IF
817  CALL dscal( n, remax, work( 1 + (nb+k)*n ), 1 )
818  END DO
819  CALL dlacpy( 'F', n, nb-iv+1,
820  $ work( 1 + (nb+iv)*n ), n,
821  $ vr( 1, ki2 ), ldvr )
822  iv = nb
823  ELSE
824  iv = iv - 1
825  END IF
826  END IF ! blocked back-transform
827 *
828  is = is - 1
829  IF( ip.NE.0 )
830  $ is = is - 1
831  140 CONTINUE
832  END IF
833 
834  IF( leftv ) THEN
835 *
836 * ============================================================
837 * Compute left eigenvectors.
838 *
839 * IV is index of column in current block.
840 * For complex left vector, uses IV for real part and IV+1 for complex part.
841 * Non-blocked version always uses IV=1;
842 * blocked version starts with IV=1, goes up to NB-1 or NB.
843 * (Note the "0-th" column is used for 1-norms computed above.)
844  iv = 1
845  ip = 0
846  is = 1
847  DO 260 ki = 1, n
848  IF( ip.EQ.1 ) THEN
849 * previous iteration (ki-1) was first of conjugate pair,
850 * so this ki is second of conjugate pair; skip to end of loop
851  ip = -1
852  GO TO 260
853  ELSE IF( ki.EQ.n ) THEN
854 * last column, so this ki must be real eigenvalue
855  ip = 0
856  ELSE IF( t( ki+1, ki ).EQ.zero ) THEN
857 * zero on sub-diagonal, so this ki is real eigenvalue
858  ip = 0
859  ELSE
860 * non-zero on sub-diagonal, so this ki is first of conjugate pair
861  ip = 1
862  END IF
863 *
864  IF( somev ) THEN
865  IF( .NOT.SELECT( ki ) )
866  $ GO TO 260
867  END IF
868 *
869 * Compute the KI-th eigenvalue (WR,WI).
870 *
871  wr = t( ki, ki )
872  wi = zero
873  IF( ip.NE.0 )
874  $ wi = sqrt( abs( t( ki, ki+1 ) ) )*
875  $ sqrt( abs( t( ki+1, ki ) ) )
876  smin = max( ulp*( abs( wr )+abs( wi ) ), smlnum )
877 *
878  IF( ip.EQ.0 ) THEN
879 *
880 * --------------------------------------------------------
881 * Real left eigenvector
882 *
883  work( ki + iv*n ) = one
884 *
885 * Form right-hand side.
886 *
887  DO 160 k = ki + 1, n
888  work( k + iv*n ) = -t( ki, k )
889  160 CONTINUE
890 *
891 * Solve transposed quasi-triangular system:
892 * [ T(KI+1:N,KI+1:N) - WR ]**T * X = SCALE*WORK
893 *
894  vmax = one
895  vcrit = bignum
896 *
897  jnxt = ki + 1
898  DO 170 j = ki + 1, n
899  IF( j.LT.jnxt )
900  $ GO TO 170
901  j1 = j
902  j2 = j
903  jnxt = j + 1
904  IF( j.LT.n ) THEN
905  IF( t( j+1, j ).NE.zero ) THEN
906  j2 = j + 1
907  jnxt = j + 2
908  END IF
909  END IF
910 *
911  IF( j1.EQ.j2 ) THEN
912 *
913 * 1-by-1 diagonal block
914 *
915 * Scale if necessary to avoid overflow when forming
916 * the right-hand side.
917 *
918  IF( work( j ).GT.vcrit ) THEN
919  rec = one / vmax
920  CALL dscal( n-ki+1, rec, work( ki+iv*n ), 1 )
921  vmax = one
922  vcrit = bignum
923  END IF
924 *
925  work( j+iv*n ) = work( j+iv*n ) -
926  $ ddot( j-ki-1, t( ki+1, j ), 1,
927  $ work( ki+1+iv*n ), 1 )
928 *
929 * Solve [ T(J,J) - WR ]**T * X = WORK
930 *
931  CALL dlaln2( .false., 1, 1, smin, one, t( j, j ),
932  $ ldt, one, one, work( j+iv*n ), n, wr,
933  $ zero, x, 2, scale, xnorm, ierr )
934 *
935 * Scale if necessary
936 *
937  IF( scale.NE.one )
938  $ CALL dscal( n-ki+1, scale, work( ki+iv*n ), 1 )
939  work( j+iv*n ) = x( 1, 1 )
940  vmax = max( abs( work( j+iv*n ) ), vmax )
941  vcrit = bignum / vmax
942 *
943  ELSE
944 *
945 * 2-by-2 diagonal block
946 *
947 * Scale if necessary to avoid overflow when forming
948 * the right-hand side.
949 *
950  beta = max( work( j ), work( j+1 ) )
951  IF( beta.GT.vcrit ) THEN
952  rec = one / vmax
953  CALL dscal( n-ki+1, rec, work( ki+iv*n ), 1 )
954  vmax = one
955  vcrit = bignum
956  END IF
957 *
958  work( j+iv*n ) = work( j+iv*n ) -
959  $ ddot( j-ki-1, t( ki+1, j ), 1,
960  $ work( ki+1+iv*n ), 1 )
961 *
962  work( j+1+iv*n ) = work( j+1+iv*n ) -
963  $ ddot( j-ki-1, t( ki+1, j+1 ), 1,
964  $ work( ki+1+iv*n ), 1 )
965 *
966 * Solve
967 * [ T(J,J)-WR T(J,J+1) ]**T * X = SCALE*( WORK1 )
968 * [ T(J+1,J) T(J+1,J+1)-WR ] ( WORK2 )
969 *
970  CALL dlaln2( .true., 2, 1, smin, one, t( j, j ),
971  $ ldt, one, one, work( j+iv*n ), n, wr,
972  $ zero, x, 2, scale, xnorm, ierr )
973 *
974 * Scale if necessary
975 *
976  IF( scale.NE.one )
977  $ CALL dscal( n-ki+1, scale, work( ki+iv*n ), 1 )
978  work( j +iv*n ) = x( 1, 1 )
979  work( j+1+iv*n ) = x( 2, 1 )
980 *
981  vmax = max( abs( work( j +iv*n ) ),
982  $ abs( work( j+1+iv*n ) ), vmax )
983  vcrit = bignum / vmax
984 *
985  END IF
986  170 CONTINUE
987 *
988 * Copy the vector x or Q*x to VL and normalize.
989 *
990  IF( .NOT.over ) THEN
991 * ------------------------------
992 * no back-transform: copy x to VL and normalize.
993  CALL dcopy( n-ki+1, work( ki + iv*n ), 1,
994  $ vl( ki, is ), 1 )
995 *
996  ii = idamax( n-ki+1, vl( ki, is ), 1 ) + ki - 1
997  remax = one / abs( vl( ii, is ) )
998  CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
999 *
1000  DO 180 k = 1, ki - 1
1001  vl( k, is ) = zero
1002  180 CONTINUE
1003 *
1004  ELSE IF( nb.EQ.1 ) THEN
1005 * ------------------------------
1006 * version 1: back-transform each vector with GEMV, Q*x.
1007  IF( ki.LT.n )
1008  $ CALL dgemv( 'N', n, n-ki, one,
1009  $ vl( 1, ki+1 ), ldvl,
1010  $ work( ki+1 + iv*n ), 1,
1011  $ work( ki + iv*n ), vl( 1, ki ), 1 )
1012 *
1013  ii = idamax( n, vl( 1, ki ), 1 )
1014  remax = one / abs( vl( ii, ki ) )
1015  CALL dscal( n, remax, vl( 1, ki ), 1 )
1016 *
1017  ELSE
1018 * ------------------------------
1019 * version 2: back-transform block of vectors with GEMM
1020 * zero out above vector
1021 * could go from KI-NV+1 to KI-1
1022  DO k = 1, ki - 1
1023  work( k + iv*n ) = zero
1024  END DO
1025  iscomplex( iv ) = ip
1026 * back-transform and normalization is done below
1027  END IF
1028  ELSE
1029 *
1030 * --------------------------------------------------------
1031 * Complex left eigenvector.
1032 *
1033 * Initial solve:
1034 * [ ( T(KI,KI) T(KI,KI+1) )**T - (WR - I* WI) ]*X = 0.
1035 * [ ( T(KI+1,KI) T(KI+1,KI+1) ) ]
1036 *
1037  IF( abs( t( ki, ki+1 ) ).GE.abs( t( ki+1, ki ) ) ) THEN
1038  work( ki + (iv )*n ) = wi / t( ki, ki+1 )
1039  work( ki+1 + (iv+1)*n ) = one
1040  ELSE
1041  work( ki + (iv )*n ) = one
1042  work( ki+1 + (iv+1)*n ) = -wi / t( ki+1, ki )
1043  END IF
1044  work( ki+1 + (iv )*n ) = zero
1045  work( ki + (iv+1)*n ) = zero
1046 *
1047 * Form right-hand side.
1048 *
1049  DO 190 k = ki + 2, n
1050  work( k+(iv )*n ) = -work( ki +(iv )*n )*t(ki, k)
1051  work( k+(iv+1)*n ) = -work( ki+1+(iv+1)*n )*t(ki+1,k)
1052  190 CONTINUE
1053 *
1054 * Solve transposed quasi-triangular system:
1055 * [ T(KI+2:N,KI+2:N)**T - (WR-i*WI) ]*X = WORK1+i*WORK2
1056 *
1057  vmax = one
1058  vcrit = bignum
1059 *
1060  jnxt = ki + 2
1061  DO 200 j = ki + 2, n
1062  IF( j.LT.jnxt )
1063  $ GO TO 200
1064  j1 = j
1065  j2 = j
1066  jnxt = j + 1
1067  IF( j.LT.n ) THEN
1068  IF( t( j+1, j ).NE.zero ) THEN
1069  j2 = j + 1
1070  jnxt = j + 2
1071  END IF
1072  END IF
1073 *
1074  IF( j1.EQ.j2 ) THEN
1075 *
1076 * 1-by-1 diagonal block
1077 *
1078 * Scale if necessary to avoid overflow when
1079 * forming the right-hand side elements.
1080 *
1081  IF( work( j ).GT.vcrit ) THEN
1082  rec = one / vmax
1083  CALL dscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1084  CALL dscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1085  vmax = one
1086  vcrit = bignum
1087  END IF
1088 *
1089  work( j+(iv )*n ) = work( j+(iv)*n ) -
1090  $ ddot( j-ki-2, t( ki+2, j ), 1,
1091  $ work( ki+2+(iv)*n ), 1 )
1092  work( j+(iv+1)*n ) = work( j+(iv+1)*n ) -
1093  $ ddot( j-ki-2, t( ki+2, j ), 1,
1094  $ work( ki+2+(iv+1)*n ), 1 )
1095 *
1096 * Solve [ T(J,J)-(WR-i*WI) ]*(X11+i*X12)= WK+I*WK2
1097 *
1098  CALL dlaln2( .false., 1, 2, smin, one, t( j, j ),
1099  $ ldt, one, one, work( j+iv*n ), n, wr,
1100  $ -wi, x, 2, scale, xnorm, ierr )
1101 *
1102 * Scale if necessary
1103 *
1104  IF( scale.NE.one ) THEN
1105  CALL dscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1106  CALL dscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1107  END IF
1108  work( j+(iv )*n ) = x( 1, 1 )
1109  work( j+(iv+1)*n ) = x( 1, 2 )
1110  vmax = max( abs( work( j+(iv )*n ) ),
1111  $ abs( work( j+(iv+1)*n ) ), vmax )
1112  vcrit = bignum / vmax
1113 *
1114  ELSE
1115 *
1116 * 2-by-2 diagonal block
1117 *
1118 * Scale if necessary to avoid overflow when forming
1119 * the right-hand side elements.
1120 *
1121  beta = max( work( j ), work( j+1 ) )
1122  IF( beta.GT.vcrit ) THEN
1123  rec = one / vmax
1124  CALL dscal( n-ki+1, rec, work(ki+(iv )*n), 1 )
1125  CALL dscal( n-ki+1, rec, work(ki+(iv+1)*n), 1 )
1126  vmax = one
1127  vcrit = bignum
1128  END IF
1129 *
1130  work( j +(iv )*n ) = work( j+(iv)*n ) -
1131  $ ddot( j-ki-2, t( ki+2, j ), 1,
1132  $ work( ki+2+(iv)*n ), 1 )
1133 *
1134  work( j +(iv+1)*n ) = work( j+(iv+1)*n ) -
1135  $ ddot( j-ki-2, t( ki+2, j ), 1,
1136  $ work( ki+2+(iv+1)*n ), 1 )
1137 *
1138  work( j+1+(iv )*n ) = work( j+1+(iv)*n ) -
1139  $ ddot( j-ki-2, t( ki+2, j+1 ), 1,
1140  $ work( ki+2+(iv)*n ), 1 )
1141 *
1142  work( j+1+(iv+1)*n ) = work( j+1+(iv+1)*n ) -
1143  $ ddot( j-ki-2, t( ki+2, j+1 ), 1,
1144  $ work( ki+2+(iv+1)*n ), 1 )
1145 *
1146 * Solve 2-by-2 complex linear equation
1147 * [ (T(j,j) T(j,j+1) )**T - (wr-i*wi)*I ]*X = SCALE*B
1148 * [ (T(j+1,j) T(j+1,j+1)) ]
1149 *
1150  CALL dlaln2( .true., 2, 2, smin, one, t( j, j ),
1151  $ ldt, one, one, work( j+iv*n ), n, wr,
1152  $ -wi, x, 2, scale, xnorm, ierr )
1153 *
1154 * Scale if necessary
1155 *
1156  IF( scale.NE.one ) THEN
1157  CALL dscal( n-ki+1, scale, work(ki+(iv )*n), 1)
1158  CALL dscal( n-ki+1, scale, work(ki+(iv+1)*n), 1)
1159  END IF
1160  work( j +(iv )*n ) = x( 1, 1 )
1161  work( j +(iv+1)*n ) = x( 1, 2 )
1162  work( j+1+(iv )*n ) = x( 2, 1 )
1163  work( j+1+(iv+1)*n ) = x( 2, 2 )
1164  vmax = max( abs( x( 1, 1 ) ), abs( x( 1, 2 ) ),
1165  $ abs( x( 2, 1 ) ), abs( x( 2, 2 ) ),
1166  $ vmax )
1167  vcrit = bignum / vmax
1168 *
1169  END IF
1170  200 CONTINUE
1171 *
1172 * Copy the vector x or Q*x to VL and normalize.
1173 *
1174  IF( .NOT.over ) THEN
1175 * ------------------------------
1176 * no back-transform: copy x to VL and normalize.
1177  CALL dcopy( n-ki+1, work( ki + (iv )*n ), 1,
1178  $ vl( ki, is ), 1 )
1179  CALL dcopy( n-ki+1, work( ki + (iv+1)*n ), 1,
1180  $ vl( ki, is+1 ), 1 )
1181 *
1182  emax = zero
1183  DO 220 k = ki, n
1184  emax = max( emax, abs( vl( k, is ) )+
1185  $ abs( vl( k, is+1 ) ) )
1186  220 CONTINUE
1187  remax = one / emax
1188  CALL dscal( n-ki+1, remax, vl( ki, is ), 1 )
1189  CALL dscal( n-ki+1, remax, vl( ki, is+1 ), 1 )
1190 *
1191  DO 230 k = 1, ki - 1
1192  vl( k, is ) = zero
1193  vl( k, is+1 ) = zero
1194  230 CONTINUE
1195 *
1196  ELSE IF( nb.EQ.1 ) THEN
1197 * ------------------------------
1198 * version 1: back-transform each vector with GEMV, Q*x.
1199  IF( ki.LT.n-1 ) THEN
1200  CALL dgemv( 'N', n, n-ki-1, one,
1201  $ vl( 1, ki+2 ), ldvl,
1202  $ work( ki+2 + (iv)*n ), 1,
1203  $ work( ki + (iv)*n ),
1204  $ vl( 1, ki ), 1 )
1205  CALL dgemv( 'N', n, n-ki-1, one,
1206  $ vl( 1, ki+2 ), ldvl,
1207  $ work( ki+2 + (iv+1)*n ), 1,
1208  $ work( ki+1 + (iv+1)*n ),
1209  $ vl( 1, ki+1 ), 1 )
1210  ELSE
1211  CALL dscal( n, work(ki+ (iv )*n), vl(1, ki ), 1)
1212  CALL dscal( n, work(ki+1+(iv+1)*n), vl(1, ki+1), 1)
1213  END IF
1214 *
1215  emax = zero
1216  DO 240 k = 1, n
1217  emax = max( emax, abs( vl( k, ki ) )+
1218  $ abs( vl( k, ki+1 ) ) )
1219  240 CONTINUE
1220  remax = one / emax
1221  CALL dscal( n, remax, vl( 1, ki ), 1 )
1222  CALL dscal( n, remax, vl( 1, ki+1 ), 1 )
1223 *
1224  ELSE
1225 * ------------------------------
1226 * version 2: back-transform block of vectors with GEMM
1227 * zero out above vector
1228 * could go from KI-NV+1 to KI-1
1229  DO k = 1, ki - 1
1230  work( k + (iv )*n ) = zero
1231  work( k + (iv+1)*n ) = zero
1232  END DO
1233  iscomplex( iv ) = ip
1234  iscomplex( iv+1 ) = -ip
1235  iv = iv + 1
1236 * back-transform and normalization is done below
1237  END IF
1238  END IF
1239 
1240  IF( nb.GT.1 ) THEN
1241 * --------------------------------------------------------
1242 * Blocked version of back-transform
1243 * For complex case, KI2 includes both vectors (KI and KI+1)
1244  IF( ip.EQ.0 ) THEN
1245  ki2 = ki
1246  ELSE
1247  ki2 = ki + 1
1248  END IF
1249 
1250 * Columns 1:IV of work are valid vectors.
1251 * When the number of vectors stored reaches NB-1 or NB,
1252 * or if this was last vector, do the GEMM
1253  IF( (iv.GE.nb-1) .OR. (ki2.EQ.n) ) THEN
1254  CALL dgemm( 'N', 'N', n, iv, n-ki2+iv, one,
1255  $ vl( 1, ki2-iv+1 ), ldvl,
1256  $ work( ki2-iv+1 + (1)*n ), n,
1257  $ zero,
1258  $ work( 1 + (nb+1)*n ), n )
1259 * normalize vectors
1260  DO k = 1, iv
1261  IF( iscomplex(k).EQ.0) THEN
1262 * real eigenvector
1263  ii = idamax( n, work( 1 + (nb+k)*n ), 1 )
1264  remax = one / abs( work( ii + (nb+k)*n ) )
1265  ELSE IF( iscomplex(k).EQ.1) THEN
1266 * first eigenvector of conjugate pair
1267  emax = zero
1268  DO ii = 1, n
1269  emax = max( emax,
1270  $ abs( work( ii + (nb+k )*n ) )+
1271  $ abs( work( ii + (nb+k+1)*n ) ) )
1272  END DO
1273  remax = one / emax
1274 * else if ISCOMPLEX(K).EQ.-1
1275 * second eigenvector of conjugate pair
1276 * reuse same REMAX as previous K
1277  END IF
1278  CALL dscal( n, remax, work( 1 + (nb+k)*n ), 1 )
1279  END DO
1280  CALL dlacpy( 'F', n, iv,
1281  $ work( 1 + (nb+1)*n ), n,
1282  $ vl( 1, ki2-iv+1 ), ldvl )
1283  iv = 1
1284  ELSE
1285  iv = iv + 1
1286  END IF
1287  END IF ! blocked back-transform
1288 *
1289  is = is + 1
1290  IF( ip.NE.0 )
1291  $ is = is + 1
1292  260 CONTINUE
1293  END IF
1294 *
1295  RETURN
1296 *
1297 * End of DTREVC3
1298 *
1299  END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
dtrevc3
subroutine dtrevc3(SIDE, HOWMNY, SELECT, N, T, LDT, VL, LDVL, VR, LDVR, MM, M, WORK, LWORK, INFO)
DTREVC3
Definition: dtrevc3.f:237
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
dgemm
subroutine dgemm(TRANSA, TRANSB, M, N, K, ALPHA, A, LDA, B, LDB, BETA, C, LDC)
DGEMM
Definition: dgemm.f:187
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
daxpy
subroutine daxpy(N, DA, DX, INCX, DY, INCY)
DAXPY
Definition: daxpy.f:89
dscal
subroutine dscal(N, DA, DX, INCX)
DSCAL
Definition: dscal.f:79
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
dlaset
subroutine dlaset(UPLO, M, N, ALPHA, BETA, A, LDA)
DLASET initializes the off-diagonal elements and the diagonal elements of a matrix to given values...
Definition: dlaset.f:110
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
dcopy
subroutine dcopy(N, DX, INCX, DY, INCY)
DCOPY
Definition: dcopy.f:82