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