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