Annotation of rpl/lapack/lapack/dgges3.f, revision 1.3

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>