![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
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