LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
cgeev.f
Go to the documentation of this file.
00001 *> \brief <b> CGEEV 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 CGEEV + dependencies 
00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/cgeev.f"> 
00011 *> [TGZ]</a> 
00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/cgeev.f"> 
00013 *> [ZIP]</a> 
00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/cgeev.f"> 
00015 *> [TXT]</a>
00016 *> \endhtmlonly 
00017 *
00018 *  Definition:
00019 *  ===========
00020 *
00021 *       SUBROUTINE CGEEV( 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 *       REAL               RWORK( * )
00030 *       COMPLEX            A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
00031 *      $                   W( * ), WORK( * )
00032 *       ..
00033 *  
00034 *
00035 *> \par Purpose:
00036 *  =============
00037 *>
00038 *> \verbatim
00039 *>
00040 *> CGEEV 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 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 array, dimension (N)
00093 *>          W contains the computed eigenvalues.
00094 *> \endverbatim
00095 *>
00096 *> \param[out] VL
00097 *> \verbatim
00098 *>          VL is COMPLEX 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 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 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 REAL 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 complexGEeigen
00175 *
00176 *  =====================================================================
00177       SUBROUTINE CGEEV( 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       REAL               RWORK( * )
00191       COMPLEX            A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
00192      $                   W( * ), WORK( * )
00193 *     ..
00194 *
00195 *  =====================================================================
00196 *
00197 *     .. Parameters ..
00198       REAL               ZERO, ONE
00199       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
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       REAL               ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
00207       COMPLEX            TMP
00208 *     ..
00209 *     .. Local Arrays ..
00210       LOGICAL            SELECT( 1 )
00211       REAL               DUM( 1 )
00212 *     ..
00213 *     .. External Subroutines ..
00214       EXTERNAL           CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY, CLASCL,
00215      $                   CSCAL, CSSCAL, CTREVC, CUNGHR, SLABAD, XERBLA
00216 *     ..
00217 *     .. External Functions ..
00218       LOGICAL            LSAME
00219       INTEGER            ILAENV, ISAMAX
00220       REAL               CLANGE, SCNRM2, SLAMCH
00221       EXTERNAL           LSAME, ILAENV, ISAMAX, CLANGE, SCNRM2, SLAMCH
00222 *     ..
00223 *     .. Intrinsic Functions ..
00224       INTRINSIC          AIMAG, CMPLX, CONJG, MAX, REAL, 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 *
00249 *     Compute workspace
00250 *      (Note: Comments in the code beginning "Workspace:" describe the
00251 *       minimal amount of workspace needed at that point in the code,
00252 *       as well as the preferred amount for good performance.
00253 *       CWorkspace refers to complex workspace, and RWorkspace to real
00254 *       workspace. NB refers to the optimal block size for the
00255 *       immediately following subroutine, as returned by ILAENV.
00256 *       HSWORK refers to the workspace preferred by CHSEQR, as
00257 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
00258 *       the worst case.)
00259 *
00260       IF( INFO.EQ.0 ) THEN
00261          IF( N.EQ.0 ) THEN
00262             MINWRK = 1
00263             MAXWRK = 1
00264          ELSE
00265             MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 )
00266             MINWRK = 2*N
00267             IF( WANTVL ) THEN
00268                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR',
00269      $                       ' ', N, 1, N, -1 ) )
00270                CALL CHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
00271      $                WORK, -1, INFO )
00272             ELSE IF( WANTVR ) THEN
00273                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR',
00274      $                       ' ', N, 1, N, -1 ) )
00275                CALL CHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
00276      $                WORK, -1, INFO )
00277             ELSE
00278                CALL CHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
00279      $                WORK, -1, INFO )
00280             END IF
00281             HSWORK = WORK( 1 )
00282             MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )
00283          END IF
00284          WORK( 1 ) = MAXWRK
00285 *
00286          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00287             INFO = -12
00288          END IF
00289       END IF
00290 *
00291       IF( INFO.NE.0 ) THEN
00292          CALL XERBLA( 'CGEEV ', -INFO )
00293          RETURN
00294       ELSE IF( LQUERY ) THEN
00295          RETURN
00296       END IF
00297 *
00298 *     Quick return if possible
00299 *
00300       IF( N.EQ.0 )
00301      $   RETURN
00302 *
00303 *     Get machine constants
00304 *
00305       EPS = SLAMCH( 'P' )
00306       SMLNUM = SLAMCH( 'S' )
00307       BIGNUM = ONE / SMLNUM
00308       CALL SLABAD( SMLNUM, BIGNUM )
00309       SMLNUM = SQRT( SMLNUM ) / EPS
00310       BIGNUM = ONE / SMLNUM
00311 *
00312 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00313 *
00314       ANRM = CLANGE( 'M', N, N, A, LDA, DUM )
00315       SCALEA = .FALSE.
00316       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00317          SCALEA = .TRUE.
00318          CSCALE = SMLNUM
00319       ELSE IF( ANRM.GT.BIGNUM ) THEN
00320          SCALEA = .TRUE.
00321          CSCALE = BIGNUM
00322       END IF
00323       IF( SCALEA )
00324      $   CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
00325 *
00326 *     Balance the matrix
00327 *     (CWorkspace: none)
00328 *     (RWorkspace: need N)
00329 *
00330       IBAL = 1
00331       CALL CGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
00332 *
00333 *     Reduce to upper Hessenberg form
00334 *     (CWorkspace: need 2*N, prefer N+N*NB)
00335 *     (RWorkspace: none)
00336 *
00337       ITAU = 1
00338       IWRK = ITAU + N
00339       CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
00340      $             LWORK-IWRK+1, IERR )
00341 *
00342       IF( WANTVL ) THEN
00343 *
00344 *        Want left eigenvectors
00345 *        Copy Householder vectors to VL
00346 *
00347          SIDE = 'L'
00348          CALL CLACPY( 'L', N, N, A, LDA, VL, LDVL )
00349 *
00350 *        Generate unitary matrix in VL
00351 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00352 *        (RWorkspace: none)
00353 *
00354          CALL CUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
00355      $                LWORK-IWRK+1, IERR )
00356 *
00357 *        Perform QR iteration, accumulating Schur vectors in VL
00358 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00359 *        (RWorkspace: none)
00360 *
00361          IWRK = ITAU
00362          CALL CHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
00363      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00364 *
00365          IF( WANTVR ) THEN
00366 *
00367 *           Want left and right eigenvectors
00368 *           Copy Schur vectors to VR
00369 *
00370             SIDE = 'B'
00371             CALL CLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
00372          END IF
00373 *
00374       ELSE IF( WANTVR ) THEN
00375 *
00376 *        Want right eigenvectors
00377 *        Copy Householder vectors to VR
00378 *
00379          SIDE = 'R'
00380          CALL CLACPY( 'L', N, N, A, LDA, VR, LDVR )
00381 *
00382 *        Generate unitary matrix in VR
00383 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00384 *        (RWorkspace: none)
00385 *
00386          CALL CUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
00387      $                LWORK-IWRK+1, IERR )
00388 *
00389 *        Perform QR iteration, accumulating Schur vectors in VR
00390 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00391 *        (RWorkspace: none)
00392 *
00393          IWRK = ITAU
00394          CALL CHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
00395      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00396 *
00397       ELSE
00398 *
00399 *        Compute eigenvalues only
00400 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00401 *        (RWorkspace: none)
00402 *
00403          IWRK = ITAU
00404          CALL CHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
00405      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00406       END IF
00407 *
00408 *     If INFO > 0 from CHSEQR, then quit
00409 *
00410       IF( INFO.GT.0 )
00411      $   GO TO 50
00412 *
00413       IF( WANTVL .OR. WANTVR ) THEN
00414 *
00415 *        Compute left and/or right eigenvectors
00416 *        (CWorkspace: need 2*N)
00417 *        (RWorkspace: need 2*N)
00418 *
00419          IRWORK = IBAL + N
00420          CALL CTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
00421      $                N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
00422       END IF
00423 *
00424       IF( WANTVL ) THEN
00425 *
00426 *        Undo balancing of left eigenvectors
00427 *        (CWorkspace: none)
00428 *        (RWorkspace: need N)
00429 *
00430          CALL CGEBAK( 'B', 'L', N, ILO, IHI, RWORK( IBAL ), N, VL, LDVL,
00431      $                IERR )
00432 *
00433 *        Normalize left eigenvectors and make largest component real
00434 *
00435          DO 20 I = 1, N
00436             SCL = ONE / SCNRM2( N, VL( 1, I ), 1 )
00437             CALL CSSCAL( N, SCL, VL( 1, I ), 1 )
00438             DO 10 K = 1, N
00439                RWORK( IRWORK+K-1 ) = REAL( VL( K, I ) )**2 +
00440      $                               AIMAG( VL( K, I ) )**2
00441    10       CONTINUE
00442             K = ISAMAX( N, RWORK( IRWORK ), 1 )
00443             TMP = CONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
00444             CALL CSCAL( N, TMP, VL( 1, I ), 1 )
00445             VL( K, I ) = CMPLX( REAL( VL( K, I ) ), ZERO )
00446    20    CONTINUE
00447       END IF
00448 *
00449       IF( WANTVR ) THEN
00450 *
00451 *        Undo balancing of right eigenvectors
00452 *        (CWorkspace: none)
00453 *        (RWorkspace: need N)
00454 *
00455          CALL CGEBAK( 'B', 'R', N, ILO, IHI, RWORK( IBAL ), N, VR, LDVR,
00456      $                IERR )
00457 *
00458 *        Normalize right eigenvectors and make largest component real
00459 *
00460          DO 40 I = 1, N
00461             SCL = ONE / SCNRM2( N, VR( 1, I ), 1 )
00462             CALL CSSCAL( N, SCL, VR( 1, I ), 1 )
00463             DO 30 K = 1, N
00464                RWORK( IRWORK+K-1 ) = REAL( VR( K, I ) )**2 +
00465      $                               AIMAG( VR( K, I ) )**2
00466    30       CONTINUE
00467             K = ISAMAX( N, RWORK( IRWORK ), 1 )
00468             TMP = CONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
00469             CALL CSCAL( N, TMP, VR( 1, I ), 1 )
00470             VR( K, I ) = CMPLX( REAL( VR( K, I ) ), ZERO )
00471    40    CONTINUE
00472       END IF
00473 *
00474 *     Undo scaling if necessary
00475 *
00476    50 CONTINUE
00477       IF( SCALEA ) THEN
00478          CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
00479      $                MAX( N-INFO, 1 ), IERR )
00480          IF( INFO.GT.0 ) THEN
00481             CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
00482          END IF
00483       END IF
00484 *
00485       WORK( 1 ) = MAXWRK
00486       RETURN
00487 *
00488 *     End of CGEEV
00489 *
00490       END
 All Files Functions