LAPACK: SRC/cgeev.f Source File

Go to the documentation of this file.
00001       SUBROUTINE CGEEV( JOBVL, JOBVR, N, A, LDA, W, VL, LDVL, VR, LDVR,
00002      $                  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, LDVL, LDVR, LWORK, N
00012 *     ..
00013 *     .. Array Arguments ..
00014       REAL               RWORK( * )
00015       COMPLEX            A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
00016      $                   W( * ), WORK( * )
00017 *     ..
00018 *
00019 *  Purpose
00020 *  =======
00021 *
00022 *  CGEEV computes for an N-by-N complex nonsymmetric matrix A, the
00023 *  eigenvalues and, optionally, the left and/or right eigenvectors.
00024 *
00025 *  The right eigenvector v(j) of A satisfies
00026 *                   A * v(j) = lambda(j) * v(j)
00027 *  where lambda(j) is its eigenvalue.
00028 *  The left eigenvector u(j) of A satisfies
00029 *                u(j)**H * A = lambda(j) * u(j)**H
00030 *  where u(j)**H denotes the conjugate transpose of u(j).
00031 *
00032 *  The computed eigenvectors are normalized to have Euclidean norm
00033 *  equal to 1 and largest component real.
00034 *
00035 *  Arguments
00036 *  =========
00037 *
00038 *  JOBVL   (input) CHARACTER*1
00039 *          = 'N': left eigenvectors of A are not computed;
00040 *          = 'V': left eigenvectors of are computed.
00041 *
00042 *  JOBVR   (input) CHARACTER*1
00043 *          = 'N': right eigenvectors of A are not computed;
00044 *          = 'V': right eigenvectors of A are computed.
00045 *
00046 *  N       (input) INTEGER
00047 *          The order of the matrix A. N >= 0.
00048 *
00049 *  A       (input/output) COMPLEX array, dimension (LDA,N)
00050 *          On entry, the N-by-N matrix A.
00051 *          On exit, A has been overwritten.
00052 *
00053 *  LDA     (input) INTEGER
00054 *          The leading dimension of the array A.  LDA >= max(1,N).
00055 *
00056 *  W       (output) COMPLEX array, dimension (N)
00057 *          W contains the computed eigenvalues.
00058 *
00059 *  VL      (output) COMPLEX array, dimension (LDVL,N)
00060 *          If JOBVL = 'V', the left eigenvectors u(j) are stored one
00061 *          after another in the columns of VL, in the same order
00062 *          as their eigenvalues.
00063 *          If JOBVL = 'N', VL is not referenced.
00064 *          u(j) = VL(:,j), the j-th column of VL.
00065 *
00066 *  LDVL    (input) INTEGER
00067 *          The leading dimension of the array VL.  LDVL >= 1; if
00068 *          JOBVL = 'V', LDVL >= N.
00069 *
00070 *  VR      (output) COMPLEX array, dimension (LDVR,N)
00071 *          If JOBVR = 'V', the right eigenvectors v(j) are stored one
00072 *          after another in the columns of VR, in the same order
00073 *          as their eigenvalues.
00074 *          If JOBVR = 'N', VR is not referenced.
00075 *          v(j) = VR(:,j), the j-th column of VR.
00076 *
00077 *  LDVR    (input) INTEGER
00078 *          The leading dimension of the array VR.  LDVR >= 1; if
00079 *          JOBVR = 'V', LDVR >= N.
00080 *
00081 *  WORK    (workspace/output) COMPLEX array, dimension (MAX(1,LWORK))
00082 *          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
00083 *
00084 *  LWORK   (input) INTEGER
00085 *          The dimension of the array WORK.  LWORK >= max(1,2*N).
00086 *          For good performance, LWORK must generally be larger.
00087 *
00088 *          If LWORK = -1, then a workspace query is assumed; the routine
00089 *          only calculates the optimal size of the WORK array, returns
00090 *          this value as the first entry of the WORK array, and no error
00091 *          message related to LWORK is issued by XERBLA.
00092 *
00093 *  RWORK   (workspace) REAL array, dimension (2*N)
00094 *
00095 *  INFO    (output) INTEGER
00096 *          = 0:  successful exit
00097 *          < 0:  if INFO = -i, the i-th argument had an illegal value.
00098 *          > 0:  if INFO = i, the QR algorithm failed to compute all the
00099 *                eigenvalues, and no eigenvectors have been computed;
00100 *                elements and i+1:N of W contain eigenvalues which have
00101 *                converged.
00102 *
00103 *  =====================================================================
00104 *
00105 *     .. Parameters ..
00106       REAL               ZERO, ONE
00107       PARAMETER          ( ZERO = 0.0E0, ONE = 1.0E0 )
00108 *     ..
00109 *     .. Local Scalars ..
00110       LOGICAL            LQUERY, SCALEA, WANTVL, WANTVR
00111       CHARACTER          SIDE
00112       INTEGER            HSWORK, I, IBAL, IERR, IHI, ILO, IRWORK, ITAU,
00113      $                   IWRK, K, MAXWRK, MINWRK, NOUT
00114       REAL               ANRM, BIGNUM, CSCALE, EPS, SCL, SMLNUM
00115       COMPLEX            TMP
00116 *     ..
00117 *     .. Local Arrays ..
00118       LOGICAL            SELECT( 1 )
00119       REAL               DUM( 1 )
00120 *     ..
00121 *     .. External Subroutines ..
00122       EXTERNAL           CGEBAK, CGEBAL, CGEHRD, CHSEQR, CLACPY, CLASCL,
00123      $                   CSCAL, CSSCAL, CTREVC, CUNGHR, SLABAD, XERBLA
00124 *     ..
00125 *     .. External Functions ..
00126       LOGICAL            LSAME
00127       INTEGER            ILAENV, ISAMAX
00128       REAL               CLANGE, SCNRM2, SLAMCH
00129       EXTERNAL           LSAME, ILAENV, ISAMAX, CLANGE, SCNRM2, SLAMCH
00130 *     ..
00131 *     .. Intrinsic Functions ..
00132       INTRINSIC          AIMAG, CMPLX, CONJG, MAX, REAL, SQRT
00133 *     ..
00134 *     .. Executable Statements ..
00135 *
00136 *     Test the input arguments
00137 *
00138       INFO = 0
00139       LQUERY = ( LWORK.EQ.-1 )
00140       WANTVL = LSAME( JOBVL, 'V' )
00141       WANTVR = LSAME( JOBVR, 'V' )
00142       IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN
00143          INFO = -1
00144       ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN
00145          INFO = -2
00146       ELSE IF( N.LT.0 ) THEN
00147          INFO = -3
00148       ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
00149          INFO = -5
00150       ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN
00151          INFO = -8
00152       ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN
00153          INFO = -10
00154       END IF
00155 
00156 *
00157 *     Compute workspace
00158 *      (Note: Comments in the code beginning "Workspace:" describe the
00159 *       minimal amount of workspace needed at that point in the code,
00160 *       as well as the preferred amount for good performance.
00161 *       CWorkspace refers to complex workspace, and RWorkspace to real
00162 *       workspace. NB refers to the optimal block size for the
00163 *       immediately following subroutine, as returned by ILAENV.
00164 *       HSWORK refers to the workspace preferred by CHSEQR, as
00165 *       calculated below. HSWORK is computed assuming ILO=1 and IHI=N,
00166 *       the worst case.)
00167 *
00168       IF( INFO.EQ.0 ) THEN
00169          IF( N.EQ.0 ) THEN
00170             MINWRK = 1
00171             MAXWRK = 1
00172          ELSE
00173             MAXWRK = N + N*ILAENV( 1, 'CGEHRD', ' ', N, 1, N, 0 )
00174             MINWRK = 2*N
00175             IF( WANTVL ) THEN
00176                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR',
00177      $                       ' ', N, 1, N, -1 ) )
00178                CALL CHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VL, LDVL,
00179      $                WORK, -1, INFO )
00180             ELSE IF( WANTVR ) THEN
00181                MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'CUNGHR',
00182      $                       ' ', N, 1, N, -1 ) )
00183                CALL CHSEQR( 'S', 'V', N, 1, N, A, LDA, W, VR, LDVR,
00184      $                WORK, -1, INFO )
00185             ELSE
00186                CALL CHSEQR( 'E', 'N', N, 1, N, A, LDA, W, VR, LDVR,
00187      $                WORK, -1, INFO )
00188             END IF
00189             HSWORK = WORK( 1 )
00190             MAXWRK = MAX( MAXWRK, HSWORK, MINWRK )
00191          END IF
00192          WORK( 1 ) = MAXWRK
00193 *
00194          IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN
00195             INFO = -12
00196          END IF
00197       END IF
00198 *
00199       IF( INFO.NE.0 ) THEN
00200          CALL XERBLA( 'CGEEV ', -INFO )
00201          RETURN
00202       ELSE IF( LQUERY ) THEN
00203          RETURN
00204       END IF
00205 *
00206 *     Quick return if possible
00207 *
00208       IF( N.EQ.0 )
00209      $   RETURN
00210 *
00211 *     Get machine constants
00212 *
00213       EPS = SLAMCH( 'P' )
00214       SMLNUM = SLAMCH( 'S' )
00215       BIGNUM = ONE / SMLNUM
00216       CALL SLABAD( SMLNUM, BIGNUM )
00217       SMLNUM = SQRT( SMLNUM ) / EPS
00218       BIGNUM = ONE / SMLNUM
00219 *
00220 *     Scale A if max element outside range [SMLNUM,BIGNUM]
00221 *
00222       ANRM = CLANGE( 'M', N, N, A, LDA, DUM )
00223       SCALEA = .FALSE.
00224       IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN
00225          SCALEA = .TRUE.
00226          CSCALE = SMLNUM
00227       ELSE IF( ANRM.GT.BIGNUM ) THEN
00228          SCALEA = .TRUE.
00229          CSCALE = BIGNUM
00230       END IF
00231       IF( SCALEA )
00232      $   CALL CLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR )
00233 *
00234 *     Balance the matrix
00235 *     (CWorkspace: none)
00236 *     (RWorkspace: need N)
00237 *
00238       IBAL = 1
00239       CALL CGEBAL( 'B', N, A, LDA, ILO, IHI, RWORK( IBAL ), IERR )
00240 *
00241 *     Reduce to upper Hessenberg form
00242 *     (CWorkspace: need 2*N, prefer N+N*NB)
00243 *     (RWorkspace: none)
00244 *
00245       ITAU = 1
00246       IWRK = ITAU + N
00247       CALL CGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ),
00248      $             LWORK-IWRK+1, IERR )
00249 *
00250       IF( WANTVL ) THEN
00251 *
00252 *        Want left eigenvectors
00253 *        Copy Householder vectors to VL
00254 *
00255          SIDE = 'L'
00256          CALL CLACPY( 'L', N, N, A, LDA, VL, LDVL )
00257 *
00258 *        Generate unitary matrix in VL
00259 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00260 *        (RWorkspace: none)
00261 *
00262          CALL CUNGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ),
00263      $                LWORK-IWRK+1, IERR )
00264 *
00265 *        Perform QR iteration, accumulating Schur vectors in VL
00266 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00267 *        (RWorkspace: none)
00268 *
00269          IWRK = ITAU
00270          CALL CHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VL, LDVL,
00271      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00272 *
00273          IF( WANTVR ) THEN
00274 *
00275 *           Want left and right eigenvectors
00276 *           Copy Schur vectors to VR
00277 *
00278             SIDE = 'B'
00279             CALL CLACPY( 'F', N, N, VL, LDVL, VR, LDVR )
00280          END IF
00281 *
00282       ELSE IF( WANTVR ) THEN
00283 *
00284 *        Want right eigenvectors
00285 *        Copy Householder vectors to VR
00286 *
00287          SIDE = 'R'
00288          CALL CLACPY( 'L', N, N, A, LDA, VR, LDVR )
00289 *
00290 *        Generate unitary matrix in VR
00291 *        (CWorkspace: need 2*N-1, prefer N+(N-1)*NB)
00292 *        (RWorkspace: none)
00293 *
00294          CALL CUNGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ),
00295      $                LWORK-IWRK+1, IERR )
00296 *
00297 *        Perform QR iteration, accumulating Schur vectors in VR
00298 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00299 *        (RWorkspace: none)
00300 *
00301          IWRK = ITAU
00302          CALL CHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, W, VR, LDVR,
00303      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00304 *
00305       ELSE
00306 *
00307 *        Compute eigenvalues only
00308 *        (CWorkspace: need 1, prefer HSWORK (see comments) )
00309 *        (RWorkspace: none)
00310 *
00311          IWRK = ITAU
00312          CALL CHSEQR( 'E', 'N', N, ILO, IHI, A, LDA, W, VR, LDVR,
00313      $                WORK( IWRK ), LWORK-IWRK+1, INFO )
00314       END IF
00315 *
00316 *     If INFO > 0 from CHSEQR, then quit
00317 *
00318       IF( INFO.GT.0 )
00319      $   GO TO 50
00320 *
00321       IF( WANTVL .OR. WANTVR ) THEN
00322 *
00323 *        Compute left and/or right eigenvectors
00324 *        (CWorkspace: need 2*N)
00325 *        (RWorkspace: need 2*N)
00326 *
00327          IRWORK = IBAL + N
00328          CALL CTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR,
00329      $                N, NOUT, WORK( IWRK ), RWORK( IRWORK ), IERR )
00330       END IF
00331 *
00332       IF( WANTVL ) THEN
00333 *
00334 *        Undo balancing of left eigenvectors
00335 *        (CWorkspace: none)
00336 *        (RWorkspace: need N)
00337 *
00338          CALL CGEBAK( 'B', 'L', N, ILO, IHI, RWORK( IBAL ), N, VL, LDVL,
00339      $                IERR )
00340 *
00341 *        Normalize left eigenvectors and make largest component real
00342 *
00343          DO 20 I = 1, N
00344             SCL = ONE / SCNRM2( N, VL( 1, I ), 1 )
00345             CALL CSSCAL( N, SCL, VL( 1, I ), 1 )
00346             DO 10 K = 1, N
00347                RWORK( IRWORK+K-1 ) = REAL( VL( K, I ) )**2 +
00348      $                               AIMAG( VL( K, I ) )**2
00349    10       CONTINUE
00350             K = ISAMAX( N, RWORK( IRWORK ), 1 )
00351             TMP = CONJG( VL( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
00352             CALL CSCAL( N, TMP, VL( 1, I ), 1 )
00353             VL( K, I ) = CMPLX( REAL( VL( K, I ) ), ZERO )
00354    20    CONTINUE
00355       END IF
00356 *
00357       IF( WANTVR ) THEN
00358 *
00359 *        Undo balancing of right eigenvectors
00360 *        (CWorkspace: none)
00361 *        (RWorkspace: need N)
00362 *
00363          CALL CGEBAK( 'B', 'R', N, ILO, IHI, RWORK( IBAL ), N, VR, LDVR,
00364      $                IERR )
00365 *
00366 *        Normalize right eigenvectors and make largest component real
00367 *
00368          DO 40 I = 1, N
00369             SCL = ONE / SCNRM2( N, VR( 1, I ), 1 )
00370             CALL CSSCAL( N, SCL, VR( 1, I ), 1 )
00371             DO 30 K = 1, N
00372                RWORK( IRWORK+K-1 ) = REAL( VR( K, I ) )**2 +
00373      $                               AIMAG( VR( K, I ) )**2
00374    30       CONTINUE
00375             K = ISAMAX( N, RWORK( IRWORK ), 1 )
00376             TMP = CONJG( VR( K, I ) ) / SQRT( RWORK( IRWORK+K-1 ) )
00377             CALL CSCAL( N, TMP, VR( 1, I ), 1 )
00378             VR( K, I ) = CMPLX( REAL( VR( K, I ) ), ZERO )
00379    40    CONTINUE
00380       END IF
00381 *
00382 *     Undo scaling if necessary
00383 *
00384    50 CONTINUE
00385       IF( SCALEA ) THEN
00386          CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, W( INFO+1 ),
00387      $                MAX( N-INFO, 1 ), IERR )
00388          IF( INFO.GT.0 ) THEN
00389             CALL CLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, W, N, IERR )
00390          END IF
00391       END IF
00392 *
00393       WORK( 1 ) = MAXWRK
00394       RETURN
00395 *
00396 *     End of CGEEV
00397 *
00398       END