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