Annotation of rpl/lapack/lapack/dggesx.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DGGESX( JOBVSL, JOBVSR, SORT, SELCTG, SENSE, N, A, LDA,
! 2: $ B, LDB, SDIM, ALPHAR, ALPHAI, BETA, VSL, LDVSL,
! 3: $ VSR, LDVSR, RCONDE, RCONDV, WORK, LWORK, IWORK,
! 4: $ LIWORK, BWORK, INFO )
! 5: *
! 6: * -- LAPACK driver routine (version 3.2.1) --
! 7: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 8: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 9: * -- April 2009 --
! 10: *
! 11: * .. Scalar Arguments ..
! 12: CHARACTER JOBVSL, JOBVSR, SENSE, SORT
! 13: INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LIWORK, LWORK, N,
! 14: $ SDIM
! 15: * ..
! 16: * .. Array Arguments ..
! 17: LOGICAL BWORK( * )
! 18: INTEGER IWORK( * )
! 19: DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
! 20: $ B( LDB, * ), BETA( * ), RCONDE( 2 ),
! 21: $ RCONDV( 2 ), VSL( LDVSL, * ), VSR( LDVSR, * ),
! 22: $ WORK( * )
! 23: * ..
! 24: * .. Function Arguments ..
! 25: LOGICAL SELCTG
! 26: EXTERNAL SELCTG
! 27: * ..
! 28: *
! 29: * Purpose
! 30: * =======
! 31: *
! 32: * DGGESX computes for a pair of N-by-N real nonsymmetric matrices
! 33: * (A,B), the generalized eigenvalues, the real Schur form (S,T), and,
! 34: * optionally, the left and/or right matrices of Schur vectors (VSL and
! 35: * VSR). This gives the generalized Schur factorization
! 36: *
! 37: * (A,B) = ( (VSL) S (VSR)**T, (VSL) T (VSR)**T )
! 38: *
! 39: * Optionally, it also orders the eigenvalues so that a selected cluster
! 40: * of eigenvalues appears in the leading diagonal blocks of the upper
! 41: * quasi-triangular matrix S and the upper triangular matrix T; computes
! 42: * a reciprocal condition number for the average of the selected
! 43: * eigenvalues (RCONDE); and computes a reciprocal condition number for
! 44: * the right and left deflating subspaces corresponding to the selected
! 45: * eigenvalues (RCONDV). The leading columns of VSL and VSR then form
! 46: * an orthonormal basis for the corresponding left and right eigenspaces
! 47: * (deflating subspaces).
! 48: *
! 49: * A generalized eigenvalue for a pair of matrices (A,B) is a scalar w
! 50: * or a ratio alpha/beta = w, such that A - w*B is singular. It is
! 51: * usually represented as the pair (alpha,beta), as there is a
! 52: * reasonable interpretation for beta=0 or for both being zero.
! 53: *
! 54: * A pair of matrices (S,T) is in generalized real Schur form if T is
! 55: * upper triangular with non-negative diagonal and S is block upper
! 56: * triangular with 1-by-1 and 2-by-2 blocks. 1-by-1 blocks correspond
! 57: * to real generalized eigenvalues, while 2-by-2 blocks of S will be
! 58: * "standardized" by making the corresponding elements of T have the
! 59: * form:
! 60: * [ a 0 ]
! 61: * [ 0 b ]
! 62: *
! 63: * and the pair of corresponding 2-by-2 blocks in S and T will have a
! 64: * complex conjugate pair of generalized eigenvalues.
! 65: *
! 66: *
! 67: * Arguments
! 68: * =========
! 69: *
! 70: * JOBVSL (input) CHARACTER*1
! 71: * = 'N': do not compute the left Schur vectors;
! 72: * = 'V': compute the left Schur vectors.
! 73: *
! 74: * JOBVSR (input) CHARACTER*1
! 75: * = 'N': do not compute the right Schur vectors;
! 76: * = 'V': compute the right Schur vectors.
! 77: *
! 78: * SORT (input) CHARACTER*1
! 79: * Specifies whether or not to order the eigenvalues on the
! 80: * diagonal of the generalized Schur form.
! 81: * = 'N': Eigenvalues are not ordered;
! 82: * = 'S': Eigenvalues are ordered (see SELCTG).
! 83: *
! 84: * SELCTG (external procedure) LOGICAL FUNCTION of three DOUBLE PRECISION arguments
! 85: * SELCTG must be declared EXTERNAL in the calling subroutine.
! 86: * If SORT = 'N', SELCTG is not referenced.
! 87: * If SORT = 'S', SELCTG is used to select eigenvalues to sort
! 88: * to the top left of the Schur form.
! 89: * An eigenvalue (ALPHAR(j)+ALPHAI(j))/BETA(j) is selected if
! 90: * SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) is true; i.e. if either
! 91: * one of a complex conjugate pair of eigenvalues is selected,
! 92: * then both complex eigenvalues are selected.
! 93: * Note that a selected complex eigenvalue may no longer satisfy
! 94: * SELCTG(ALPHAR(j),ALPHAI(j),BETA(j)) = .TRUE. after ordering,
! 95: * since ordering may change the value of complex eigenvalues
! 96: * (especially if the eigenvalue is ill-conditioned), in this
! 97: * case INFO is set to N+3.
! 98: *
! 99: * SENSE (input) CHARACTER*1
! 100: * Determines which reciprocal condition numbers are computed.
! 101: * = 'N' : None are computed;
! 102: * = 'E' : Computed for average of selected eigenvalues only;
! 103: * = 'V' : Computed for selected deflating subspaces only;
! 104: * = 'B' : Computed for both.
! 105: * If SENSE = 'E', 'V', or 'B', SORT must equal 'S'.
! 106: *
! 107: * N (input) INTEGER
! 108: * The order of the matrices A, B, VSL, and VSR. N >= 0.
! 109: *
! 110: * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
! 111: * On entry, the first of the pair of matrices.
! 112: * On exit, A has been overwritten by its generalized Schur
! 113: * form S.
! 114: *
! 115: * LDA (input) INTEGER
! 116: * The leading dimension of A. LDA >= max(1,N).
! 117: *
! 118: * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
! 119: * On entry, the second of the pair of matrices.
! 120: * On exit, B has been overwritten by its generalized Schur
! 121: * form T.
! 122: *
! 123: * LDB (input) INTEGER
! 124: * The leading dimension of B. LDB >= max(1,N).
! 125: *
! 126: * SDIM (output) INTEGER
! 127: * If SORT = 'N', SDIM = 0.
! 128: * If SORT = 'S', SDIM = number of eigenvalues (after sorting)
! 129: * for which SELCTG is true. (Complex conjugate pairs for which
! 130: * SELCTG is true for either eigenvalue count as 2.)
! 131: *
! 132: * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
! 133: * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
! 134: * BETA (output) DOUBLE PRECISION array, dimension (N)
! 135: * On exit, (ALPHAR(j) + ALPHAI(j)*i)/BETA(j), j=1,...,N, will
! 136: * be the generalized eigenvalues. ALPHAR(j) + ALPHAI(j)*i
! 137: * and BETA(j),j=1,...,N are the diagonals of the complex Schur
! 138: * form (S,T) that would result if the 2-by-2 diagonal blocks of
! 139: * the real Schur form of (A,B) were further reduced to
! 140: * triangular form using 2-by-2 complex unitary transformations.
! 141: * If ALPHAI(j) is zero, then the j-th eigenvalue is real; if
! 142: * positive, then the j-th and (j+1)-st eigenvalues are a
! 143: * complex conjugate pair, with ALPHAI(j+1) negative.
! 144: *
! 145: * Note: the quotients ALPHAR(j)/BETA(j) and ALPHAI(j)/BETA(j)
! 146: * may easily over- or underflow, and BETA(j) may even be zero.
! 147: * Thus, the user should avoid naively computing the ratio.
! 148: * However, ALPHAR and ALPHAI will be always less than and
! 149: * usually comparable with norm(A) in magnitude, and BETA always
! 150: * less than and usually comparable with norm(B).
! 151: *
! 152: * VSL (output) DOUBLE PRECISION array, dimension (LDVSL,N)
! 153: * If JOBVSL = 'V', VSL will contain the left Schur vectors.
! 154: * Not referenced if JOBVSL = 'N'.
! 155: *
! 156: * LDVSL (input) INTEGER
! 157: * The leading dimension of the matrix VSL. LDVSL >=1, and
! 158: * if JOBVSL = 'V', LDVSL >= N.
! 159: *
! 160: * VSR (output) DOUBLE PRECISION array, dimension (LDVSR,N)
! 161: * If JOBVSR = 'V', VSR will contain the right Schur vectors.
! 162: * Not referenced if JOBVSR = 'N'.
! 163: *
! 164: * LDVSR (input) INTEGER
! 165: * The leading dimension of the matrix VSR. LDVSR >= 1, and
! 166: * if JOBVSR = 'V', LDVSR >= N.
! 167: *
! 168: * RCONDE (output) DOUBLE PRECISION array, dimension ( 2 )
! 169: * If SENSE = 'E' or 'B', RCONDE(1) and RCONDE(2) contain the
! 170: * reciprocal condition numbers for the average of the selected
! 171: * eigenvalues.
! 172: * Not referenced if SENSE = 'N' or 'V'.
! 173: *
! 174: * RCONDV (output) DOUBLE PRECISION array, dimension ( 2 )
! 175: * If SENSE = 'V' or 'B', RCONDV(1) and RCONDV(2) contain the
! 176: * reciprocal condition numbers for the selected deflating
! 177: * subspaces.
! 178: * Not referenced if SENSE = 'N' or 'E'.
! 179: *
! 180: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 181: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 182: *
! 183: * LWORK (input) INTEGER
! 184: * The dimension of the array WORK.
! 185: * If N = 0, LWORK >= 1, else if SENSE = 'E', 'V', or 'B',
! 186: * LWORK >= max( 8*N, 6*N+16, 2*SDIM*(N-SDIM) ), else
! 187: * LWORK >= max( 8*N, 6*N+16 ).
! 188: * Note that 2*SDIM*(N-SDIM) <= N*N/2.
! 189: * Note also that an error is only returned if
! 190: * LWORK < max( 8*N, 6*N+16), but if SENSE = 'E' or 'V' or 'B'
! 191: * this may not be large enough.
! 192: *
! 193: * If LWORK = -1, then a workspace query is assumed; the routine
! 194: * only calculates the bound on the optimal size of the WORK
! 195: * array and the minimum size of the IWORK array, returns these
! 196: * values as the first entries of the WORK and IWORK arrays, and
! 197: * no error message related to LWORK or LIWORK is issued by
! 198: * XERBLA.
! 199: *
! 200: * IWORK (workspace) INTEGER array, dimension (MAX(1,LIWORK))
! 201: * On exit, if INFO = 0, IWORK(1) returns the minimum LIWORK.
! 202: *
! 203: * LIWORK (input) INTEGER
! 204: * The dimension of the array IWORK.
! 205: * If SENSE = 'N' or N = 0, LIWORK >= 1, otherwise
! 206: * LIWORK >= N+6.
! 207: *
! 208: * If LIWORK = -1, then a workspace query is assumed; the
! 209: * routine only calculates the bound on the optimal size of the
! 210: * WORK array and the minimum size of the IWORK array, returns
! 211: * these values as the first entries of the WORK and IWORK
! 212: * arrays, and no error message related to LWORK or LIWORK is
! 213: * issued by XERBLA.
! 214: *
! 215: * BWORK (workspace) LOGICAL array, dimension (N)
! 216: * Not referenced if SORT = 'N'.
! 217: *
! 218: * INFO (output) INTEGER
! 219: * = 0: successful exit
! 220: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 221: * = 1,...,N:
! 222: * The QZ iteration failed. (A,B) are not in Schur
! 223: * form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
! 224: * be correct for j=INFO+1,...,N.
! 225: * > N: =N+1: other than QZ iteration failed in DHGEQZ
! 226: * =N+2: after reordering, roundoff changed values of
! 227: * some complex eigenvalues so that leading
! 228: * eigenvalues in the Generalized Schur form no
! 229: * longer satisfy SELCTG=.TRUE. This could also
! 230: * be caused due to scaling.
! 231: * =N+3: reordering failed in DTGSEN.
! 232: *
! 233: * Further Details
! 234: * ===============
! 235: *
! 236: * An approximate (asymptotic) bound on the average absolute error of
! 237: * the selected eigenvalues is
! 238: *
! 239: * EPS * norm((A, B)) / RCONDE( 1 ).
! 240: *
! 241: * An approximate (asymptotic) bound on the maximum angular error in
! 242: * the computed deflating subspaces is
! 243: *
! 244: * EPS * norm((A, B)) / RCONDV( 2 ).
! 245: *
! 246: * See LAPACK User's Guide, section 4.11 for more information.
! 247: *
! 248: * =====================================================================
! 249: *
! 250: * .. Parameters ..
! 251: DOUBLE PRECISION ZERO, ONE
! 252: PARAMETER ( ZERO = 0.0D+0, ONE = 1.0D+0 )
! 253: * ..
! 254: * .. Local Scalars ..
! 255: LOGICAL CURSL, ILASCL, ILBSCL, ILVSL, ILVSR, LASTSL,
! 256: $ LQUERY, LST2SL, WANTSB, WANTSE, WANTSN, WANTST,
! 257: $ WANTSV
! 258: INTEGER I, ICOLS, IERR, IHI, IJOB, IJOBVL, IJOBVR,
! 259: $ ILEFT, ILO, IP, IRIGHT, IROWS, ITAU, IWRK,
! 260: $ LIWMIN, LWRK, MAXWRK, MINWRK
! 261: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS, PL,
! 262: $ PR, SAFMAX, SAFMIN, SMLNUM
! 263: * ..
! 264: * .. Local Arrays ..
! 265: DOUBLE PRECISION DIF( 2 )
! 266: * ..
! 267: * .. External Subroutines ..
! 268: EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLABAD,
! 269: $ DLACPY, DLASCL, DLASET, DORGQR, DORMQR, DTGSEN,
! 270: $ XERBLA
! 271: * ..
! 272: * .. External Functions ..
! 273: LOGICAL LSAME
! 274: INTEGER ILAENV
! 275: DOUBLE PRECISION DLAMCH, DLANGE
! 276: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
! 277: * ..
! 278: * .. Intrinsic Functions ..
! 279: INTRINSIC ABS, MAX, SQRT
! 280: * ..
! 281: * .. Executable Statements ..
! 282: *
! 283: * Decode the input arguments
! 284: *
! 285: IF( LSAME( JOBVSL, 'N' ) ) THEN
! 286: IJOBVL = 1
! 287: ILVSL = .FALSE.
! 288: ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
! 289: IJOBVL = 2
! 290: ILVSL = .TRUE.
! 291: ELSE
! 292: IJOBVL = -1
! 293: ILVSL = .FALSE.
! 294: END IF
! 295: *
! 296: IF( LSAME( JOBVSR, 'N' ) ) THEN
! 297: IJOBVR = 1
! 298: ILVSR = .FALSE.
! 299: ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
! 300: IJOBVR = 2
! 301: ILVSR = .TRUE.
! 302: ELSE
! 303: IJOBVR = -1
! 304: ILVSR = .FALSE.
! 305: END IF
! 306: *
! 307: WANTST = LSAME( SORT, 'S' )
! 308: WANTSN = LSAME( SENSE, 'N' )
! 309: WANTSE = LSAME( SENSE, 'E' )
! 310: WANTSV = LSAME( SENSE, 'V' )
! 311: WANTSB = LSAME( SENSE, 'B' )
! 312: LQUERY = ( LWORK.EQ.-1 .OR. LIWORK.EQ.-1 )
! 313: IF( WANTSN ) THEN
! 314: IJOB = 0
! 315: ELSE IF( WANTSE ) THEN
! 316: IJOB = 1
! 317: ELSE IF( WANTSV ) THEN
! 318: IJOB = 2
! 319: ELSE IF( WANTSB ) THEN
! 320: IJOB = 4
! 321: END IF
! 322: *
! 323: * Test the input arguments
! 324: *
! 325: INFO = 0
! 326: IF( IJOBVL.LE.0 ) THEN
! 327: INFO = -1
! 328: ELSE IF( IJOBVR.LE.0 ) THEN
! 329: INFO = -2
! 330: ELSE IF( ( .NOT.WANTST ) .AND. ( .NOT.LSAME( SORT, 'N' ) ) ) THEN
! 331: INFO = -3
! 332: ELSE IF( .NOT.( WANTSN .OR. WANTSE .OR. WANTSV .OR. WANTSB ) .OR.
! 333: $ ( .NOT.WANTST .AND. .NOT.WANTSN ) ) THEN
! 334: INFO = -5
! 335: ELSE IF( N.LT.0 ) THEN
! 336: INFO = -6
! 337: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 338: INFO = -8
! 339: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 340: INFO = -10
! 341: ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
! 342: INFO = -16
! 343: ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
! 344: INFO = -18
! 345: END IF
! 346: *
! 347: * Compute workspace
! 348: * (Note: Comments in the code beginning "Workspace:" describe the
! 349: * minimal amount of workspace needed at that point in the code,
! 350: * as well as the preferred amount for good performance.
! 351: * NB refers to the optimal block size for the immediately
! 352: * following subroutine, as returned by ILAENV.)
! 353: *
! 354: IF( INFO.EQ.0 ) THEN
! 355: IF( N.GT.0) THEN
! 356: MINWRK = MAX( 8*N, 6*N + 16 )
! 357: MAXWRK = MINWRK - N +
! 358: $ N*ILAENV( 1, 'DGEQRF', ' ', N, 1, N, 0 )
! 359: MAXWRK = MAX( MAXWRK, MINWRK - N +
! 360: $ N*ILAENV( 1, 'DORMQR', ' ', N, 1, N, -1 ) )
! 361: IF( ILVSL ) THEN
! 362: MAXWRK = MAX( MAXWRK, MINWRK - N +
! 363: $ N*ILAENV( 1, 'DORGQR', ' ', N, 1, N, -1 ) )
! 364: END IF
! 365: LWRK = MAXWRK
! 366: IF( IJOB.GE.1 )
! 367: $ LWRK = MAX( LWRK, N*N/2 )
! 368: ELSE
! 369: MINWRK = 1
! 370: MAXWRK = 1
! 371: LWRK = 1
! 372: END IF
! 373: WORK( 1 ) = LWRK
! 374: IF( WANTSN .OR. N.EQ.0 ) THEN
! 375: LIWMIN = 1
! 376: ELSE
! 377: LIWMIN = N + 6
! 378: END IF
! 379: IWORK( 1 ) = LIWMIN
! 380: *
! 381: IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
! 382: INFO = -22
! 383: ELSE IF( LIWORK.LT.LIWMIN .AND. .NOT.LQUERY ) THEN
! 384: INFO = -24
! 385: END IF
! 386: END IF
! 387: *
! 388: IF( INFO.NE.0 ) THEN
! 389: CALL XERBLA( 'DGGESX', -INFO )
! 390: RETURN
! 391: ELSE IF (LQUERY) THEN
! 392: RETURN
! 393: END IF
! 394: *
! 395: * Quick return if possible
! 396: *
! 397: IF( N.EQ.0 ) THEN
! 398: SDIM = 0
! 399: RETURN
! 400: END IF
! 401: *
! 402: * Get machine constants
! 403: *
! 404: EPS = DLAMCH( 'P' )
! 405: SAFMIN = DLAMCH( 'S' )
! 406: SAFMAX = ONE / SAFMIN
! 407: CALL DLABAD( SAFMIN, SAFMAX )
! 408: SMLNUM = SQRT( SAFMIN ) / EPS
! 409: BIGNUM = ONE / SMLNUM
! 410: *
! 411: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 412: *
! 413: ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
! 414: ILASCL = .FALSE.
! 415: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 416: ANRMTO = SMLNUM
! 417: ILASCL = .TRUE.
! 418: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 419: ANRMTO = BIGNUM
! 420: ILASCL = .TRUE.
! 421: END IF
! 422: IF( ILASCL )
! 423: $ CALL DLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
! 424: *
! 425: * Scale B if max element outside range [SMLNUM,BIGNUM]
! 426: *
! 427: BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
! 428: ILBSCL = .FALSE.
! 429: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 430: BNRMTO = SMLNUM
! 431: ILBSCL = .TRUE.
! 432: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 433: BNRMTO = BIGNUM
! 434: ILBSCL = .TRUE.
! 435: END IF
! 436: IF( ILBSCL )
! 437: $ CALL DLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
! 438: *
! 439: * Permute the matrix to make it more nearly triangular
! 440: * (Workspace: need 6*N + 2*N for permutation parameters)
! 441: *
! 442: ILEFT = 1
! 443: IRIGHT = N + 1
! 444: IWRK = IRIGHT + N
! 445: CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
! 446: $ WORK( IRIGHT ), WORK( IWRK ), IERR )
! 447: *
! 448: * Reduce B to triangular form (QR decomposition of B)
! 449: * (Workspace: need N, prefer N*NB)
! 450: *
! 451: IROWS = IHI + 1 - ILO
! 452: ICOLS = N + 1 - ILO
! 453: ITAU = IWRK
! 454: IWRK = ITAU + IROWS
! 455: CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
! 456: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
! 457: *
! 458: * Apply the orthogonal transformation to matrix A
! 459: * (Workspace: need N, prefer N*NB)
! 460: *
! 461: CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
! 462: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
! 463: $ LWORK+1-IWRK, IERR )
! 464: *
! 465: * Initialize VSL
! 466: * (Workspace: need N, prefer N*NB)
! 467: *
! 468: IF( ILVSL ) THEN
! 469: CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
! 470: IF( IROWS.GT.1 ) THEN
! 471: CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
! 472: $ VSL( ILO+1, ILO ), LDVSL )
! 473: END IF
! 474: CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
! 475: $ WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
! 476: END IF
! 477: *
! 478: * Initialize VSR
! 479: *
! 480: IF( ILVSR )
! 481: $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
! 482: *
! 483: * Reduce to generalized Hessenberg form
! 484: * (Workspace: none needed)
! 485: *
! 486: CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
! 487: $ LDVSL, VSR, LDVSR, IERR )
! 488: *
! 489: SDIM = 0
! 490: *
! 491: * Perform QZ algorithm, computing Schur vectors if desired
! 492: * (Workspace: need N)
! 493: *
! 494: IWRK = ITAU
! 495: CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
! 496: $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
! 497: $ WORK( IWRK ), LWORK+1-IWRK, IERR )
! 498: IF( IERR.NE.0 ) THEN
! 499: IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
! 500: INFO = IERR
! 501: ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
! 502: INFO = IERR - N
! 503: ELSE
! 504: INFO = N + 1
! 505: END IF
! 506: GO TO 60
! 507: END IF
! 508: *
! 509: * Sort eigenvalues ALPHA/BETA and compute the reciprocal of
! 510: * condition number(s)
! 511: * (Workspace: If IJOB >= 1, need MAX( 8*(N+1), 2*SDIM*(N-SDIM) )
! 512: * otherwise, need 8*(N+1) )
! 513: *
! 514: IF( WANTST ) THEN
! 515: *
! 516: * Undo scaling on eigenvalues before SELCTGing
! 517: *
! 518: IF( ILASCL ) THEN
! 519: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N,
! 520: $ IERR )
! 521: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N,
! 522: $ IERR )
! 523: END IF
! 524: IF( ILBSCL )
! 525: $ CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
! 526: *
! 527: * Select eigenvalues
! 528: *
! 529: DO 10 I = 1, N
! 530: BWORK( I ) = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
! 531: 10 CONTINUE
! 532: *
! 533: * Reorder eigenvalues, transform Generalized Schur vectors, and
! 534: * compute reciprocal condition numbers
! 535: *
! 536: CALL DTGSEN( IJOB, ILVSL, ILVSR, BWORK, N, A, LDA, B, LDB,
! 537: $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
! 538: $ SDIM, PL, PR, DIF, WORK( IWRK ), LWORK-IWRK+1,
! 539: $ IWORK, LIWORK, IERR )
! 540: *
! 541: IF( IJOB.GE.1 )
! 542: $ MAXWRK = MAX( MAXWRK, 2*SDIM*( N-SDIM ) )
! 543: IF( IERR.EQ.-22 ) THEN
! 544: *
! 545: * not enough real workspace
! 546: *
! 547: INFO = -22
! 548: ELSE
! 549: IF( IJOB.EQ.1 .OR. IJOB.EQ.4 ) THEN
! 550: RCONDE( 1 ) = PL
! 551: RCONDE( 2 ) = PR
! 552: END IF
! 553: IF( IJOB.EQ.2 .OR. IJOB.EQ.4 ) THEN
! 554: RCONDV( 1 ) = DIF( 1 )
! 555: RCONDV( 2 ) = DIF( 2 )
! 556: END IF
! 557: IF( IERR.EQ.1 )
! 558: $ INFO = N + 3
! 559: END IF
! 560: *
! 561: END IF
! 562: *
! 563: * Apply permutation to VSL and VSR
! 564: * (Workspace: none needed)
! 565: *
! 566: IF( ILVSL )
! 567: $ CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
! 568: $ WORK( IRIGHT ), N, VSL, LDVSL, IERR )
! 569: *
! 570: IF( ILVSR )
! 571: $ CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
! 572: $ WORK( IRIGHT ), N, VSR, LDVSR, IERR )
! 573: *
! 574: * Check if unscaling would cause over/underflow, if so, rescale
! 575: * (ALPHAR(I),ALPHAI(I),BETA(I)) so BETA(I) is on the order of
! 576: * B(I,I) and ALPHAR(I) and ALPHAI(I) are on the order of A(I,I)
! 577: *
! 578: IF( ILASCL ) THEN
! 579: DO 20 I = 1, N
! 580: IF( ALPHAI( I ).NE.ZERO ) THEN
! 581: IF( ( ALPHAR( I ) / SAFMAX ).GT.( ANRMTO / ANRM ) .OR.
! 582: $ ( SAFMIN / ALPHAR( I ) ).GT.( ANRM / ANRMTO ) ) THEN
! 583: WORK( 1 ) = ABS( A( I, I ) / ALPHAR( I ) )
! 584: BETA( I ) = BETA( I )*WORK( 1 )
! 585: ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
! 586: ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
! 587: ELSE IF( ( ALPHAI( I ) / SAFMAX ).GT.
! 588: $ ( ANRMTO / ANRM ) .OR.
! 589: $ ( SAFMIN / ALPHAI( I ) ).GT.( ANRM / ANRMTO ) )
! 590: $ THEN
! 591: WORK( 1 ) = ABS( A( I, I+1 ) / ALPHAI( I ) )
! 592: BETA( I ) = BETA( I )*WORK( 1 )
! 593: ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
! 594: ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
! 595: END IF
! 596: END IF
! 597: 20 CONTINUE
! 598: END IF
! 599: *
! 600: IF( ILBSCL ) THEN
! 601: DO 30 I = 1, N
! 602: IF( ALPHAI( I ).NE.ZERO ) THEN
! 603: IF( ( BETA( I ) / SAFMAX ).GT.( BNRMTO / BNRM ) .OR.
! 604: $ ( SAFMIN / BETA( I ) ).GT.( BNRM / BNRMTO ) ) THEN
! 605: WORK( 1 ) = ABS( B( I, I ) / BETA( I ) )
! 606: BETA( I ) = BETA( I )*WORK( 1 )
! 607: ALPHAR( I ) = ALPHAR( I )*WORK( 1 )
! 608: ALPHAI( I ) = ALPHAI( I )*WORK( 1 )
! 609: END IF
! 610: END IF
! 611: 30 CONTINUE
! 612: END IF
! 613: *
! 614: * Undo scaling
! 615: *
! 616: IF( ILASCL ) THEN
! 617: CALL DLASCL( 'H', 0, 0, ANRMTO, ANRM, N, N, A, LDA, IERR )
! 618: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAR, N, IERR )
! 619: CALL DLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHAI, N, IERR )
! 620: END IF
! 621: *
! 622: IF( ILBSCL ) THEN
! 623: CALL DLASCL( 'U', 0, 0, BNRMTO, BNRM, N, N, B, LDB, IERR )
! 624: CALL DLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
! 625: END IF
! 626: *
! 627: IF( WANTST ) THEN
! 628: *
! 629: * Check if reordering is correct
! 630: *
! 631: LASTSL = .TRUE.
! 632: LST2SL = .TRUE.
! 633: SDIM = 0
! 634: IP = 0
! 635: DO 50 I = 1, N
! 636: CURSL = SELCTG( ALPHAR( I ), ALPHAI( I ), BETA( I ) )
! 637: IF( ALPHAI( I ).EQ.ZERO ) THEN
! 638: IF( CURSL )
! 639: $ SDIM = SDIM + 1
! 640: IP = 0
! 641: IF( CURSL .AND. .NOT.LASTSL )
! 642: $ INFO = N + 2
! 643: ELSE
! 644: IF( IP.EQ.1 ) THEN
! 645: *
! 646: * Last eigenvalue of conjugate pair
! 647: *
! 648: CURSL = CURSL .OR. LASTSL
! 649: LASTSL = CURSL
! 650: IF( CURSL )
! 651: $ SDIM = SDIM + 2
! 652: IP = -1
! 653: IF( CURSL .AND. .NOT.LST2SL )
! 654: $ INFO = N + 2
! 655: ELSE
! 656: *
! 657: * First eigenvalue of conjugate pair
! 658: *
! 659: IP = 1
! 660: END IF
! 661: END IF
! 662: LST2SL = LASTSL
! 663: LASTSL = CURSL
! 664: 50 CONTINUE
! 665: *
! 666: END IF
! 667: *
! 668: 60 CONTINUE
! 669: *
! 670: WORK( 1 ) = MAXWRK
! 671: IWORK( 1 ) = LIWMIN
! 672: *
! 673: RETURN
! 674: *
! 675: * End of DGGESX
! 676: *
! 677: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dggesx.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack