![]() |
LAPACK
3.4.1
LAPACK: Linear Algebra PACKage
|
00001 *> \brief <b> DGEEVX 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 DGEEVX + dependencies 00010 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.tgz?format=tgz&filename=/lapack/lapack_routine/dgeevx.f"> 00011 *> [TGZ]</a> 00012 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.zip?format=zip&filename=/lapack/lapack_routine/dgeevx.f"> 00013 *> [ZIP]</a> 00014 *> <a href="http://www.netlib.org/cgi-bin/netlibfiles.txt?format=txt&filename=/lapack/lapack_routine/dgeevx.f"> 00015 *> [TXT]</a> 00016 *> \endhtmlonly 00017 * 00018 * Definition: 00019 * =========== 00020 * 00021 * SUBROUTINE DGEEVX( 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 * DOUBLE PRECISION ABNRM 00029 * .. 00030 * .. Array Arguments .. 00031 * INTEGER IWORK( * ) 00032 * DOUBLE PRECISION 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 *> DGEEVX 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 DOUBLE PRECISION 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 DOUBLE PRECISION array, dimension (N) 00150 *> \endverbatim 00151 *> 00152 *> \param[out] WI 00153 *> \verbatim 00154 *> WI is DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 DOUBLE PRECISION 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 doubleGEeigen 00300 * 00301 * ===================================================================== 00302 SUBROUTINE DGEEVX( 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 DOUBLE PRECISION ABNRM 00315 * .. 00316 * .. Array Arguments .. 00317 INTEGER IWORK( * ) 00318 DOUBLE PRECISION A( LDA, * ), RCONDE( * ), RCONDV( * ), 00319 $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ), 00320 $ WI( * ), WORK( * ), WR( * ) 00321 * .. 00322 * 00323 * ===================================================================== 00324 * 00325 * .. Parameters .. 00326 DOUBLE PRECISION ZERO, ONE 00327 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 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 DOUBLE PRECISION ANRM, BIGNUM, CS, CSCALE, EPS, R, SCL, SMLNUM, 00336 $ SN 00337 * .. 00338 * .. Local Arrays .. 00339 LOGICAL SELECT( 1 ) 00340 DOUBLE PRECISION DUM( 1 ) 00341 * .. 00342 * .. External Subroutines .. 00343 EXTERNAL DGEBAK, DGEBAL, DGEHRD, DHSEQR, DLABAD, DLACPY, 00344 $ DLARTG, DLASCL, DORGHR, DROT, DSCAL, DTREVC, 00345 $ DTRSNA, XERBLA 00346 * .. 00347 * .. External Functions .. 00348 LOGICAL LSAME 00349 INTEGER IDAMAX, ILAENV 00350 DOUBLE PRECISION DLAMCH, DLANGE, DLAPY2, DNRM2 00351 EXTERNAL LSAME, IDAMAX, ILAENV, DLAMCH, DLANGE, DLAPY2, 00352 $ DNRM2 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, 00370 $ 'S' ) .OR. LSAME( BALANC, 'P' ) .OR. LSAME( BALANC, 'B' ) ) ) 00371 $ THEN 00372 INFO = -1 00373 ELSE IF( ( .NOT.WANTVL ) .AND. ( .NOT.LSAME( JOBVL, 'N' ) ) ) THEN 00374 INFO = -2 00375 ELSE IF( ( .NOT.WANTVR ) .AND. ( .NOT.LSAME( JOBVR, 'N' ) ) ) THEN 00376 INFO = -3 00377 ELSE IF( .NOT.( WNTSNN .OR. WNTSNE .OR. WNTSNB .OR. WNTSNV ) .OR. 00378 $ ( ( WNTSNE .OR. WNTSNB ) .AND. .NOT.( WANTVL .AND. 00379 $ WANTVR ) ) ) THEN 00380 INFO = -4 00381 ELSE IF( N.LT.0 ) THEN 00382 INFO = -5 00383 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00384 INFO = -7 00385 ELSE IF( LDVL.LT.1 .OR. ( WANTVL .AND. LDVL.LT.N ) ) THEN 00386 INFO = -11 00387 ELSE IF( LDVR.LT.1 .OR. ( WANTVR .AND. LDVR.LT.N ) ) THEN 00388 INFO = -13 00389 END IF 00390 * 00391 * Compute workspace 00392 * (Note: Comments in the code beginning "Workspace:" describe the 00393 * minimal amount of workspace needed at that point in the code, 00394 * as well as the preferred amount for good performance. 00395 * NB refers to the optimal block size for the immediately 00396 * following subroutine, as returned by ILAENV. 00397 * HSWORK refers to the workspace preferred by DHSEQR, as 00398 * calculated below. HSWORK is computed assuming ILO=1 and IHI=N, 00399 * the worst case.) 00400 * 00401 IF( INFO.EQ.0 ) THEN 00402 IF( N.EQ.0 ) THEN 00403 MINWRK = 1 00404 MAXWRK = 1 00405 ELSE 00406 MAXWRK = N + N*ILAENV( 1, 'DGEHRD', ' ', N, 1, N, 0 ) 00407 * 00408 IF( WANTVL ) THEN 00409 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VL, LDVL, 00410 $ WORK, -1, INFO ) 00411 ELSE IF( WANTVR ) THEN 00412 CALL DHSEQR( 'S', 'V', N, 1, N, A, LDA, WR, WI, VR, LDVR, 00413 $ WORK, -1, INFO ) 00414 ELSE 00415 IF( WNTSNN ) THEN 00416 CALL DHSEQR( 'E', 'N', N, 1, N, A, LDA, WR, WI, VR, 00417 $ LDVR, WORK, -1, INFO ) 00418 ELSE 00419 CALL DHSEQR( 'S', 'N', N, 1, N, A, LDA, WR, WI, VR, 00420 $ LDVR, WORK, -1, INFO ) 00421 END IF 00422 END IF 00423 HSWORK = WORK( 1 ) 00424 * 00425 IF( ( .NOT.WANTVL ) .AND. ( .NOT.WANTVR ) ) THEN 00426 MINWRK = 2*N 00427 IF( .NOT.WNTSNN ) 00428 $ MINWRK = MAX( MINWRK, N*N+6*N ) 00429 MAXWRK = MAX( MAXWRK, HSWORK ) 00430 IF( .NOT.WNTSNN ) 00431 $ MAXWRK = MAX( MAXWRK, N*N + 6*N ) 00432 ELSE 00433 MINWRK = 3*N 00434 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) ) 00435 $ MINWRK = MAX( MINWRK, N*N + 6*N ) 00436 MAXWRK = MAX( MAXWRK, HSWORK ) 00437 MAXWRK = MAX( MAXWRK, N + ( N - 1 )*ILAENV( 1, 'DORGHR', 00438 $ ' ', N, 1, N, -1 ) ) 00439 IF( ( .NOT.WNTSNN ) .AND. ( .NOT.WNTSNE ) ) 00440 $ MAXWRK = MAX( MAXWRK, N*N + 6*N ) 00441 MAXWRK = MAX( MAXWRK, 3*N ) 00442 END IF 00443 MAXWRK = MAX( MAXWRK, MINWRK ) 00444 END IF 00445 WORK( 1 ) = MAXWRK 00446 * 00447 IF( LWORK.LT.MINWRK .AND. .NOT.LQUERY ) THEN 00448 INFO = -21 00449 END IF 00450 END IF 00451 * 00452 IF( INFO.NE.0 ) THEN 00453 CALL XERBLA( 'DGEEVX', -INFO ) 00454 RETURN 00455 ELSE IF( LQUERY ) THEN 00456 RETURN 00457 END IF 00458 * 00459 * Quick return if possible 00460 * 00461 IF( N.EQ.0 ) 00462 $ RETURN 00463 * 00464 * Get machine constants 00465 * 00466 EPS = DLAMCH( 'P' ) 00467 SMLNUM = DLAMCH( 'S' ) 00468 BIGNUM = ONE / SMLNUM 00469 CALL DLABAD( SMLNUM, BIGNUM ) 00470 SMLNUM = SQRT( SMLNUM ) / EPS 00471 BIGNUM = ONE / SMLNUM 00472 * 00473 * Scale A if max element outside range [SMLNUM,BIGNUM] 00474 * 00475 ICOND = 0 00476 ANRM = DLANGE( 'M', N, N, A, LDA, DUM ) 00477 SCALEA = .FALSE. 00478 IF( ANRM.GT.ZERO .AND. ANRM.LT.SMLNUM ) THEN 00479 SCALEA = .TRUE. 00480 CSCALE = SMLNUM 00481 ELSE IF( ANRM.GT.BIGNUM ) THEN 00482 SCALEA = .TRUE. 00483 CSCALE = BIGNUM 00484 END IF 00485 IF( SCALEA ) 00486 $ CALL DLASCL( 'G', 0, 0, ANRM, CSCALE, N, N, A, LDA, IERR ) 00487 * 00488 * Balance the matrix and compute ABNRM 00489 * 00490 CALL DGEBAL( BALANC, N, A, LDA, ILO, IHI, SCALE, IERR ) 00491 ABNRM = DLANGE( '1', N, N, A, LDA, DUM ) 00492 IF( SCALEA ) THEN 00493 DUM( 1 ) = ABNRM 00494 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, 1, 1, DUM, 1, IERR ) 00495 ABNRM = DUM( 1 ) 00496 END IF 00497 * 00498 * Reduce to upper Hessenberg form 00499 * (Workspace: need 2*N, prefer N+N*NB) 00500 * 00501 ITAU = 1 00502 IWRK = ITAU + N 00503 CALL DGEHRD( N, ILO, IHI, A, LDA, WORK( ITAU ), WORK( IWRK ), 00504 $ LWORK-IWRK+1, IERR ) 00505 * 00506 IF( WANTVL ) THEN 00507 * 00508 * Want left eigenvectors 00509 * Copy Householder vectors to VL 00510 * 00511 SIDE = 'L' 00512 CALL DLACPY( 'L', N, N, A, LDA, VL, LDVL ) 00513 * 00514 * Generate orthogonal matrix in VL 00515 * (Workspace: need 2*N-1, prefer N+(N-1)*NB) 00516 * 00517 CALL DORGHR( N, ILO, IHI, VL, LDVL, WORK( ITAU ), WORK( IWRK ), 00518 $ LWORK-IWRK+1, IERR ) 00519 * 00520 * Perform QR iteration, accumulating Schur vectors in VL 00521 * (Workspace: need 1, prefer HSWORK (see comments) ) 00522 * 00523 IWRK = ITAU 00524 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VL, LDVL, 00525 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00526 * 00527 IF( WANTVR ) THEN 00528 * 00529 * Want left and right eigenvectors 00530 * Copy Schur vectors to VR 00531 * 00532 SIDE = 'B' 00533 CALL DLACPY( 'F', N, N, VL, LDVL, VR, LDVR ) 00534 END IF 00535 * 00536 ELSE IF( WANTVR ) THEN 00537 * 00538 * Want right eigenvectors 00539 * Copy Householder vectors to VR 00540 * 00541 SIDE = 'R' 00542 CALL DLACPY( 'L', N, N, A, LDA, VR, LDVR ) 00543 * 00544 * Generate orthogonal matrix in VR 00545 * (Workspace: need 2*N-1, prefer N+(N-1)*NB) 00546 * 00547 CALL DORGHR( N, ILO, IHI, VR, LDVR, WORK( ITAU ), WORK( IWRK ), 00548 $ LWORK-IWRK+1, IERR ) 00549 * 00550 * Perform QR iteration, accumulating Schur vectors in VR 00551 * (Workspace: need 1, prefer HSWORK (see comments) ) 00552 * 00553 IWRK = ITAU 00554 CALL DHSEQR( 'S', 'V', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 00555 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00556 * 00557 ELSE 00558 * 00559 * Compute eigenvalues only 00560 * If condition numbers desired, compute Schur form 00561 * 00562 IF( WNTSNN ) THEN 00563 JOB = 'E' 00564 ELSE 00565 JOB = 'S' 00566 END IF 00567 * 00568 * (Workspace: need 1, prefer HSWORK (see comments) ) 00569 * 00570 IWRK = ITAU 00571 CALL DHSEQR( JOB, 'N', N, ILO, IHI, A, LDA, WR, WI, VR, LDVR, 00572 $ WORK( IWRK ), LWORK-IWRK+1, INFO ) 00573 END IF 00574 * 00575 * If INFO > 0 from DHSEQR, then quit 00576 * 00577 IF( INFO.GT.0 ) 00578 $ GO TO 50 00579 * 00580 IF( WANTVL .OR. WANTVR ) THEN 00581 * 00582 * Compute left and/or right eigenvectors 00583 * (Workspace: need 3*N) 00584 * 00585 CALL DTREVC( SIDE, 'B', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 00586 $ N, NOUT, WORK( IWRK ), IERR ) 00587 END IF 00588 * 00589 * Compute condition numbers if desired 00590 * (Workspace: need N*N+6*N unless SENSE = 'E') 00591 * 00592 IF( .NOT.WNTSNN ) THEN 00593 CALL DTRSNA( SENSE, 'A', SELECT, N, A, LDA, VL, LDVL, VR, LDVR, 00594 $ RCONDE, RCONDV, N, NOUT, WORK( IWRK ), N, IWORK, 00595 $ ICOND ) 00596 END IF 00597 * 00598 IF( WANTVL ) THEN 00599 * 00600 * Undo balancing of left eigenvectors 00601 * 00602 CALL DGEBAK( BALANC, 'L', N, ILO, IHI, SCALE, N, VL, LDVL, 00603 $ IERR ) 00604 * 00605 * Normalize left eigenvectors and make largest component real 00606 * 00607 DO 20 I = 1, N 00608 IF( WI( I ).EQ.ZERO ) THEN 00609 SCL = ONE / DNRM2( N, VL( 1, I ), 1 ) 00610 CALL DSCAL( N, SCL, VL( 1, I ), 1 ) 00611 ELSE IF( WI( I ).GT.ZERO ) THEN 00612 SCL = ONE / DLAPY2( DNRM2( N, VL( 1, I ), 1 ), 00613 $ DNRM2( N, VL( 1, I+1 ), 1 ) ) 00614 CALL DSCAL( N, SCL, VL( 1, I ), 1 ) 00615 CALL DSCAL( N, SCL, VL( 1, I+1 ), 1 ) 00616 DO 10 K = 1, N 00617 WORK( K ) = VL( K, I )**2 + VL( K, I+1 )**2 00618 10 CONTINUE 00619 K = IDAMAX( N, WORK, 1 ) 00620 CALL DLARTG( VL( K, I ), VL( K, I+1 ), CS, SN, R ) 00621 CALL DROT( N, VL( 1, I ), 1, VL( 1, I+1 ), 1, CS, SN ) 00622 VL( K, I+1 ) = ZERO 00623 END IF 00624 20 CONTINUE 00625 END IF 00626 * 00627 IF( WANTVR ) THEN 00628 * 00629 * Undo balancing of right eigenvectors 00630 * 00631 CALL DGEBAK( BALANC, 'R', N, ILO, IHI, SCALE, N, VR, LDVR, 00632 $ IERR ) 00633 * 00634 * Normalize right eigenvectors and make largest component real 00635 * 00636 DO 40 I = 1, N 00637 IF( WI( I ).EQ.ZERO ) THEN 00638 SCL = ONE / DNRM2( N, VR( 1, I ), 1 ) 00639 CALL DSCAL( N, SCL, VR( 1, I ), 1 ) 00640 ELSE IF( WI( I ).GT.ZERO ) THEN 00641 SCL = ONE / DLAPY2( DNRM2( N, VR( 1, I ), 1 ), 00642 $ DNRM2( N, VR( 1, I+1 ), 1 ) ) 00643 CALL DSCAL( N, SCL, VR( 1, I ), 1 ) 00644 CALL DSCAL( N, SCL, VR( 1, I+1 ), 1 ) 00645 DO 30 K = 1, N 00646 WORK( K ) = VR( K, I )**2 + VR( K, I+1 )**2 00647 30 CONTINUE 00648 K = IDAMAX( N, WORK, 1 ) 00649 CALL DLARTG( VR( K, I ), VR( K, I+1 ), CS, SN, R ) 00650 CALL DROT( N, VR( 1, I ), 1, VR( 1, I+1 ), 1, CS, SN ) 00651 VR( K, I+1 ) = ZERO 00652 END IF 00653 40 CONTINUE 00654 END IF 00655 * 00656 * Undo scaling if necessary 00657 * 00658 50 CONTINUE 00659 IF( SCALEA ) THEN 00660 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WR( INFO+1 ), 00661 $ MAX( N-INFO, 1 ), IERR ) 00662 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N-INFO, 1, WI( INFO+1 ), 00663 $ MAX( N-INFO, 1 ), IERR ) 00664 IF( INFO.EQ.0 ) THEN 00665 IF( ( WNTSNV .OR. WNTSNB ) .AND. ICOND.EQ.0 ) 00666 $ CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, N, 1, RCONDV, N, 00667 $ IERR ) 00668 ELSE 00669 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WR, N, 00670 $ IERR ) 00671 CALL DLASCL( 'G', 0, 0, CSCALE, ANRM, ILO-1, 1, WI, N, 00672 $ IERR ) 00673 END IF 00674 END IF 00675 * 00676 WORK( 1 ) = MAXWRK 00677 RETURN 00678 * 00679 * End of DGEEVX 00680 * 00681 END