![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> SGEEVX 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 SGEEVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/sgeevx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/sgeevx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/sgeevx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE SGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, 00022 * VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, 00023 * RCONDE, RCONDV, WORK, LWORK, IWORK, INFO ) 00024 * 00025 * .. Scalar Arguments .. 00026 * CHARACTER BALANC, JOBVL, JOBVR, SENSE 00027 * INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N 00028 * REAL ABNRM 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IWORK( * ) 00032 * REAL A( LDA, * ), RCONDE( * ), RCONDV( * ), 00033 * $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ), 00034 * $ WI( * ), WORK( * ), WR( * ) 00035 * .. 00036 * 00037 * 00038 *> \par Purpose: 00039 * ============= 00040 *> 00041 *> \verbatim 00042 *> 00043 *> SGEEVX computes for an N-by-N real nonsymmetric matrix A, the 00044 *> eigenvalues and, optionally, the left and/or right eigenvectors. 00045 *> 00046 *> Optionally also, it computes a balancing transformation to improve 00047 *> the conditioning of the eigenvalues and eigenvectors (ILO, IHI, 00048 *> SCALE, and ABNRM), reciprocal condition numbers for the eigenvalues 00049 *> (RCONDE), and reciprocal condition numbers for the right 00050 *> eigenvectors (RCONDV). 00051 *> 00052 *> The right eigenvector v(j) of A satisfies 00053 *> A * v(j) = lambda(j) * v(j) 00054 *> where lambda(j) is its eigenvalue. 00055 *> The left eigenvector u(j) of A satisfies 00056 *> u(j)**T * A = lambda(j) * u(j)**T 00057 *> where u(j)**T denotes the transpose of u(j). 00058 *> 00059 *> The computed eigenvectors are normalized to have Euclidean norm 00060 *> equal to 1 and largest component real. 00061 *> 00062 *> Balancing a matrix means permuting the rows and columns to make it 00063 *> more nearly upper triangular, and applying a diagonal similarity 00064 *> transformation D * A * D**(-1), where D is a diagonal matrix, to 00065 *> make its rows and columns closer in norm and the condition numbers 00066 *> of its eigenvalues and eigenvectors smaller. The computed 00067 *> reciprocal condition numbers correspond to the balanced matrix. 00068 *> Permuting rows and columns will not change the condition numbers 00069 *> (in exact arithmetic) but diagonal scaling will. For further 00070 *> explanation of balancing, see section 4.10.2 of the LAPACK 00071 *> Users' Guide. 00072 *> \endverbatim 00073 * 00074 * Arguments: 00075 * ========== 00076 * 00077 *> \param[in] BALANC 00078 *> \verbatim 00079 *> BALANC is CHARACTER*1 00080 *> Indicates how the input matrix should be diagonally scaled 00081 *> and/or permuted to improve the conditioning of its 00082 *> eigenvalues. 00083 *> = 'N': Do not diagonally scale or permute; 00084 *> = 'P': Perform permutations to make the matrix more nearly 00085 *> upper triangular. Do not diagonally scale; 00086 *> = 'S': Diagonally scale the matrix, i.e. replace A by 00087 *> D*A*D**(-1), where D is a diagonal matrix chosen 00088 *> to make the rows and columns of A more equal in 00089 *> norm. Do not permute; 00090 *> = 'B': Both diagonally scale and permute A. 00091 *> 00092 *> Computed reciprocal condition numbers will be for the matrix 00093 *> after balancing and/or permuting. Permuting does not change 00094 *> condition numbers (in exact arithmetic), but balancing does. 00095 *> \endverbatim 00096 *> 00097 *> \param[in] JOBVL 00098 *> \verbatim 00099 *> JOBVL is CHARACTER*1 00100 *> = 'N': left eigenvectors of A are not computed; 00101 *> = 'V': left eigenvectors of A are computed. 00102 *> If SENSE = 'E' or 'B', JOBVL must = 'V'. 00103 *> \endverbatim 00104 *> 00105 *> \param[in] JOBVR 00106 *> \verbatim 00107 *> JOBVR is CHARACTER*1 00108 *> = 'N': right eigenvectors of A are not computed; 00109 *> = 'V': right eigenvectors of A are computed. 00110 *> If SENSE = 'E' or 'B', JOBVR must = 'V'. 00111 *> \endverbatim 00112 *> 00113 *> \param[in] SENSE 00114 *> \verbatim 00115 *> SENSE is CHARACTER*1 00116 *> Determines which reciprocal condition numbers are computed. 00117 *> = 'N': None are computed; 00118 *> = 'E': Computed for eigenvalues only; 00119 *> = 'V': Computed for right eigenvectors only; 00120 *> = 'B': Computed for eigenvalues and right eigenvectors. 00121 *> 00122 *> If SENSE = 'E' or 'B', both left and right eigenvectors 00123 *> must also be computed (JOBVL = 'V' and JOBVR = 'V'). 00124 *> \endverbatim 00125 *> 00126 *> \param[in] N 00127 *> \verbatim 00128 *> N is INTEGER 00129 *> The order of the matrix A. N >= 0. 00130 *> \endverbatim 00131 *> 00132 *> \param[in,out] A 00133 *> \verbatim 00134 *> A is REAL array, dimension (LDA,N) 00135 *> On entry, the N-by-N matrix A. 00136 *> On exit, A has been overwritten. If JOBVL = 'V' or 00137 *> JOBVR = 'V', A contains the real Schur form of the balanced 00138 *> version of the input matrix A. 00139 *> \endverbatim 00140 *> 00141 *> \param[in] LDA 00142 *> \verbatim 00143 *> LDA is INTEGER 00144 *> The leading dimension of the array A. LDA >= max(1,N). 00145 *> \endverbatim 00146 *> 00147 *> \param[out] WR 00148 *> \verbatim 00149 *> WR is REAL array, dimension (N) 00150 *> \endverbatim 00151 *> 00152 *> \param[out] WI 00153 *> \verbatim 00154 *> WI is REAL array, dimension (N) 00155 *> WR and WI contain the real and imaginary parts, 00156 *> respectively, of the computed eigenvalues. Complex 00157 *> conjugate pairs of eigenvalues will appear consecutively 00158 *> with the eigenvalue having the positive imaginary part 00159 *> first. 00160 *> \endverbatim 00161 *> 00162 *> \param[out] VL 00163 *> \verbatim 00164 *> VL is REAL array, dimension (LDVL,N) 00165 *> If JOBVL = 'V', the left eigenvectors u(j) are stored one 00166 *> after another in the columns of VL, in the same order 00167 *> as their eigenvalues. 00168 *> If JOBVL = 'N', VL is not referenced. 00169 *> If the j-th eigenvalue is real, then u(j) = VL(:,j), 00170 *> the j-th column of VL. 00171 *> If the j-th and (j+1)-st eigenvalues form a complex 00172 *> conjugate pair, then u(j) = VL(:,j) + i*VL(:,j+1) and 00173 *> u(j+1) = VL(:,j) - i*VL(:,j+1). 00174 *> \endverbatim 00175 *> 00176 *> \param[in] LDVL 00177 *> \verbatim 00178 *> LDVL is INTEGER 00179 *> The leading dimension of the array VL. LDVL >= 1; if 00180 *> JOBVL = 'V', LDVL >= N. 00181 *> \endverbatim 00182 *> 00183 *> \param[out] VR 00184 *> \verbatim 00185 *> VR is REAL array, dimension (LDVR,N) 00186 *> If JOBVR = 'V', the right eigenvectors v(j) are stored one 00187 *> after another in the columns of VR, in the same order 00188 *> as their eigenvalues. 00189 *> If JOBVR = 'N', VR is not referenced. 00190 *> If the j-th eigenvalue is real, then v(j) = VR(:,j), 00191 *> the j-th column of VR. 00192 *> If the j-th and (j+1)-st eigenvalues form a complex 00193 *> conjugate pair, then v(j) = VR(:,j) + i*VR(:,j+1) and 00194 *> v(j+1) = VR(:,j) - i*VR(:,j+1). 00195 *> \endverbatim 00196 *> 00197 *> \param[in] LDVR 00198 *> \verbatim 00199 *> LDVR is INTEGER 00200 *> The leading dimension of the array VR. LDVR >= 1, and if 00201 *> JOBVR = 'V', LDVR >= N. 00202 *> \endverbatim 00203 *> 00204 *> \param[out] ILO 00205 *> \verbatim 00206 *> ILO is INTEGER 00207 *> \endverbatim 00208 *> 00209 *> \param[out] IHI 00210 *> \verbatim 00211 *> IHI is INTEGER 00212 *> ILO and IHI are integer values determined when A was 00213 *> balanced. The balanced A(i,j) = 0 if I > J and 00214 *> J = 1,...,ILO-1 or I = IHI+1,...,N. 00215 *> \endverbatim 00216 *> 00217 *> \param[out] SCALE 00218 *> \verbatim 00219 *> SCALE is REAL array, dimension (N) 00220 *> Details of the permutations and scaling factors applied 00221 *> when balancing A. If P(j) is the index of the row and column 00222 *> interchanged with row and column j, and D(j) is the scaling 00223 *> factor applied to row and column j, then 00224 *> SCALE(J) = P(J), for J = 1,...,ILO-1 00225 *> = D(J), for J = ILO,...,IHI 00226 *> = P(J) for J = IHI+1,...,N. 00227 *> The order in which the interchanges are made is N to IHI+1, 00228 *> then 1 to ILO-1. 00229 *> \endverbatim 00230 *> 00231 *> \param[out] ABNRM 00232 *> \verbatim 00233 *> ABNRM is REAL 00234 *> The one-norm of the balanced matrix (the maximum 00235 *> of the sum of absolute values of elements of any column). 00236 *> \endverbatim 00237 *> 00238 *> \param[out] RCONDE 00239 *> \verbatim 00240 *> RCONDE is REAL array, dimension (N) 00241 *> RCONDE(j) is the reciprocal condition number of the j-th 00242 *> eigenvalue. 00243 *> \endverbatim 00244 *> 00245 *> \param[out] RCONDV 00246 *> \verbatim 00247 *> RCONDV is REAL array, dimension (N) 00248 *> RCONDV(j) is the reciprocal condition number of the j-th 00249 *> right eigenvector. 00250 *> \endverbatim 00251 *> 00252 *> \param[out] WORK 00253 *> \verbatim 00254 *> WORK is REAL array, dimension (MAX(1,LWORK)) 00255 *> On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00256 *> \endverbatim 00257 *> 00258 *> \param[in] LWORK 00259 *> \verbatim 00260 *> LWORK is INTEGER 00261 *> The dimension of the array WORK. If SENSE = 'N' or 'E', 00262 *> LWORK >= max(1,2*N), and if JOBVL = 'V' or JOBVR = 'V', 00263 *> LWORK >= 3*N. If SENSE = 'V' or 'B', LWORK >= N*(N+6). 00264 *> For good performance, LWORK must generally be larger. 00265 *> 00266 *> If LWORK = -1, then a workspace query is assumed; the routine 00267 *> only calculates the optimal size of the WORK array, returns 00268 *> this value as the first entry of the WORK array, and no error 00269 *> message related to LWORK is issued by XERBLA. 00270 *> \endverbatim 00271 *> 00272 *> \param[out] IWORK 00273 *> \verbatim 00274 *> IWORK is INTEGER array, dimension (2*N-2) 00275 *> If SENSE = 'N' or 'E', not referenced. 00276 *> \endverbatim 00277 *> 00278 *> \param[out] INFO 00279 *> \verbatim 00280 *> INFO is INTEGER 00281 *> = 0: successful exit 00282 *> < 0: if INFO = -i, the i-th argument had an illegal value. 00283 *> > 0: if INFO = i, the QR algorithm failed to compute all the 00284 *> eigenvalues, and no eigenvectors or condition numbers 00285 *> have been computed; elements 1:ILO-1 and i+1:N of WR 00286 *> and WI contain eigenvalues which have converged. 00287 *> \endverbatim 00288 * 00289 * Authors: 00290 * ======== 00291 * 00292 *> \author Univ. of Tennessee 00293 *> \author Univ. of California Berkeley 00294 *> \author Univ. of Colorado Denver 00295 *> \author NAG Ltd. 00296 * 00297 *> \date November 2011 00298 * 00299 *> \ingroup realGEeigen 00300 * 00301 * ===================================================================== 00302 SUBROUTINE SGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI, 00303 $ VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, 00304 $ RCONDE, RCONDV, WORK, LWORK, IWORK, INFO ) 00305 * 00306 * -- LAPACK driver routine (version 3.4.0) -- 00307 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00308 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00309 * November 2011 00310 * 00311 * .. Scalar Arguments .. 00312 CHARACTER BALANC, JOBVL, JOBVR, SENSE 00313 INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N 00314 REAL ABNRM 00315 * .. 00316 * .. Array Arguments .. 00317 INTEGER IWORK( * ) 00318 REAL A( LDA, * ), RCONDE( * ), RCONDV( * ), 00319 $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ), 00320 $ WI( * ), WORK( * ), WR( * ) 00321 * .. 00322 * 00323 * ===================================================================== 00324 * 00325 * .. Parameters .. 00326 REAL ZERO, ONE 00327 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00328 * .. 00329 * .. Local Scalars .. 00330 LOGICAL LQUERY, SCALEA, WANTVL, WANTVR, WNTSNB, WNTSNE, 00331 $ WNTSNN, WNTSNV 00332 CHARACTER JOB, SIDE 00333 INTEGER HSWORK, I, ICOND, IERR, ITAU, IWRK, K, MAXWRK, 00334 $ MINWRK, NOUT 00335 REAL ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM, 00336 $ SN 00337 * .. 00338 * .. Local Arrays .. 00339 LOGICAL SELECT( 1 ) 00340 REAL DUM( 1 ) 00341 * .. 00342 * .. External Subroutines .. 00343 EXTERNAL SGEBAK, SGEBAL, SGEHRD, SHSEQR, SLABAD, SLACPY, 00344 $ SLARTG, SLASCL, SORGHR, SROT, SSCAL, STREVC, 00345 $ STRSNA, XERBLA 00346 * .. 00347 * .. External Functions .. 00348 LOGICAL LSAME 00349 INTEGER ILAENV, ISAMAX 00350 REAL SLAMCH, SLANGE, SLAPY2, SNRM2 00351 EXTERNAL LSAME, ILAENV, ISAMAX, SLAMCH, SLANGE, SLAPY2, 00352 $ SNRM2 00353 * .. 00354 * .. Intrinsic Functions .. 00355 INTRINSIC MAX, SQRT 00356 * .. 00357 * .. Executable Statements .. 00358 * 00359 * Test the input arguments 00360 * 00361 INFO = 0 00362 LQUERY = ( LWORK.EQ.-1 ) 00363 WANTVL = LSAME( JOBVL, 'V' ) 00364 WANTVR = LSAME( JOBVR, 'V' ) 00365 WNTSNN = LSAME( SENSE, 'N' ) 00366 WNTSNE = LSAME( SENSE, 'E' ) 00367 WNTSNV = LSAME( SENSE, 'V' ) 00368 WNTSNB = LSAME( SENSE, 'B' ) 00369 IF( .NOT.( LSAME( BALANC, 'N' ) .OR. LSAME( BALANC, 'S' ) .OR. 00370 $ LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) THEN 00371 INFO = -1 00372 ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN 00373 INFO = -2 00374 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN 00375 INFO = -3 00376 ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR. 00377 $ ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND. 00378 $ WANTVR ) ) ) THEN 00379 INFO = -4 00380 ELSE IF( N.LT.0 ) THEN 00381 INFO = -5 00382 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00383 INFO = -7 00384 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN 00385 INFO = -11 00386 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN 00387 INFO = -13 00388 END IF 00389 * 00390 * Compute workspace 00391 * (Note: Comments in the code beginning "Workspace:" describe the 00392 * minimal amount of workspace needed at that point in the code, 00393 * as well as the preferred amount for good performance. 00394 * NB refers to the optimal block size for the immediately 00395 * following subroutine, as returned by ILAENV. 00396 * HSWORK refers to the workspace preferred by SHSEQR, as 00397 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 00398 * the worst case.) 00399 * 00400 IF( INFO.EQ.0 ) THEN 00401 IF( N.EQ.0 ) THEN 00402 MINWRK = 1 00403 MAXWRK = 1 00404 ELSE 00405 MAXWRK = N + N*ILAENV( 1, 'SGEHRD', ' ', N, 1, N, 0 ) 00406 * 00407 IF( WANTVL ) THEN 00408 CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL, 00409 $ WORK, -1, INFO ) 00410 ELSE IF( WANTVR ) THEN 00411 CALL SHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR, 00412 $ WORK, -1, INFO ) 00413 ELSE 00414 IF( WNTSNN ) THEN 00415 CALL SHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, 00416 $ LDVR, WORK, -1, INFO ) 00417 ELSE 00418 CALL SHSEQR( 'S', 'N', N, 1, N, A, LDA, WR, WI, VR, 00419 $ LDVR, WORK, -1, INFO ) 00420 END IF 00421 END IF 00422 HSWORK = WORK( 1 ) 00423 * 00424 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN 00425 MINWRK = 2*N 00426 IF( .NOT.WNTSNN ) 00427 $ MINWRK = MAX( MINWRK, N*N+6*N ) 00428 MAXWRK = MAX( MAXWRK, HSWORK ) 00429 IF( .NOT.WNTSNN ) 00430 $ MAXWRK = MAX( MAXWRK, N*N + 6*N ) 00431 ELSE 00432 MINWRK = 3*N 00433 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) ) 00434 $ MINWRK = MAX( MINWRK, N*N + 6*N ) 00435 MAXWRK = MAX( MAXWRK, HSWORK ) 00436 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'SORGHR', 00437 $ ' ', N, 1, N, -1 ) ) 00438 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) ) 00439 $ MAXWRK = MAX( MAXWRK, N*N + 6*N ) 00440 MAXWRK = MAX( MAXWRK, 3*N ) 00441 END IF 00442 MAXWRK = MAX( MAXWRK, MINWRK ) 00443 END IF 00444 WORK( 1 ) = MAXWRK 00445 * 00446 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00447 INFO = -21 00448 END IF 00449 END IF 00450 * 00451 IF( INFO.NE.0 ) THEN 00452 CALL XERBLA( 'SGEEVX', -INFO ) 00453 RETURN 00454 ELSE IF( LQUERY ) THEN 00455 RETURN 00456 END IF 00457 * 00458 * Quick return if possible 00459 * 00460 IF( N.EQ.0 ) 00461 $ RETURN 00462 * 00463 * Get machine constants 00464 * 00465 EPS = SLAMCH( 'P' ) 00466 SMLNUM = SLAMCH( 'S' ) 00467 BIGNUM = ONE / SMLNUM 00468 CALL SLABAD( SMLNUM, BIGNUM ) 00469 SMLNUM = SQRT( SMLNUM ) / EPS 00470 BIGNUM = ONE / SMLNUM 00471 * 00472 * Scale A if max element outside range [SMLNUM,BIGNUM] 00473 * 00474 ICOND = 0 00475 ANRM = SLANGE( 'M', N, N, A, LDA, DUM ) 00476 SCALEA = .FALSE. 00477 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00478 SCALEA = .TRUE. 00479 CSCALE = SMLNUM 00480 ELSE IF( ANRM.GT.BIGNUM ) THEN 00481 SCALEA = .TRUE. 00482 CSCALE = BIGNUM 00483 END IF 00484 IF( SCALEA ) 00485 $ CALL SLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 00486 * 00487 * Balance the matrix and compute ABNRM 00488 * 00489 CALL SGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR ) 00490 ABNRM = SLANGE( '1', N, N, A, LDA, DUM ) 00491 IF( SCALEA ) THEN 00492 DUM( 1 ) = ABNRM 00493 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR ) 00494 ABNRM = DUM( 1 ) 00495 END IF 00496 * 00497 * Reduce to upper Hessenberg form 00498 * (Workspace: need 2*N, prefer N+N*NB) 00499 * 00500 ITAU = 1 00501 IWRK = ITAU + N 00502 CALL SGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 00503 $ LWORK-IWRK+1, IERR ) 00504 * 00505 IF( WANTVL ) THEN 00506 * 00507 * Want left eigenvectors 00508 * Copy Householder vectors to VL 00509 * 00510 SIDE = 'L' 00511 CALL SLACPY( 'L', N, N, A, LDA, VL, LDVL ) 00512 * 00513 * Generate orthogonal matrix in VL 00514 * (Workspace: need 2*N-1, prefer N+(N-1)*NB) 00515 * 00516 CALL SORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ), 00517 $ LWORK-IWRK+1, IERR ) 00518 * 00519 * Perform QR iteration, accumulating Schur vectors in VL 00520 * (Workspace: need 1, prefer HSWORK (see comments) ) 00521 * 00522 IWRK = ITAU 00523 CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL, 00524 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00525 * 00526 IF( WANTVR ) THEN 00527 * 00528 * Want left and right eigenvectors 00529 * Copy Schur vectors to VR 00530 * 00531 SIDE = 'B' 00532 CALL SLACPY( 'F', N, N, VL, LDVL, VR, LDVR ) 00533 END IF 00534 * 00535 ELSE IF( WANTVR ) THEN 00536 * 00537 * Want right eigenvectors 00538 * Copy Householder vectors to VR 00539 * 00540 SIDE = 'R' 00541 CALL SLACPY( 'L', N, N, A, LDA, VR, LDVR ) 00542 * 00543 * Generate orthogonal matrix in VR 00544 * (Workspace: need 2*N-1, prefer N+(N-1)*NB) 00545 * 00546 CALL SORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ), 00547 $ LWORK-IWRK+1, IERR ) 00548 * 00549 * Perform QR iteration, accumulating Schur vectors in VR 00550 * (Workspace: need 1, prefer HSWORK (see comments) ) 00551 * 00552 IWRK = ITAU 00553 CALL SHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 00554 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00555 * 00556 ELSE 00557 * 00558 * Compute eigenvalues only 00559 * If condition numbers desired, compute Schur form 00560 * 00561 IF( WNTSNN ) THEN 00562 JOB = 'E' 00563 ELSE 00564 JOB = 'S' 00565 END IF 00566 * 00567 * (Workspace: need 1, prefer HSWORK (see comments) ) 00568 * 00569 IWRK = ITAU 00570 CALL SHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 00571 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00572 END IF 00573 * 00574 * If INFO > 0 from SHSEQR, then quit 00575 * 00576 IF( INFO.GT.0 ) 00577 $ GO TO 50 00578 * 00579 IF( WANTVL .OR. WANTVR ) THEN 00580 * 00581 * Compute left and/or right eigenvectors 00582 * (Workspace: need 3*N) 00583 * 00584 CALL STREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 00585 $ N, NOUT, WORK( IWRK ), IERR ) 00586 END IF 00587 * 00588 * Compute condition numbers if desired 00589 * (Workspace: need N*N+6*N unless SENSE = 'E') 00590 * 00591 IF( .NOT.WNTSNN ) THEN 00592 CALL STRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 00593 $ RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, IWORK, 00594 $ ICOND ) 00595 END IF 00596 * 00597 IF( WANTVL ) THEN 00598 * 00599 * Undo balancing of left eigenvectors 00600 * 00601 CALL SGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL, 00602 $ IERR ) 00603 * 00604 * Normalize left eigenvectors and make largest component real 00605 * 00606 DO 20 I = 1, N 00607 IF( WI( I ).EQ.ZERO ) THEN 00608 SCL = ONE / SNRM2( N, VL( 1, I ), 1 ) 00609 CALL SSCAL( N, SCL, VL( 1, I ), 1 ) 00610 ELSE IF( WI( I ).GT.ZERO ) THEN 00611 SCL = ONE / SLAPY2( SNRM2( N, VL( 1, I ), 1 ), 00612 $ SNRM2( N, VL( 1, I+1 ), 1 ) ) 00613 CALL SSCAL( N, SCL, VL( 1, I ), 1 ) 00614 CALL SSCAL( N, SCL, VL( 1, I+1 ), 1 ) 00615 DO 10 K = 1, N 00616 WORK( K ) = VL( K, I )**2 + VL( K, I+1 )**2 00617 10 CONTINUE 00618 K = ISAMAX( N, WORK, 1 ) 00619 CALL SLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R ) 00620 CALL SROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN ) 00621 VL( K, I+1 ) = ZERO 00622 END IF 00623 20 CONTINUE 00624 END IF 00625 * 00626 IF( WANTVR ) THEN 00627 * 00628 * Undo balancing of right eigenvectors 00629 * 00630 CALL SGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR, 00631 $ IERR ) 00632 * 00633 * Normalize right eigenvectors and make largest component real 00634 * 00635 DO 40 I = 1, N 00636 IF( WI( I ).EQ.ZERO ) THEN 00637 SCL = ONE / SNRM2( N, VR( 1, I ), 1 ) 00638 CALL SSCAL( N, SCL, VR( 1, I ), 1 ) 00639 ELSE IF( WI( I ).GT.ZERO ) THEN 00640 SCL = ONE / SLAPY2( SNRM2( N, VR( 1, I ), 1 ), 00641 $ SNRM2( N, VR( 1, I+1 ), 1 ) ) 00642 CALL SSCAL( N, SCL, VR( 1, I ), 1 ) 00643 CALL SSCAL( N, SCL, VR( 1, I+1 ), 1 ) 00644 DO 30 K = 1, N 00645 WORK( K ) = VR( K, I )**2 + VR( K, I+1 )**2 00646 30 CONTINUE 00647 K = ISAMAX( N, WORK, 1 ) 00648 CALL SLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R ) 00649 CALL SROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN ) 00650 VR( K, I+1 ) = ZERO 00651 END IF 00652 40 CONTINUE 00653 END IF 00654 * 00655 * Undo scaling if necessary 00656 * 00657 50 CONTINUE 00658 IF( SCALEA ) THEN 00659 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ), 00660 $ MAX( N-INFO, 1 ), IERR ) 00661 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ), 00662 $ MAX( N-INFO, 1 ), IERR ) 00663 IF( INFO.EQ.0 ) THEN 00664 IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 ) 00665 $ CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N, 00666 $ IERR ) 00667 ELSE 00668 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N, 00669 $ IERR ) 00670 CALL SLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N, 00671 $ IERR ) 00672 END IF 00673 END IF 00674 * 00675 WORK( 1 ) = MAXWRK 00676 RETURN 00677 * 00678 * End of SGEEVX 00679 * 00680 END