LAPACK  3.4.2
LAPACK: Linear Algebra PACKage
 All Files Functions Groups
dtgevc.f
Go to the documentation of this file.
1 *> \brief \b DTGEVC
2 *
3 * =========== DOCUMENTATION ===========
4 *
5 * Online html documentation available at
6 * http://www.netlib.org/lapack/explore-html/
7 *
8 *> \htmlonly
9 *> Download DTGEVC + dependencies
10 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dtgevc.f">
11 *> [TGZ]</a>
12 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dtgevc.f">
13 *> [ZIP]</a>
14 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dtgevc.f">
15 *> [TXT]</a>
16 *> \endhtmlonly
17 *
18 * Definition:
19 * ===========
20 *
21 * SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
22 * LDVL, VR, LDVR, MM, M, WORK, INFO )
23 *
24 * .. Scalar Arguments ..
25 * CHARACTER HOWMNY, SIDE
26 * INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
27 * ..
28 * .. Array Arguments ..
29 * LOGICAL SELECT( * )
30 * DOUBLE PRECISION P( LDP, * ), S( LDS, * ), VL( LDVL, * ),
31 * $ VR( LDVR, * ), WORK( * )
32 * ..
33 *
34 *
35 *
36 *> \par Purpose:
37 * =============
38 *>
39 *> \verbatim
40 *>
41 *> DTGEVC computes some or all of the right and/or left eigenvectors of
42 *> a pair of real matrices (S,P), where S is a quasi-triangular matrix
43 *> and P is upper triangular. Matrix pairs of this type are produced by
44 *> the generalized Schur factorization of a matrix pair (A,B):
45 *>
46 *> A = Q*S*Z**T, B = Q*P*Z**T
47 *>
48 *> as computed by DGGHRD + DHGEQZ.
49 *>
50 *> The right eigenvector x and the left eigenvector y of (S,P)
51 *> corresponding to an eigenvalue w are defined by:
52 *>
53 *> S*x = w*P*x, (y**H)*S = w*(y**H)*P,
54 *>
55 *> where y**H denotes the conjugate tranpose of y.
56 *> The eigenvalues are not input to this routine, but are computed
57 *> directly from the diagonal blocks of S and P.
58 *>
59 *> This routine returns the matrices X and/or Y of right and left
60 *> eigenvectors of (S,P), or the products Z*X and/or Q*Y,
61 *> where Z and Q are input matrices.
62 *> If Q and Z are the orthogonal factors from the generalized Schur
63 *> factorization of a matrix pair (A,B), then Z*X and Q*Y
64 *> are the matrices of right and left eigenvectors of (A,B).
65 *>
66 *> \endverbatim
67 *
68 * Arguments:
69 * ==========
70 *
71 *> \param[in] SIDE
72 *> \verbatim
73 *> SIDE is CHARACTER*1
74 *> = 'R': compute right eigenvectors only;
75 *> = 'L': compute left eigenvectors only;
76 *> = 'B': compute both right and left eigenvectors.
77 *> \endverbatim
78 *>
79 *> \param[in] HOWMNY
80 *> \verbatim
81 *> HOWMNY is CHARACTER*1
82 *> = 'A': compute all right and/or left eigenvectors;
83 *> = 'B': compute all right and/or left eigenvectors,
84 *> backtransformed by the matrices in VR and/or VL;
85 *> = 'S': compute selected right and/or left eigenvectors,
86 *> specified by the logical array SELECT.
87 *> \endverbatim
88 *>
89 *> \param[in] SELECT
90 *> \verbatim
91 *> SELECT is LOGICAL array, dimension (N)
92 *> If HOWMNY='S', SELECT specifies the eigenvectors to be
93 *> computed. If w(j) is a real eigenvalue, the corresponding
94 *> real eigenvector is computed if SELECT(j) is .TRUE..
95 *> If w(j) and w(j+1) are the real and imaginary parts of a
96 *> complex eigenvalue, the corresponding complex eigenvector
97 *> is computed if either SELECT(j) or SELECT(j+1) is .TRUE.,
98 *> and on exit SELECT(j) is set to .TRUE. and SELECT(j+1) is
99 *> set to .FALSE..
100 *> Not referenced if HOWMNY = 'A' or 'B'.
101 *> \endverbatim
102 *>
103 *> \param[in] N
104 *> \verbatim
105 *> N is INTEGER
106 *> The order of the matrices S and P. N >= 0.
107 *> \endverbatim
108 *>
109 *> \param[in] S
110 *> \verbatim
111 *> S is DOUBLE PRECISION array, dimension (LDS,N)
112 *> The upper quasi-triangular matrix S from a generalized Schur
113 *> factorization, as computed by DHGEQZ.
114 *> \endverbatim
115 *>
116 *> \param[in] LDS
117 *> \verbatim
118 *> LDS is INTEGER
119 *> The leading dimension of array S. LDS >= max(1,N).
120 *> \endverbatim
121 *>
122 *> \param[in] P
123 *> \verbatim
124 *> P is DOUBLE PRECISION array, dimension (LDP,N)
125 *> The upper triangular matrix P from a generalized Schur
126 *> factorization, as computed by DHGEQZ.
127 *> 2-by-2 diagonal blocks of P corresponding to 2-by-2 blocks
128 *> of S must be in positive diagonal form.
129 *> \endverbatim
130 *>
131 *> \param[in] LDP
132 *> \verbatim
133 *> LDP is INTEGER
134 *> The leading dimension of array P. LDP >= max(1,N).
135 *> \endverbatim
136 *>
137 *> \param[in,out] VL
138 *> \verbatim
139 *> VL is DOUBLE PRECISION array, dimension (LDVL,MM)
140 *> On entry, if SIDE = 'L' or 'B' and HOWMNY = 'B', VL must
141 *> contain an N-by-N matrix Q (usually the orthogonal matrix Q
142 *> of left Schur vectors returned by DHGEQZ).
143 *> On exit, if SIDE = 'L' or 'B', VL contains:
144 *> if HOWMNY = 'A', the matrix Y of left eigenvectors of (S,P);
145 *> if HOWMNY = 'B', the matrix Q*Y;
146 *> if HOWMNY = 'S', the left eigenvectors of (S,P) specified by
147 *> SELECT, stored consecutively in the columns of
148 *> VL, in the same order as their eigenvalues.
149 *>
150 *> A complex eigenvector corresponding to a complex eigenvalue
151 *> is stored in two consecutive columns, the first holding the
152 *> real part, and the second the imaginary part.
153 *>
154 *> Not referenced if SIDE = 'R'.
155 *> \endverbatim
156 *>
157 *> \param[in] LDVL
158 *> \verbatim
159 *> LDVL is INTEGER
160 *> The leading dimension of array VL. LDVL >= 1, and if
161 *> SIDE = 'L' or 'B', LDVL >= N.
162 *> \endverbatim
163 *>
164 *> \param[in,out] VR
165 *> \verbatim
166 *> VR is DOUBLE PRECISION array, dimension (LDVR,MM)
167 *> On entry, if SIDE = 'R' or 'B' and HOWMNY = 'B', VR must
168 *> contain an N-by-N matrix Z (usually the orthogonal matrix Z
169 *> of right Schur vectors returned by DHGEQZ).
170 *>
171 *> On exit, if SIDE = 'R' or 'B', VR contains:
172 *> if HOWMNY = 'A', the matrix X of right eigenvectors of (S,P);
173 *> if HOWMNY = 'B' or 'b', the matrix Z*X;
174 *> if HOWMNY = 'S' or 's', the right eigenvectors of (S,P)
175 *> specified by SELECT, stored consecutively in the
176 *> columns of VR, in the same order as their
177 *> eigenvalues.
178 *>
179 *> A complex eigenvector corresponding to a complex eigenvalue
180 *> is stored in two consecutive columns, the first holding the
181 *> real part and the second the imaginary part.
182 *>
183 *> Not referenced if SIDE = 'L'.
184 *> \endverbatim
185 *>
186 *> \param[in] LDVR
187 *> \verbatim
188 *> LDVR is INTEGER
189 *> The leading dimension of the array VR. LDVR >= 1, and if
190 *> SIDE = 'R' or 'B', LDVR >= N.
191 *> \endverbatim
192 *>
193 *> \param[in] MM
194 *> \verbatim
195 *> MM is INTEGER
196 *> The number of columns in the arrays VL and/or VR. MM >= M.
197 *> \endverbatim
198 *>
199 *> \param[out] M
200 *> \verbatim
201 *> M is INTEGER
202 *> The number of columns in the arrays VL and/or VR actually
203 *> used to store the eigenvectors. If HOWMNY = 'A' or 'B', M
204 *> is set to N. Each selected real eigenvector occupies one
205 *> column and each selected complex eigenvector occupies two
206 *> columns.
207 *> \endverbatim
208 *>
209 *> \param[out] WORK
210 *> \verbatim
211 *> WORK is DOUBLE PRECISION array, dimension (6*N)
212 *> \endverbatim
213 *>
214 *> \param[out] INFO
215 *> \verbatim
216 *> INFO is INTEGER
217 *> = 0: successful exit.
218 *> < 0: if INFO = -i, the i-th argument had an illegal value.
219 *> > 0: the 2-by-2 block (INFO:INFO+1) does not have a complex
220 *> eigenvalue.
221 *> \endverbatim
222 *
223 * Authors:
224 * ========
225 *
226 *> \author Univ. of Tennessee
227 *> \author Univ. of California Berkeley
228 *> \author Univ. of Colorado Denver
229 *> \author NAG Ltd.
230 *
231 *> \date November 2011
232 *
233 *> \ingroup doubleGEcomputational
234 *
235 *> \par Further Details:
236 * =====================
237 *>
238 *> \verbatim
239 *>
240 *> Allocation of workspace:
241 *> ---------- -- ---------
242 *>
243 *> WORK( j ) = 1-norm of j-th column of A, above the diagonal
244 *> WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
245 *> WORK( 2*N+1:3*N ) = real part of eigenvector
246 *> WORK( 3*N+1:4*N ) = imaginary part of eigenvector
247 *> WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
248 *> WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
249 *>
250 *> Rowwise vs. columnwise solution methods:
251 *> ------- -- ---------- -------- -------
252 *>
253 *> Finding a generalized eigenvector consists basically of solving the
254 *> singular triangular system
255 *>
256 *> (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left)
257 *>
258 *> Consider finding the i-th right eigenvector (assume all eigenvalues
259 *> are real). The equation to be solved is:
260 *> n i
261 *> 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1
262 *> k=j k=j
263 *>
264 *> where C = (A - w B) (The components v(i+1:n) are 0.)
265 *>
266 *> The "rowwise" method is:
267 *>
268 *> (1) v(i) := 1
269 *> for j = i-1,. . .,1:
270 *> i
271 *> (2) compute s = - sum C(j,k) v(k) and
272 *> k=j+1
273 *>
274 *> (3) v(j) := s / C(j,j)
275 *>
276 *> Step 2 is sometimes called the "dot product" step, since it is an
277 *> inner product between the j-th row and the portion of the eigenvector
278 *> that has been computed so far.
279 *>
280 *> The "columnwise" method consists basically in doing the sums
281 *> for all the rows in parallel. As each v(j) is computed, the
282 *> contribution of v(j) times the j-th column of C is added to the
283 *> partial sums. Since FORTRAN arrays are stored columnwise, this has
284 *> the advantage that at each step, the elements of C that are accessed
285 *> are adjacent to one another, whereas with the rowwise method, the
286 *> elements accessed at a step are spaced LDS (and LDP) words apart.
287 *>
288 *> When finding left eigenvectors, the matrix in question is the
289 *> transpose of the one in storage, so the rowwise method then
290 *> actually accesses columns of A and B at each step, and so is the
291 *> preferred method.
292 *> \endverbatim
293 *>
294 * =====================================================================
295  SUBROUTINE dtgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
296  $ ldvl, vr, ldvr, mm, m, work, info )
297 *
298 * -- LAPACK computational routine (version 3.4.0) --
299 * -- LAPACK is a software package provided by Univ. of Tennessee, --
300 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
301 * November 2011
302 *
303 * .. Scalar Arguments ..
304  CHARACTER howmny, side
305  INTEGER info, ldp, lds, ldvl, ldvr, m, mm, n
306 * ..
307 * .. Array Arguments ..
308  LOGICAL select( * )
309  DOUBLE PRECISION p( ldp, * ), s( lds, * ), vl( ldvl, * ),
310  $ vr( ldvr, * ), work( * )
311 * ..
312 *
313 *
314 * =====================================================================
315 *
316 * .. Parameters ..
317  DOUBLE PRECISION zero, one, safety
318  parameter( zero = 0.0d+0, one = 1.0d+0,
319  $ safety = 1.0d+2 )
320 * ..
321 * .. Local Scalars ..
322  LOGICAL compl, compr, il2by2, ilabad, ilall, ilback,
323  $ ilbbad, ilcomp, ilcplx, lsa, lsb
324  INTEGER i, ibeg, ieig, iend, ihwmny, iinfo, im, iside,
325  $ j, ja, jc, je, jr, jw, na, nw
326  DOUBLE PRECISION acoef, acoefa, anorm, ascale, bcoefa, bcoefi,
327  $ bcoefr, big, bignum, bnorm, bscale, cim2a,
328  $ cim2b, cimaga, cimagb, cre2a, cre2b, creala,
329  $ crealb, dmin, safmin, salfar, sbeta, scale,
330  $ small, temp, temp2, temp2i, temp2r, ulp, xmax,
331  $ xscale
332 * ..
333 * .. Local Arrays ..
334  DOUBLE PRECISION bdiag( 2 ), sum( 2, 2 ), sums( 2, 2 ),
335  $ sump( 2, 2 )
336 * ..
337 * .. External Functions ..
338  LOGICAL lsame
339  DOUBLE PRECISION dlamch
340  EXTERNAL lsame, dlamch
341 * ..
342 * .. External Subroutines ..
343  EXTERNAL dgemv, dlabad, dlacpy, dlag2, dlaln2, xerbla
344 * ..
345 * .. Intrinsic Functions ..
346  INTRINSIC abs, max, min
347 * ..
348 * .. Executable Statements ..
349 *
350 * Decode and Test the input parameters
351 *
352  IF( lsame( howmny, 'A' ) ) THEN
353  ihwmny = 1
354  ilall = .true.
355  ilback = .false.
356  ELSE IF( lsame( howmny, 'S' ) ) THEN
357  ihwmny = 2
358  ilall = .false.
359  ilback = .false.
360  ELSE IF( lsame( howmny, 'B' ) ) THEN
361  ihwmny = 3
362  ilall = .true.
363  ilback = .true.
364  ELSE
365  ihwmny = -1
366  ilall = .true.
367  END IF
368 *
369  IF( lsame( side, 'R' ) ) THEN
370  iside = 1
371  compl = .false.
372  compr = .true.
373  ELSE IF( lsame( side, 'L' ) ) THEN
374  iside = 2
375  compl = .true.
376  compr = .false.
377  ELSE IF( lsame( side, 'B' ) ) THEN
378  iside = 3
379  compl = .true.
380  compr = .true.
381  ELSE
382  iside = -1
383  END IF
384 *
385  info = 0
386  IF( iside.LT.0 ) THEN
387  info = -1
388  ELSE IF( ihwmny.LT.0 ) THEN
389  info = -2
390  ELSE IF( n.LT.0 ) THEN
391  info = -4
392  ELSE IF( lds.LT.max( 1, n ) ) THEN
393  info = -6
394  ELSE IF( ldp.LT.max( 1, n ) ) THEN
395  info = -8
396  END IF
397  IF( info.NE.0 ) THEN
398  CALL xerbla( 'DTGEVC', -info )
399  return
400  END IF
401 *
402 * Count the number of eigenvectors to be computed
403 *
404  IF( .NOT.ilall ) THEN
405  im = 0
406  ilcplx = .false.
407  DO 10 j = 1, n
408  IF( ilcplx ) THEN
409  ilcplx = .false.
410  go to 10
411  END IF
412  IF( j.LT.n ) THEN
413  IF( s( j+1, j ).NE.zero )
414  $ ilcplx = .true.
415  END IF
416  IF( ilcplx ) THEN
417  IF( SELECT( j ) .OR. SELECT( j+1 ) )
418  $ im = im + 2
419  ELSE
420  IF( SELECT( j ) )
421  $ im = im + 1
422  END IF
423  10 continue
424  ELSE
425  im = n
426  END IF
427 *
428 * Check 2-by-2 diagonal blocks of A, B
429 *
430  ilabad = .false.
431  ilbbad = .false.
432  DO 20 j = 1, n - 1
433  IF( s( j+1, j ).NE.zero ) THEN
434  IF( p( j, j ).EQ.zero .OR. p( j+1, j+1 ).EQ.zero .OR.
435  $ p( j, j+1 ).NE.zero )ilbbad = .true.
436  IF( j.LT.n-1 ) THEN
437  IF( s( j+2, j+1 ).NE.zero )
438  $ ilabad = .true.
439  END IF
440  END IF
441  20 continue
442 *
443  IF( ilabad ) THEN
444  info = -5
445  ELSE IF( ilbbad ) THEN
446  info = -7
447  ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
448  info = -10
449  ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
450  info = -12
451  ELSE IF( mm.LT.im ) THEN
452  info = -13
453  END IF
454  IF( info.NE.0 ) THEN
455  CALL xerbla( 'DTGEVC', -info )
456  return
457  END IF
458 *
459 * Quick return if possible
460 *
461  m = im
462  IF( n.EQ.0 )
463  $ return
464 *
465 * Machine Constants
466 *
467  safmin = dlamch( 'Safe minimum' )
468  big = one / safmin
469  CALL dlabad( safmin, big )
470  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
471  small = safmin*n / ulp
472  big = one / small
473  bignum = one / ( safmin*n )
474 *
475 * Compute the 1-norm of each column of the strictly upper triangular
476 * part (i.e., excluding all elements belonging to the diagonal
477 * blocks) of A and B to check for possible overflow in the
478 * triangular solver.
479 *
480  anorm = abs( s( 1, 1 ) )
481  IF( n.GT.1 )
482  $ anorm = anorm + abs( s( 2, 1 ) )
483  bnorm = abs( p( 1, 1 ) )
484  work( 1 ) = zero
485  work( n+1 ) = zero
486 *
487  DO 50 j = 2, n
488  temp = zero
489  temp2 = zero
490  IF( s( j, j-1 ).EQ.zero ) THEN
491  iend = j - 1
492  ELSE
493  iend = j - 2
494  END IF
495  DO 30 i = 1, iend
496  temp = temp + abs( s( i, j ) )
497  temp2 = temp2 + abs( p( i, j ) )
498  30 continue
499  work( j ) = temp
500  work( n+j ) = temp2
501  DO 40 i = iend + 1, min( j+1, n )
502  temp = temp + abs( s( i, j ) )
503  temp2 = temp2 + abs( p( i, j ) )
504  40 continue
505  anorm = max( anorm, temp )
506  bnorm = max( bnorm, temp2 )
507  50 continue
508 *
509  ascale = one / max( anorm, safmin )
510  bscale = one / max( bnorm, safmin )
511 *
512 * Left eigenvectors
513 *
514  IF( compl ) THEN
515  ieig = 0
516 *
517 * Main loop over eigenvalues
518 *
519  ilcplx = .false.
520  DO 220 je = 1, n
521 *
522 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
523 * (b) this would be the second of a complex pair.
524 * Check for complex eigenvalue, so as to be sure of which
525 * entry(-ies) of SELECT to look at.
526 *
527  IF( ilcplx ) THEN
528  ilcplx = .false.
529  go to 220
530  END IF
531  nw = 1
532  IF( je.LT.n ) THEN
533  IF( s( je+1, je ).NE.zero ) THEN
534  ilcplx = .true.
535  nw = 2
536  END IF
537  END IF
538  IF( ilall ) THEN
539  ilcomp = .true.
540  ELSE IF( ilcplx ) THEN
541  ilcomp = SELECT( je ) .OR. SELECT( je+1 )
542  ELSE
543  ilcomp = SELECT( je )
544  END IF
545  IF( .NOT.ilcomp )
546  $ go to 220
547 *
548 * Decide if (a) singular pencil, (b) real eigenvalue, or
549 * (c) complex eigenvalue.
550 *
551  IF( .NOT.ilcplx ) THEN
552  IF( abs( s( je, je ) ).LE.safmin .AND.
553  $ abs( p( je, je ) ).LE.safmin ) THEN
554 *
555 * Singular matrix pencil -- return unit eigenvector
556 *
557  ieig = ieig + 1
558  DO 60 jr = 1, n
559  vl( jr, ieig ) = zero
560  60 continue
561  vl( ieig, ieig ) = one
562  go to 220
563  END IF
564  END IF
565 *
566 * Clear vector
567 *
568  DO 70 jr = 1, nw*n
569  work( 2*n+jr ) = zero
570  70 continue
571 * T
572 * Compute coefficients in ( a A - b B ) y = 0
573 * a is ACOEF
574 * b is BCOEFR + i*BCOEFI
575 *
576  IF( .NOT.ilcplx ) THEN
577 *
578 * Real eigenvalue
579 *
580  temp = one / max( abs( s( je, je ) )*ascale,
581  $ abs( p( je, je ) )*bscale, safmin )
582  salfar = ( temp*s( je, je ) )*ascale
583  sbeta = ( temp*p( je, je ) )*bscale
584  acoef = sbeta*ascale
585  bcoefr = salfar*bscale
586  bcoefi = zero
587 *
588 * Scale to avoid underflow
589 *
590  scale = one
591  lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
592  lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
593  $ small
594  IF( lsa )
595  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
596  IF( lsb )
597  $ scale = max( scale, ( small / abs( salfar ) )*
598  $ min( bnorm, big ) )
599  IF( lsa .OR. lsb ) THEN
600  scale = min( scale, one /
601  $ ( safmin*max( one, abs( acoef ),
602  $ abs( bcoefr ) ) ) )
603  IF( lsa ) THEN
604  acoef = ascale*( scale*sbeta )
605  ELSE
606  acoef = scale*acoef
607  END IF
608  IF( lsb ) THEN
609  bcoefr = bscale*( scale*salfar )
610  ELSE
611  bcoefr = scale*bcoefr
612  END IF
613  END IF
614  acoefa = abs( acoef )
615  bcoefa = abs( bcoefr )
616 *
617 * First component is 1
618 *
619  work( 2*n+je ) = one
620  xmax = one
621  ELSE
622 *
623 * Complex eigenvalue
624 *
625  CALL dlag2( s( je, je ), lds, p( je, je ), ldp,
626  $ safmin*safety, acoef, temp, bcoefr, temp2,
627  $ bcoefi )
628  bcoefi = -bcoefi
629  IF( bcoefi.EQ.zero ) THEN
630  info = je
631  return
632  END IF
633 *
634 * Scale to avoid over/underflow
635 *
636  acoefa = abs( acoef )
637  bcoefa = abs( bcoefr ) + abs( bcoefi )
638  scale = one
639  IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
640  $ scale = ( safmin / ulp ) / acoefa
641  IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
642  $ scale = max( scale, ( safmin / ulp ) / bcoefa )
643  IF( safmin*acoefa.GT.ascale )
644  $ scale = ascale / ( safmin*acoefa )
645  IF( safmin*bcoefa.GT.bscale )
646  $ scale = min( scale, bscale / ( safmin*bcoefa ) )
647  IF( scale.NE.one ) THEN
648  acoef = scale*acoef
649  acoefa = abs( acoef )
650  bcoefr = scale*bcoefr
651  bcoefi = scale*bcoefi
652  bcoefa = abs( bcoefr ) + abs( bcoefi )
653  END IF
654 *
655 * Compute first two components of eigenvector
656 *
657  temp = acoef*s( je+1, je )
658  temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
659  temp2i = -bcoefi*p( je, je )
660  IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) ) THEN
661  work( 2*n+je ) = one
662  work( 3*n+je ) = zero
663  work( 2*n+je+1 ) = -temp2r / temp
664  work( 3*n+je+1 ) = -temp2i / temp
665  ELSE
666  work( 2*n+je+1 ) = one
667  work( 3*n+je+1 ) = zero
668  temp = acoef*s( je, je+1 )
669  work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
670  $ s( je+1, je+1 ) ) / temp
671  work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
672  END IF
673  xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
674  $ abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
675  END IF
676 *
677  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
678 *
679 * T
680 * Triangular solve of (a A - b B) y = 0
681 *
682 * T
683 * (rowwise in (a A - b B) , or columnwise in (a A - b B) )
684 *
685  il2by2 = .false.
686 *
687  DO 160 j = je + nw, n
688  IF( il2by2 ) THEN
689  il2by2 = .false.
690  go to 160
691  END IF
692 *
693  na = 1
694  bdiag( 1 ) = p( j, j )
695  IF( j.LT.n ) THEN
696  IF( s( j+1, j ).NE.zero ) THEN
697  il2by2 = .true.
698  bdiag( 2 ) = p( j+1, j+1 )
699  na = 2
700  END IF
701  END IF
702 *
703 * Check whether scaling is necessary for dot products
704 *
705  xscale = one / max( one, xmax )
706  temp = max( work( j ), work( n+j ),
707  $ acoefa*work( j )+bcoefa*work( n+j ) )
708  IF( il2by2 )
709  $ temp = max( temp, work( j+1 ), work( n+j+1 ),
710  $ acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
711  IF( temp.GT.bignum*xscale ) THEN
712  DO 90 jw = 0, nw - 1
713  DO 80 jr = je, j - 1
714  work( ( jw+2 )*n+jr ) = xscale*
715  $ work( ( jw+2 )*n+jr )
716  80 continue
717  90 continue
718  xmax = xmax*xscale
719  END IF
720 *
721 * Compute dot products
722 *
723 * j-1
724 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
725 * k=je
726 *
727 * To reduce the op count, this is done as
728 *
729 * _ j-1 _ j-1
730 * a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) )
731 * k=je k=je
732 *
733 * which may cause underflow problems if A or B are close
734 * to underflow. (E.g., less than SMALL.)
735 *
736 *
737  DO 120 jw = 1, nw
738  DO 110 ja = 1, na
739  sums( ja, jw ) = zero
740  sump( ja, jw ) = zero
741 *
742  DO 100 jr = je, j - 1
743  sums( ja, jw ) = sums( ja, jw ) +
744  $ s( jr, j+ja-1 )*
745  $ work( ( jw+1 )*n+jr )
746  sump( ja, jw ) = sump( ja, jw ) +
747  $ p( jr, j+ja-1 )*
748  $ work( ( jw+1 )*n+jr )
749  100 continue
750  110 continue
751  120 continue
752 *
753  DO 130 ja = 1, na
754  IF( ilcplx ) THEN
755  sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
756  $ bcoefr*sump( ja, 1 ) -
757  $ bcoefi*sump( ja, 2 )
758  sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
759  $ bcoefr*sump( ja, 2 ) +
760  $ bcoefi*sump( ja, 1 )
761  ELSE
762  sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
763  $ bcoefr*sump( ja, 1 )
764  END IF
765  130 continue
766 *
767 * T
768 * Solve ( a A - b B ) y = SUM(,)
769 * with scaling and perturbation of the denominator
770 *
771  CALL dlaln2( .true., na, nw, dmin, acoef, s( j, j ), lds,
772  $ bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
773  $ bcoefi, work( 2*n+j ), n, scale, temp,
774  $ iinfo )
775  IF( scale.LT.one ) THEN
776  DO 150 jw = 0, nw - 1
777  DO 140 jr = je, j - 1
778  work( ( jw+2 )*n+jr ) = scale*
779  $ work( ( jw+2 )*n+jr )
780  140 continue
781  150 continue
782  xmax = scale*xmax
783  END IF
784  xmax = max( xmax, temp )
785  160 continue
786 *
787 * Copy eigenvector to VL, back transforming if
788 * HOWMNY='B'.
789 *
790  ieig = ieig + 1
791  IF( ilback ) THEN
792  DO 170 jw = 0, nw - 1
793  CALL dgemv( 'N', n, n+1-je, one, vl( 1, je ), ldvl,
794  $ work( ( jw+2 )*n+je ), 1, zero,
795  $ work( ( jw+4 )*n+1 ), 1 )
796  170 continue
797  CALL dlacpy( ' ', n, nw, work( 4*n+1 ), n, vl( 1, je ),
798  $ ldvl )
799  ibeg = 1
800  ELSE
801  CALL dlacpy( ' ', n, nw, work( 2*n+1 ), n, vl( 1, ieig ),
802  $ ldvl )
803  ibeg = je
804  END IF
805 *
806 * Scale eigenvector
807 *
808  xmax = zero
809  IF( ilcplx ) THEN
810  DO 180 j = ibeg, n
811  xmax = max( xmax, abs( vl( j, ieig ) )+
812  $ abs( vl( j, ieig+1 ) ) )
813  180 continue
814  ELSE
815  DO 190 j = ibeg, n
816  xmax = max( xmax, abs( vl( j, ieig ) ) )
817  190 continue
818  END IF
819 *
820  IF( xmax.GT.safmin ) THEN
821  xscale = one / xmax
822 *
823  DO 210 jw = 0, nw - 1
824  DO 200 jr = ibeg, n
825  vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
826  200 continue
827  210 continue
828  END IF
829  ieig = ieig + nw - 1
830 *
831  220 continue
832  END IF
833 *
834 * Right eigenvectors
835 *
836  IF( compr ) THEN
837  ieig = im + 1
838 *
839 * Main loop over eigenvalues
840 *
841  ilcplx = .false.
842  DO 500 je = n, 1, -1
843 *
844 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
845 * (b) this would be the second of a complex pair.
846 * Check for complex eigenvalue, so as to be sure of which
847 * entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
848 * or SELECT(JE-1).
849 * If this is a complex pair, the 2-by-2 diagonal block
850 * corresponding to the eigenvalue is in rows/columns JE-1:JE
851 *
852  IF( ilcplx ) THEN
853  ilcplx = .false.
854  go to 500
855  END IF
856  nw = 1
857  IF( je.GT.1 ) THEN
858  IF( s( je, je-1 ).NE.zero ) THEN
859  ilcplx = .true.
860  nw = 2
861  END IF
862  END IF
863  IF( ilall ) THEN
864  ilcomp = .true.
865  ELSE IF( ilcplx ) THEN
866  ilcomp = SELECT( je ) .OR. SELECT( je-1 )
867  ELSE
868  ilcomp = SELECT( je )
869  END IF
870  IF( .NOT.ilcomp )
871  $ go to 500
872 *
873 * Decide if (a) singular pencil, (b) real eigenvalue, or
874 * (c) complex eigenvalue.
875 *
876  IF( .NOT.ilcplx ) THEN
877  IF( abs( s( je, je ) ).LE.safmin .AND.
878  $ abs( p( je, je ) ).LE.safmin ) THEN
879 *
880 * Singular matrix pencil -- unit eigenvector
881 *
882  ieig = ieig - 1
883  DO 230 jr = 1, n
884  vr( jr, ieig ) = zero
885  230 continue
886  vr( ieig, ieig ) = one
887  go to 500
888  END IF
889  END IF
890 *
891 * Clear vector
892 *
893  DO 250 jw = 0, nw - 1
894  DO 240 jr = 1, n
895  work( ( jw+2 )*n+jr ) = zero
896  240 continue
897  250 continue
898 *
899 * Compute coefficients in ( a A - b B ) x = 0
900 * a is ACOEF
901 * b is BCOEFR + i*BCOEFI
902 *
903  IF( .NOT.ilcplx ) THEN
904 *
905 * Real eigenvalue
906 *
907  temp = one / max( abs( s( je, je ) )*ascale,
908  $ abs( p( je, je ) )*bscale, safmin )
909  salfar = ( temp*s( je, je ) )*ascale
910  sbeta = ( temp*p( je, je ) )*bscale
911  acoef = sbeta*ascale
912  bcoefr = salfar*bscale
913  bcoefi = zero
914 *
915 * Scale to avoid underflow
916 *
917  scale = one
918  lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
919  lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
920  $ small
921  IF( lsa )
922  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
923  IF( lsb )
924  $ scale = max( scale, ( small / abs( salfar ) )*
925  $ min( bnorm, big ) )
926  IF( lsa .OR. lsb ) THEN
927  scale = min( scale, one /
928  $ ( safmin*max( one, abs( acoef ),
929  $ abs( bcoefr ) ) ) )
930  IF( lsa ) THEN
931  acoef = ascale*( scale*sbeta )
932  ELSE
933  acoef = scale*acoef
934  END IF
935  IF( lsb ) THEN
936  bcoefr = bscale*( scale*salfar )
937  ELSE
938  bcoefr = scale*bcoefr
939  END IF
940  END IF
941  acoefa = abs( acoef )
942  bcoefa = abs( bcoefr )
943 *
944 * First component is 1
945 *
946  work( 2*n+je ) = one
947  xmax = one
948 *
949 * Compute contribution from column JE of A and B to sum
950 * (See "Further Details", above.)
951 *
952  DO 260 jr = 1, je - 1
953  work( 2*n+jr ) = bcoefr*p( jr, je ) -
954  $ acoef*s( jr, je )
955  260 continue
956  ELSE
957 *
958 * Complex eigenvalue
959 *
960  CALL dlag2( s( je-1, je-1 ), lds, p( je-1, je-1 ), ldp,
961  $ safmin*safety, acoef, temp, bcoefr, temp2,
962  $ bcoefi )
963  IF( bcoefi.EQ.zero ) THEN
964  info = je - 1
965  return
966  END IF
967 *
968 * Scale to avoid over/underflow
969 *
970  acoefa = abs( acoef )
971  bcoefa = abs( bcoefr ) + abs( bcoefi )
972  scale = one
973  IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
974  $ scale = ( safmin / ulp ) / acoefa
975  IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
976  $ scale = max( scale, ( safmin / ulp ) / bcoefa )
977  IF( safmin*acoefa.GT.ascale )
978  $ scale = ascale / ( safmin*acoefa )
979  IF( safmin*bcoefa.GT.bscale )
980  $ scale = min( scale, bscale / ( safmin*bcoefa ) )
981  IF( scale.NE.one ) THEN
982  acoef = scale*acoef
983  acoefa = abs( acoef )
984  bcoefr = scale*bcoefr
985  bcoefi = scale*bcoefi
986  bcoefa = abs( bcoefr ) + abs( bcoefi )
987  END IF
988 *
989 * Compute first two components of eigenvector
990 * and contribution to sums
991 *
992  temp = acoef*s( je, je-1 )
993  temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
994  temp2i = -bcoefi*p( je, je )
995  IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) ) THEN
996  work( 2*n+je ) = one
997  work( 3*n+je ) = zero
998  work( 2*n+je-1 ) = -temp2r / temp
999  work( 3*n+je-1 ) = -temp2i / temp
1000  ELSE
1001  work( 2*n+je-1 ) = one
1002  work( 3*n+je-1 ) = zero
1003  temp = acoef*s( je-1, je )
1004  work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
1005  $ s( je-1, je-1 ) ) / temp
1006  work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
1007  END IF
1008 *
1009  xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
1010  $ abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
1011 *
1012 * Compute contribution from columns JE and JE-1
1013 * of A and B to the sums.
1014 *
1015  creala = acoef*work( 2*n+je-1 )
1016  cimaga = acoef*work( 3*n+je-1 )
1017  crealb = bcoefr*work( 2*n+je-1 ) -
1018  $ bcoefi*work( 3*n+je-1 )
1019  cimagb = bcoefi*work( 2*n+je-1 ) +
1020  $ bcoefr*work( 3*n+je-1 )
1021  cre2a = acoef*work( 2*n+je )
1022  cim2a = acoef*work( 3*n+je )
1023  cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
1024  cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
1025  DO 270 jr = 1, je - 2
1026  work( 2*n+jr ) = -creala*s( jr, je-1 ) +
1027  $ crealb*p( jr, je-1 ) -
1028  $ cre2a*s( jr, je ) + cre2b*p( jr, je )
1029  work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
1030  $ cimagb*p( jr, je-1 ) -
1031  $ cim2a*s( jr, je ) + cim2b*p( jr, je )
1032  270 continue
1033  END IF
1034 *
1035  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
1036 *
1037 * Columnwise triangular solve of (a A - b B) x = 0
1038 *
1039  il2by2 = .false.
1040  DO 370 j = je - nw, 1, -1
1041 *
1042 * If a 2-by-2 block, is in position j-1:j, wait until
1043 * next iteration to process it (when it will be j:j+1)
1044 *
1045  IF( .NOT.il2by2 .AND. j.GT.1 ) THEN
1046  IF( s( j, j-1 ).NE.zero ) THEN
1047  il2by2 = .true.
1048  go to 370
1049  END IF
1050  END IF
1051  bdiag( 1 ) = p( j, j )
1052  IF( il2by2 ) THEN
1053  na = 2
1054  bdiag( 2 ) = p( j+1, j+1 )
1055  ELSE
1056  na = 1
1057  END IF
1058 *
1059 * Compute x(j) (and x(j+1), if 2-by-2 block)
1060 *
1061  CALL dlaln2( .false., na, nw, dmin, acoef, s( j, j ),
1062  $ lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
1063  $ n, bcoefr, bcoefi, sum, 2, scale, temp,
1064  $ iinfo )
1065  IF( scale.LT.one ) THEN
1066 *
1067  DO 290 jw = 0, nw - 1
1068  DO 280 jr = 1, je
1069  work( ( jw+2 )*n+jr ) = scale*
1070  $ work( ( jw+2 )*n+jr )
1071  280 continue
1072  290 continue
1073  END IF
1074  xmax = max( scale*xmax, temp )
1075 *
1076  DO 310 jw = 1, nw
1077  DO 300 ja = 1, na
1078  work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
1079  300 continue
1080  310 continue
1081 *
1082 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1083 *
1084  IF( j.GT.1 ) THEN
1085 *
1086 * Check whether scaling is necessary for sum.
1087 *
1088  xscale = one / max( one, xmax )
1089  temp = acoefa*work( j ) + bcoefa*work( n+j )
1090  IF( il2by2 )
1091  $ temp = max( temp, acoefa*work( j+1 )+bcoefa*
1092  $ work( n+j+1 ) )
1093  temp = max( temp, acoefa, bcoefa )
1094  IF( temp.GT.bignum*xscale ) THEN
1095 *
1096  DO 330 jw = 0, nw - 1
1097  DO 320 jr = 1, je
1098  work( ( jw+2 )*n+jr ) = xscale*
1099  $ work( ( jw+2 )*n+jr )
1100  320 continue
1101  330 continue
1102  xmax = xmax*xscale
1103  END IF
1104 *
1105 * Compute the contributions of the off-diagonals of
1106 * column j (and j+1, if 2-by-2 block) of A and B to the
1107 * sums.
1108 *
1109 *
1110  DO 360 ja = 1, na
1111  IF( ilcplx ) THEN
1112  creala = acoef*work( 2*n+j+ja-1 )
1113  cimaga = acoef*work( 3*n+j+ja-1 )
1114  crealb = bcoefr*work( 2*n+j+ja-1 ) -
1115  $ bcoefi*work( 3*n+j+ja-1 )
1116  cimagb = bcoefi*work( 2*n+j+ja-1 ) +
1117  $ bcoefr*work( 3*n+j+ja-1 )
1118  DO 340 jr = 1, j - 1
1119  work( 2*n+jr ) = work( 2*n+jr ) -
1120  $ creala*s( jr, j+ja-1 ) +
1121  $ crealb*p( jr, j+ja-1 )
1122  work( 3*n+jr ) = work( 3*n+jr ) -
1123  $ cimaga*s( jr, j+ja-1 ) +
1124  $ cimagb*p( jr, j+ja-1 )
1125  340 continue
1126  ELSE
1127  creala = acoef*work( 2*n+j+ja-1 )
1128  crealb = bcoefr*work( 2*n+j+ja-1 )
1129  DO 350 jr = 1, j - 1
1130  work( 2*n+jr ) = work( 2*n+jr ) -
1131  $ creala*s( jr, j+ja-1 ) +
1132  $ crealb*p( jr, j+ja-1 )
1133  350 continue
1134  END IF
1135  360 continue
1136  END IF
1137 *
1138  il2by2 = .false.
1139  370 continue
1140 *
1141 * Copy eigenvector to VR, back transforming if
1142 * HOWMNY='B'.
1143 *
1144  ieig = ieig - nw
1145  IF( ilback ) THEN
1146 *
1147  DO 410 jw = 0, nw - 1
1148  DO 380 jr = 1, n
1149  work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
1150  $ vr( jr, 1 )
1151  380 continue
1152 *
1153 * A series of compiler directives to defeat
1154 * vectorization for the next loop
1155 *
1156 *
1157  DO 400 jc = 2, je
1158  DO 390 jr = 1, n
1159  work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
1160  $ work( ( jw+2 )*n+jc )*vr( jr, jc )
1161  390 continue
1162  400 continue
1163  410 continue
1164 *
1165  DO 430 jw = 0, nw - 1
1166  DO 420 jr = 1, n
1167  vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
1168  420 continue
1169  430 continue
1170 *
1171  iend = n
1172  ELSE
1173  DO 450 jw = 0, nw - 1
1174  DO 440 jr = 1, n
1175  vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
1176  440 continue
1177  450 continue
1178 *
1179  iend = je
1180  END IF
1181 *
1182 * Scale eigenvector
1183 *
1184  xmax = zero
1185  IF( ilcplx ) THEN
1186  DO 460 j = 1, iend
1187  xmax = max( xmax, abs( vr( j, ieig ) )+
1188  $ abs( vr( j, ieig+1 ) ) )
1189  460 continue
1190  ELSE
1191  DO 470 j = 1, iend
1192  xmax = max( xmax, abs( vr( j, ieig ) ) )
1193  470 continue
1194  END IF
1195 *
1196  IF( xmax.GT.safmin ) THEN
1197  xscale = one / xmax
1198  DO 490 jw = 0, nw - 1
1199  DO 480 jr = 1, iend
1200  vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
1201  480 continue
1202  490 continue
1203  END IF
1204  500 continue
1205  END IF
1206 *
1207  return
1208 *
1209 * End of DTGEVC
1210 *
1211  END