Annotation of rpl/lapack/lapack/dgegs.f, revision 1.1
1.1 ! bertrand 1: SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
! 2: $ ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
! 3: $ LWORK, INFO )
! 4: *
! 5: * -- LAPACK driver routine (version 3.2) --
! 6: * -- LAPACK is a software package provided by Univ. of Tennessee, --
! 7: * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
! 8: * November 2006
! 9: *
! 10: * .. Scalar Arguments ..
! 11: CHARACTER JOBVSL, JOBVSR
! 12: INTEGER INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
! 13: * ..
! 14: * .. Array Arguments ..
! 15: DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
! 16: $ B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
! 17: $ VSR( LDVSR, * ), WORK( * )
! 18: * ..
! 19: *
! 20: * Purpose
! 21: * =======
! 22: *
! 23: * This routine is deprecated and has been replaced by routine DGGES.
! 24: *
! 25: * DGEGS computes the eigenvalues, real Schur form, and, optionally,
! 26: * left and or/right Schur vectors of a real matrix pair (A,B).
! 27: * Given two square matrices A and B, the generalized real Schur
! 28: * factorization has the form
! 29: *
! 30: * A = Q*S*Z**T, B = Q*T*Z**T
! 31: *
! 32: * where Q and Z are orthogonal matrices, T is upper triangular, and S
! 33: * is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
! 34: * blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
! 35: * of eigenvalues of (A,B). The columns of Q are the left Schur vectors
! 36: * and the columns of Z are the right Schur vectors.
! 37: *
! 38: * If only the eigenvalues of (A,B) are needed, the driver routine
! 39: * DGEGV should be used instead. See DGEGV for a description of the
! 40: * eigenvalues of the generalized nonsymmetric eigenvalue problem
! 41: * (GNEP).
! 42: *
! 43: * Arguments
! 44: * =========
! 45: *
! 46: * JOBVSL (input) CHARACTER*1
! 47: * = 'N': do not compute the left Schur vectors;
! 48: * = 'V': compute the left Schur vectors (returned in VSL).
! 49: *
! 50: * JOBVSR (input) CHARACTER*1
! 51: * = 'N': do not compute the right Schur vectors;
! 52: * = 'V': compute the right Schur vectors (returned in VSR).
! 53: *
! 54: * N (input) INTEGER
! 55: * The order of the matrices A, B, VSL, and VSR. N >= 0.
! 56: *
! 57: * A (input/output) DOUBLE PRECISION array, dimension (LDA, N)
! 58: * On entry, the matrix A.
! 59: * On exit, the upper quasi-triangular matrix S from the
! 60: * generalized real Schur factorization.
! 61: *
! 62: * LDA (input) INTEGER
! 63: * The leading dimension of A. LDA >= max(1,N).
! 64: *
! 65: * B (input/output) DOUBLE PRECISION array, dimension (LDB, N)
! 66: * On entry, the matrix B.
! 67: * On exit, the upper triangular matrix T from the generalized
! 68: * real Schur factorization.
! 69: *
! 70: * LDB (input) INTEGER
! 71: * The leading dimension of B. LDB >= max(1,N).
! 72: *
! 73: * ALPHAR (output) DOUBLE PRECISION array, dimension (N)
! 74: * The real parts of each scalar alpha defining an eigenvalue
! 75: * of GNEP.
! 76: *
! 77: * ALPHAI (output) DOUBLE PRECISION array, dimension (N)
! 78: * The imaginary parts of each scalar alpha defining an
! 79: * eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th
! 80: * eigenvalue is real; if positive, then the j-th and (j+1)-st
! 81: * eigenvalues are a complex conjugate pair, with
! 82: * ALPHAI(j+1) = -ALPHAI(j).
! 83: *
! 84: * BETA (output) DOUBLE PRECISION array, dimension (N)
! 85: * The scalars beta that define the eigenvalues of GNEP.
! 86: * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
! 87: * beta = BETA(j) represent the j-th eigenvalue of the matrix
! 88: * pair (A,B), in one of the forms lambda = alpha/beta or
! 89: * mu = beta/alpha. Since either lambda or mu may overflow,
! 90: * they should not, in general, be computed.
! 91: *
! 92: * VSL (output) DOUBLE PRECISION array, dimension (LDVSL,N)
! 93: * If JOBVSL = 'V', the matrix of left Schur vectors Q.
! 94: * Not referenced if JOBVSL = 'N'.
! 95: *
! 96: * LDVSL (input) INTEGER
! 97: * The leading dimension of the matrix VSL. LDVSL >=1, and
! 98: * if JOBVSL = 'V', LDVSL >= N.
! 99: *
! 100: * VSR (output) DOUBLE PRECISION array, dimension (LDVSR,N)
! 101: * If JOBVSR = 'V', the matrix of right Schur vectors Z.
! 102: * Not referenced if JOBVSR = 'N'.
! 103: *
! 104: * LDVSR (input) INTEGER
! 105: * The leading dimension of the matrix VSR. LDVSR >= 1, and
! 106: * if JOBVSR = 'V', LDVSR >= N.
! 107: *
! 108: * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
! 109: * On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
! 110: *
! 111: * LWORK (input) INTEGER
! 112: * The dimension of the array WORK. LWORK >= max(1,4*N).
! 113: * For good performance, LWORK must generally be larger.
! 114: * To compute the optimal value of LWORK, call ILAENV to get
! 115: * blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute:
! 116: * NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR
! 117: * The optimal LWORK is 2*N + N*(NB+1).
! 118: *
! 119: * If LWORK = -1, then a workspace query is assumed; the routine
! 120: * only calculates the optimal size of the WORK array, returns
! 121: * this value as the first entry of the WORK array, and no error
! 122: * message related to LWORK is issued by XERBLA.
! 123: *
! 124: * INFO (output) INTEGER
! 125: * = 0: successful exit
! 126: * < 0: if INFO = -i, the i-th argument had an illegal value.
! 127: * = 1,...,N:
! 128: * The QZ iteration failed. (A,B) are not in Schur
! 129: * form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
! 130: * be correct for j=INFO+1,...,N.
! 131: * > N: errors that usually indicate LAPACK problems:
! 132: * =N+1: error return from DGGBAL
! 133: * =N+2: error return from DGEQRF
! 134: * =N+3: error return from DORMQR
! 135: * =N+4: error return from DORGQR
! 136: * =N+5: error return from DGGHRD
! 137: * =N+6: error return from DHGEQZ (other than failed
! 138: * iteration)
! 139: * =N+7: error return from DGGBAK (computing VSL)
! 140: * =N+8: error return from DGGBAK (computing VSR)
! 141: * =N+9: error return from DLASCL (various places)
! 142: *
! 143: * =====================================================================
! 144: *
! 145: * .. Parameters ..
! 146: DOUBLE PRECISION ZERO, ONE
! 147: PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 )
! 148: * ..
! 149: * .. Local Scalars ..
! 150: LOGICAL ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
! 151: INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
! 152: $ IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
! 153: $ LWKOPT, NB, NB1, NB2, NB3
! 154: DOUBLE PRECISION ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
! 155: $ SAFMIN, SMLNUM
! 156: * ..
! 157: * .. External Subroutines ..
! 158: EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
! 159: $ DLASCL, DLASET, DORGQR, DORMQR, XERBLA
! 160: * ..
! 161: * .. External Functions ..
! 162: LOGICAL LSAME
! 163: INTEGER ILAENV
! 164: DOUBLE PRECISION DLAMCH, DLANGE
! 165: EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE
! 166: * ..
! 167: * .. Intrinsic Functions ..
! 168: INTRINSIC INT, MAX
! 169: * ..
! 170: * .. Executable Statements ..
! 171: *
! 172: * Decode the input arguments
! 173: *
! 174: IF( LSAME( JOBVSL, 'N' ) ) THEN
! 175: IJOBVL = 1
! 176: ILVSL = .FALSE.
! 177: ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
! 178: IJOBVL = 2
! 179: ILVSL = .TRUE.
! 180: ELSE
! 181: IJOBVL = -1
! 182: ILVSL = .FALSE.
! 183: END IF
! 184: *
! 185: IF( LSAME( JOBVSR, 'N' ) ) THEN
! 186: IJOBVR = 1
! 187: ILVSR = .FALSE.
! 188: ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
! 189: IJOBVR = 2
! 190: ILVSR = .TRUE.
! 191: ELSE
! 192: IJOBVR = -1
! 193: ILVSR = .FALSE.
! 194: END IF
! 195: *
! 196: * Test the input arguments
! 197: *
! 198: LWKMIN = MAX( 4*N, 1 )
! 199: LWKOPT = LWKMIN
! 200: WORK( 1 ) = LWKOPT
! 201: LQUERY = ( LWORK.EQ.-1 )
! 202: INFO = 0
! 203: IF( IJOBVL.LE.0 ) THEN
! 204: INFO = -1
! 205: ELSE IF( IJOBVR.LE.0 ) THEN
! 206: INFO = -2
! 207: ELSE IF( N.LT.0 ) THEN
! 208: INFO = -3
! 209: ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
! 210: INFO = -5
! 211: ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
! 212: INFO = -7
! 213: ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
! 214: INFO = -12
! 215: ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
! 216: INFO = -14
! 217: ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
! 218: INFO = -16
! 219: END IF
! 220: *
! 221: IF( INFO.EQ.0 ) THEN
! 222: NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
! 223: NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
! 224: NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
! 225: NB = MAX( NB1, NB2, NB3 )
! 226: LOPT = 2*N + N*( NB+1 )
! 227: WORK( 1 ) = LOPT
! 228: END IF
! 229: *
! 230: IF( INFO.NE.0 ) THEN
! 231: CALL XERBLA( 'DGEGS ', -INFO )
! 232: RETURN
! 233: ELSE IF( LQUERY ) THEN
! 234: RETURN
! 235: END IF
! 236: *
! 237: * Quick return if possible
! 238: *
! 239: IF( N.EQ.0 )
! 240: $ RETURN
! 241: *
! 242: * Get machine constants
! 243: *
! 244: EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
! 245: SAFMIN = DLAMCH( 'S' )
! 246: SMLNUM = N*SAFMIN / EPS
! 247: BIGNUM = ONE / SMLNUM
! 248: *
! 249: * Scale A if max element outside range [SMLNUM,BIGNUM]
! 250: *
! 251: ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
! 252: ILASCL = .FALSE.
! 253: IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
! 254: ANRMTO = SMLNUM
! 255: ILASCL = .TRUE.
! 256: ELSE IF( ANRM.GT.BIGNUM ) THEN
! 257: ANRMTO = BIGNUM
! 258: ILASCL = .TRUE.
! 259: END IF
! 260: *
! 261: IF( ILASCL ) THEN
! 262: CALL DLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
! 263: IF( IINFO.NE.0 ) THEN
! 264: INFO = N + 9
! 265: RETURN
! 266: END IF
! 267: END IF
! 268: *
! 269: * Scale B if max element outside range [SMLNUM,BIGNUM]
! 270: *
! 271: BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
! 272: ILBSCL = .FALSE.
! 273: IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
! 274: BNRMTO = SMLNUM
! 275: ILBSCL = .TRUE.
! 276: ELSE IF( BNRM.GT.BIGNUM ) THEN
! 277: BNRMTO = BIGNUM
! 278: ILBSCL = .TRUE.
! 279: END IF
! 280: *
! 281: IF( ILBSCL ) THEN
! 282: CALL DLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
! 283: IF( IINFO.NE.0 ) THEN
! 284: INFO = N + 9
! 285: RETURN
! 286: END IF
! 287: END IF
! 288: *
! 289: * Permute the matrix to make it more nearly triangular
! 290: * Workspace layout: (2*N words -- "work..." not actually used)
! 291: * left_permutation, right_permutation, work...
! 292: *
! 293: ILEFT = 1
! 294: IRIGHT = N + 1
! 295: IWORK = IRIGHT + N
! 296: CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
! 297: $ WORK( IRIGHT ), WORK( IWORK ), IINFO )
! 298: IF( IINFO.NE.0 ) THEN
! 299: INFO = N + 1
! 300: GO TO 10
! 301: END IF
! 302: *
! 303: * Reduce B to triangular form, and initialize VSL and/or VSR
! 304: * Workspace layout: ("work..." must have at least N words)
! 305: * left_permutation, right_permutation, tau, work...
! 306: *
! 307: IROWS = IHI + 1 - ILO
! 308: ICOLS = N + 1 - ILO
! 309: ITAU = IWORK
! 310: IWORK = ITAU + IROWS
! 311: CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
! 312: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
! 313: IF( IINFO.GE.0 )
! 314: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
! 315: IF( IINFO.NE.0 ) THEN
! 316: INFO = N + 2
! 317: GO TO 10
! 318: END IF
! 319: *
! 320: CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
! 321: $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
! 322: $ LWORK+1-IWORK, IINFO )
! 323: IF( IINFO.GE.0 )
! 324: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
! 325: IF( IINFO.NE.0 ) THEN
! 326: INFO = N + 3
! 327: GO TO 10
! 328: END IF
! 329: *
! 330: IF( ILVSL ) THEN
! 331: CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
! 332: CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
! 333: $ VSL( ILO+1, ILO ), LDVSL )
! 334: CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
! 335: $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
! 336: $ IINFO )
! 337: IF( IINFO.GE.0 )
! 338: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
! 339: IF( IINFO.NE.0 ) THEN
! 340: INFO = N + 4
! 341: GO TO 10
! 342: END IF
! 343: END IF
! 344: *
! 345: IF( ILVSR )
! 346: $ CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
! 347: *
! 348: * Reduce to generalized Hessenberg form
! 349: *
! 350: CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
! 351: $ LDVSL, VSR, LDVSR, IINFO )
! 352: IF( IINFO.NE.0 ) THEN
! 353: INFO = N + 5
! 354: GO TO 10
! 355: END IF
! 356: *
! 357: * Perform QZ algorithm, computing Schur vectors if desired
! 358: * Workspace layout: ("work..." must have at least 1 word)
! 359: * left_permutation, right_permutation, work...
! 360: *
! 361: IWORK = ITAU
! 362: CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
! 363: $ ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
! 364: $ WORK( IWORK ), LWORK+1-IWORK, IINFO )
! 365: IF( IINFO.GE.0 )
! 366: $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
! 367: IF( IINFO.NE.0 ) THEN
! 368: IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
! 369: INFO = IINFO
! 370: ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
! 371: INFO = IINFO - N
! 372: ELSE
! 373: INFO = N + 6
! 374: END IF
! 375: GO TO 10
! 376: END IF
! 377: *
! 378: * Apply permutation to VSL and VSR
! 379: *
! 380: IF( ILVSL ) THEN
! 381: CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
! 382: $ WORK( IRIGHT ), N, VSL, LDVSL, IINFO )
! 383: IF( IINFO.NE.0 ) THEN
! 384: INFO = N + 7
! 385: GO TO 10
! 386: END IF
! 387: END IF
! 388: IF( ILVSR ) THEN
! 389: CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
! 390: $ WORK( IRIGHT ), N, VSR, LDVSR, IINFO )
! 391: IF( IINFO.NE.0 ) THEN
! 392: INFO = N + 8
! 393: GO TO 10
! 394: END IF
! 395: END IF
! 396: *
! 397: * Undo scaling
! 398: *
! 399: IF( ILASCL ) THEN
! 400: CALL DLASCL( 'H', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
! 401: IF( IINFO.NE.0 ) THEN
! 402: INFO = N + 9
! 403: RETURN
! 404: END IF
! 405: CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAR, N,
! 406: $ IINFO )
! 407: IF( IINFO.NE.0 ) THEN
! 408: INFO = N + 9
! 409: RETURN
! 410: END IF
! 411: CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAI, N,
! 412: $ IINFO )
! 413: IF( IINFO.NE.0 ) THEN
! 414: INFO = N + 9
! 415: RETURN
! 416: END IF
! 417: END IF
! 418: *
! 419: IF( ILBSCL ) THEN
! 420: CALL DLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
! 421: IF( IINFO.NE.0 ) THEN
! 422: INFO = N + 9
! 423: RETURN
! 424: END IF
! 425: CALL DLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
! 426: IF( IINFO.NE.0 ) THEN
! 427: INFO = N + 9
! 428: RETURN
! 429: END IF
! 430: END IF
! 431: *
! 432: 10 CONTINUE
! 433: WORK( 1 ) = LWKOPT
! 434: *
! 435: RETURN
! 436: *
! 437: * End of DGEGS
! 438: *
! 439: END
CVSweb interface <joel.bertrand@systella.fr>
![[BACK]](/cvsweb/icons/back.gif){kind=icon}
Return to dgegs.f CVS log ![[TXT]](/cvsweb/icons/text.gif){kind=icon}
![[DIR]](/cvsweb/icons/dir.gif){kind=icon}
Up to [local] / rpl / lapack / lapack