LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
zgeev.f
Go to the documentation of this file.
00001 *> \brief <b> ZGEEV 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 ZGEEV + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/zgeev.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/zgeev.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/zgeev.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
00022 *                         WORK, LWORK, RWORK, INFO )
00023 * 
00024 *       .. Scalar Arguments ..
00025 *       CHARACTER          JOBVL, JOBVR
00026 *       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
00027 *       ..
00028 *       .. Array Arguments ..
00029 *       DOUBLE PRECISION   RWORK( * )
00030 *       COMPLEX*16         A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
00031 *      $                   W( * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> ZGEEV computes for an N-by-N complex nonsymmetric matrix A, the
00041 *> eigenvalues and, optionally, the left and/or right eigenvectors.
00042 *>
00043 *> The right eigenvector v(j) of A satisfies
00044 *>                  A * v(j) = lambda(j) * v(j)
00045 *> where lambda(j) is its eigenvalue.
00046 *> The left eigenvector u(j) of A satisfies
00047 *>               u(j)**H * A = lambda(j) * u(j)**H
00048 *> where u(j)**H denotes the conjugate transpose of u(j).
00049 *>
00050 *> The computed eigenvectors are normalized to have Euclidean norm
00051 *> equal to 1 and largest component real.
00052 *> \endverbatim
00053 *
00054 *  Arguments:
00055 *  ==========
00056 *
00057 *> \param[in] JOBVL
00058 *> \verbatim
00059 *>          JOBVL is CHARACTER*1
00060 *>          = 'N': left eigenvectors of A are not computed;
00061 *>          = 'V': left eigenvectors of are computed.
00062 *> \endverbatim
00063 *>
00064 *> \param[in] JOBVR
00065 *> \verbatim
00066 *>          JOBVR is CHARACTER*1
00067 *>          = 'N': right eigenvectors of A are not computed;
00068 *>          = 'V': right eigenvectors of A are computed.
00069 *> \endverbatim
00070 *>
00071 *> \param[in] N
00072 *> \verbatim
00073 *>          N is INTEGER
00074 *>          The order of the matrix A. N >= 0.
00075 *> \endverbatim
00076 *>
00077 *> \param[in,out] A
00078 *> \verbatim
00079 *>          A is COMPLEX*16 array, dimension (LDA,N)
00080 *>          On entry, the N-by-N matrix A.
00081 *>          On exit, A has been overwritten.
00082 *> \endverbatim
00083 *>
00084 *> \param[in] LDA
00085 *> \verbatim
00086 *>          LDA is INTEGER
00087 *>          The leading dimension of the array A.  LDA >= max(1,N).
00088 *> \endverbatim
00089 *>
00090 *> \param[out] W
00091 *> \verbatim
00092 *>          W is COMPLEX*16 array, dimension (N)
00093 *>          W contains the computed eigenvalues.
00094 *> \endverbatim
00095 *>
00096 *> \param[out] VL
00097 *> \verbatim
00098 *>          VL is COMPLEX*16 array, dimension (LDVL,N)
00099 *>          If JOBVL = 'V', the left eigenvectors u(j) are stored one
00100 *>          after another in the columns of VL, in the same order
00101 *>          as their eigenvalues.
00102 *>          If JOBVL = 'N', VL is not referenced.
00103 *>          u(j) = VL(:,j), the j-th column of VL.
00104 *> \endverbatim
00105 *>
00106 *> \param[in] LDVL
00107 *> \verbatim
00108 *>          LDVL is INTEGER
00109 *>          The leading dimension of the array VL.  LDVL >= 1; if
00110 *>          JOBVL = 'V', LDVL >= N.
00111 *> \endverbatim
00112 *>
00113 *> \param[out] VR
00114 *> \verbatim
00115 *>          VR is COMPLEX*16 array, dimension (LDVR,N)
00116 *>          If JOBVR = 'V', the right eigenvectors v(j) are stored one
00117 *>          after another in the columns of VR, in the same order
00118 *>          as their eigenvalues.
00119 *>          If JOBVR = 'N', VR is not referenced.
00120 *>          v(j) = VR(:,j), the j-th column of VR.
00121 *> \endverbatim
00122 *>
00123 *> \param[in] LDVR
00124 *> \verbatim
00125 *>          LDVR is INTEGER
00126 *>          The leading dimension of the array VR.  LDVR >= 1; if
00127 *>          JOBVR = 'V', LDVR >= N.
00128 *> \endverbatim
00129 *>
00130 *> \param[out] WORK
00131 *> \verbatim
00132 *>          WORK is COMPLEX*16 array, dimension (MAX(1,LWORK))
00133 *>          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00134 *> \endverbatim
00135 *>
00136 *> \param[in] LWORK
00137 *> \verbatim
00138 *>          LWORK is INTEGER
00139 *>          The dimension of the array WORK.  LWORK >= max(1,2*N).
00140 *>          For good performance, LWORK must generally be larger.
00141 *>
00142 *>          If LWORK = -1, then a workspace query is assumed; the routine
00143 *>          only calculates the optimal size of the WORK array, returns
00144 *>          this value as the first entry of the WORK array, and no error
00145 *>          message related to LWORK is issued by XERBLA.
00146 *> \endverbatim
00147 *>
00148 *> \param[out] RWORK
00149 *> \verbatim
00150 *>          RWORK is DOUBLE PRECISION array, dimension (2*N)
00151 *> \endverbatim
00152 *>
00153 *> \param[out] INFO
00154 *> \verbatim
00155 *>          INFO is INTEGER
00156 *>          = 0:  successful exit
00157 *>          < 0:  if INFO = -i, the i-th argument had an illegal value.
00158 *>          > 0:  if INFO = i, the QR algorithm failed to compute all the
00159 *>                eigenvalues, and no eigenvectors have been computed;
00160 *>                elements and i+1:N of W contain eigenvalues which have
00161 *>                converged.
00162 *> \endverbatim
00163 *
00164 *  Authors:
00165 *  ========
00166 *
00167 *> \author Univ. of Tennessee 
00168 *> \author Univ. of California Berkeley 
00169 *> \author Univ. of Colorado Denver 
00170 *> \author NAG Ltd. 
00171 *
00172 *> \date November 2011
00173 *
00174 *> \ingroup complex16GEeigen
00175 *
00176 *  =====================================================================
00177       SUBROUTINE ZGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
00178      $                  WORK, LWORK, RWORK, INFO )
00179 *
00180 *  -- LAPACK driver routine (version 3.4.0) --
00181 *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
00182 *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
00183 *     November 2011
00184 *
00185 *     .. Scalar Arguments ..
00186       CHARACTER          JOBVL, JOBVR
00187       INTEGER            INFO, LDA, LDVL, LDVR, LWORK, N
00188 *     ..
00189 *     .. Array Arguments ..
00190       DOUBLE PRECISION   RWORK( * )
00191       COMPLEX*16         A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
00192      $                   W( * ), WORK( * )
00193 *     ..
00194 *
00195 *  =====================================================================
00196 *
00197 *     .. Parameters ..
00198       DOUBLE PRECISION   ZERO, ONE
00199       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
00200 *     ..
00201 *     .. Local Scalars ..
00202       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
00203       CHARACTER          SIDE
00204       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
00205      $                   IWRK, K, MAXWRK, MINWRK, NOUT
00206       DOUBLE PRECISION   ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
00207       COMPLEX*16         TMP
00208 *     ..
00209 *     .. Local Arrays ..
00210       LOGICAL            SELECT( 1 )
00211       DOUBLE PRECISION   DUM( 1 )
00212 *     ..
00213 *     .. External Subroutines ..
00214       EXTERNAL           DLABAD, XERBLA, ZDSCAL, ZGEBAK, ZGEBAL, ZGEHRD,
00215      $                   ZHSEQR, ZLACPY, ZLASCL, ZSCAL, ZTREVC, ZUNGHR
00216 *     ..
00217 *     .. External Functions ..
00218       LOGICAL            LSAME
00219       INTEGER            IDAMAX, ILAENV
00220       DOUBLE PRECISION   DLAMCH, DZNRM2, ZLANGE
00221       EXTERNAL           LSAME, IDAMAX, ILAENV, DLAMCH, DZNRM2, ZLANGE
00222 *     ..
00223 *     .. Intrinsic Functions ..
00224       INTRINSIC          DBLE, DCMPLX, DCONJG, DIMAG, MAX, SQRT
00225 *     ..
00226 *     .. Executable Statements ..
00227 *
00228 *     Test the input arguments
00229 *
00230       INFO = 0
00231       LQUERY = ( LWORK.EQ.-1 )
00232       WANTVL = LSAME( JOBVL, 'V' )
00233       WANTVR = LSAME( JOBVR, 'V' )
00234       IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
00235          INFO = -1
00236       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
00237          INFO = -2
00238       ELSE IF( N.LT.0 ) THEN
00239          INFO = -3
00240       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00241          INFO = -5
00242       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
00243          INFO = -8
00244       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
00245          INFO = -10
00246       END IF
00247 *
00248 *     Compute workspace
00249 *      (Note: Comments in the code beginning "Workspace:" describe the
00250 *       minimal amount of workspace needed at that point in the code,
00251 *       as well as the preferred amount for good performance.
00252 *       CWorkspace refers to complex workspace, and RWorkspace to real
00253 *       workspace. NB refers to the optimal block size for the
00254 *       immediately following subroutine, as returned by ILAENV.
00255 *       HSWORK refers to the workspace preferred by ZHSEQR, as
00256 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
00257 *       the worst case.)
00258 *
00259       IF( INFO.EQ.0 ) THEN
00260          IF( N.EQ.0 ) THEN
00261             MINWRK = 1
00262             MAXWRK = 1
00263          ELSE
00264             MAXWRK = N + N*ILAENV( 1, 'ZGEHRD', ' ', N, 1, N, 0 )
00265             MINWRK = 2*N
00266             IF( WANTVL ) THEN
00267                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
00268      $                       ' ', N, 1, N, -1 ) )
00269                CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
00270      $                WORK, -1, INFO )
00271             ELSE IF( WANTVR ) THEN
00272                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'ZUNGHR',
00273      $                       ' ', N, 1, N, -1 ) )
00274                CALL ZHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
00275      $                WORK, -1, INFO )
00276             ELSE
00277                CALL ZHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
00278      $                WORK, -1, INFO )
00279             END IF
00280             HSWORK = WORK( 1 )
00281             MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )
00282          END IF
00283          WORK( 1 ) = MAXWRK
00284 *
00285          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00286             INFO = -12
00287          END IF
00288       END IF
00289 *
00290       IF( INFO.NE.0 ) THEN
00291          CALL XERBLA( 'ZGEEV ', -INFO )
00292          RETURN
00293       ELSE IF( LQUERY ) THEN
00294          RETURN
00295       END IF
00296 *
00297 *     Quick return if possible
00298 *
00299       IF( N.EQ.0 )
00300      $   RETURN
00301 *
00302 *     Get machine constants
00303 *
00304       EPS = DLAMCH( 'P' )
00305       SMLNUM = DLAMCH( 'S' )
00306       BIGNUM = ONE / SMLNUM
00307       CALL DLABAD( SMLNUM, BIGNUM )
00308       SMLNUM = SQRT( SMLNUM ) / EPS
00309       BIGNUM = ONE / SMLNUM
00310 *
00311 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00312 *
00313       ANRM = ZLANGE( 'M', N, N, A, LDA, DUM )
00314       SCALEA = .FALSE.
00315       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00316          SCALEA = .TRUE.
00317          CSCALE = SMLNUM
00318       ELSE IF( ANRM.GT.BIGNUM ) THEN
00319          SCALEA = .TRUE.
00320          CSCALE = BIGNUM
00321       END IF
00322       IF( SCALEA )
00323      $   CALL ZLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
00324 *
00325 *     Balance the matrix
00326 *     (CWorkspace: none)
00327 *     (RWorkspace: need N)
00328 *
00329       IBAL = 1
00330       CALL ZGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
00331 *
00332 *     Reduce to upper Hessenberg form
00333 *     (CWorkspace: need 2*N, prefer N+N*NB)
00334 *     (RWorkspace: none)
00335 *
00336       ITAU = 1
00337       IWRK = ITAU + N
00338       CALL ZGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
00339      $             LWORK-IWRK+1, IERR )
00340 *
00341       IF( WANTVL ) THEN
00342 *
00343 *        Want left eigenvectors
00344 *        Copy Householder vectors to VL
00345 *
00346          SIDE = 'L'
00347          CALL ZLACPY( 'L', N, N, A, LDA, VL, LDVL )
00348 *
00349 *        Generate unitary matrix in VL
00350 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00351 *        (RWorkspace: none)
00352 *
00353          CALL ZUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
00354      $                LWORK-IWRK+1, IERR )
00355 *
00356 *        Perform QR iteration, accumulating Schur vectors in VL
00357 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00358 *        (RWorkspace: none)
00359 *
00360          IWRK = ITAU
00361          CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
00362      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00363 *
00364          IF( WANTVR ) THEN
00365 *
00366 *           Want left and right eigenvectors
00367 *           Copy Schur vectors to VR
00368 *
00369             SIDE = 'B'
00370             CALL ZLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
00371          END IF
00372 *
00373       ELSE IF( WANTVR ) THEN
00374 *
00375 *        Want right eigenvectors
00376 *        Copy Householder vectors to VR
00377 *
00378          SIDE = 'R'
00379          CALL ZLACPY( 'L', N, N, A, LDA, VR, LDVR )
00380 *
00381 *        Generate unitary matrix in VR
00382 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00383 *        (RWorkspace: none)
00384 *
00385          CALL ZUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
00386      $                LWORK-IWRK+1, IERR )
00387 *
00388 *        Perform QR iteration, accumulating Schur vectors in VR
00389 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00390 *        (RWorkspace: none)
00391 *
00392          IWRK = ITAU
00393          CALL ZHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
00394      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00395 *
00396       ELSE
00397 *
00398 *        Compute eigenvalues only
00399 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00400 *        (RWorkspace: none)
00401 *
00402          IWRK = ITAU
00403          CALL ZHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
00404      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00405       END IF
00406 *
00407 *     If INFO > 0 from ZHSEQR, then quit
00408 *
00409       IF( INFO.GT.0 )
00410      $   GO TO 50
00411 *
00412       IF( WANTVL .OR. WANTVR ) THEN
00413 *
00414 *        Compute left and/or right eigenvectors
00415 *        (CWorkspace: need 2*N)
00416 *        (RWorkspace: need 2*N)
00417 *
00418          IRWORK = IBAL + N
00419          CALL ZTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
00420      $                N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
00421       END IF
00422 *
00423       IF( WANTVL ) THEN
00424 *
00425 *        Undo balancing of left eigenvectors
00426 *        (CWorkspace: none)
00427 *        (RWorkspace: need N)
00428 *
00429          CALL ZGEBAK( 'B', 'L', N, ILO, IHI, RWORK( IBAL ), N, VL, LDVL,
00430      $                IERR )
00431 *
00432 *        Normalize left eigenvectors and make largest component real
00433 *
00434          DO 20 I = 1, N
00435             SCL = ONE / DZNRM2( N, VL( 1, I ), 1 )
00436             CALL ZDSCAL( N, SCL, VL( 1, I ), 1 )
00437             DO 10 K = 1, N
00438                RWORK( IRWORK+K-1 ) = DBLE( VL( K, I ) )**2 +
00439      $                               DIMAG( VL( K, I ) )**2
00440    10       CONTINUE
00441             K = IDAMAX( N, RWORK( IRWORK ), 1 )
00442             TMP = DCONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
00443             CALL ZSCAL( N, TMP, VL( 1, I ), 1 )
00444             VL( K, I ) = DCMPLX( DBLE( VL( K, I ) ), ZERO )
00445    20    CONTINUE
00446       END IF
00447 *
00448       IF( WANTVR ) THEN
00449 *
00450 *        Undo balancing of right eigenvectors
00451 *        (CWorkspace: none)
00452 *        (RWorkspace: need N)
00453 *
00454          CALL ZGEBAK( 'B', 'R', N, ILO, IHI, RWORK( IBAL ), N, VR, LDVR,
00455      $                IERR )
00456 *
00457 *        Normalize right eigenvectors and make largest component real
00458 *
00459          DO 40 I = 1, N
00460             SCL = ONE / DZNRM2( N, VR( 1, I ), 1 )
00461             CALL ZDSCAL( N, SCL, VR( 1, I ), 1 )
00462             DO 30 K = 1, N
00463                RWORK( IRWORK+K-1 ) = DBLE( VR( K, I ) )**2 +
00464      $                               DIMAG( VR( K, I ) )**2
00465    30       CONTINUE
00466             K = IDAMAX( N, RWORK( IRWORK ), 1 )
00467             TMP = DCONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
00468             CALL ZSCAL( N, TMP, VR( 1, I ), 1 )
00469             VR( K, I ) = DCMPLX( DBLE( VR( K, I ) ), ZERO )
00470    40    CONTINUE
00471       END IF
00472 *
00473 *     Undo scaling if necessary
00474 *
00475    50 CONTINUE
00476       IF( SCALEA ) THEN
00477          CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
00478      $                MAX( N-INFO, 1 ), IERR )
00479          IF( INFO.GT.0 ) THEN
00480             CALL ZLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
00481          END IF
00482       END IF
00483 *
00484       WORK( 1 ) = MAXWRK
00485       RETURN
00486 *
00487 *     End of ZGEEV
00488 *
00489       END
 All Files Functions