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.2) --
5: * -- LAPACK is a software package provided by Univ. of Tennessee, --
6: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
7: * November 2006
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: * A series of compiler directives to defeat vectorization
635: * for the next loop
636: *
637: *$PL$ CMCHAR=' '
638: CDIR$ NEXTSCALAR
639: C$DIR SCALAR
640: CDIR$ NEXT SCALAR
641: CVD$L NOVECTOR
642: CDEC$ NOVECTOR
643: CVD$ NOVECTOR
644: *VDIR NOVECTOR
645: *VOCL LOOP,SCALAR
646: CIBM PREFER SCALAR
647: *$PL$ CMCHAR='*'
648: *
649: DO 120 JW = 1, NW
650: *
651: *$PL$ CMCHAR=' '
652: CDIR$ NEXTSCALAR
653: C$DIR SCALAR
654: CDIR$ NEXT SCALAR
655: CVD$L NOVECTOR
656: CDEC$ NOVECTOR
657: CVD$ NOVECTOR
658: *VDIR NOVECTOR
659: *VOCL LOOP,SCALAR
660: CIBM PREFER SCALAR
661: *$PL$ CMCHAR='*'
662: *
663: DO 110 JA = 1, NA
664: SUMS( JA, JW ) = ZERO
665: SUMP( JA, JW ) = ZERO
666: *
667: DO 100 JR = JE, J - 1
668: SUMS( JA, JW ) = SUMS( JA, JW ) +
669: $ S( JR, J+JA-1 )*
670: $ WORK( ( JW+1 )*N+JR )
671: SUMP( JA, JW ) = SUMP( JA, JW ) +
672: $ P( JR, J+JA-1 )*
673: $ WORK( ( JW+1 )*N+JR )
674: 100 CONTINUE
675: 110 CONTINUE
676: 120 CONTINUE
677: *
678: *$PL$ CMCHAR=' '
679: CDIR$ NEXTSCALAR
680: C$DIR SCALAR
681: CDIR$ NEXT SCALAR
682: CVD$L NOVECTOR
683: CDEC$ NOVECTOR
684: CVD$ NOVECTOR
685: *VDIR NOVECTOR
686: *VOCL LOOP,SCALAR
687: CIBM PREFER SCALAR
688: *$PL$ CMCHAR='*'
689: *
690: DO 130 JA = 1, NA
691: IF( ILCPLX ) THEN
692: SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
693: $ BCOEFR*SUMP( JA, 1 ) -
694: $ BCOEFI*SUMP( JA, 2 )
695: SUM( JA, 2 ) = -ACOEF*SUMS( JA, 2 ) +
696: $ BCOEFR*SUMP( JA, 2 ) +
697: $ BCOEFI*SUMP( JA, 1 )
698: ELSE
699: SUM( JA, 1 ) = -ACOEF*SUMS( JA, 1 ) +
700: $ BCOEFR*SUMP( JA, 1 )
701: END IF
702: 130 CONTINUE
703: *
704: * T
705: * Solve ( a A - b B ) y = SUM(,)
706: * with scaling and perturbation of the denominator
707: *
708: CALL DLALN2( .TRUE., NA, NW, DMIN, ACOEF, S( J, J ), LDS,
709: $ BDIAG( 1 ), BDIAG( 2 ), SUM, 2, BCOEFR,
710: $ BCOEFI, WORK( 2*N+J ), N, SCALE, TEMP,
711: $ IINFO )
712: IF( SCALE.LT.ONE ) THEN
713: DO 150 JW = 0, NW - 1
714: DO 140 JR = JE, J - 1
715: WORK( ( JW+2 )*N+JR ) = SCALE*
716: $ WORK( ( JW+2 )*N+JR )
717: 140 CONTINUE
718: 150 CONTINUE
719: XMAX = SCALE*XMAX
720: END IF
721: XMAX = MAX( XMAX, TEMP )
722: 160 CONTINUE
723: *
724: * Copy eigenvector to VL, back transforming if
725: * HOWMNY='B'.
726: *
727: IEIG = IEIG + 1
728: IF( ILBACK ) THEN
729: DO 170 JW = 0, NW - 1
730: CALL DGEMV( 'N', N, N+1-JE, ONE, VL( 1, JE ), LDVL,
731: $ WORK( ( JW+2 )*N+JE ), 1, ZERO,
732: $ WORK( ( JW+4 )*N+1 ), 1 )
733: 170 CONTINUE
734: CALL DLACPY( ' ', N, NW, WORK( 4*N+1 ), N, VL( 1, JE ),
735: $ LDVL )
736: IBEG = 1
737: ELSE
738: CALL DLACPY( ' ', N, NW, WORK( 2*N+1 ), N, VL( 1, IEIG ),
739: $ LDVL )
740: IBEG = JE
741: END IF
742: *
743: * Scale eigenvector
744: *
745: XMAX = ZERO
746: IF( ILCPLX ) THEN
747: DO 180 J = IBEG, N
748: XMAX = MAX( XMAX, ABS( VL( J, IEIG ) )+
749: $ ABS( VL( J, IEIG+1 ) ) )
750: 180 CONTINUE
751: ELSE
752: DO 190 J = IBEG, N
753: XMAX = MAX( XMAX, ABS( VL( J, IEIG ) ) )
754: 190 CONTINUE
755: END IF
756: *
757: IF( XMAX.GT.SAFMIN ) THEN
758: XSCALE = ONE / XMAX
759: *
760: DO 210 JW = 0, NW - 1
761: DO 200 JR = IBEG, N
762: VL( JR, IEIG+JW ) = XSCALE*VL( JR, IEIG+JW )
763: 200 CONTINUE
764: 210 CONTINUE
765: END IF
766: IEIG = IEIG + NW - 1
767: *
768: 220 CONTINUE
769: END IF
770: *
771: * Right eigenvectors
772: *
773: IF( COMPR ) THEN
774: IEIG = IM + 1
775: *
776: * Main loop over eigenvalues
777: *
778: ILCPLX = .FALSE.
779: DO 500 JE = N, 1, -1
780: *
781: * Skip this iteration if (a) HOWMNY='S' and SELECT=.FALSE., or
782: * (b) this would be the second of a complex pair.
783: * Check for complex eigenvalue, so as to be sure of which
784: * entry(-ies) of SELECT to look at -- if complex, SELECT(JE)
785: * or SELECT(JE-1).
786: * If this is a complex pair, the 2-by-2 diagonal block
787: * corresponding to the eigenvalue is in rows/columns JE-1:JE
788: *
789: IF( ILCPLX ) THEN
790: ILCPLX = .FALSE.
791: GO TO 500
792: END IF
793: NW = 1
794: IF( JE.GT.1 ) THEN
795: IF( S( JE, JE-1 ).NE.ZERO ) THEN
796: ILCPLX = .TRUE.
797: NW = 2
798: END IF
799: END IF
800: IF( ILALL ) THEN
801: ILCOMP = .TRUE.
802: ELSE IF( ILCPLX ) THEN
803: ILCOMP = SELECT( JE ) .OR. SELECT( JE-1 )
804: ELSE
805: ILCOMP = SELECT( JE )
806: END IF
807: IF( .NOT.ILCOMP )
808: $ GO TO 500
809: *
810: * Decide if (a) singular pencil, (b) real eigenvalue, or
811: * (c) complex eigenvalue.
812: *
813: IF( .NOT.ILCPLX ) THEN
814: IF( ABS( S( JE, JE ) ).LE.SAFMIN .AND.
815: $ ABS( P( JE, JE ) ).LE.SAFMIN ) THEN
816: *
817: * Singular matrix pencil -- unit eigenvector
818: *
819: IEIG = IEIG - 1
820: DO 230 JR = 1, N
821: VR( JR, IEIG ) = ZERO
822: 230 CONTINUE
823: VR( IEIG, IEIG ) = ONE
824: GO TO 500
825: END IF
826: END IF
827: *
828: * Clear vector
829: *
830: DO 250 JW = 0, NW - 1
831: DO 240 JR = 1, N
832: WORK( ( JW+2 )*N+JR ) = ZERO
833: 240 CONTINUE
834: 250 CONTINUE
835: *
836: * Compute coefficients in ( a A - b B ) x = 0
837: * a is ACOEF
838: * b is BCOEFR + i*BCOEFI
839: *
840: IF( .NOT.ILCPLX ) THEN
841: *
842: * Real eigenvalue
843: *
844: TEMP = ONE / MAX( ABS( S( JE, JE ) )*ASCALE,
845: $ ABS( P( JE, JE ) )*BSCALE, SAFMIN )
846: SALFAR = ( TEMP*S( JE, JE ) )*ASCALE
847: SBETA = ( TEMP*P( JE, JE ) )*BSCALE
848: ACOEF = SBETA*ASCALE
849: BCOEFR = SALFAR*BSCALE
850: BCOEFI = ZERO
851: *
852: * Scale to avoid underflow
853: *
854: SCALE = ONE
855: LSA = ABS( SBETA ).GE.SAFMIN .AND. ABS( ACOEF ).LT.SMALL
856: LSB = ABS( SALFAR ).GE.SAFMIN .AND. ABS( BCOEFR ).LT.
857: $ SMALL
858: IF( LSA )
859: $ SCALE = ( SMALL / ABS( SBETA ) )*MIN( ANORM, BIG )
860: IF( LSB )
861: $ SCALE = MAX( SCALE, ( SMALL / ABS( SALFAR ) )*
862: $ MIN( BNORM, BIG ) )
863: IF( LSA .OR. LSB ) THEN
864: SCALE = MIN( SCALE, ONE /
865: $ ( SAFMIN*MAX( ONE, ABS( ACOEF ),
866: $ ABS( BCOEFR ) ) ) )
867: IF( LSA ) THEN
868: ACOEF = ASCALE*( SCALE*SBETA )
869: ELSE
870: ACOEF = SCALE*ACOEF
871: END IF
872: IF( LSB ) THEN
873: BCOEFR = BSCALE*( SCALE*SALFAR )
874: ELSE
875: BCOEFR = SCALE*BCOEFR
876: END IF
877: END IF
878: ACOEFA = ABS( ACOEF )
879: BCOEFA = ABS( BCOEFR )
880: *
881: * First component is 1
882: *
883: WORK( 2*N+JE ) = ONE
884: XMAX = ONE
885: *
886: * Compute contribution from column JE of A and B to sum
887: * (See "Further Details", above.)
888: *
889: DO 260 JR = 1, JE - 1
890: WORK( 2*N+JR ) = BCOEFR*P( JR, JE ) -
891: $ ACOEF*S( JR, JE )
892: 260 CONTINUE
893: ELSE
894: *
895: * Complex eigenvalue
896: *
897: CALL DLAG2( S( JE-1, JE-1 ), LDS, P( JE-1, JE-1 ), LDP,
898: $ SAFMIN*SAFETY, ACOEF, TEMP, BCOEFR, TEMP2,
899: $ BCOEFI )
900: IF( BCOEFI.EQ.ZERO ) THEN
901: INFO = JE - 1
902: RETURN
903: END IF
904: *
905: * Scale to avoid over/underflow
906: *
907: ACOEFA = ABS( ACOEF )
908: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
909: SCALE = ONE
910: IF( ACOEFA*ULP.LT.SAFMIN .AND. ACOEFA.GE.SAFMIN )
911: $ SCALE = ( SAFMIN / ULP ) / ACOEFA
912: IF( BCOEFA*ULP.LT.SAFMIN .AND. BCOEFA.GE.SAFMIN )
913: $ SCALE = MAX( SCALE, ( SAFMIN / ULP ) / BCOEFA )
914: IF( SAFMIN*ACOEFA.GT.ASCALE )
915: $ SCALE = ASCALE / ( SAFMIN*ACOEFA )
916: IF( SAFMIN*BCOEFA.GT.BSCALE )
917: $ SCALE = MIN( SCALE, BSCALE / ( SAFMIN*BCOEFA ) )
918: IF( SCALE.NE.ONE ) THEN
919: ACOEF = SCALE*ACOEF
920: ACOEFA = ABS( ACOEF )
921: BCOEFR = SCALE*BCOEFR
922: BCOEFI = SCALE*BCOEFI
923: BCOEFA = ABS( BCOEFR ) + ABS( BCOEFI )
924: END IF
925: *
926: * Compute first two components of eigenvector
927: * and contribution to sums
928: *
929: TEMP = ACOEF*S( JE, JE-1 )
930: TEMP2R = ACOEF*S( JE, JE ) - BCOEFR*P( JE, JE )
931: TEMP2I = -BCOEFI*P( JE, JE )
932: IF( ABS( TEMP ).GE.ABS( TEMP2R )+ABS( TEMP2I ) ) THEN
933: WORK( 2*N+JE ) = ONE
934: WORK( 3*N+JE ) = ZERO
935: WORK( 2*N+JE-1 ) = -TEMP2R / TEMP
936: WORK( 3*N+JE-1 ) = -TEMP2I / TEMP
937: ELSE
938: WORK( 2*N+JE-1 ) = ONE
939: WORK( 3*N+JE-1 ) = ZERO
940: TEMP = ACOEF*S( JE-1, JE )
941: WORK( 2*N+JE ) = ( BCOEFR*P( JE-1, JE-1 )-ACOEF*
942: $ S( JE-1, JE-1 ) ) / TEMP
943: WORK( 3*N+JE ) = BCOEFI*P( JE-1, JE-1 ) / TEMP
944: END IF
945: *
946: XMAX = MAX( ABS( WORK( 2*N+JE ) )+ABS( WORK( 3*N+JE ) ),
947: $ ABS( WORK( 2*N+JE-1 ) )+ABS( WORK( 3*N+JE-1 ) ) )
948: *
949: * Compute contribution from columns JE and JE-1
950: * of A and B to the sums.
951: *
952: CREALA = ACOEF*WORK( 2*N+JE-1 )
953: CIMAGA = ACOEF*WORK( 3*N+JE-1 )
954: CREALB = BCOEFR*WORK( 2*N+JE-1 ) -
955: $ BCOEFI*WORK( 3*N+JE-1 )
956: CIMAGB = BCOEFI*WORK( 2*N+JE-1 ) +
957: $ BCOEFR*WORK( 3*N+JE-1 )
958: CRE2A = ACOEF*WORK( 2*N+JE )
959: CIM2A = ACOEF*WORK( 3*N+JE )
960: CRE2B = BCOEFR*WORK( 2*N+JE ) - BCOEFI*WORK( 3*N+JE )
961: CIM2B = BCOEFI*WORK( 2*N+JE ) + BCOEFR*WORK( 3*N+JE )
962: DO 270 JR = 1, JE - 2
963: WORK( 2*N+JR ) = -CREALA*S( JR, JE-1 ) +
964: $ CREALB*P( JR, JE-1 ) -
965: $ CRE2A*S( JR, JE ) + CRE2B*P( JR, JE )
966: WORK( 3*N+JR ) = -CIMAGA*S( JR, JE-1 ) +
967: $ CIMAGB*P( JR, JE-1 ) -
968: $ CIM2A*S( JR, JE ) + CIM2B*P( JR, JE )
969: 270 CONTINUE
970: END IF
971: *
972: DMIN = MAX( ULP*ACOEFA*ANORM, ULP*BCOEFA*BNORM, SAFMIN )
973: *
974: * Columnwise triangular solve of (a A - b B) x = 0
975: *
976: IL2BY2 = .FALSE.
977: DO 370 J = JE - NW, 1, -1
978: *
979: * If a 2-by-2 block, is in position j-1:j, wait until
980: * next iteration to process it (when it will be j:j+1)
981: *
982: IF( .NOT.IL2BY2 .AND. J.GT.1 ) THEN
983: IF( S( J, J-1 ).NE.ZERO ) THEN
984: IL2BY2 = .TRUE.
985: GO TO 370
986: END IF
987: END IF
988: BDIAG( 1 ) = P( J, J )
989: IF( IL2BY2 ) THEN
990: NA = 2
991: BDIAG( 2 ) = P( J+1, J+1 )
992: ELSE
993: NA = 1
994: END IF
995: *
996: * Compute x(j) (and x(j+1), if 2-by-2 block)
997: *
998: CALL DLALN2( .FALSE., NA, NW, DMIN, ACOEF, S( J, J ),
999: $ LDS, BDIAG( 1 ), BDIAG( 2 ), WORK( 2*N+J ),
1000: $ N, BCOEFR, BCOEFI, SUM, 2, SCALE, TEMP,
1001: $ IINFO )
1002: IF( SCALE.LT.ONE ) THEN
1003: *
1004: DO 290 JW = 0, NW - 1
1005: DO 280 JR = 1, JE
1006: WORK( ( JW+2 )*N+JR ) = SCALE*
1007: $ WORK( ( JW+2 )*N+JR )
1008: 280 CONTINUE
1009: 290 CONTINUE
1010: END IF
1011: XMAX = MAX( SCALE*XMAX, TEMP )
1012: *
1013: DO 310 JW = 1, NW
1014: DO 300 JA = 1, NA
1015: WORK( ( JW+1 )*N+J+JA-1 ) = SUM( JA, JW )
1016: 300 CONTINUE
1017: 310 CONTINUE
1018: *
1019: * w = w + x(j)*(a S(*,j) - b P(*,j) ) with scaling
1020: *
1021: IF( J.GT.1 ) THEN
1022: *
1023: * Check whether scaling is necessary for sum.
1024: *
1025: XSCALE = ONE / MAX( ONE, XMAX )
1026: TEMP = ACOEFA*WORK( J ) + BCOEFA*WORK( N+J )
1027: IF( IL2BY2 )
1028: $ TEMP = MAX( TEMP, ACOEFA*WORK( J+1 )+BCOEFA*
1029: $ WORK( N+J+1 ) )
1030: TEMP = MAX( TEMP, ACOEFA, BCOEFA )
1031: IF( TEMP.GT.BIGNUM*XSCALE ) THEN
1032: *
1033: DO 330 JW = 0, NW - 1
1034: DO 320 JR = 1, JE
1035: WORK( ( JW+2 )*N+JR ) = XSCALE*
1036: $ WORK( ( JW+2 )*N+JR )
1037: 320 CONTINUE
1038: 330 CONTINUE
1039: XMAX = XMAX*XSCALE
1040: END IF
1041: *
1042: * Compute the contributions of the off-diagonals of
1043: * column j (and j+1, if 2-by-2 block) of A and B to the
1044: * sums.
1045: *
1046: *
1047: DO 360 JA = 1, NA
1048: IF( ILCPLX ) THEN
1049: CREALA = ACOEF*WORK( 2*N+J+JA-1 )
1050: CIMAGA = ACOEF*WORK( 3*N+J+JA-1 )
1051: CREALB = BCOEFR*WORK( 2*N+J+JA-1 ) -
1052: $ BCOEFI*WORK( 3*N+J+JA-1 )
1053: CIMAGB = BCOEFI*WORK( 2*N+J+JA-1 ) +
1054: $ BCOEFR*WORK( 3*N+J+JA-1 )
1055: DO 340 JR = 1, J - 1
1056: WORK( 2*N+JR ) = WORK( 2*N+JR ) -
1057: $ CREALA*S( JR, J+JA-1 ) +
1058: $ CREALB*P( JR, J+JA-1 )
1059: WORK( 3*N+JR ) = WORK( 3*N+JR ) -
1060: $ CIMAGA*S( JR, J+JA-1 ) +
1061: $ CIMAGB*P( JR, J+JA-1 )
1062: 340 CONTINUE
1063: ELSE
1064: CREALA = ACOEF*WORK( 2*N+J+JA-1 )
1065: CREALB = BCOEFR*WORK( 2*N+J+JA-1 )
1066: DO 350 JR = 1, J - 1
1067: WORK( 2*N+JR ) = WORK( 2*N+JR ) -
1068: $ CREALA*S( JR, J+JA-1 ) +
1069: $ CREALB*P( JR, J+JA-1 )
1070: 350 CONTINUE
1071: END IF
1072: 360 CONTINUE
1073: END IF
1074: *
1075: IL2BY2 = .FALSE.
1076: 370 CONTINUE
1077: *
1078: * Copy eigenvector to VR, back transforming if
1079: * HOWMNY='B'.
1080: *
1081: IEIG = IEIG - NW
1082: IF( ILBACK ) THEN
1083: *
1084: DO 410 JW = 0, NW - 1
1085: DO 380 JR = 1, N
1086: WORK( ( JW+4 )*N+JR ) = WORK( ( JW+2 )*N+1 )*
1087: $ VR( JR, 1 )
1088: 380 CONTINUE
1089: *
1090: * A series of compiler directives to defeat
1091: * vectorization for the next loop
1092: *
1093: *
1094: DO 400 JC = 2, JE
1095: DO 390 JR = 1, N
1096: WORK( ( JW+4 )*N+JR ) = WORK( ( JW+4 )*N+JR ) +
1097: $ WORK( ( JW+2 )*N+JC )*VR( JR, JC )
1098: 390 CONTINUE
1099: 400 CONTINUE
1100: 410 CONTINUE
1101: *
1102: DO 430 JW = 0, NW - 1
1103: DO 420 JR = 1, N
1104: VR( JR, IEIG+JW ) = WORK( ( JW+4 )*N+JR )
1105: 420 CONTINUE
1106: 430 CONTINUE
1107: *
1108: IEND = N
1109: ELSE
1110: DO 450 JW = 0, NW - 1
1111: DO 440 JR = 1, N
1112: VR( JR, IEIG+JW ) = WORK( ( JW+2 )*N+JR )
1113: 440 CONTINUE
1114: 450 CONTINUE
1115: *
1116: IEND = JE
1117: END IF
1118: *
1119: * Scale eigenvector
1120: *
1121: XMAX = ZERO
1122: IF( ILCPLX ) THEN
1123: DO 460 J = 1, IEND
1124: XMAX = MAX( XMAX, ABS( VR( J, IEIG ) )+
1125: $ ABS( VR( J, IEIG+1 ) ) )
1126: 460 CONTINUE
1127: ELSE
1128: DO 470 J = 1, IEND
1129: XMAX = MAX( XMAX, ABS( VR( J, IEIG ) ) )
1130: 470 CONTINUE
1131: END IF
1132: *
1133: IF( XMAX.GT.SAFMIN ) THEN
1134: XSCALE = ONE / XMAX
1135: DO 490 JW = 0, NW - 1
1136: DO 480 JR = 1, IEND
1137: VR( JR, IEIG+JW ) = XSCALE*VR( JR, IEIG+JW )
1138: 480 CONTINUE
1139: 490 CONTINUE
1140: END IF
1141: 500 CONTINUE
1142: END IF
1143: *
1144: RETURN
1145: *
1146: * End of DTGEVC
1147: *
1148: END
CVSweb interface <joel.bertrand@systella.fr>
![[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