![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dtgevc.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack|
Return to dtgevc.f CVS log | Up to [local] / rpl / lapack / lapack |
1.1 bertrand 1: SUBROUTINE DTGEVC( SIDE, HOWMNY, SELECT, N, S, LDS, P, LDP, VL,
2: $ LDVL, VR, LDVR, MM, M, WORK, INFO )
3: *
1.8 ! bertrand 4: * -- LAPACK routine (version 3.3.1) --
1.1 bertrand 5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.8 ! bertrand 7: * -- April 2011 --
1.1 bertrand 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