LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zgegs.f
Go to the documentation of this file.
00001 *> \brief <b> ZGEEVX computes the eigenvalues and, optionally, the left and/or right eigenvectors for GE matrices</b>
00002 *
00003 *  =========== DOCUMENTATION ===========
00004 *
00005 * Online html documentation available at 
00006 *            http://www.netlib.org/lapack/explore-html/ 
00007 *
00008 *> \htmlonly
00009 *> Download ZGEGS + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgegs.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgegs.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgegs.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA,
00022 *                         VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK,
00023 *                         INFO )
00024 * 
00025 *       .. Scalar Arguments ..
00026 *       CHARACTER          JOBVSL, JOBVSR
00027 *       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
00028 *       ..
00029 *       .. Array Arguments ..
00030 *       DOUBLE PRECISION   RWORK( * )
00031 *       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
00032 *      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
00033 *      $                   WORK( * )
00034 *       ..
00035 *  
00036 *
00037 *> \par Purpose:
00038 *  =============
00039 *>
00040 *> \verbatim
00041 *>
00042 *> This routine is deprecated and has been replaced by routine ZGGES.
00043 *>
00044 *> ZGEGS computes the eigenvalues, Schur form, and, optionally, the
00045 *> left and or/right Schur vectors of a complex matrix pair (A,B).
00046 *> Given two square matrices A and B, the generalized Schur
00047 *> factorization has the form
00048 *> 
00049 *>    A = Q*S*Z**H,  B = Q*T*Z**H
00050 *> 
00051 *> where Q and Z are unitary matrices and S and T are upper triangular.
00052 *> The columns of Q are the left Schur vectors
00053 *> and the columns of Z are the right Schur vectors.
00054 *> 
00055 *> If only the eigenvalues of (A,B) are needed, the driver routine
00056 *> ZGEGV should be used instead.  See ZGEGV for a description of the
00057 *> eigenvalues of the generalized nonsymmetric eigenvalue problem
00058 *> (GNEP).
00059 *> \endverbatim
00060 *
00061 *  Arguments:
00062 *  ==========
00063 *
00064 *> \param[in] JOBVSL
00065 *> \verbatim
00066 *>          JOBVSL is CHARACTER*1
00067 *>          = 'N':  do not compute the left Schur vectors;
00068 *>          = 'V':  compute the left Schur vectors (returned in VSL).
00069 *> \endverbatim
00070 *>
00071 *> \param[in] JOBVSR
00072 *> \verbatim
00073 *>          JOBVSR is CHARACTER*1
00074 *>          = 'N':  do not compute the right Schur vectors;
00075 *>          = 'V':  compute the right Schur vectors (returned in VSR).
00076 *> \endverbatim
00077 *>
00078 *> \param[in] N
00079 *> \verbatim
00080 *>          N is INTEGER
00081 *>          The order of the matrices A, B, VSL, and VSR.  N >= 0.
00082 *> \endverbatim
00083 *>
00084 *> \param[in,out] A
00085 *> \verbatim
00086 *>          A is COMPLEX*16 array, dimension (LDA, N)
00087 *>          On entry, the matrix A.
00088 *>          On exit, the upper triangular matrix S from the generalized
00089 *>          Schur factorization.
00090 *> \endverbatim
00091 *>
00092 *> \param[in] LDA
00093 *> \verbatim
00094 *>          LDA is INTEGER
00095 *>          The leading dimension of A.  LDA >= max(1,N).
00096 *> \endverbatim
00097 *>
00098 *> \param[in,out] B
00099 *> \verbatim
00100 *>          B is COMPLEX*16 array, dimension (LDB, N)
00101 *>          On entry, the matrix B.
00102 *>          On exit, the upper triangular matrix T from the generalized
00103 *>          Schur factorization.
00104 *> \endverbatim
00105 *>
00106 *> \param[in] LDB
00107 *> \verbatim
00108 *>          LDB is INTEGER
00109 *>          The leading dimension of B.  LDB >= max(1,N).
00110 *> \endverbatim
00111 *>
00112 *> \param[out] ALPHA
00113 *> \verbatim
00114 *>          ALPHA is COMPLEX*16 array, dimension (N)
00115 *>          The complex scalars alpha that define the eigenvalues of
00116 *>          GNEP.  ALPHA(j) = S(j,j), the diagonal element of the Schur
00117 *>          form of A.
00118 *> \endverbatim
00119 *>
00120 *> \param[out] BETA
00121 *> \verbatim
00122 *>          BETA is COMPLEX*16 array, dimension (N)
00123 *>          The non-negative real scalars beta that define the
00124 *>          eigenvalues of GNEP.  BETA(j) = T(j,j), the diagonal element
00125 *>          of the triangular factor T.
00126 *>
00127 *>          Together, the quantities alpha = ALPHA(j) and beta = BETA(j)
00128 *>          represent the j-th eigenvalue of the matrix pair (A,B), in
00129 *>          one of the forms lambda = alpha/beta or mu = beta/alpha.
00130 *>          Since either lambda or mu may overflow, they should not,
00131 *>          in general, be computed.
00132 *> \endverbatim
00133 *>
00134 *> \param[out] VSL
00135 *> \verbatim
00136 *>          VSL is COMPLEX*16 array, dimension (LDVSL,N)
00137 *>          If JOBVSL = 'V', the matrix of left Schur vectors Q.
00138 *>          Not referenced if JOBVSL = 'N'.
00139 *> \endverbatim
00140 *>
00141 *> \param[in] LDVSL
00142 *> \verbatim
00143 *>          LDVSL is INTEGER
00144 *>          The leading dimension of the matrix VSL. LDVSL >= 1, and
00145 *>          if JOBVSL = 'V', LDVSL >= N.
00146 *> \endverbatim
00147 *>
00148 *> \param[out] VSR
00149 *> \verbatim
00150 *>          VSR is COMPLEX*16 array, dimension (LDVSR,N)
00151 *>          If JOBVSR = 'V', the matrix of right Schur vectors Z.
00152 *>          Not referenced if JOBVSR = 'N'.
00153 *> \endverbatim
00154 *>
00155 *> \param[in] LDVSR
00156 *> \verbatim
00157 *>          LDVSR is INTEGER
00158 *>          The leading dimension of the matrix VSR. LDVSR >= 1, and
00159 *>          if JOBVSR = 'V', LDVSR >= N.
00160 *> \endverbatim
00161 *>
00162 *> \param[out] WORK
00163 *> \verbatim
00164 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
00165 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00166 *> \endverbatim
00167 *>
00168 *> \param[in] LWORK
00169 *> \verbatim
00170 *>          LWORK is INTEGER
00171 *>          The dimension of the array WORK.  LWORK >= max(1,2*N).
00172 *>          For good performance, LWORK must generally be larger.
00173 *>          To compute the optimal value of LWORK, call ILAENV to get
00174 *>          blocksizes (for ZGEQRF, ZUNMQR, and CUNGQR.)  Then compute:
00175 *>          NB  -- MAX of the blocksizes for ZGEQRF, ZUNMQR, and CUNGQR;
00176 *>          the optimal LWORK is N*(NB+1).
00177 *>
00178 *>          If LWORK = -1, then a workspace query is assumed; the routine
00179 *>          only calculates the optimal size of the WORK array, returns
00180 *>          this value as the first entry of the WORK array, and no error
00181 *>          message related to LWORK is issued by XERBLA.
00182 *> \endverbatim
00183 *>
00184 *> \param[out] RWORK
00185 *> \verbatim
00186 *>          RWORK is DOUBLE PRECISION array, dimension (3*N)
00187 *> \endverbatim
00188 *>
00189 *> \param[out] INFO
00190 *> \verbatim
00191 *>          INFO is INTEGER
00192 *>          = 0:  successful exit
00193 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00194 *>          =1,...,N:
00195 *>                The QZ iteration failed.  (A,B) are not in Schur
00196 *>                form, but ALPHA(j) and BETA(j) should be correct for
00197 *>                j=INFO+1,...,N.
00198 *>          > N:  errors that usually indicate LAPACK problems:
00199 *>                =N+1: error return from ZGGBAL
00200 *>                =N+2: error return from ZGEQRF
00201 *>                =N+3: error return from ZUNMQR
00202 *>                =N+4: error return from ZUNGQR
00203 *>                =N+5: error return from ZGGHRD
00204 *>                =N+6: error return from ZHGEQZ (other than failed
00205 *>                                               iteration)
00206 *>                =N+7: error return from ZGGBAK (computing VSL)
00207 *>                =N+8: error return from ZGGBAK (computing VSR)
00208 *>                =N+9: error return from ZLASCL (various places)
00209 *> \endverbatim
00210 *
00211 *  Authors:
00212 *  ========
00213 *
00214 *> \author Univ. of Tennessee 
00215 *> \author Univ. of California Berkeley 
00216 *> \author Univ. of Colorado Denver 
00217 *> \author NAG Ltd. 
00218 *
00219 *> \date November 2011
00220 *
00221 *> \ingroup complex16GEeigen
00222 *
00223 *  =====================================================================
00224       SUBROUTINE ZGEGS( JOBVSL, JOBVSR, N, A, LDA, B, LDB, ALPHA, BETA,
00225      $                  VSL, LDVSL, VSR, LDVSR, WORK, LWORK, RWORK,
00226      $                  INFO )
00227 *
00228 *  -- LAPACK driver routine (version 3.4.0) --
00229 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00230 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00231 *     November 2011
00232 *
00233 *     .. Scalar Arguments ..
00234       CHARACTER          JOBVSL, JOBVSR
00235       INTEGER            INFO, LDA, LDB, LDVSL, LDVSR, LWORK, N
00236 *     ..
00237 *     .. Array Arguments ..
00238       DOUBLE PRECISION   RWORK( * )
00239       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
00240      $                   BETA( * ), VSL( LDVSL, * ), VSR( LDVSR, * ),
00241      $                   WORK( * )
00242 *     ..
00243 *
00244 *  =====================================================================
00245 *
00246 *     .. Parameters ..
00247       DOUBLE PRECISION   ZERO, ONE
00248       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00249       COMPLEX*16         CZERO, CONE
00250       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
00251      $                   CONE = ( 1.0D0, 0.0D0 ) )
00252 *     ..
00253 *     .. Local Scalars ..
00254       LOGICAL            ILASCL, ILBSCL, ILVSL, ILVSR, LQUERY
00255       INTEGER            ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO,
00256      $                   IRIGHT, IROWS, IRWORK, ITAU, IWORK, LOPT,
00257      $                   LWKMIN, LWKOPT, NB, NB1, NB2, NB3
00258       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
00259      $                   SAFMIN, SMLNUM
00260 *     ..
00261 *     .. External Subroutines ..
00262       EXTERNAL           XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD, ZHGEQZ,
00263      $                   ZLACPY, ZLASCL, ZLASET, ZUNGQR, ZUNMQR
00264 *     ..
00265 *     .. External Functions ..
00266       LOGICAL            LSAME
00267       INTEGER            ILAENV
00268       DOUBLE PRECISION   DLAMCH, ZLANGE
00269       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
00270 *     ..
00271 *     .. Intrinsic Functions ..
00272       INTRINSIC          INT, MAX
00273 *     ..
00274 *     .. Executable Statements ..
00275 *
00276 *     Decode the input arguments
00277 *
00278       IF( LSAME( JOBVSL, 'N' ) ) THEN
00279          IJOBVL = 1
00280          ILVSL = .FALSE.
00281       ELSE IF( LSAME( JOBVSL, 'V' ) ) THEN
00282          IJOBVL = 2
00283          ILVSL = .TRUE.
00284       ELSE
00285          IJOBVL = -1
00286          ILVSL = .FALSE.
00287       END IF
00288 *
00289       IF( LSAME( JOBVSR, 'N' ) ) THEN
00290          IJOBVR = 1
00291          ILVSR = .FALSE.
00292       ELSE IF( LSAME( JOBVSR, 'V' ) ) THEN
00293          IJOBVR = 2
00294          ILVSR = .TRUE.
00295       ELSE
00296          IJOBVR = -1
00297          ILVSR = .FALSE.
00298       END IF
00299 *
00300 *     Test the input arguments
00301 *
00302       LWKMIN = MAX( 2*N, 1 )
00303       LWKOPT = LWKMIN
00304       WORK( 1 ) = LWKOPT
00305       LQUERY = ( LWORK.EQ.-1 )
00306       INFO = 0
00307       IF( IJOBVL.LE.0 ) THEN
00308          INFO = -1
00309       ELSE IF( IJOBVR.LE.0 ) THEN
00310          INFO = -2
00311       ELSE IF( N.LT.0 ) THEN
00312          INFO = -3
00313       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00314          INFO = -5
00315       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00316          INFO = -7
00317       ELSE IF( LDVSL.LT.1 .OR. ( ILVSL .AND. LDVSL.LT.N ) ) THEN
00318          INFO = -11
00319       ELSE IF( LDVSR.LT.1 .OR. ( ILVSR .AND. LDVSR.LT.N ) ) THEN
00320          INFO = -13
00321       ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
00322          INFO = -15
00323       END IF
00324 *
00325       IF( INFO.EQ.0 ) THEN
00326          NB1 = ILAENV( 1, 'ZGEQRF', ' ', N, N, -1, -1 )
00327          NB2 = ILAENV( 1, 'ZUNMQR', ' ', N, N, N, -1 )
00328          NB3 = ILAENV( 1, 'ZUNGQR', ' ', N, N, N, -1 )
00329          NB = MAX( NB1, NB2, NB3 )
00330          LOPT = N*( NB+1 )
00331          WORK( 1 ) = LOPT
00332       END IF
00333 *
00334       IF( INFO.NE.0 ) THEN
00335          CALL XERBLA( 'ZGEGS ', -INFO )
00336          RETURN
00337       ELSE IF( LQUERY ) THEN
00338          RETURN
00339       END IF
00340 *
00341 *     Quick return if possible
00342 *
00343       IF( N.EQ.0 )
00344      $   RETURN
00345 *
00346 *     Get machine constants
00347 *
00348       EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
00349       SAFMIN = DLAMCH( 'S' )
00350       SMLNUM = N*SAFMIN / EPS
00351       BIGNUM = ONE / SMLNUM
00352 *
00353 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00354 *
00355       ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
00356       ILASCL = .FALSE.
00357       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00358          ANRMTO = SMLNUM
00359          ILASCL = .TRUE.
00360       ELSE IF( ANRM.GT.BIGNUM ) THEN
00361          ANRMTO = BIGNUM
00362          ILASCL = .TRUE.
00363       END IF
00364 *
00365       IF( ILASCL ) THEN
00366          CALL ZLASCL( 'G', -1, -1, ANRM, ANRMTO, N, N, A, LDA, IINFO )
00367          IF( IINFO.NE.0 ) THEN
00368             INFO = N + 9
00369             RETURN
00370          END IF
00371       END IF
00372 *
00373 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00374 *
00375       BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
00376       ILBSCL = .FALSE.
00377       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00378          BNRMTO = SMLNUM
00379          ILBSCL = .TRUE.
00380       ELSE IF( BNRM.GT.BIGNUM ) THEN
00381          BNRMTO = BIGNUM
00382          ILBSCL = .TRUE.
00383       END IF
00384 *
00385       IF( ILBSCL ) THEN
00386          CALL ZLASCL( 'G', -1, -1, BNRM, BNRMTO, N, N, B, LDB, IINFO )
00387          IF( IINFO.NE.0 ) THEN
00388             INFO = N + 9
00389             RETURN
00390          END IF
00391       END IF
00392 *
00393 *     Permute the matrix to make it more nearly triangular
00394 *
00395       ILEFT = 1
00396       IRIGHT = N + 1
00397       IRWORK = IRIGHT + N
00398       IWORK = 1
00399       CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
00400      $             RWORK( IRIGHT ), RWORK( IRWORK ), IINFO )
00401       IF( IINFO.NE.0 ) THEN
00402          INFO = N + 1
00403          GO TO 10
00404       END IF
00405 *
00406 *     Reduce B to triangular form, and initialize VSL and/or VSR
00407 *
00408       IROWS = IHI + 1 - ILO
00409       ICOLS = N + 1 - ILO
00410       ITAU = IWORK
00411       IWORK = ITAU + IROWS
00412       CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00413      $             WORK( IWORK ), LWORK+1-IWORK, IINFO )
00414       IF( IINFO.GE.0 )
00415      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00416       IF( IINFO.NE.0 ) THEN
00417          INFO = N + 2
00418          GO TO 10
00419       END IF
00420 *
00421       CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00422      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ),
00423      $             LWORK+1-IWORK, IINFO )
00424       IF( IINFO.GE.0 )
00425      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00426       IF( IINFO.NE.0 ) THEN
00427          INFO = N + 3
00428          GO TO 10
00429       END IF
00430 *
00431       IF( ILVSL ) THEN
00432          CALL ZLASET( 'Full', N, N, CZERO, CONE, VSL, LDVSL )
00433          CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00434      $                VSL( ILO+1, ILO ), LDVSL )
00435          CALL ZUNGQR( IROWS, IROWS, IROWS, VSL( ILO, ILO ), LDVSL,
00436      $                WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK,
00437      $                IINFO )
00438          IF( IINFO.GE.0 )
00439      $      LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00440          IF( IINFO.NE.0 ) THEN
00441             INFO = N + 4
00442             GO TO 10
00443          END IF
00444       END IF
00445 *
00446       IF( ILVSR )
00447      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, VSR, LDVSR )
00448 *
00449 *     Reduce to generalized Hessenberg form
00450 *
00451       CALL ZGGHRD( JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB, VSL,
00452      $             LDVSL, VSR, LDVSR, IINFO )
00453       IF( IINFO.NE.0 ) THEN
00454          INFO = N + 5
00455          GO TO 10
00456       END IF
00457 *
00458 *     Perform QZ algorithm, computing Schur vectors if desired
00459 *
00460       IWORK = ITAU
00461       CALL ZHGEQZ( 'S', JOBVSL, JOBVSR, N, ILO, IHI, A, LDA, B, LDB,
00462      $             ALPHA, BETA, VSL, LDVSL, VSR, LDVSR, WORK( IWORK ),
00463      $             LWORK+1-IWORK, RWORK( IRWORK ), IINFO )
00464       IF( IINFO.GE.0 )
00465      $   LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 )
00466       IF( IINFO.NE.0 ) THEN
00467          IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN
00468             INFO = IINFO
00469          ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN
00470             INFO = IINFO - N
00471          ELSE
00472             INFO = N + 6
00473          END IF
00474          GO TO 10
00475       END IF
00476 *
00477 *     Apply permutation to VSL and VSR
00478 *
00479       IF( ILVSL ) THEN
00480          CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
00481      $                RWORK( IRIGHT ), N, VSL, LDVSL, IINFO )
00482          IF( IINFO.NE.0 ) THEN
00483             INFO = N + 7
00484             GO TO 10
00485          END IF
00486       END IF
00487       IF( ILVSR ) THEN
00488          CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
00489      $                RWORK( IRIGHT ), N, VSR, LDVSR, IINFO )
00490          IF( IINFO.NE.0 ) THEN
00491             INFO = N + 8
00492             GO TO 10
00493          END IF
00494       END IF
00495 *
00496 *     Undo scaling
00497 *
00498       IF( ILASCL ) THEN
00499          CALL ZLASCL( 'U', -1, -1, ANRMTO, ANRM, N, N, A, LDA, IINFO )
00500          IF( IINFO.NE.0 ) THEN
00501             INFO = N + 9
00502             RETURN
00503          END IF
00504          CALL ZLASCL( 'G', -1, -1, ANRMTO, ANRM, N, 1, ALPHA, N, IINFO )
00505          IF( IINFO.NE.0 ) THEN
00506             INFO = N + 9
00507             RETURN
00508          END IF
00509       END IF
00510 *
00511       IF( ILBSCL ) THEN
00512          CALL ZLASCL( 'U', -1, -1, BNRMTO, BNRM, N, N, B, LDB, IINFO )
00513          IF( IINFO.NE.0 ) THEN
00514             INFO = N + 9
00515             RETURN
00516          END IF
00517          CALL ZLASCL( 'G', -1, -1, BNRMTO, BNRM, N, 1, BETA, N, IINFO )
00518          IF( IINFO.NE.0 ) THEN
00519             INFO = N + 9
00520             RETURN
00521          END IF
00522       END IF
00523 *
00524    10 CONTINUE
00525       WORK( 1 ) = LWKOPT
00526 *
00527       RETURN
00528 *
00529 *     End of ZGEGS
00530 *
00531       END
 All Files Functions