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