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