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