LAPACK: SRC/dgegs.f Source File

Go to the documentation of this file.
00001       SUBROUTINE DGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHAR,
00002      $                  ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR, WORK,
00003      $                  LWORK, INFO )
00004 *
00005 *  -- LAPACK driver routine (version 3.2) --
00006 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00007 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00008 *     November 2006
00009 *
00010 *     .. Scalar Arguments ..
00011       CHARACTER          JOBVSL, JOBVSR
00012       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
00013 *     ..
00014 *     .. Array Arguments ..
00015       DOUBLE PRECISION   A( LDA, * ), ALPHAI( * ), ALPHAR( * ),
00016      $                   B( LDB, * ), BETA( * ), VSL( LDVSL, * ),
00017      $                   VSR( LDVSR, * ), WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  This routine is deprecated and has been replaced by routine DGGES.
00024 *
00025 *  DGEGS computes the eigenvalues, real Schur form, and, optionally,
00026 *  left and or/right Schur vectors of a real matrix pair (A,B).
00027 *  Given two square matrices A and B, the generalized real Schur
00028 *  factorization has the form
00029 *
00030 *    A = Q*S*Z**T,  B = Q*T*Z**T
00031 *
00032 *  where Q and Z are orthogonal matrices, T is upper triangular, and S
00033 *  is an upper quasi-triangular matrix with 1-by-1 and 2-by-2 diagonal
00034 *  blocks, the 2-by-2 blocks corresponding to complex conjugate pairs
00035 *  of eigenvalues of (A,B).  The columns of Q are the left Schur vectors
00036 *  and the columns of Z are the right Schur vectors.
00037 *
00038 *  If only the eigenvalues of (A,B) are needed, the driver routine
00039 *  DGEGV should be used instead.  See DGEGV for a description of the
00040 *  eigenvalues of the generalized nonsymmetric eigenvalue problem
00041 *  (GNEP).
00042 *
00043 *  Arguments
00044 *  =========
00045 *
00046 *  JOBVSL  (input) CHARACTER*1
00047 *          = 'N':  do not compute the left Schur vectors;
00048 *          = 'V':  compute the left Schur vectors (returned in VSL).
00049 *
00050 *  JOBVSR  (input) CHARACTER*1
00051 *          = 'N':  do not compute the right Schur vectors;
00052 *          = 'V':  compute the right Schur vectors (returned in VSR).
00053 *
00054 *  N       (input) INTEGER
00055 *          The order of the matrices A, B, VSL, and VSR.  N >= 0.
00056 *
00057 *  A       (input/output) DOUBLE PRECISION array, dimension (LDA, N)
00058 *          On entry, the matrix A.
00059 *          On exit, the upper quasi-triangular matrix S from the
00060 *          generalized real Schur factorization.
00061 *
00062 *  LDA     (input) INTEGER
00063 *          The leading dimension of A.  LDA >= max(1,N).
00064 *
00065 *  B       (input/output) DOUBLE PRECISION array, dimension (LDB, N)
00066 *          On entry, the matrix B.
00067 *          On exit, the upper triangular matrix T from the generalized
00068 *          real Schur factorization.
00069 *
00070 *  LDB     (input) INTEGER
00071 *          The leading dimension of B.  LDB >= max(1,N).
00072 *
00073 *  ALPHAR  (output) DOUBLE PRECISION array, dimension (N)
00074 *          The real parts of each scalar alpha defining an eigenvalue
00075 *          of GNEP.
00076 *
00077 *  ALPHAI  (output) DOUBLE PRECISION array, dimension (N)
00078 *          The imaginary parts of each scalar alpha defining an
00079 *          eigenvalue of GNEP.  If ALPHAI(j) is zero, then the j-th
00080 *          eigenvalue is real; if positive, then the j-th and (j+1)-st
00081 *          eigenvalues are a complex conjugate pair, with
00082 *          ALPHAI(j+1) = -ALPHAI(j).
00083 *
00084 *  BETA    (output) DOUBLE PRECISION array, dimension (N)
00085 *          The scalars beta that define the eigenvalues of GNEP.
00086 *          Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and
00087 *          beta = BETA(j) represent the j-th eigenvalue of the matrix
00088 *          pair (A,B), in one of the forms lambda = alpha/beta or
00089 *          mu = beta/alpha.  Since either lambda or mu may overflow,
00090 *          they should not, in general, be computed.
00091 *
00092 *  VSL     (output) DOUBLE PRECISION array, dimension (LDVSL,N)
00093 *          If JOBVSL = 'V', the matrix of left Schur vectors Q.
00094 *          Not referenced if JOBVSL = 'N'.
00095 *
00096 *  LDVSL   (input) INTEGER
00097 *          The leading dimension of the matrix VSL. LDVSL >=1, and
00098 *          if JOBVSL = 'V', LDVSL >= N.
00099 *
00100 *  VSR     (output) DOUBLE PRECISION array, dimension (LDVSR,N)
00101 *          If JOBVSR = 'V', the matrix of right Schur vectors Z.
00102 *          Not referenced if JOBVSR = 'N'.
00103 *
00104 *  LDVSR   (input) INTEGER
00105 *          The leading dimension of the matrix VSR. LDVSR >= 1, and
00106 *          if JOBVSR = 'V', LDVSR >= N.
00107 *
00108 *  WORK    (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK))
00109 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00110 *
00111 *  LWORK   (input) INTEGER
00112 *          The dimension of the array WORK.  LWORK >= max(1,4*N).
00113 *          For good performance, LWORK must generally be larger.
00114 *          To compute the optimal value of LWORK, call ILAENV to get
00115 *          blocksizes (for DGEQRF, DORMQR, and DORGQR.)  Then compute:
00116 *          NB  -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR
00117 *          The optimal LWORK is  2*N + N*(NB+1).
00118 *
00119 *          If LWORK = -1, then a workspace query is assumed; the routine
00120 *          only calculates the optimal size of the WORK array, returns
00121 *          this value as the first entry of the WORK array, and no error
00122 *          message related to LWORK is issued by XERBLA.
00123 *
00124 *  INFO    (output) INTEGER
00125 *          = 0:  successful exit
00126 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00127 *          = 1,...,N:
00128 *                The QZ iteration failed.  (A,B) are not in Schur
00129 *                form, but ALPHAR(j), ALPHAI(j), and BETA(j) should
00130 *                be correct for j=INFO+1,...,N.
00131 *          > N:  errors that usually indicate LAPACK problems:
00132 *                =N+1: error return from DGGBAL
00133 *                =N+2: error return from DGEQRF
00134 *                =N+3: error return from DORMQR
00135 *                =N+4: error return from DORGQR
00136 *                =N+5: error return from DGGHRD
00137 *                =N+6: error return from DHGEQZ (other than failed
00138 *                                                iteration)
00139 *                =N+7: error return from DGGBAK (computing VSL)
00140 *                =N+8: error return from DGGBAK (computing VSR)
00141 *                =N+9: error return from DLASCL (various places)
00142 *
00143 *  =====================================================================
00144 *
00145 *     .. Parameters ..
00146       DOUBLE PRECISION   ZERO, ONE
00147       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00148 *     ..
00149 *     .. Local Scalars ..
00150       LOGICAL            ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
00151       INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
00152      $                   IRIGHT, IROWS, ITAU, IWORK, LOPT, LWKMIN,
00153      $                   LWKOPT, NB, NB1, NB2, NB3
00154       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
00155      $                   SAFMIN, SMLNUM
00156 *     ..
00157 *     .. External Subroutines ..
00158       EXTERNAL           DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY,
00159      $                   DLASCL, DLASET, DORGQR, DORMQR, XERBLA
00160 *     ..
00161 *     .. External Functions ..
00162       LOGICAL            LSAME
00163       INTEGER            ILAENV
00164       DOUBLE PRECISION   DLAMCH, DLANGE
00165       EXTERNAL           LSAME, ILAENV, DLAMCH, DLANGE
00166 *     ..
00167 *     .. Intrinsic Functions ..
00168       INTRINSIC          INT, MAX
00169 *     ..
00170 *     .. Executable Statements ..
00171 *
00172 *     Decode the input arguments
00173 *
00174       IF( LSAME( JOBVSL, 'N' ) ) THEN
00175          IJOBVL = 1
00176          ILVSL = .FALSE.
00177       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
00178          IJOBVL = 2
00179          ILVSL = .TRUE.
00180       ELSE
00181          IJOBVL = -1
00182          ILVSL = .FALSE.
00183       END IF
00184 *
00185       IF( LSAME( JOBVSR, 'N' ) ) THEN
00186          IJOBVR = 1
00187          ILVSR = .FALSE.
00188       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
00189          IJOBVR = 2
00190          ILVSR = .TRUE.
00191       ELSE
00192          IJOBVR = -1
00193          ILVSR = .FALSE.
00194       END IF
00195 *
00196 *     Test the input arguments
00197 *
00198       LWKMIN = MAX( 4*N, 1 )
00199       LWKOPT = LWKMIN
00200       WORK( 1 ) = LWKOPT
00201       LQUERY = ( LWORK.EQ.-1 )
00202       INFO = 0
00203       IF( IJOBVL.LE.0 ) THEN
00204          INFO = -1
00205       ELSE IF( IJOBVR.LE.0 ) THEN
00206          INFO = -2
00207       ELSE IF( N.LT.0 ) THEN
00208          INFO = -3
00209       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00210          INFO = -5
00211       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00212          INFO = -7
00213       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
00214          INFO = -12
00215       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
00216          INFO = -14
00217       ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00218          INFO = -16
00219       END IF
00220 *
00221       IF( INFO.EQ.0 ) THEN
00222          NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 )
00223          NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 )
00224          NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 )
00225          NB = MAX( NB1, NB2, NB3 )
00226          LOPT = 2*N + N*( NB+1 )
00227          WORK( 1 ) = LOPT
00228       END IF
00229 *
00230       IF( INFO.NE.0 ) THEN
00231          CALL XERBLA( 'DGEGS ', -INFO )
00232          RETURN
00233       ELSE IF( LQUERY ) THEN
00234          RETURN
00235       END IF
00236 *
00237 *     Quick return if possible
00238 *
00239       IF( N.EQ.0 )
00240      $   RETURN
00241 *
00242 *     Get machine constants
00243 *
00244       EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
00245       SAFMIN = DLAMCH( 'S' )
00246       SMLNUM = N*SAFMIN / EPS
00247       BIGNUM = ONE / SMLNUM
00248 *
00249 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00250 *
00251       ANRM = DLANGE( 'M', N, N, A, LDA, WORK )
00252       ILASCL = .FALSE.
00253       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00254          ANRMTO = SMLNUM
00255          ILASCL = .TRUE.
00256       ELSE IF( ANRM.GT.BIGNUM ) THEN
00257          ANRMTO = BIGNUM
00258          ILASCL = .TRUE.
00259       END IF
00260 *
00261       IF( ILASCL ) THEN
00262          CALL DLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
00263          IF( IINFO.NE.0 ) THEN
00264             INFO = N + 9
00265             RETURN
00266          END IF
00267       END IF
00268 *
00269 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00270 *
00271       BNRM = DLANGE( 'M', N, N, B, LDB, WORK )
00272       ILBSCL = .FALSE.
00273       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00274          BNRMTO = SMLNUM
00275          ILBSCL = .TRUE.
00276       ELSE IF( BNRM.GT.BIGNUM ) THEN
00277          BNRMTO = BIGNUM
00278          ILBSCL = .TRUE.
00279       END IF
00280 *
00281       IF( ILBSCL ) THEN
00282          CALL DLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
00283          IF( IINFO.NE.0 ) THEN
00284             INFO = N + 9
00285             RETURN
00286          END IF
00287       END IF
00288 *
00289 *     Permute the matrix to make it more nearly triangular
00290 *     Workspace layout:  (2*N words -- "work..." not actually used)
00291 *        left_permutation, right_permutation, work...
00292 *
00293       ILEFT = 1
00294       IRIGHT = N + 1
00295       IWORK = IRIGHT + N
00296       CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ),
00297      $             WORK( IRIGHT ), WORK( IWORK ), IINFO )
00298       IF( IINFO.NE.0 ) THEN
00299          INFO = N + 1
00300          GO TO 10
00301       END IF
00302 *
00303 *     Reduce B to triangular form, and initialize VSL and/or VSR
00304 *     Workspace layout:  ("work..." must have at least N words)
00305 *        left_permutation, right_permutation, tau, work...
00306 *
00307       IROWS = IHI + 1 - ILO
00308       ICOLS = N + 1 - ILO
00309       ITAU = IWORK
00310       IWORK = ITAU + IROWS
00311       CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00312      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
00313       IF( IINFO.GE.0 )
00314      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00315       IF( IINFO.NE.0 ) THEN
00316          INFO = N + 2
00317          GO TO 10
00318       END IF
00319 *
00320       CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00321      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
00322      $             LWORK+1-IWORK, IINFO )
00323       IF( IINFO.GE.0 )
00324      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00325       IF( IINFO.NE.0 ) THEN
00326          INFO = N + 3
00327          GO TO 10
00328       END IF
00329 *
00330       IF( ILVSL ) THEN
00331          CALL DLASET( 'Full', N, N, ZERO, ONE, VSL, LDVSL )
00332          CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00333      $                VSL( ILO+1, ILO ), LDVSL )
00334          CALL DORGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
00335      $                WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
00336      $                IINFO )
00337          IF( IINFO.GE.0 )
00338      $      LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00339          IF( IINFO.NE.0 ) THEN
00340             INFO = N + 4
00341             GO TO 10
00342          END IF
00343       END IF
00344 *
00345       IF( ILVSR )
00346      $   CALL DLASET( 'Full', N, N, ZERO, ONE, VSR, LDVSR )
00347 *
00348 *     Reduce to generalized Hessenberg form
00349 *
00350       CALL DGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
00351      $             LDVSL, VSR, LDVSR, IINFO )
00352       IF( IINFO.NE.0 ) THEN
00353          INFO = N + 5
00354          GO TO 10
00355       END IF
00356 *
00357 *     Perform QZ algorithm, computing Schur vectors if desired
00358 *     Workspace layout:  ("work..." must have at least 1 word)
00359 *        left_permutation, right_permutation, work...
00360 *
00361       IWORK = ITAU
00362       CALL DHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
00363      $             ALPHAR, ALPHAI, BETA, VSL, LDVSL, VSR, LDVSR,
00364      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
00365       IF( IINFO.GE.0 )
00366      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00367       IF( IINFO.NE.0 ) THEN
00368          IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
00369             INFO = IINFO
00370          ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
00371             INFO = IINFO - N
00372          ELSE
00373             INFO = N + 6
00374          END IF
00375          GO TO 10
00376       END IF
00377 *
00378 *     Apply permutation to VSL and VSR
00379 *
00380       IF( ILVSL ) THEN
00381          CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ),
00382      $                WORK( IRIGHT ), N, VSL, LDVSL, IINFO )
00383          IF( IINFO.NE.0 ) THEN
00384             INFO = N + 7
00385             GO TO 10
00386          END IF
00387       END IF
00388       IF( ILVSR ) THEN
00389          CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ),
00390      $                WORK( IRIGHT ), N, VSR, LDVSR, IINFO )
00391          IF( IINFO.NE.0 ) THEN
00392             INFO = N + 8
00393             GO TO 10
00394          END IF
00395       END IF
00396 *
00397 *     Undo scaling
00398 *
00399       IF( ILASCL ) THEN
00400          CALL DLASCL( 'H', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
00401          IF( IINFO.NE.0 ) THEN
00402             INFO = N + 9
00403             RETURN
00404          END IF
00405          CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAR, N,
00406      $                IINFO )
00407          IF( IINFO.NE.0 ) THEN
00408             INFO = N + 9
00409             RETURN
00410          END IF
00411          CALL DLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHAI, N,
00412      $                IINFO )
00413          IF( IINFO.NE.0 ) THEN
00414             INFO = N + 9
00415             RETURN
00416          END IF
00417       END IF
00418 *
00419       IF( ILBSCL ) THEN
00420          CALL DLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
00421          IF( IINFO.NE.0 ) THEN
00422             INFO = N + 9
00423             RETURN
00424          END IF
00425          CALL DLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
00426          IF( IINFO.NE.0 ) THEN
00427             INFO = N + 9
00428             RETURN
00429          END IF
00430       END IF
00431 *
00432    10 CONTINUE
00433       WORK( 1 ) = LWKOPT
00434 *
00435       RETURN
00436 *
00437 *     End of DGEGS
00438 *
00439       END