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