LAPACK  3.11.0
LAPACK: Linear Algebra PACKage
dtgevc.f
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 *> \ingroup tgevc
232 *
233 *> \par Further Details:
234 * =====================
235 *>
236 *> \verbatim
237 *>
238 *> Allocation of workspace:
239 *> ---------- -- ---------
240 *>
241 *> WORK( j ) = 1-norm of j-th column of A, above the diagonal
242 *> WORK( N+j ) = 1-norm of j-th column of B, above the diagonal
243 *> WORK( 2*N+1:3*N ) = real part of eigenvector
244 *> WORK( 3*N+1:4*N ) = imaginary part of eigenvector
245 *> WORK( 4*N+1:5*N ) = real part of back-transformed eigenvector
246 *> WORK( 5*N+1:6*N ) = imaginary part of back-transformed eigenvector
247 *>
248 *> Rowwise vs. columnwise solution methods:
249 *> ------- -- ---------- -------- -------
250 *>
251 *> Finding a generalized eigenvector consists basically of solving the
252 *> singular triangular system
253 *>
254 *> (A - w B) x = 0 (for right) or: (A - w B)**H y = 0 (for left)
255 *>
256 *> Consider finding the i-th right eigenvector (assume all eigenvalues
257 *> are real). The equation to be solved is:
258 *> n i
259 *> 0 = sum C(j,k) v(k) = sum C(j,k) v(k) for j = i,. . .,1
260 *> k=j k=j
261 *>
262 *> where C = (A - w B) (The components v(i+1:n) are 0.)
263 *>
264 *> The "rowwise" method is:
265 *>
266 *> (1) v(i) := 1
267 *> for j = i-1,. . .,1:
268 *> i
269 *> (2) compute s = - sum C(j,k) v(k) and
270 *> k=j+1
271 *>
272 *> (3) v(j) := s / C(j,j)
273 *>
274 *> Step 2 is sometimes called the "dot product" step, since it is an
275 *> inner product between the j-th row and the portion of the eigenvector
276 *> that has been computed so far.
277 *>
278 *> The "columnwise" method consists basically in doing the sums
279 *> for all the rows in parallel. As each v(j) is computed, the
280 *> contribution of v(j) times the j-th column of C is added to the
281 *> partial sums. Since FORTRAN arrays are stored columnwise, this has
282 *> the advantage that at each step, the elements of C that are accessed
283 *> are adjacent to one another, whereas with the rowwise method, the
284 *> elements accessed at a step are spaced LDS (and LDP) words apart.
285 *>
286 *> When finding left eigenvectors, the matrix in question is the
287 *> transpose of the one in storage, so the rowwise method then
288 *> actually accesses columns of A and B at each step, and so is the
289 *> preferred method.
290 *> \endverbatim
291 *>
292 * =====================================================================
293  SUBROUTINE dtgevc( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
294  $ LDVL, VR, LDVR, MM, M, WORK, INFO )
295 *
296 * -- LAPACK computational routine --
297 * -- LAPACK is a software package provided by Univ. of Tennessee, --
298 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
299 *
300 * .. Scalar Arguments ..
301  CHARACTER HOWMNY, SIDE
302  INTEGER INFO, LDP, LDS, LDVL, LDVR, M, MM, N
303 * ..
304 * .. Array Arguments ..
305  LOGICAL SELECT( * )
306  DOUBLE PRECISION P( ldp, * ), S( lds, * ), VL( ldvl, * ),
307  $ vr( ldvr, * ), work( * )
308 * ..
309 *
310 *
311 * =====================================================================
312 *
313 * .. Parameters ..
314  DOUBLE PRECISION ZERO, ONE, SAFETY
315  parameter( zero = 0.0d+0, one = 1.0d+0,
316  $ safety = 1.0d+2 )
317 * ..
318 * .. Local Scalars ..
319  LOGICAL COMPL, COMPR, IL2BY2, ILABAD, ILALL, ILBACK,
320  $ ilbbad, ilcomp, ilcplx, lsa, lsb
321  INTEGER I, IBEG, IEIG, IEND, IHWMNY, IINFO, IM, ISIDE,
322  $ j, ja, jc, je, jr, jw, na, nw
323  DOUBLE PRECISION ACOEF, ACOEFA, ANORM, ASCALE, BCOEFA, BCOEFI,
324  $ bcoefr, big, bignum, bnorm, bscale, cim2a,
325  $ cim2b, cimaga, cimagb, cre2a, cre2b, creala,
326  $ crealb, dmin, safmin, salfar, sbeta, scale,
327  $ small, temp, temp2, temp2i, temp2r, ulp, xmax,
328  $ xscale
329 * ..
330 * .. Local Arrays ..
331  DOUBLE PRECISION BDIAG( 2 ), SUM( 2, 2 ), SUMS( 2, 2 ),
332  $ sump( 2, 2 )
333 * ..
334 * .. External Functions ..
335  LOGICAL LSAME
336  DOUBLE PRECISION DLAMCH
337  EXTERNAL lsame, dlamch
338 * ..
339 * .. External Subroutines ..
340  EXTERNAL dgemv, dlabad, dlacpy, dlag2, dlaln2, xerbla
341 * ..
342 * .. Intrinsic Functions ..
343  INTRINSIC abs, max, min
344 * ..
345 * .. Executable Statements ..
346 *
347 * Decode and Test the input parameters
348 *
349  IF( lsame( howmny, 'A' ) ) THEN
350  ihwmny = 1
351  ilall = .true.
352  ilback = .false.
353  ELSE IF( lsame( howmny, 'S' ) ) THEN
354  ihwmny = 2
355  ilall = .false.
356  ilback = .false.
357  ELSE IF( lsame( howmny, 'B' ) ) THEN
358  ihwmny = 3
359  ilall = .true.
360  ilback = .true.
361  ELSE
362  ihwmny = -1
363  ilall = .true.
364  END IF
365 *
366  IF( lsame( side, 'R' ) ) THEN
367  iside = 1
368  compl = .false.
369  compr = .true.
370  ELSE IF( lsame( side, 'L' ) ) THEN
371  iside = 2
372  compl = .true.
373  compr = .false.
374  ELSE IF( lsame( side, 'B' ) ) THEN
375  iside = 3
376  compl = .true.
377  compr = .true.
378  ELSE
379  iside = -1
380  END IF
381 *
382  info = 0
383  IF( iside.LT.0 ) THEN
384  info = -1
385  ELSE IF( ihwmny.LT.0 ) THEN
386  info = -2
387  ELSE IF( n.LT.0 ) THEN
388  info = -4
389  ELSE IF( lds.LT.max( 1, n ) ) THEN
390  info = -6
391  ELSE IF( ldp.LT.max( 1, n ) ) THEN
392  info = -8
393  END IF
394  IF( info.NE.0 ) THEN
395  CALL xerbla( 'DTGEVC', -info )
396  RETURN
397  END IF
398 *
399 * Count the number of eigenvectors to be computed
400 *
401  IF( .NOT.ilall ) THEN
402  im = 0
403  ilcplx = .false.
404  DO 10 j = 1, n
405  IF( ilcplx ) THEN
406  ilcplx = .false.
407  GO TO 10
408  END IF
409  IF( j.LT.n ) THEN
410  IF( s( j+1, j ).NE.zero )
411  $ ilcplx = .true.
412  END IF
413  IF( ilcplx ) THEN
414  IF( SELECT( j ) .OR. SELECT( j+1 ) )
415  $ im = im + 2
416  ELSE
417  IF( SELECT( j ) )
418  $ im = im + 1
419  END IF
420  10 CONTINUE
421  ELSE
422  im = n
423  END IF
424 *
425 * Check 2-by-2 diagonal blocks of A, B
426 *
427  ilabad = .false.
428  ilbbad = .false.
429  DO 20 j = 1, n - 1
430  IF( s( j+1, j ).NE.zero ) THEN
431  IF( p( j, j ).EQ.zero .OR. p( j+1, j+1 ).EQ.zero .OR.
432  $ p( j, j+1 ).NE.zero )ilbbad = .true.
433  IF( j.LT.n-1 ) THEN
434  IF( s( j+2, j+1 ).NE.zero )
435  $ ilabad = .true.
436  END IF
437  END IF
438  20 CONTINUE
439 *
440  IF( ilabad ) THEN
441  info = -5
442  ELSE IF( ilbbad ) THEN
443  info = -7
444  ELSE IF( compl .AND. ldvl.LT.n .OR. ldvl.LT.1 ) THEN
445  info = -10
446  ELSE IF( compr .AND. ldvr.LT.n .OR. ldvr.LT.1 ) THEN
447  info = -12
448  ELSE IF( mm.LT.im ) THEN
449  info = -13
450  END IF
451  IF( info.NE.0 ) THEN
452  CALL xerbla( 'DTGEVC', -info )
453  RETURN
454  END IF
455 *
456 * Quick return if possible
457 *
458  m = im
459  IF( n.EQ.0 )
460  $ RETURN
461 *
462 * Machine Constants
463 *
464  safmin = dlamch( 'Safe minimum' )
465  big = one / safmin
466  CALL dlabad( safmin, big )
467  ulp = dlamch( 'Epsilon' )*dlamch( 'Base' )
468  small = safmin*n / ulp
469  big = one / small
470  bignum = one / ( safmin*n )
471 *
472 * Compute the 1-norm of each column of the strictly upper triangular
473 * part (i.e., excluding all elements belonging to the diagonal
474 * blocks) of A and B to check for possible overflow in the
475 * triangular solver.
476 *
477  anorm = abs( s( 1, 1 ) )
478  IF( n.GT.1 )
479  $ anorm = anorm + abs( s( 2, 1 ) )
480  bnorm = abs( p( 1, 1 ) )
481  work( 1 ) = zero
482  work( n+1 ) = zero
483 *
484  DO 50 j = 2, n
485  temp = zero
486  temp2 = zero
487  IF( s( j, j-1 ).EQ.zero ) THEN
488  iend = j - 1
489  ELSE
490  iend = j - 2
491  END IF
492  DO 30 i = 1, iend
493  temp = temp + abs( s( i, j ) )
494  temp2 = temp2 + abs( p( i, j ) )
495  30 CONTINUE
496  work( j ) = temp
497  work( n+j ) = temp2
498  DO 40 i = iend + 1, min( j+1, n )
499  temp = temp + abs( s( i, j ) )
500  temp2 = temp2 + abs( p( i, j ) )
501  40 CONTINUE
502  anorm = max( anorm, temp )
503  bnorm = max( bnorm, temp2 )
504  50 CONTINUE
505 *
506  ascale = one / max( anorm, safmin )
507  bscale = one / max( bnorm, safmin )
508 *
509 * Left eigenvectors
510 *
511  IF( compl ) THEN
512  ieig = 0
513 *
514 * Main loop over eigenvalues
515 *
516  ilcplx = .false.
517  DO 220 je = 1, n
518 *
519 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
520 * (b) this would be the second of a complex pair.
521 * Check for complex eigenvalue, so as to be sure of which
522 * entry(-ies) of SELECT to look at.
523 *
524  IF( ilcplx ) THEN
525  ilcplx = .false.
526  GO TO 220
527  END IF
528  nw = 1
529  IF( je.LT.n ) THEN
530  IF( s( je+1, je ).NE.zero ) THEN
531  ilcplx = .true.
532  nw = 2
533  END IF
534  END IF
535  IF( ilall ) THEN
536  ilcomp = .true.
537  ELSE IF( ilcplx ) THEN
538  ilcomp = SELECT( je ) .OR. SELECT( je+1 )
539  ELSE
540  ilcomp = SELECT( je )
541  END IF
542  IF( .NOT.ilcomp )
543  $ GO TO 220
544 *
545 * Decide if (a) singular pencil, (b) real eigenvalue, or
546 * (c) complex eigenvalue.
547 *
548  IF( .NOT.ilcplx ) THEN
549  IF( abs( s( je, je ) ).LE.safmin .AND.
550  $ abs( p( je, je ) ).LE.safmin ) THEN
551 *
552 * Singular matrix pencil -- return unit eigenvector
553 *
554  ieig = ieig + 1
555  DO 60 jr = 1, n
556  vl( jr, ieig ) = zero
557  60 CONTINUE
558  vl( ieig, ieig ) = one
559  GO TO 220
560  END IF
561  END IF
562 *
563 * Clear vector
564 *
565  DO 70 jr = 1, nw*n
566  work( 2*n+jr ) = zero
567  70 CONTINUE
568 * T
569 * Compute coefficients in ( a A - b B ) y = 0
570 * a is ACOEF
571 * b is BCOEFR + i*BCOEFI
572 *
573  IF( .NOT.ilcplx ) THEN
574 *
575 * Real eigenvalue
576 *
577  temp = one / max( abs( s( je, je ) )*ascale,
578  $ abs( p( je, je ) )*bscale, safmin )
579  salfar = ( temp*s( je, je ) )*ascale
580  sbeta = ( temp*p( je, je ) )*bscale
581  acoef = sbeta*ascale
582  bcoefr = salfar*bscale
583  bcoefi = zero
584 *
585 * Scale to avoid underflow
586 *
587  scale = one
588  lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
589  lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
590  $ small
591  IF( lsa )
592  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
593  IF( lsb )
594  $ scale = max( scale, ( small / abs( salfar ) )*
595  $ min( bnorm, big ) )
596  IF( lsa .OR. lsb ) THEN
597  scale = min( scale, one /
598  $ ( safmin*max( one, abs( acoef ),
599  $ abs( bcoefr ) ) ) )
600  IF( lsa ) THEN
601  acoef = ascale*( scale*sbeta )
602  ELSE
603  acoef = scale*acoef
604  END IF
605  IF( lsb ) THEN
606  bcoefr = bscale*( scale*salfar )
607  ELSE
608  bcoefr = scale*bcoefr
609  END IF
610  END IF
611  acoefa = abs( acoef )
612  bcoefa = abs( bcoefr )
613 *
614 * First component is 1
615 *
616  work( 2*n+je ) = one
617  xmax = one
618  ELSE
619 *
620 * Complex eigenvalue
621 *
622  CALL dlag2( s( je, je ), lds, p( je, je ), ldp,
623  $ safmin*safety, acoef, temp, bcoefr, temp2,
624  $ bcoefi )
625  bcoefi = -bcoefi
626  IF( bcoefi.EQ.zero ) THEN
627  info = je
628  RETURN
629  END IF
630 *
631 * Scale to avoid over/underflow
632 *
633  acoefa = abs( acoef )
634  bcoefa = abs( bcoefr ) + abs( bcoefi )
635  scale = one
636  IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
637  $ scale = ( safmin / ulp ) / acoefa
638  IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
639  $ scale = max( scale, ( safmin / ulp ) / bcoefa )
640  IF( safmin*acoefa.GT.ascale )
641  $ scale = ascale / ( safmin*acoefa )
642  IF( safmin*bcoefa.GT.bscale )
643  $ scale = min( scale, bscale / ( safmin*bcoefa ) )
644  IF( scale.NE.one ) THEN
645  acoef = scale*acoef
646  acoefa = abs( acoef )
647  bcoefr = scale*bcoefr
648  bcoefi = scale*bcoefi
649  bcoefa = abs( bcoefr ) + abs( bcoefi )
650  END IF
651 *
652 * Compute first two components of eigenvector
653 *
654  temp = acoef*s( je+1, je )
655  temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
656  temp2i = -bcoefi*p( je, je )
657  IF( abs( temp ).GT.abs( temp2r )+abs( temp2i ) ) THEN
658  work( 2*n+je ) = one
659  work( 3*n+je ) = zero
660  work( 2*n+je+1 ) = -temp2r / temp
661  work( 3*n+je+1 ) = -temp2i / temp
662  ELSE
663  work( 2*n+je+1 ) = one
664  work( 3*n+je+1 ) = zero
665  temp = acoef*s( je, je+1 )
666  work( 2*n+je ) = ( bcoefr*p( je+1, je+1 )-acoef*
667  $ s( je+1, je+1 ) ) / temp
668  work( 3*n+je ) = bcoefi*p( je+1, je+1 ) / temp
669  END IF
670  xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
671  $ abs( work( 2*n+je+1 ) )+abs( work( 3*n+je+1 ) ) )
672  END IF
673 *
674  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
675 *
676 * T
677 * Triangular solve of (a A - b B) y = 0
678 *
679 * T
680 * (rowwise in (a A - b B) , or columnwise in (a A - b B) )
681 *
682  il2by2 = .false.
683 *
684  DO 160 j = je + nw, n
685  IF( il2by2 ) THEN
686  il2by2 = .false.
687  GO TO 160
688  END IF
689 *
690  na = 1
691  bdiag( 1 ) = p( j, j )
692  IF( j.LT.n ) THEN
693  IF( s( j+1, j ).NE.zero ) THEN
694  il2by2 = .true.
695  bdiag( 2 ) = p( j+1, j+1 )
696  na = 2
697  END IF
698  END IF
699 *
700 * Check whether scaling is necessary for dot products
701 *
702  xscale = one / max( one, xmax )
703  temp = max( work( j ), work( n+j ),
704  $ acoefa*work( j )+bcoefa*work( n+j ) )
705  IF( il2by2 )
706  $ temp = max( temp, work( j+1 ), work( n+j+1 ),
707  $ acoefa*work( j+1 )+bcoefa*work( n+j+1 ) )
708  IF( temp.GT.bignum*xscale ) THEN
709  DO 90 jw = 0, nw - 1
710  DO 80 jr = je, j - 1
711  work( ( jw+2 )*n+jr ) = xscale*
712  $ work( ( jw+2 )*n+jr )
713  80 CONTINUE
714  90 CONTINUE
715  xmax = xmax*xscale
716  END IF
717 *
718 * Compute dot products
719 *
720 * j-1
721 * SUM = sum conjg( a*S(k,j) - b*P(k,j) )*x(k)
722 * k=je
723 *
724 * To reduce the op count, this is done as
725 *
726 * _ j-1 _ j-1
727 * a*conjg( sum S(k,j)*x(k) ) - b*conjg( sum P(k,j)*x(k) )
728 * k=je k=je
729 *
730 * which may cause underflow problems if A or B are close
731 * to underflow. (E.g., less than SMALL.)
732 *
733 *
734  DO 120 jw = 1, nw
735  DO 110 ja = 1, na
736  sums( ja, jw ) = zero
737  sump( ja, jw ) = zero
738 *
739  DO 100 jr = je, j - 1
740  sums( ja, jw ) = sums( ja, jw ) +
741  $ s( jr, j+ja-1 )*
742  $ work( ( jw+1 )*n+jr )
743  sump( ja, jw ) = sump( ja, jw ) +
744  $ p( jr, j+ja-1 )*
745  $ work( ( jw+1 )*n+jr )
746  100 CONTINUE
747  110 CONTINUE
748  120 CONTINUE
749 *
750  DO 130 ja = 1, na
751  IF( ilcplx ) THEN
752  sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
753  $ bcoefr*sump( ja, 1 ) -
754  $ bcoefi*sump( ja, 2 )
755  sum( ja, 2 ) = -acoef*sums( ja, 2 ) +
756  $ bcoefr*sump( ja, 2 ) +
757  $ bcoefi*sump( ja, 1 )
758  ELSE
759  sum( ja, 1 ) = -acoef*sums( ja, 1 ) +
760  $ bcoefr*sump( ja, 1 )
761  END IF
762  130 CONTINUE
763 *
764 * T
765 * Solve ( a A - b B ) y = SUM(,)
766 * with scaling and perturbation of the denominator
767 *
768  CALL dlaln2( .true., na, nw, dmin, acoef, s( j, j ), lds,
769  $ bdiag( 1 ), bdiag( 2 ), sum, 2, bcoefr,
770  $ bcoefi, work( 2*n+j ), n, scale, temp,
771  $ iinfo )
772  IF( scale.LT.one ) THEN
773  DO 150 jw = 0, nw - 1
774  DO 140 jr = je, j - 1
775  work( ( jw+2 )*n+jr ) = scale*
776  $ work( ( jw+2 )*n+jr )
777  140 CONTINUE
778  150 CONTINUE
779  xmax = scale*xmax
780  END IF
781  xmax = max( xmax, temp )
782  160 CONTINUE
783 *
784 * Copy eigenvector to VL, back transforming if
785 * HOWMNY='B'.
786 *
787  ieig = ieig + 1
788  IF( ilback ) THEN
789  DO 170 jw = 0, nw - 1
790  CALL dgemv( 'N', n, n+1-je, one, vl( 1, je ), ldvl,
791  $ work( ( jw+2 )*n+je ), 1, zero,
792  $ work( ( jw+4 )*n+1 ), 1 )
793  170 CONTINUE
794  CALL dlacpy( ' ', n, nw, work( 4*n+1 ), n, vl( 1, je ),
795  $ ldvl )
796  ibeg = 1
797  ELSE
798  CALL dlacpy( ' ', n, nw, work( 2*n+1 ), n, vl( 1, ieig ),
799  $ ldvl )
800  ibeg = je
801  END IF
802 *
803 * Scale eigenvector
804 *
805  xmax = zero
806  IF( ilcplx ) THEN
807  DO 180 j = ibeg, n
808  xmax = max( xmax, abs( vl( j, ieig ) )+
809  $ abs( vl( j, ieig+1 ) ) )
810  180 CONTINUE
811  ELSE
812  DO 190 j = ibeg, n
813  xmax = max( xmax, abs( vl( j, ieig ) ) )
814  190 CONTINUE
815  END IF
816 *
817  IF( xmax.GT.safmin ) THEN
818  xscale = one / xmax
819 *
820  DO 210 jw = 0, nw - 1
821  DO 200 jr = ibeg, n
822  vl( jr, ieig+jw ) = xscale*vl( jr, ieig+jw )
823  200 CONTINUE
824  210 CONTINUE
825  END IF
826  ieig = ieig + nw - 1
827 *
828  220 CONTINUE
829  END IF
830 *
831 * Right eigenvectors
832 *
833  IF( compr ) THEN
834  ieig = im + 1
835 *
836 * Main loop over eigenvalues
837 *
838  ilcplx = .false.
839  DO 500 je = n, 1, -1
840 *
841 * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
842 * (b) this would be the second of a complex pair.
843 * Check for complex eigenvalue, so as to be sure of which
844 * entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
845 * or SELECT(JE-1).
846 * If this is a complex pair, the 2-by-2 diagonal block
847 * corresponding to the eigenvalue is in rows/columns JE-1:JE
848 *
849  IF( ilcplx ) THEN
850  ilcplx = .false.
851  GO TO 500
852  END IF
853  nw = 1
854  IF( je.GT.1 ) THEN
855  IF( s( je, je-1 ).NE.zero ) THEN
856  ilcplx = .true.
857  nw = 2
858  END IF
859  END IF
860  IF( ilall ) THEN
861  ilcomp = .true.
862  ELSE IF( ilcplx ) THEN
863  ilcomp = SELECT( je ) .OR. SELECT( je-1 )
864  ELSE
865  ilcomp = SELECT( je )
866  END IF
867  IF( .NOT.ilcomp )
868  $ GO TO 500
869 *
870 * Decide if (a) singular pencil, (b) real eigenvalue, or
871 * (c) complex eigenvalue.
872 *
873  IF( .NOT.ilcplx ) THEN
874  IF( abs( s( je, je ) ).LE.safmin .AND.
875  $ abs( p( je, je ) ).LE.safmin ) THEN
876 *
877 * Singular matrix pencil -- unit eigenvector
878 *
879  ieig = ieig - 1
880  DO 230 jr = 1, n
881  vr( jr, ieig ) = zero
882  230 CONTINUE
883  vr( ieig, ieig ) = one
884  GO TO 500
885  END IF
886  END IF
887 *
888 * Clear vector
889 *
890  DO 250 jw = 0, nw - 1
891  DO 240 jr = 1, n
892  work( ( jw+2 )*n+jr ) = zero
893  240 CONTINUE
894  250 CONTINUE
895 *
896 * Compute coefficients in ( a A - b B ) x = 0
897 * a is ACOEF
898 * b is BCOEFR + i*BCOEFI
899 *
900  IF( .NOT.ilcplx ) THEN
901 *
902 * Real eigenvalue
903 *
904  temp = one / max( abs( s( je, je ) )*ascale,
905  $ abs( p( je, je ) )*bscale, safmin )
906  salfar = ( temp*s( je, je ) )*ascale
907  sbeta = ( temp*p( je, je ) )*bscale
908  acoef = sbeta*ascale
909  bcoefr = salfar*bscale
910  bcoefi = zero
911 *
912 * Scale to avoid underflow
913 *
914  scale = one
915  lsa = abs( sbeta ).GE.safmin .AND. abs( acoef ).LT.small
916  lsb = abs( salfar ).GE.safmin .AND. abs( bcoefr ).LT.
917  $ small
918  IF( lsa )
919  $ scale = ( small / abs( sbeta ) )*min( anorm, big )
920  IF( lsb )
921  $ scale = max( scale, ( small / abs( salfar ) )*
922  $ min( bnorm, big ) )
923  IF( lsa .OR. lsb ) THEN
924  scale = min( scale, one /
925  $ ( safmin*max( one, abs( acoef ),
926  $ abs( bcoefr ) ) ) )
927  IF( lsa ) THEN
928  acoef = ascale*( scale*sbeta )
929  ELSE
930  acoef = scale*acoef
931  END IF
932  IF( lsb ) THEN
933  bcoefr = bscale*( scale*salfar )
934  ELSE
935  bcoefr = scale*bcoefr
936  END IF
937  END IF
938  acoefa = abs( acoef )
939  bcoefa = abs( bcoefr )
940 *
941 * First component is 1
942 *
943  work( 2*n+je ) = one
944  xmax = one
945 *
946 * Compute contribution from column JE of A and B to sum
947 * (See "Further Details", above.)
948 *
949  DO 260 jr = 1, je - 1
950  work( 2*n+jr ) = bcoefr*p( jr, je ) -
951  $ acoef*s( jr, je )
952  260 CONTINUE
953  ELSE
954 *
955 * Complex eigenvalue
956 *
957  CALL dlag2( s( je-1, je-1 ), lds, p( je-1, je-1 ), ldp,
958  $ safmin*safety, acoef, temp, bcoefr, temp2,
959  $ bcoefi )
960  IF( bcoefi.EQ.zero ) THEN
961  info = je - 1
962  RETURN
963  END IF
964 *
965 * Scale to avoid over/underflow
966 *
967  acoefa = abs( acoef )
968  bcoefa = abs( bcoefr ) + abs( bcoefi )
969  scale = one
970  IF( acoefa*ulp.LT.safmin .AND. acoefa.GE.safmin )
971  $ scale = ( safmin / ulp ) / acoefa
972  IF( bcoefa*ulp.LT.safmin .AND. bcoefa.GE.safmin )
973  $ scale = max( scale, ( safmin / ulp ) / bcoefa )
974  IF( safmin*acoefa.GT.ascale )
975  $ scale = ascale / ( safmin*acoefa )
976  IF( safmin*bcoefa.GT.bscale )
977  $ scale = min( scale, bscale / ( safmin*bcoefa ) )
978  IF( scale.NE.one ) THEN
979  acoef = scale*acoef
980  acoefa = abs( acoef )
981  bcoefr = scale*bcoefr
982  bcoefi = scale*bcoefi
983  bcoefa = abs( bcoefr ) + abs( bcoefi )
984  END IF
985 *
986 * Compute first two components of eigenvector
987 * and contribution to sums
988 *
989  temp = acoef*s( je, je-1 )
990  temp2r = acoef*s( je, je ) - bcoefr*p( je, je )
991  temp2i = -bcoefi*p( je, je )
992  IF( abs( temp ).GE.abs( temp2r )+abs( temp2i ) ) THEN
993  work( 2*n+je ) = one
994  work( 3*n+je ) = zero
995  work( 2*n+je-1 ) = -temp2r / temp
996  work( 3*n+je-1 ) = -temp2i / temp
997  ELSE
998  work( 2*n+je-1 ) = one
999  work( 3*n+je-1 ) = zero
1000  temp = acoef*s( je-1, je )
1001  work( 2*n+je ) = ( bcoefr*p( je-1, je-1 )-acoef*
1002  $ s( je-1, je-1 ) ) / temp
1003  work( 3*n+je ) = bcoefi*p( je-1, je-1 ) / temp
1004  END IF
1005 *
1006  xmax = max( abs( work( 2*n+je ) )+abs( work( 3*n+je ) ),
1007  $ abs( work( 2*n+je-1 ) )+abs( work( 3*n+je-1 ) ) )
1008 *
1009 * Compute contribution from columns JE and JE-1
1010 * of A and B to the sums.
1011 *
1012  creala = acoef*work( 2*n+je-1 )
1013  cimaga = acoef*work( 3*n+je-1 )
1014  crealb = bcoefr*work( 2*n+je-1 ) -
1015  $ bcoefi*work( 3*n+je-1 )
1016  cimagb = bcoefi*work( 2*n+je-1 ) +
1017  $ bcoefr*work( 3*n+je-1 )
1018  cre2a = acoef*work( 2*n+je )
1019  cim2a = acoef*work( 3*n+je )
1020  cre2b = bcoefr*work( 2*n+je ) - bcoefi*work( 3*n+je )
1021  cim2b = bcoefi*work( 2*n+je ) + bcoefr*work( 3*n+je )
1022  DO 270 jr = 1, je - 2
1023  work( 2*n+jr ) = -creala*s( jr, je-1 ) +
1024  $ crealb*p( jr, je-1 ) -
1025  $ cre2a*s( jr, je ) + cre2b*p( jr, je )
1026  work( 3*n+jr ) = -cimaga*s( jr, je-1 ) +
1027  $ cimagb*p( jr, je-1 ) -
1028  $ cim2a*s( jr, je ) + cim2b*p( jr, je )
1029  270 CONTINUE
1030  END IF
1031 *
1032  dmin = max( ulp*acoefa*anorm, ulp*bcoefa*bnorm, safmin )
1033 *
1034 * Columnwise triangular solve of (a A - b B) x = 0
1035 *
1036  il2by2 = .false.
1037  DO 370 j = je - nw, 1, -1
1038 *
1039 * If a 2-by-2 block, is in position j-1:j, wait until
1040 * next iteration to process it (when it will be j:j+1)
1041 *
1042  IF( .NOT.il2by2 .AND. j.GT.1 ) THEN
1043  IF( s( j, j-1 ).NE.zero ) THEN
1044  il2by2 = .true.
1045  GO TO 370
1046  END IF
1047  END IF
1048  bdiag( 1 ) = p( j, j )
1049  IF( il2by2 ) THEN
1050  na = 2
1051  bdiag( 2 ) = p( j+1, j+1 )
1052  ELSE
1053  na = 1
1054  END IF
1055 *
1056 * Compute x(j) (and x(j+1), if 2-by-2 block)
1057 *
1058  CALL dlaln2( .false., na, nw, dmin, acoef, s( j, j ),
1059  $ lds, bdiag( 1 ), bdiag( 2 ), work( 2*n+j ),
1060  $ n, bcoefr, bcoefi, sum, 2, scale, temp,
1061  $ iinfo )
1062  IF( scale.LT.one ) THEN
1063 *
1064  DO 290 jw = 0, nw - 1
1065  DO 280 jr = 1, je
1066  work( ( jw+2 )*n+jr ) = scale*
1067  $ work( ( jw+2 )*n+jr )
1068  280 CONTINUE
1069  290 CONTINUE
1070  END IF
1071  xmax = max( scale*xmax, temp )
1072 *
1073  DO 310 jw = 1, nw
1074  DO 300 ja = 1, na
1075  work( ( jw+1 )*n+j+ja-1 ) = sum( ja, jw )
1076  300 CONTINUE
1077  310 CONTINUE
1078 *
1079 * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1080 *
1081  IF( j.GT.1 ) THEN
1082 *
1083 * Check whether scaling is necessary for sum.
1084 *
1085  xscale = one / max( one, xmax )
1086  temp = acoefa*work( j ) + bcoefa*work( n+j )
1087  IF( il2by2 )
1088  $ temp = max( temp, acoefa*work( j+1 )+bcoefa*
1089  $ work( n+j+1 ) )
1090  temp = max( temp, acoefa, bcoefa )
1091  IF( temp.GT.bignum*xscale ) THEN
1092 *
1093  DO 330 jw = 0, nw - 1
1094  DO 320 jr = 1, je
1095  work( ( jw+2 )*n+jr ) = xscale*
1096  $ work( ( jw+2 )*n+jr )
1097  320 CONTINUE
1098  330 CONTINUE
1099  xmax = xmax*xscale
1100  END IF
1101 *
1102 * Compute the contributions of the off-diagonals of
1103 * column j (and j+1, if 2-by-2 block) of A and B to the
1104 * sums.
1105 *
1106 *
1107  DO 360 ja = 1, na
1108  IF( ilcplx ) THEN
1109  creala = acoef*work( 2*n+j+ja-1 )
1110  cimaga = acoef*work( 3*n+j+ja-1 )
1111  crealb = bcoefr*work( 2*n+j+ja-1 ) -
1112  $ bcoefi*work( 3*n+j+ja-1 )
1113  cimagb = bcoefi*work( 2*n+j+ja-1 ) +
1114  $ bcoefr*work( 3*n+j+ja-1 )
1115  DO 340 jr = 1, j - 1
1116  work( 2*n+jr ) = work( 2*n+jr ) -
1117  $ creala*s( jr, j+ja-1 ) +
1118  $ crealb*p( jr, j+ja-1 )
1119  work( 3*n+jr ) = work( 3*n+jr ) -
1120  $ cimaga*s( jr, j+ja-1 ) +
1121  $ cimagb*p( jr, j+ja-1 )
1122  340 CONTINUE
1123  ELSE
1124  creala = acoef*work( 2*n+j+ja-1 )
1125  crealb = bcoefr*work( 2*n+j+ja-1 )
1126  DO 350 jr = 1, j - 1
1127  work( 2*n+jr ) = work( 2*n+jr ) -
1128  $ creala*s( jr, j+ja-1 ) +
1129  $ crealb*p( jr, j+ja-1 )
1130  350 CONTINUE
1131  END IF
1132  360 CONTINUE
1133  END IF
1134 *
1135  il2by2 = .false.
1136  370 CONTINUE
1137 *
1138 * Copy eigenvector to VR, back transforming if
1139 * HOWMNY='B'.
1140 *
1141  ieig = ieig - nw
1142  IF( ilback ) THEN
1143 *
1144  DO 410 jw = 0, nw - 1
1145  DO 380 jr = 1, n
1146  work( ( jw+4 )*n+jr ) = work( ( jw+2 )*n+1 )*
1147  $ vr( jr, 1 )
1148  380 CONTINUE
1149 *
1150 * A series of compiler directives to defeat
1151 * vectorization for the next loop
1152 *
1153 *
1154  DO 400 jc = 2, je
1155  DO 390 jr = 1, n
1156  work( ( jw+4 )*n+jr ) = work( ( jw+4 )*n+jr ) +
1157  $ work( ( jw+2 )*n+jc )*vr( jr, jc )
1158  390 CONTINUE
1159  400 CONTINUE
1160  410 CONTINUE
1161 *
1162  DO 430 jw = 0, nw - 1
1163  DO 420 jr = 1, n
1164  vr( jr, ieig+jw ) = work( ( jw+4 )*n+jr )
1165  420 CONTINUE
1166  430 CONTINUE
1167 *
1168  iend = n
1169  ELSE
1170  DO 450 jw = 0, nw - 1
1171  DO 440 jr = 1, n
1172  vr( jr, ieig+jw ) = work( ( jw+2 )*n+jr )
1173  440 CONTINUE
1174  450 CONTINUE
1175 *
1176  iend = je
1177  END IF
1178 *
1179 * Scale eigenvector
1180 *
1181  xmax = zero
1182  IF( ilcplx ) THEN
1183  DO 460 j = 1, iend
1184  xmax = max( xmax, abs( vr( j, ieig ) )+
1185  $ abs( vr( j, ieig+1 ) ) )
1186  460 CONTINUE
1187  ELSE
1188  DO 470 j = 1, iend
1189  xmax = max( xmax, abs( vr( j, ieig ) ) )
1190  470 CONTINUE
1191  END IF
1192 *
1193  IF( xmax.GT.safmin ) THEN
1194  xscale = one / xmax
1195  DO 490 jw = 0, nw - 1
1196  DO 480 jr = 1, iend
1197  vr( jr, ieig+jw ) = xscale*vr( jr, ieig+jw )
1198  480 CONTINUE
1199  490 CONTINUE
1200  END IF
1201  500 CONTINUE
1202  END IF
1203 *
1204  RETURN
1205 *
1206 * End of DTGEVC
1207 *
1208  END
xerbla
subroutine xerbla(SRNAME, INFO)
XERBLA
Definition: xerbla.f:60
dlacpy
subroutine dlacpy(UPLO, M, N, A, LDA, B, LDB)
DLACPY copies all or part of one two-dimensional array to another.
Definition: dlacpy.f:103
dgemv
subroutine dgemv(TRANS, M, N, ALPHA, A, LDA, X, INCX, BETA, Y, INCY)
DGEMV
Definition: dgemv.f:156
dlag2
subroutine dlag2(A, LDA, B, LDB, SAFMIN, SCALE1, SCALE2, WR1, WR2, WI)
DLAG2 computes the eigenvalues of a 2-by-2 generalized eigenvalue problem, with scaling as necessary ...
Definition: dlag2.f:156
dlaln2
subroutine dlaln2(LTRANS, NA, NW, SMIN, CA, A, LDA, D1, D2, B, LDB, WR, WI, X, LDX, SCALE, XNORM, INFO)
DLALN2 solves a 1-by-1 or 2-by-2 linear system of equations of the specified form.
Definition: dlaln2.f:218
dlabad
subroutine dlabad(SMALL, LARGE)
DLABAD
Definition: dlabad.f:74
dtgevc
subroutine dtgevc(SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL, LDVL, VR, LDVR, MM, M, WORK, INFO)
DTGEVC
Definition: dtgevc.f:295