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