LAPACK: SRC/zggev.f Source File

Go to the documentation of this file.
00001       SUBROUTINE ZGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
00002      $                  VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
00003 *
00004 *  -- LAPACK driver routine (version 3.2) --
00005 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00006 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00007 *     November 2006
00008 *
00009 *     .. Scalar Arguments ..
00010       CHARACTER          JOBVL, JOBVR
00011       INTEGER            INFO, LDA, LDB, LDVL, LDVR, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       DOUBLE PRECISION   RWORK( * )
00015       COMPLEX*16         A( LDA, * ), ALPHA( * ), B( LDB, * ),
00016      $                   BETA( * ), VL( LDVL, * ), VR( LDVR, * ),
00017      $                   WORK( * )
00018 *     ..
00019 *
00020 *  Purpose
00021 *  =======
00022 *
00023 *  ZGGEV computes for a pair of N-by-N complex nonsymmetric matrices
00024 *  (A,B), the generalized eigenvalues, and optionally, the left and/or
00025 *  right generalized eigenvectors.
00026 *
00027 *  A generalized eigenvalue for a pair of matrices (A,B) is a scalar
00028 *  lambda or a ratio alpha/beta = lambda, such that A - lambda*B is
00029 *  singular. It is usually represented as the pair (alpha,beta), as
00030 *  there is a reasonable interpretation for beta=0, and even for both
00031 *  being zero.
00032 *
00033 *  The right generalized eigenvector v(j) corresponding to the
00034 *  generalized eigenvalue lambda(j) of (A,B) satisfies
00035 *
00036 *               A * v(j) = lambda(j) * B * v(j).
00037 *
00038 *  The left generalized eigenvector u(j) corresponding to the
00039 *  generalized eigenvalues lambda(j) of (A,B) satisfies
00040 *
00041 *               u(j)**H * A = lambda(j) * u(j)**H * B
00042 *
00043 *  where u(j)**H is the conjugate-transpose of u(j).
00044 *
00045 *  Arguments
00046 *  =========
00047 *
00048 *  JOBVL   (input) CHARACTER*1
00049 *          = 'N':  do not compute the left generalized eigenvectors;
00050 *          = 'V':  compute the left generalized eigenvectors.
00051 *
00052 *  JOBVR   (input) CHARACTER*1
00053 *          = 'N':  do not compute the right generalized eigenvectors;
00054 *          = 'V':  compute the right generalized eigenvectors.
00055 *
00056 *  N       (input) INTEGER
00057 *          The order of the matrices A, B, VL, and VR.  N >= 0.
00058 *
00059 *  A       (input/output) COMPLEX*16 array, dimension (LDA, N)
00060 *          On entry, the matrix A in the pair (A,B).
00061 *          On exit, A has been overwritten.
00062 *
00063 *  LDA     (input) INTEGER
00064 *          The leading dimension of A.  LDA >= max(1,N).
00065 *
00066 *  B       (input/output) COMPLEX*16 array, dimension (LDB, N)
00067 *          On entry, the matrix B in the pair (A,B).
00068 *          On exit, B has been overwritten.
00069 *
00070 *  LDB     (input) INTEGER
00071 *          The leading dimension of B.  LDB >= max(1,N).
00072 *
00073 *  ALPHA   (output) COMPLEX*16 array, dimension (N)
00074 *  BETA    (output) COMPLEX*16 array, dimension (N)
00075 *          On exit, ALPHA(j)/BETA(j), j=1,...,N, will be the
00076 *          generalized eigenvalues.
00077 *
00078 *          Note: the quotients ALPHA(j)/BETA(j) may easily over- or
00079 *          underflow, and BETA(j) may even be zero.  Thus, the user
00080 *          should avoid naively computing the ratio alpha/beta.
00081 *          However, ALPHA will be always less than and usually
00082 *          comparable with norm(A) in magnitude, and BETA always less
00083 *          than and usually comparable with norm(B).
00084 *
00085 *  VL      (output) COMPLEX*16 array, dimension (LDVL,N)
00086 *          If JOBVL = 'V', the left generalized eigenvectors u(j) are
00087 *          stored one after another in the columns of VL, in the same
00088 *          order as their eigenvalues.
00089 *          Each eigenvector is scaled so the largest component has
00090 *          abs(real part) + abs(imag. part) = 1.
00091 *          Not referenced if JOBVL = 'N'.
00092 *
00093 *  LDVL    (input) INTEGER
00094 *          The leading dimension of the matrix VL. LDVL >= 1, and
00095 *          if JOBVL = 'V', LDVL >= N.
00096 *
00097 *  VR      (output) COMPLEX*16 array, dimension (LDVR,N)
00098 *          If JOBVR = 'V', the right generalized eigenvectors v(j) are
00099 *          stored one after another in the columns of VR, in the same
00100 *          order as their eigenvalues.
00101 *          Each eigenvector is scaled so the largest component has
00102 *          abs(real part) + abs(imag. part) = 1.
00103 *          Not referenced if JOBVR = 'N'.
00104 *
00105 *  LDVR    (input) INTEGER
00106 *          The leading dimension of the matrix VR. LDVR >= 1, and
00107 *          if JOBVR = 'V', LDVR >= N.
00108 *
00109 *  WORK    (workspace/output) COMPLEX*16 array, dimension (MAX(1,LWORK))
00110 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00111 *
00112 *  LWORK   (input) INTEGER
00113 *          The dimension of the array WORK.  LWORK >= max(1,2*N).
00114 *          For good performance, LWORK must generally be larger.
00115 *
00116 *          If LWORK = -1, then a workspace query is assumed; the routine
00117 *          only calculates the optimal size of the WORK array, returns
00118 *          this value as the first entry of the WORK array, and no error
00119 *          message related to LWORK is issued by XERBLA.
00120 *
00121 *  RWORK   (workspace/output) DOUBLE PRECISION array, dimension (8*N)
00122 *
00123 *  INFO    (output) INTEGER
00124 *          = 0:  successful exit
00125 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00126 *          =1,...,N:
00127 *                The QZ iteration failed.  No eigenvectors have been
00128 *                calculated, but ALPHA(j) and BETA(j) should be
00129 *                correct for j=INFO+1,...,N.
00130 *          > N:  =N+1: other then QZ iteration failed in DHGEQZ,
00131 *                =N+2: error return from DTGEVC.
00132 *
00133 *  =====================================================================
00134 *
00135 *     .. Parameters ..
00136       DOUBLE PRECISION   ZERO, ONE
00137       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00138       COMPLEX*16         CZERO, CONE
00139       PARAMETER          ( CZERO = ( 0.0D0, 0.0D0 ),
00140      $                   CONE = ( 1.0D0, 0.0D0 ) )
00141 *     ..
00142 *     .. Local Scalars ..
00143       LOGICAL            ILASCL, ILBSCL, ILV, ILVL, ILVR, LQUERY
00144       CHARACTER          CHTEMP
00145       INTEGER            ICOLS, IERR, IHI, IJOBVL, IJOBVR, ILEFT, ILO,
00146      $                   IN, IRIGHT, IROWS, IRWRK, ITAU, IWRK, JC, JR,
00147      $                   LWKMIN, LWKOPT
00148       DOUBLE PRECISION   ANRM, ANRMTO, BIGNUM, BNRM, BNRMTO, EPS,
00149      $                   SMLNUM, TEMP
00150       COMPLEX*16         X
00151 *     ..
00152 *     .. Local Arrays ..
00153       LOGICAL            LDUMMA( 1 )
00154 *     ..
00155 *     .. External Subroutines ..
00156       EXTERNAL           DLABAD, XERBLA, ZGEQRF, ZGGBAK, ZGGBAL, ZGGHRD,
00157      $                   ZHGEQZ, ZLACPY, ZLASCL, ZLASET, ZTGEVC, ZUNGQR,
00158      $                   ZUNMQR
00159 *     ..
00160 *     .. External Functions ..
00161       LOGICAL            LSAME
00162       INTEGER            ILAENV
00163       DOUBLE PRECISION   DLAMCH, ZLANGE
00164       EXTERNAL           LSAME, ILAENV, DLAMCH, ZLANGE
00165 *     ..
00166 *     .. Intrinsic Functions ..
00167       INTRINSIC          ABS, DBLE, DIMAG, MAX, SQRT
00168 *     ..
00169 *     .. Statement Functions ..
00170       DOUBLE PRECISION   ABS1
00171 *     ..
00172 *     .. Statement Function definitions ..
00173       ABS1( X ) = ABS( DBLE( X ) ) + ABS( DIMAG( X ) )
00174 *     ..
00175 *     .. Executable Statements ..
00176 *
00177 *     Decode the input arguments
00178 *
00179       IF( LSAME( JOBVL, 'N' ) ) THEN
00180          IJOBVL = 1
00181          ILVL = .FALSE.
00182       ELSE IF( LSAME( JOBVL, 'V' ) ) THEN
00183          IJOBVL = 2
00184          ILVL = .TRUE.
00185       ELSE
00186          IJOBVL = -1
00187          ILVL = .FALSE.
00188       END IF
00189 *
00190       IF( LSAME( JOBVR, 'N' ) ) THEN
00191          IJOBVR = 1
00192          ILVR = .FALSE.
00193       ELSE IF( LSAME( JOBVR, 'V' ) ) THEN
00194          IJOBVR = 2
00195          ILVR = .TRUE.
00196       ELSE
00197          IJOBVR = -1
00198          ILVR = .FALSE.
00199       END IF
00200       ILV = ILVL .OR. ILVR
00201 *
00202 *     Test the input arguments
00203 *
00204       INFO = 0
00205       LQUERY = ( LWORK.EQ.-1 )
00206       IF( IJOBVL.LE.0 ) THEN
00207          INFO = -1
00208       ELSE IF( IJOBVR.LE.0 ) THEN
00209          INFO = -2
00210       ELSE IF( N.LT.0 ) THEN
00211          INFO = -3
00212       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00213          INFO = -5
00214       ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
00215          INFO = -7
00216       ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN
00217          INFO = -11
00218       ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN
00219          INFO = -13
00220       END IF
00221 *
00222 *     Compute workspace
00223 *      (Note: Comments in the code beginning "Workspace:" describe the
00224 *       minimal amount of workspace needed at that point in the code,
00225 *       as well as the preferred amount for good performance.
00226 *       NB refers to the optimal block size for the immediately
00227 *       following subroutine, as returned by ILAENV. The workspace is
00228 *       computed assuming ILO = 1 and IHI = N, the worst case.)
00229 *
00230       IF( INFO.EQ.0 ) THEN
00231          LWKMIN = MAX( 1, 2*N )
00232          LWKOPT = MAX( 1, N + N*ILAENV( 1, 'ZGEQRF', ' ', N, 1, N, 0 ) )
00233          LWKOPT = MAX( LWKOPT, N +
00234      $                 N*ILAENV( 1, 'ZUNMQR', ' ', N, 1, N, 0 ) )
00235          IF( ILVL ) THEN
00236             LWKOPT = MAX( LWKOPT, N +
00237      $                    N*ILAENV( 1, 'ZUNGQR', ' ', N, 1, N, -1 ) )
00238          END IF
00239          WORK( 1 ) = LWKOPT
00240 *
00241          IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY )
00242      $      INFO = -15
00243       END IF
00244 *
00245       IF( INFO.NE.0 ) THEN
00246          CALL XERBLA( 'ZGGEV ', -INFO )
00247          RETURN
00248       ELSE IF( LQUERY ) THEN
00249          RETURN
00250       END IF
00251 *
00252 *     Quick return if possible
00253 *
00254       IF( N.EQ.0 )
00255      $   RETURN
00256 *
00257 *     Get machine constants
00258 *
00259       EPS = DLAMCH( 'E' )*DLAMCH( 'B' )
00260       SMLNUM = DLAMCH( 'S' )
00261       BIGNUM = ONE / SMLNUM
00262       CALL DLABAD( SMLNUM, BIGNUM )
00263       SMLNUM = SQRT( SMLNUM ) / EPS
00264       BIGNUM = ONE / SMLNUM
00265 *
00266 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00267 *
00268       ANRM = ZLANGE( 'M', N, N, A, LDA, RWORK )
00269       ILASCL = .FALSE.
00270       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00271          ANRMTO = SMLNUM
00272          ILASCL = .TRUE.
00273       ELSE IF( ANRM.GT.BIGNUM ) THEN
00274          ANRMTO = BIGNUM
00275          ILASCL = .TRUE.
00276       END IF
00277       IF( ILASCL )
00278      $   CALL ZLASCL( 'G', 0, 0, ANRM, ANRMTO, N, N, A, LDA, IERR )
00279 *
00280 *     Scale B if max element outside range [SMLNUM,BIGNUM]
00281 *
00282       BNRM = ZLANGE( 'M', N, N, B, LDB, RWORK )
00283       ILBSCL = .FALSE.
00284       IF( BNRM.GT.ZERO .AND. BNRM.LT.SMLNUM ) THEN
00285          BNRMTO = SMLNUM
00286          ILBSCL = .TRUE.
00287       ELSE IF( BNRM.GT.BIGNUM ) THEN
00288          BNRMTO = BIGNUM
00289          ILBSCL = .TRUE.
00290       END IF
00291       IF( ILBSCL )
00292      $   CALL ZLASCL( 'G', 0, 0, BNRM, BNRMTO, N, N, B, LDB, IERR )
00293 *
00294 *     Permute the matrices A, B to isolate eigenvalues if possible
00295 *     (Real Workspace: need 6*N)
00296 *
00297       ILEFT = 1
00298       IRIGHT = N + 1
00299       IRWRK = IRIGHT + N
00300       CALL ZGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, RWORK( ILEFT ),
00301      $             RWORK( IRIGHT ), RWORK( IRWRK ), IERR )
00302 *
00303 *     Reduce B to triangular form (QR decomposition of B)
00304 *     (Complex Workspace: need N, prefer N*NB)
00305 *
00306       IROWS = IHI + 1 - ILO
00307       IF( ILV ) THEN
00308          ICOLS = N + 1 - ILO
00309       ELSE
00310          ICOLS = IROWS
00311       END IF
00312       ITAU = 1
00313       IWRK = ITAU + IROWS
00314       CALL ZGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ),
00315      $             WORK( IWRK ), LWORK+1-IWRK, IERR )
00316 *
00317 *     Apply the orthogonal transformation to matrix A
00318 *     (Complex Workspace: need N, prefer N*NB)
00319 *
00320       CALL ZUNMQR( 'L', 'C', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB,
00321      $             WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWRK ),
00322      $             LWORK+1-IWRK, IERR )
00323 *
00324 *     Initialize VL
00325 *     (Complex Workspace: need N, prefer N*NB)
00326 *
00327       IF( ILVL ) THEN
00328          CALL ZLASET( 'Full', N, N, CZERO, CONE, VL, LDVL )
00329          IF( IROWS.GT.1 ) THEN
00330             CALL ZLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB,
00331      $                   VL( ILO+1, ILO ), LDVL )
00332          END IF
00333          CALL ZUNGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL,
00334      $                WORK( ITAU ), WORK( IWRK ), LWORK+1-IWRK, IERR )
00335       END IF
00336 *
00337 *     Initialize VR
00338 *
00339       IF( ILVR )
00340      $   CALL ZLASET( 'Full', N, N, CZERO, CONE, VR, LDVR )
00341 *
00342 *     Reduce to generalized Hessenberg form
00343 *
00344       IF( ILV ) THEN
00345 *
00346 *        Eigenvectors requested -- work on whole matrix.
00347 *
00348          CALL ZGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL,
00349      $                LDVL, VR, LDVR, IERR )
00350       ELSE
00351          CALL ZGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA,
00352      $                B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IERR )
00353       END IF
00354 *
00355 *     Perform QZ algorithm (Compute eigenvalues, and optionally, the
00356 *     Schur form and Schur vectors)
00357 *     (Complex Workspace: need N)
00358 *     (Real Workspace: need N)
00359 *
00360       IWRK = ITAU
00361       IF( ILV ) THEN
00362          CHTEMP = 'S'
00363       ELSE
00364          CHTEMP = 'E'
00365       END IF
00366       CALL ZHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB,
00367      $             ALPHA, BETA, VL, LDVL, VR, LDVR, WORK( IWRK ),
00368      $             LWORK+1-IWRK, RWORK( IRWRK ), IERR )
00369       IF( IERR.NE.0 ) THEN
00370          IF( IERR.GT.0 .AND. IERR.LE.N ) THEN
00371             INFO = IERR
00372          ELSE IF( IERR.GT.N .AND. IERR.LE.2*N ) THEN
00373             INFO = IERR - N
00374          ELSE
00375             INFO = N + 1
00376          END IF
00377          GO TO 70
00378       END IF
00379 *
00380 *     Compute Eigenvectors
00381 *     (Real Workspace: need 2*N)
00382 *     (Complex Workspace: need 2*N)
00383 *
00384       IF( ILV ) THEN
00385          IF( ILVL ) THEN
00386             IF( ILVR ) THEN
00387                CHTEMP = 'B'
00388             ELSE
00389                CHTEMP = 'L'
00390             END IF
00391          ELSE
00392             CHTEMP = 'R'
00393          END IF
00394 *
00395          CALL ZTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL,
00396      $                VR, LDVR, N, IN, WORK( IWRK ), RWORK( IRWRK ),
00397      $                IERR )
00398          IF( IERR.NE.0 ) THEN
00399             INFO = N + 2
00400             GO TO 70
00401          END IF
00402 *
00403 *        Undo balancing on VL and VR and normalization
00404 *        (Workspace: none needed)
00405 *
00406          IF( ILVL ) THEN
00407             CALL ZGGBAK( 'P', 'L', N, ILO, IHI, RWORK( ILEFT ),
00408      $                   RWORK( IRIGHT ), N, VL, LDVL, IERR )
00409             DO 30 JC = 1, N
00410                TEMP = ZERO
00411                DO 10 JR = 1, N
00412                   TEMP = MAX( TEMP, ABS1( VL( JR, JC ) ) )
00413    10          CONTINUE
00414                IF( TEMP.LT.SMLNUM )
00415      $            GO TO 30
00416                TEMP = ONE / TEMP
00417                DO 20 JR = 1, N
00418                   VL( JR, JC ) = VL( JR, JC )*TEMP
00419    20          CONTINUE
00420    30       CONTINUE
00421          END IF
00422          IF( ILVR ) THEN
00423             CALL ZGGBAK( 'P', 'R', N, ILO, IHI, RWORK( ILEFT ),
00424      $                   RWORK( IRIGHT ), N, VR, LDVR, IERR )
00425             DO 60 JC = 1, N
00426                TEMP = ZERO
00427                DO 40 JR = 1, N
00428                   TEMP = MAX( TEMP, ABS1( VR( JR, JC ) ) )
00429    40          CONTINUE
00430                IF( TEMP.LT.SMLNUM )
00431      $            GO TO 60
00432                TEMP = ONE / TEMP
00433                DO 50 JR = 1, N
00434                   VR( JR, JC ) = VR( JR, JC )*TEMP
00435    50          CONTINUE
00436    60       CONTINUE
00437          END IF
00438       END IF
00439 *
00440 *     Undo scaling if necessary
00441 *
00442       IF( ILASCL )
00443      $   CALL ZLASCL( 'G', 0, 0, ANRMTO, ANRM, N, 1, ALPHA, N, IERR )
00444 *
00445       IF( ILBSCL )
00446      $   CALL ZLASCL( 'G', 0, 0, BNRMTO, BNRM, N, 1, BETA, N, IERR )
00447 *
00448    70 CONTINUE
00449       WORK( 1 ) = LWKOPT
00450 *
00451       RETURN
00452 *
00453 *     End of ZGGEV
00454 *
00455       END