Annotation of rpl/lapack/lapack/dgeesx.f, revision 1.9
1.9 ! bertrand 1: *> \brief <b> DGEESX computes the eigenvalues, the Schur form, and, optionally, the matrix of Schur vectors for GE matrices</b>
! 2: *
! 3: * =========== DOCUMENTATION ===========
! 4: *
! 5: * Online html documentation available at
! 6: * http://www.netlib.org/lapack/explore-html/
! 7: *
! 8: *> \htmlonly
! 9: *> Download DGEESX + dependencies
! 10: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeesx.f">
! 11: *> [TGZ]</a>
! 12: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeesx.f">
! 13: *> [ZIP]</a>
! 14: *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeesx.f">
! 15: *> [TXT]</a>
! 16: *> \endhtmlonly
! 17: *
! 18: * Definition:
! 19: * ===========
! 20: *
! 21: * SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
! 22: * WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
! 23: * IWORK, LIWORK, BWORK, INFO )
! 24: *
! 25: * .. Scalar Arguments ..
! 26: * CHARACTER JOBVS, SENSE, SORT
! 27: * INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
! 28: * DOUBLE PRECISION RCONDE, RCONDV
! 29: * ..
! 30: * .. Array Arguments ..
! 31: * LOGICAL BWORK( * )
! 32: * INTEGER IWORK( * )
! 33: * DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
! 34: * $ WR( * )
! 35: * ..
! 36: * .. Function Arguments ..
! 37: * LOGICAL SELECT
! 38: * EXTERNAL SELECT
! 39: * ..
! 40: *
! 41: *
! 42: *> \par Purpose:
! 43: * =============
! 44: *>
! 45: *> \verbatim
! 46: *>
! 47: *> DGEESX computes for an N-by-N real nonsymmetric matrix A, the
! 48: *> eigenvalues, the real Schur form T, and, optionally, the matrix of
! 49: *> Schur vectors Z. This gives the Schur factorization A = Z*T*(Z**T).
! 50: *>
! 51: *> Optionally, it also orders the eigenvalues on the diagonal of the
! 52: *> real Schur form so that selected eigenvalues are at the top left;
! 53: *> computes a reciprocal condition number for the average of the
! 54: *> selected eigenvalues (RCONDE); and computes a reciprocal condition
! 55: *> number for the right invariant subspace corresponding to the
! 56: *> selected eigenvalues (RCONDV). The leading columns of Z form an
! 57: *> orthonormal basis for this invariant subspace.
! 58: *>
! 59: *> For further explanation of the reciprocal condition numbers RCONDE
! 60: *> and RCONDV, see Section 4.10 of the LAPACK Users' Guide (where
! 61: *> these quantities are called s and sep respectively).
! 62: *>
! 63: *> A real matrix is in real Schur form if it is upper quasi-triangular
! 64: *> with 1-by-1 and 2-by-2 blocks. 2-by-2 blocks will be standardized in
! 65: *> the form
! 66: *> [ a b ]
! 67: *> [ c a ]
! 68: *>
! 69: *> where b*c < 0. The eigenvalues of such a block are a +- sqrt(bc).
! 70: *> \endverbatim
! 71: *
! 72: * Arguments:
! 73: * ==========
! 74: *
! 75: *> \param[in] JOBVS
! 76: *> \verbatim
! 77: *> JOBVS is CHARACTER*1
! 78: *> = 'N': Schur vectors are not computed;
! 79: *> = 'V': Schur vectors are computed.
! 80: *> \endverbatim
! 81: *>
! 82: *> \param[in] SORT
! 83: *> \verbatim
! 84: *> SORT is CHARACTER*1
! 85: *> Specifies whether or not to order the eigenvalues on the
! 86: *> diagonal of the Schur form.
! 87: *> = 'N': Eigenvalues are not ordered;
! 88: *> = 'S': Eigenvalues are ordered (see SELECT).
! 89: *> \endverbatim
! 90: *>
! 91: *> \param[in] SELECT
! 92: *> \verbatim
! 93: *> SELECT is procedure) LOGICAL FUNCTION of two DOUBLE PRECISION arguments
! 94: *> SELECT must be declared EXTERNAL in the calling subroutine.
! 95: *> If SORT = 'S', SELECT is used to select eigenvalues to sort
! 96: *> to the top left of the Schur form.
! 97: *> If SORT = 'N', SELECT is not referenced.
! 98: *> An eigenvalue WR(j)+sqrt(-1)*WI(j) is selected if
! 99: *> SELECT(WR(j),WI(j)) is true; i.e., if either one of a
! 100: *> complex conjugate pair of eigenvalues is selected, then both
! 101: *> are. Note that a selected complex eigenvalue may no longer
! 102: *> satisfy SELECT(WR(j),WI(j)) = .TRUE. after ordering, since
! 103: *> ordering may change the value of complex eigenvalues
! 104: *> (especially if the eigenvalue is ill-conditioned); in this
! 105: *> case INFO may be set to N+3 (see INFO below).
! 106: *> \endverbatim
! 107: *>
! 108: *> \param[in] SENSE
! 109: *> \verbatim
! 110: *> SENSE is CHARACTER*1
! 111: *> Determines which reciprocal condition numbers are computed.
! 112: *> = 'N': None are computed;
! 113: *> = 'E': Computed for average of selected eigenvalues only;
! 114: *> = 'V': Computed for selected right invariant subspace only;
! 115: *> = 'B': Computed for both.
! 116: *> If SENSE = 'E', 'V' or 'B', SORT must equal 'S'.
! 117: *> \endverbatim
! 118: *>
! 119: *> \param[in] N
! 120: *> \verbatim
! 121: *> N is INTEGER
! 122: *> The order of the matrix A. N >= 0.
! 123: *> \endverbatim
! 124: *>
! 125: *> \param[in,out] A
! 126: *> \verbatim
! 127: *> A is DOUBLE PRECISION array, dimension (LDA, N)
! 128: *> On entry, the N-by-N matrix A.
! 129: *> On exit, A is overwritten by its real Schur form T.
! 130: *> \endverbatim
! 131: *>
! 132: *> \param[in] LDA
! 133: *> \verbatim
! 134: *> LDA is INTEGER
! 135: *> The leading dimension of the array A. LDA >= max(1,N).
! 136: *> \endverbatim
! 137: *>
! 138: *> \param[out] SDIM
! 139: *> \verbatim
! 140: *> SDIM is INTEGER
! 141: *> If SORT = 'N', SDIM = 0.
! 142: *> If SORT = 'S', SDIM = number of eigenvalues (after sorting)
! 143: *> for which SELECT is true. (Complex conjugate
! 144: *> pairs for which SELECT is true for either
! 145: *> eigenvalue count as 2.)
! 146: *> \endverbatim
! 147: *>
! 148: *> \param[out] WR
! 149: *> \verbatim
! 150: *> WR is DOUBLE PRECISION array, dimension (N)
! 151: *> \endverbatim
! 152: *>
! 153: *> \param[out] WI
! 154: *> \verbatim
! 155: *> WI is DOUBLE PRECISION array, dimension (N)
! 156: *> WR and WI contain the real and imaginary parts, respectively,
! 157: *> of the computed eigenvalues, in the same order that they
! 158: *> appear on the diagonal of the output Schur form T. Complex
! 159: *> conjugate pairs of eigenvalues appear consecutively with the
! 160: *> eigenvalue having the positive imaginary part first.
! 161: *> \endverbatim
! 162: *>
! 163: *> \param[out] VS
! 164: *> \verbatim
! 165: *> VS is DOUBLE PRECISION array, dimension (LDVS,N)
! 166: *> If JOBVS = 'V', VS contains the orthogonal matrix Z of Schur
! 167: *> vectors.
! 168: *> If JOBVS = 'N', VS is not referenced.
! 169: *> \endverbatim
! 170: *>
! 171: *> \param[in] LDVS
! 172: *> \verbatim
! 173: *> LDVS is INTEGER
! 174: *> The leading dimension of the array VS. LDVS >= 1, and if
! 175: *> JOBVS = 'V', LDVS >= N.
! 176: *> \endverbatim
! 177: *>
! 178: *> \param[out] RCONDE
! 179: *> \verbatim
! 180: *> RCONDE is DOUBLE PRECISION
! 181: *> If SENSE = 'E' or 'B', RCONDE contains the reciprocal
! 182: *> condition number for the average of the selected eigenvalues.
! 183: *> Not referenced if SENSE = 'N' or 'V'.
! 184: *> \endverbatim
! 185: *>
! 186: *> \param[out] RCONDV
! 187: *> \verbatim
! 188: *> RCONDV is DOUBLE PRECISION
! 189: *> If SENSE = 'V' or 'B', RCONDV contains the reciprocal
! 190: *> condition number for the selected right invariant subspace.
! 191: *> Not referenced if SENSE = 'N' or 'E'.
! 192: *> \endverbatim
! 193: *>
! 194: *> \param[out] WORK
! 195: *> \verbatim
! 196: *> WORK is DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 197: *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 198: *> \endverbatim
! 199: *>
! 200: *> \param[in] LWORK
! 201: *> \verbatim
! 202: *> LWORK is INTEGER
! 203: *> The dimension of the array WORK. LWORK >= max(1,3*N).
! 204: *> Also, if SENSE = 'E' or 'V' or 'B',
! 205: *> LWORK >= N+2*SDIM*(N-SDIM), where SDIM is the number of
! 206: *> selected eigenvalues computed by this routine. Note that
! 207: *> N+2*SDIM*(N-SDIM) <= N+N*N/2. Note also that an error is only
! 208: *> returned if LWORK < max(1,3*N), but if SENSE = 'E' or 'V' or
! 209: *> 'B' this may not be large enough.
! 210: *> For good performance, LWORK must generally be larger.
! 211: *>
! 212: *> If LWORK = -1, then a workspace query is assumed; the routine
! 213: *> only calculates upper bounds on the optimal sizes of the
! 214: *> arrays WORK and IWORK, returns these values as the first
! 215: *> entries of the WORK and IWORK arrays, and no error messages
! 216: *> related to LWORK or LIWORK are issued by XERBLA.
! 217: *> \endverbatim
! 218: *>
! 219: *> \param[out] IWORK
! 220: *> \verbatim
! 221: *> IWORK is INTEGER array, dimension (MAX(1,LIWORK))
! 222: *> On exit, if INFO = 0, IWORK(1) returns the optimal LIWORK.
! 223: *> \endverbatim
! 224: *>
! 225: *> \param[in] LIWORK
! 226: *> \verbatim
! 227: *> LIWORK is INTEGER
! 228: *> The dimension of the array IWORK.
! 229: *> LIWORK >= 1; if SENSE = 'V' or 'B', LIWORK >= SDIM*(N-SDIM).
! 230: *> Note that SDIM*(N-SDIM) <= N*N/4. Note also that an error is
! 231: *> only returned if LIWORK < 1, but if SENSE = 'V' or 'B' this
! 232: *> may not be large enough.
! 233: *>
! 234: *> If LIWORK = -1, then a workspace query is assumed; the
! 235: *> routine only calculates upper bounds on the optimal sizes of
! 236: *> the arrays WORK and IWORK, returns these values as the first
! 237: *> entries of the WORK and IWORK arrays, and no error messages
! 238: *> related to LWORK or LIWORK are issued by XERBLA.
! 239: *> \endverbatim
! 240: *>
! 241: *> \param[out] BWORK
! 242: *> \verbatim
! 243: *> BWORK is LOGICAL array, dimension (N)
! 244: *> Not referenced if SORT = 'N'.
! 245: *> \endverbatim
! 246: *>
! 247: *> \param[out] INFO
! 248: *> \verbatim
! 249: *> INFO is INTEGER
! 250: *> = 0: successful exit
! 251: *> < 0: if INFO = -i, the i-th argument had an illegal value.
! 252: *> > 0: if INFO = i, and i is
! 253: *> <= N: the QR algorithm failed to compute all the
! 254: *> eigenvalues; elements 1:ILO-1 and i+1:N of WR and WI
! 255: *> contain those eigenvalues which have converged; if
! 256: *> JOBVS = 'V', VS contains the transformation which
! 257: *> reduces A to its partially converged Schur form.
! 258: *> = N+1: the eigenvalues could not be reordered because some
! 259: *> eigenvalues were too close to separate (the problem
! 260: *> is very ill-conditioned);
! 261: *> = N+2: after reordering, roundoff changed values of some
! 262: *> complex eigenvalues so that leading eigenvalues in
! 263: *> the Schur form no longer satisfy SELECT=.TRUE. This
! 264: *> could also be caused by underflow due to scaling.
! 265: *> \endverbatim
! 266: *
! 267: * Authors:
! 268: * ========
! 269: *
! 270: *> \author Univ. of Tennessee
! 271: *> \author Univ. of California Berkeley
! 272: *> \author Univ. of Colorado Denver
! 273: *> \author NAG Ltd.
! 274: *
! 275: *> \date November 2011
! 276: *
! 277: *> \ingroup doubleGEeigen
! 278: *
! 279: * =====================================================================
1.1 bertrand 280: SUBROUTINE DGEESX( JOBVS, SORT, SELECT, SENSE, N, A, LDA, SDIM,
281: $ WR, WI, VS, LDVS, RCONDE, RCONDV, WORK, LWORK,
282: $ IWORK, LIWORK, BWORK, INFO )
283: *
1.9 ! bertrand 284: * -- LAPACK driver routine (version 3.4.0) --
1.1 bertrand 285: * -- LAPACK is a software package provided by Univ. of Tennessee, --
286: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
1.9 ! bertrand 287: * November 2011
1.1 bertrand 288: *
289: * .. Scalar Arguments ..
290: CHARACTER JOBVS, SENSE, SORT
291: INTEGER INFO, LDA, LDVS, LIWORK, LWORK, N, SDIM
292: DOUBLE PRECISION RCONDE, RCONDV
293: * ..
294: * .. Array Arguments ..
295: LOGICAL BWORK( * )
296: INTEGER IWORK( * )
297: DOUBLE PRECISION A( LDA, * ), VS( LDVS, * ), WI( * ), WORK( * ),
298: $ WR( * )
299: * ..
300: * .. Function Arguments ..
301: LOGICAL SELECT
302: EXTERNAL SELECT
303: * ..
304: *
305: * =====================================================================
306: *
307: * .. Parameters ..
308: DOUBLE PRECISION ZERO, ONE
309: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
310: * ..
311: * .. Local Scalars ..
312: LOGICAL CURSL, LASTSL, LQUERY, LST2SL, SCALEA, WANTSB,
313: $ WANTSE, WANTSN, WANTST, WANTSV, WANTVS
314: INTEGER HSWORK, I, I1, I2, IBAL, ICOND, IERR, IEVAL,
315: $ IHI, ILO, INXT, IP, ITAU, IWRK, LIWRK, LWRK,
316: $ MAXWRK, MINWRK
317: DOUBLE PRECISION ANRM, BIGNUM, CSCALE, EPS, SMLNUM
318: * ..
319: * .. Local Arrays ..
320: DOUBLE PRECISION DUM( 1 )
321: * ..
322: * .. External Subroutines ..
323: EXTERNAL DCOPY, DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLACPY,
324: $ DLASCL, DORGHR, DSWAP, DTRSEN, XERBLA
325: * ..
326: * .. External Functions ..
327: LOGICAL LSAME
328: INTEGER ILAENV
329: DOUBLE PRECISION DLAMCH, DLANGE
330: EXTERNAL LSAME, ILAENV, DLABAD, DLAMCH, DLANGE
331: * ..
332: * .. Intrinsic Functions ..
333: INTRINSIC MAX, SQRT
334: * ..
335: * .. Executable Statements ..
336: *
337: * Test the input arguments
338: *
339: INFO = 0
340: WANTVS = LSAME( JOBVS, 'V' )
341: WANTST = LSAME( SORT, 'S' )
342: WANTSN = LSAME( SENSE, 'N' )
343: WANTSE = LSAME( SENSE, 'E' )
344: WANTSV = LSAME( SENSE, 'V' )
345: WANTSB = LSAME( SENSE, 'B' )
346: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
1.5 bertrand 347: *
1.1 bertrand 348: IF( ( .NOT.WANTVS ) .AND. ( .NOT.LSAME( JOBVS, 'N' ) ) ) THEN
349: INFO = -1
350: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
351: INFO = -2
352: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
353: $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
354: INFO = -4
355: ELSE IF( N.LT.0 ) THEN
356: INFO = -5
357: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
358: INFO = -7
359: ELSE IF( LDVS.LT.1 .OR. ( WANTVS .AND. LDVS.LT.N ) ) THEN
360: INFO = -12
361: END IF
362: *
363: * Compute workspace
364: * (Note: Comments in the code beginning "RWorkspace:" describe the
365: * minimal amount of real workspace needed at that point in the
366: * code, as well as the preferred amount for good performance.
367: * IWorkspace refers to integer workspace.
368: * NB refers to the optimal block size for the immediately
369: * following subroutine, as returned by ILAENV.
370: * HSWORK refers to the workspace preferred by DHSEQR, as
371: * calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
372: * the worst case.
373: * If SENSE = 'E', 'V' or 'B', then the amount of workspace needed
374: * depends on SDIM, which is computed by the routine DTRSEN later
375: * in the code.)
376: *
377: IF( INFO.EQ.0 ) THEN
378: LIWRK = 1
379: IF( N.EQ.0 ) THEN
380: MINWRK = 1
381: LWRK = 1
382: ELSE
383: MAXWRK = 2*N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 )
384: MINWRK = 3*N
385: *
386: CALL DHSEQR( 'S', JOBVS, N, 1, N, A, LDA, WR, WI, VS, LDVS,
387: $ WORK, -1, IEVAL )
388: HSWORK = WORK( 1 )
389: *
390: IF( .NOT.WANTVS ) THEN
391: MAXWRK = MAX( MAXWRK, N + HSWORK )
392: ELSE
393: MAXWRK = MAX( MAXWRK, 2*N + ( N - 1 )*ILAENV( 1,
394: $ 'DORGHR', ' ', N, 1, N, -1 ) )
395: MAXWRK = MAX( MAXWRK, N + HSWORK )
396: END IF
397: LWRK = MAXWRK
398: IF( .NOT.WANTSN )
399: $ LWRK = MAX( LWRK, N + ( N*N )/2 )
400: IF( WANTSV .OR. WANTSB )
401: $ LIWRK = ( N*N )/4
402: END IF
403: IWORK( 1 ) = LIWRK
404: WORK( 1 ) = LWRK
405: *
406: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
407: INFO = -16
408: ELSE IF( LIWORK.LT.1 .AND. .NOT.LQUERY ) THEN
409: INFO = -18
410: END IF
411: END IF
412: *
413: IF( INFO.NE.0 ) THEN
414: CALL XERBLA( 'DGEESX', -INFO )
415: RETURN
1.5 bertrand 416: ELSE IF( LQUERY ) THEN
417: RETURN
1.1 bertrand 418: END IF
419: *
420: * Quick return if possible
421: *
422: IF( N.EQ.0 ) THEN
423: SDIM = 0
424: RETURN
425: END IF
426: *
427: * Get machine constants
428: *
429: EPS = DLAMCH( 'P' )
430: SMLNUM = DLAMCH( 'S' )
431: BIGNUM = ONE / SMLNUM
432: CALL DLABAD( SMLNUM, BIGNUM )
433: SMLNUM = SQRT( SMLNUM ) / EPS
434: BIGNUM = ONE / SMLNUM
435: *
436: * Scale A if max element outside range [SMLNUM,BIGNUM]
437: *
438: ANRM = DLANGE( 'M', N, N, A, LDA, DUM )
439: SCALEA = .FALSE.
440: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
441: SCALEA = .TRUE.
442: CSCALE = SMLNUM
443: ELSE IF( ANRM.GT.BIGNUM ) THEN
444: SCALEA = .TRUE.
445: CSCALE = BIGNUM
446: END IF
447: IF( SCALEA )
448: $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
449: *
450: * Permute the matrix to make it more nearly triangular
451: * (RWorkspace: need N)
452: *
453: IBAL = 1
454: CALL DGEBAL( 'P', N, A, LDA, ILO, IHI, WORK( IBAL ), IERR )
455: *
456: * Reduce to upper Hessenberg form
457: * (RWorkspace: need 3*N, prefer 2*N+N*NB)
458: *
459: ITAU = N + IBAL
460: IWRK = N + ITAU
461: CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
462: $ LWORK-IWRK+1, IERR )
463: *
464: IF( WANTVS ) THEN
465: *
466: * Copy Householder vectors to VS
467: *
468: CALL DLACPY( 'L', N, N, A, LDA, VS, LDVS )
469: *
470: * Generate orthogonal matrix in VS
471: * (RWorkspace: need 3*N-1, prefer 2*N+(N-1)*NB)
472: *
473: CALL DORGHR( N, ILO, IHI, VS, LDVS, WORK( ITAU ), WORK( IWRK ),
474: $ LWORK-IWRK+1, IERR )
475: END IF
476: *
477: SDIM = 0
478: *
479: * Perform QR iteration, accumulating Schur vectors in VS if desired
480: * (RWorkspace: need N+1, prefer N+HSWORK (see comments) )
481: *
482: IWRK = ITAU
483: CALL DHSEQR( 'S', JOBVS, N, ILO, IHI, A, LDA, WR, WI, VS, LDVS,
484: $ WORK( IWRK ), LWORK-IWRK+1, IEVAL )
485: IF( IEVAL.GT.0 )
486: $ INFO = IEVAL
487: *
488: * Sort eigenvalues if desired
489: *
490: IF( WANTST .AND. INFO.EQ.0 ) THEN
491: IF( SCALEA ) THEN
492: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WR, N, IERR )
493: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, WI, N, IERR )
494: END IF
495: DO 10 I = 1, N
496: BWORK( I ) = SELECT( WR( I ), WI( I ) )
497: 10 CONTINUE
498: *
499: * Reorder eigenvalues, transform Schur vectors, and compute
500: * reciprocal condition numbers
501: * (RWorkspace: if SENSE is not 'N', need N+2*SDIM*(N-SDIM)
502: * otherwise, need N )
503: * (IWorkspace: if SENSE is 'V' or 'B', need SDIM*(N-SDIM)
504: * otherwise, need 0 )
505: *
506: CALL DTRSEN( SENSE, JOBVS, BWORK, N, A, LDA, VS, LDVS, WR, WI,
507: $ SDIM, RCONDE, RCONDV, WORK( IWRK ), LWORK-IWRK+1,
508: $ IWORK, LIWORK, ICOND )
509: IF( .NOT.WANTSN )
510: $ MAXWRK = MAX( MAXWRK, N+2*SDIM*( N-SDIM ) )
511: IF( ICOND.EQ.-15 ) THEN
512: *
513: * Not enough real workspace
514: *
515: INFO = -16
516: ELSE IF( ICOND.EQ.-17 ) THEN
517: *
518: * Not enough integer workspace
519: *
520: INFO = -18
521: ELSE IF( ICOND.GT.0 ) THEN
522: *
523: * DTRSEN failed to reorder or to restore standard Schur form
524: *
525: INFO = ICOND + N
526: END IF
527: END IF
528: *
529: IF( WANTVS ) THEN
530: *
531: * Undo balancing
532: * (RWorkspace: need N)
533: *
534: CALL DGEBAK( 'P', 'R', N, ILO, IHI, WORK( IBAL ), N, VS, LDVS,
535: $ IERR )
536: END IF
537: *
538: IF( SCALEA ) THEN
539: *
540: * Undo scaling for the Schur form of A
541: *
542: CALL DLASCL( 'H', 0, 0, CSCALE, ANRM, N, N, A, LDA, IERR )
543: CALL DCOPY( N, A, LDA+1, WR, 1 )
544: IF( ( WANTSV .OR. WANTSB ) .AND. INFO.EQ.0 ) THEN
545: DUM( 1 ) = RCONDV
546: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR )
547: RCONDV = DUM( 1 )
548: END IF
549: IF( CSCALE.EQ.SMLNUM ) THEN
550: *
551: * If scaling back towards underflow, adjust WI if an
552: * offdiagonal element of a 2-by-2 block in the Schur form
553: * underflows.
554: *
555: IF( IEVAL.GT.0 ) THEN
556: I1 = IEVAL + 1
557: I2 = IHI - 1
558: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N,
559: $ IERR )
560: ELSE IF( WANTST ) THEN
561: I1 = 1
562: I2 = N - 1
563: ELSE
564: I1 = ILO
565: I2 = IHI - 1
566: END IF
567: INXT = I1 - 1
568: DO 20 I = I1, I2
569: IF( I.LT.INXT )
570: $ GO TO 20
571: IF( WI( I ).EQ.ZERO ) THEN
572: INXT = I + 1
573: ELSE
574: IF( A( I+1, I ).EQ.ZERO ) THEN
575: WI( I ) = ZERO
576: WI( I+1 ) = ZERO
577: ELSE IF( A( I+1, I ).NE.ZERO .AND. A( I, I+1 ).EQ.
578: $ ZERO ) THEN
579: WI( I ) = ZERO
580: WI( I+1 ) = ZERO
581: IF( I.GT.1 )
582: $ CALL DSWAP( I-1, A( 1, I ), 1, A( 1, I+1 ), 1 )
583: IF( N.GT.I+1 )
584: $ CALL DSWAP( N-I-1, A( I, I+2 ), LDA,
585: $ A( I+1, I+2 ), LDA )
586: CALL DSWAP( N, VS( 1, I ), 1, VS( 1, I+1 ), 1 )
587: A( I, I+1 ) = A( I+1, I )
588: A( I+1, I ) = ZERO
589: END IF
590: INXT = I + 2
591: END IF
592: 20 CONTINUE
593: END IF
594: CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-IEVAL, 1,
595: $ WI( IEVAL+1 ), MAX( N-IEVAL, 1 ), IERR )
596: END IF
597: *
598: IF( WANTST .AND. INFO.EQ.0 ) THEN
599: *
600: * Check if reordering successful
601: *
602: LASTSL = .TRUE.
603: LST2SL = .TRUE.
604: SDIM = 0
605: IP = 0
606: DO 30 I = 1, N
607: CURSL = SELECT( WR( I ), WI( I ) )
608: IF( WI( I ).EQ.ZERO ) THEN
609: IF( CURSL )
610: $ SDIM = SDIM + 1
611: IP = 0
612: IF( CURSL .AND. .NOT.LASTSL )
613: $ INFO = N + 2
614: ELSE
615: IF( IP.EQ.1 ) THEN
616: *
617: * Last eigenvalue of conjugate pair
618: *
619: CURSL = CURSL .OR. LASTSL
620: LASTSL = CURSL
621: IF( CURSL )
622: $ SDIM = SDIM + 2
623: IP = -1
624: IF( CURSL .AND. .NOT.LST2SL )
625: $ INFO = N + 2
626: ELSE
627: *
628: * First eigenvalue of conjugate pair
629: *
630: IP = 1
631: END IF
632: END IF
633: LST2SL = LASTSL
634: LASTSL = CURSL
635: 30 CONTINUE
636: END IF
637: *
638: WORK( 1 ) = MAXWRK
639: IF( WANTSV .OR. WANTSB ) THEN
640: IWORK( 1 ) = MAX( 1, SDIM*( N-SDIM ) )
641: ELSE
642: IWORK( 1 ) = 1
643: END IF
644: *
645: RETURN
646: *
647: * End of DGEESX
648: *
649: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgeesx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack