 |
LAPACK 3.3.1
Linear Algebra PACKage
|
|
LAPACK 3.3.1
Linear Algebra PACKage
|
00001 SUBROUTINE DGEGV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI, 00002 $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO ) 00003 * 00004 * -- LAPACK driver routine (version 3.2) -- 00005 * -- LAPACK is a software package provided by Univ. of Tennessee, -- 00006 * -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..-- 00007 * November 2006 00008 * 00009 * .. Scalar Arguments .. 00010 CHARACTER JOBVL, JOBVR 00011 INTEGER INFO, LDA, LDB, LDVL, LDVR, LWORK, N 00012 * .. 00013 * .. Array Arguments .. 00014 DOUBLE PRECISION A( LDA, * ), ALPHAI( * ), ALPHAR( * ), 00015 $ B( LDB, * ), BETA( * ), VL( LDVL, * ), 00016 $ VR( LDVR, * ), WORK( * ) 00017 * .. 00018 * 00019 * Purpose 00020 * ======= 00021 * 00022 * This routine is deprecated and has been replaced by routine DGGEV. 00023 * 00024 * DGEGV computes the eigenvalues and, optionally, the left and/or right 00025 * eigenvectors of a real matrix pair (A,B). 00026 * Given two square matrices A and B, 00027 * the generalized nonsymmetric eigenvalue problem (GNEP) is to find the 00028 * eigenvalues lambda and corresponding (non-zero) eigenvectors x such 00029 * that 00030 * 00031 * A*x = lambda*B*x. 00032 * 00033 * An alternate form is to find the eigenvalues mu and corresponding 00034 * eigenvectors y such that 00035 * 00036 * mu*A*y = B*y. 00037 * 00038 * These two forms are equivalent with mu = 1/lambda and x = y if 00039 * neither lambda nor mu is zero. In order to deal with the case that 00040 * lambda or mu is zero or small, two values alpha and beta are returned 00041 * for each eigenvalue, such that lambda = alpha/beta and 00042 * mu = beta/alpha. 00043 * 00044 * The vectors x and y in the above equations are right eigenvectors of 00045 * the matrix pair (A,B). Vectors u and v satisfying 00046 * 00047 * u**H*A = lambda*u**H*B or mu*v**H*A = v**H*B 00048 * 00049 * are left eigenvectors of (A,B). 00050 * 00051 * Note: this routine performs "full balancing" on A and B -- see 00052 * "Further Details", below. 00053 * 00054 * Arguments 00055 * ========= 00056 * 00057 * JOBVL (input) CHARACTER*1 00058 * = 'N': do not compute the left generalized eigenvectors; 00059 * = 'V': compute the left generalized eigenvectors (returned 00060 * in VL). 00061 * 00062 * JOBVR (input) CHARACTER*1 00063 * = 'N': do not compute the right generalized eigenvectors; 00064 * = 'V': compute the right generalized eigenvectors (returned 00065 * in VR). 00066 * 00067 * N (input) INTEGER 00068 * The order of the matrices A, B, VL, and VR. N >= 0. 00069 * 00070 * A (input/output) DOUBLE PRECISION array, dimension (LDA, N) 00071 * On entry, the matrix A. 00072 * If JOBVL = 'V' or JOBVR = 'V', then on exit A 00073 * contains the real Schur form of A from the generalized Schur 00074 * factorization of the pair (A,B) after balancing. 00075 * If no eigenvectors were computed, then only the diagonal 00076 * blocks from the Schur form will be correct. See DGGHRD and 00077 * DHGEQZ for details. 00078 * 00079 * LDA (input) INTEGER 00080 * The leading dimension of A. LDA >= max(1,N). 00081 * 00082 * B (input/output) DOUBLE PRECISION array, dimension (LDB, N) 00083 * On entry, the matrix B. 00084 * If JOBVL = 'V' or JOBVR = 'V', then on exit B contains the 00085 * upper triangular matrix obtained from B in the generalized 00086 * Schur factorization of the pair (A,B) after balancing. 00087 * If no eigenvectors were computed, then only those elements of 00088 * B corresponding to the diagonal blocks from the Schur form of 00089 * A will be correct. See DGGHRD and DHGEQZ for details. 00090 * 00091 * LDB (input) INTEGER 00092 * The leading dimension of B. LDB >= max(1,N). 00093 * 00094 * ALPHAR (output) DOUBLE PRECISION array, dimension (N) 00095 * The real parts of each scalar alpha defining an eigenvalue of 00096 * GNEP. 00097 * 00098 * ALPHAI (output) DOUBLE PRECISION array, dimension (N) 00099 * The imaginary parts of each scalar alpha defining an 00100 * eigenvalue of GNEP. If ALPHAI(j) is zero, then the j-th 00101 * eigenvalue is real; if positive, then the j-th and 00102 * (j+1)-st eigenvalues are a complex conjugate pair, with 00103 * ALPHAI(j+1) = -ALPHAI(j). 00104 * 00105 * BETA (output) DOUBLE PRECISION array, dimension (N) 00106 * The scalars beta that define the eigenvalues of GNEP. 00107 * 00108 * Together, the quantities alpha = (ALPHAR(j),ALPHAI(j)) and 00109 * beta = BETA(j) represent the j-th eigenvalue of the matrix 00110 * pair (A,B), in one of the forms lambda = alpha/beta or 00111 * mu = beta/alpha. Since either lambda or mu may overflow, 00112 * they should not, in general, be computed. 00113 * 00114 * VL (output) DOUBLE PRECISION array, dimension (LDVL,N) 00115 * If JOBVL = 'V', the left eigenvectors u(j) are stored 00116 * in the columns of VL, in the same order as their eigenvalues. 00117 * If the j-th eigenvalue is real, then u(j) = VL(:,j). 00118 * If the j-th and (j+1)-st eigenvalues form a complex conjugate 00119 * pair, then 00120 * u(j) = VL(:,j) + i*VL(:,j+1) 00121 * and 00122 * u(j+1) = VL(:,j) - i*VL(:,j+1). 00123 * 00124 * Each eigenvector is scaled so that its largest component has 00125 * abs(real part) + abs(imag. part) = 1, except for eigenvectors 00126 * corresponding to an eigenvalue with alpha = beta = 0, which 00127 * are set to zero. 00128 * Not referenced if JOBVL = 'N'. 00129 * 00130 * LDVL (input) INTEGER 00131 * The leading dimension of the matrix VL. LDVL >= 1, and 00132 * if JOBVL = 'V', LDVL >= N. 00133 * 00134 * VR (output) DOUBLE PRECISION array, dimension (LDVR,N) 00135 * If JOBVR = 'V', the right eigenvectors x(j) are stored 00136 * in the columns of VR, in the same order as their eigenvalues. 00137 * If the j-th eigenvalue is real, then x(j) = VR(:,j). 00138 * If the j-th and (j+1)-st eigenvalues form a complex conjugate 00139 * pair, then 00140 * x(j) = VR(:,j) + i*VR(:,j+1) 00141 * and 00142 * x(j+1) = VR(:,j) - i*VR(:,j+1). 00143 * 00144 * Each eigenvector is scaled so that its largest component has 00145 * abs(real part) + abs(imag. part) = 1, except for eigenvalues 00146 * corresponding to an eigenvalue with alpha = beta = 0, which 00147 * are set to zero. 00148 * Not referenced if JOBVR = 'N'. 00149 * 00150 * LDVR (input) INTEGER 00151 * The leading dimension of the matrix VR. LDVR >= 1, and 00152 * if JOBVR = 'V', LDVR >= N. 00153 * 00154 * WORK (workspace/output) DOUBLE PRECISION array, dimension (MAX(1,LWORK)) 00155 * On exit, if INFO = 0, WORK(1) returns the optimal LWORK. 00156 * 00157 * LWORK (input) INTEGER 00158 * The dimension of the array WORK. LWORK >= max(1,8*N). 00159 * For good performance, LWORK must generally be larger. 00160 * To compute the optimal value of LWORK, call ILAENV to get 00161 * blocksizes (for DGEQRF, DORMQR, and DORGQR.) Then compute: 00162 * NB -- MAX of the blocksizes for DGEQRF, DORMQR, and DORGQR; 00163 * The optimal LWORK is: 00164 * 2*N + MAX( 6*N, N*(NB+1) ). 00165 * 00166 * If LWORK = -1, then a workspace query is assumed; the routine 00167 * only calculates the optimal size of the WORK array, returns 00168 * this value as the first entry of the WORK array, and no error 00169 * message related to LWORK is issued by XERBLA. 00170 * 00171 * INFO (output) INTEGER 00172 * = 0: successful exit 00173 * < 0: if INFO = -i, the i-th argument had an illegal value. 00174 * = 1,...,N: 00175 * The QZ iteration failed. No eigenvectors have been 00176 * calculated, but ALPHAR(j), ALPHAI(j), and BETA(j) 00177 * should be correct for j=INFO+1,...,N. 00178 * > N: errors that usually indicate LAPACK problems: 00179 * =N+1: error return from DGGBAL 00180 * =N+2: error return from DGEQRF 00181 * =N+3: error return from DORMQR 00182 * =N+4: error return from DORGQR 00183 * =N+5: error return from DGGHRD 00184 * =N+6: error return from DHGEQZ (other than failed 00185 * iteration) 00186 * =N+7: error return from DTGEVC 00187 * =N+8: error return from DGGBAK (computing VL) 00188 * =N+9: error return from DGGBAK (computing VR) 00189 * =N+10: error return from DLASCL (various calls) 00190 * 00191 * Further Details 00192 * =============== 00193 * 00194 * Balancing 00195 * --------- 00196 * 00197 * This driver calls DGGBAL to both permute and scale rows and columns 00198 * of A and B. The permutations PL and PR are chosen so that PL*A*PR 00199 * and PL*B*R will be upper triangular except for the diagonal blocks 00200 * A(i:j,i:j) and B(i:j,i:j), with i and j as close together as 00201 * possible. The diagonal scaling matrices DL and DR are chosen so 00202 * that the pair DL*PL*A*PR*DR, DL*PL*B*PR*DR have elements close to 00203 * one (except for the elements that start out zero.) 00204 * 00205 * After the eigenvalues and eigenvectors of the balanced matrices 00206 * have been computed, DGGBAK transforms the eigenvectors back to what 00207 * they would have been (in perfect arithmetic) if they had not been 00208 * balanced. 00209 * 00210 * Contents of A and B on Exit 00211 * -------- -- - --- - -- ---- 00212 * 00213 * If any eigenvectors are computed (either JOBVL='V' or JOBVR='V' or 00214 * both), then on exit the arrays A and B will contain the real Schur 00215 * form[*] of the "balanced" versions of A and B. If no eigenvectors 00216 * are computed, then only the diagonal blocks will be correct. 00217 * 00218 * [*] See DHGEQZ, DGEGS, or read the book "Matrix Computations", 00219 * by Golub & van Loan, pub. by Johns Hopkins U. Press. 00220 * 00221 * ===================================================================== 00222 * 00223 * .. Parameters .. 00224 DOUBLE PRECISION ZERO, ONE 00225 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00226 * .. 00227 * .. Local Scalars .. 00228 LOGICAL ILIMIT, ILV, ILVL, ILVR, LQUERY 00229 CHARACTER CHTEMP 00230 INTEGER ICOLS, IHI, IINFO, IJOBVL, IJOBVR, ILEFT, ILO, 00231 $ IN, IRIGHT, IROWS, ITAU, IWORK, JC, JR, LOPT, 00232 $ LWKMIN, LWKOPT, NB, NB1, NB2, NB3 00233 DOUBLE PRECISION ABSAI, ABSAR, ABSB, ANRM, ANRM1, ANRM2, BNRM, 00234 $ BNRM1, BNRM2, EPS, ONEPLS, SAFMAX, SAFMIN, 00235 $ SALFAI, SALFAR, SBETA, SCALE, TEMP 00236 * .. 00237 * .. Local Arrays .. 00238 LOGICAL LDUMMA( 1 ) 00239 * .. 00240 * .. External Subroutines .. 00241 EXTERNAL DGEQRF, DGGBAK, DGGBAL, DGGHRD, DHGEQZ, DLACPY, 00242 $ DLASCL, DLASET, DORGQR, DORMQR, DTGEVC, XERBLA 00243 * .. 00244 * .. External Functions .. 00245 LOGICAL LSAME 00246 INTEGER ILAENV 00247 DOUBLE PRECISION DLAMCH, DLANGE 00248 EXTERNAL LSAME, ILAENV, DLAMCH, DLANGE 00249 * .. 00250 * .. Intrinsic Functions .. 00251 INTRINSIC ABS, INT, MAX 00252 * .. 00253 * .. Executable Statements .. 00254 * 00255 * Decode the input arguments 00256 * 00257 IF( LSAME( JOBVL, 'N' ) ) THEN 00258 IJOBVL = 1 00259 ILVL = .FALSE. 00260 ELSE IF( LSAME( JOBVL, 'V' ) ) THEN 00261 IJOBVL = 2 00262 ILVL = .TRUE. 00263 ELSE 00264 IJOBVL = -1 00265 ILVL = .FALSE. 00266 END IF 00267 * 00268 IF( LSAME( JOBVR, 'N' ) ) THEN 00269 IJOBVR = 1 00270 ILVR = .FALSE. 00271 ELSE IF( LSAME( JOBVR, 'V' ) ) THEN 00272 IJOBVR = 2 00273 ILVR = .TRUE. 00274 ELSE 00275 IJOBVR = -1 00276 ILVR = .FALSE. 00277 END IF 00278 ILV = ILVL .OR. ILVR 00279 * 00280 * Test the input arguments 00281 * 00282 LWKMIN = MAX( 8*N, 1 ) 00283 LWKOPT = LWKMIN 00284 WORK( 1 ) = LWKOPT 00285 LQUERY = ( LWORK.EQ.-1 ) 00286 INFO = 0 00287 IF( IJOBVL.LE.0 ) THEN 00288 INFO = -1 00289 ELSE IF( IJOBVR.LE.0 ) THEN 00290 INFO = -2 00291 ELSE IF( N.LT.0 ) THEN 00292 INFO = -3 00293 ELSE IF( LDA.LT.MAX( 1, N ) ) THEN 00294 INFO = -5 00295 ELSE IF( LDB.LT.MAX( 1, N ) ) THEN 00296 INFO = -7 00297 ELSE IF( LDVL.LT.1 .OR. ( ILVL .AND. LDVL.LT.N ) ) THEN 00298 INFO = -12 00299 ELSE IF( LDVR.LT.1 .OR. ( ILVR .AND. LDVR.LT.N ) ) THEN 00300 INFO = -14 00301 ELSE IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN 00302 INFO = -16 00303 END IF 00304 * 00305 IF( INFO.EQ.0 ) THEN 00306 NB1 = ILAENV( 1, 'DGEQRF', ' ', N, N, -1, -1 ) 00307 NB2 = ILAENV( 1, 'DORMQR', ' ', N, N, N, -1 ) 00308 NB3 = ILAENV( 1, 'DORGQR', ' ', N, N, N, -1 ) 00309 NB = MAX( NB1, NB2, NB3 ) 00310 LOPT = 2*N + MAX( 6*N, N*( NB+1 ) ) 00311 WORK( 1 ) = LOPT 00312 END IF 00313 * 00314 IF( INFO.NE.0 ) THEN 00315 CALL XERBLA( 'DGEGV ', -INFO ) 00316 RETURN 00317 ELSE IF( LQUERY ) THEN 00318 RETURN 00319 END IF 00320 * 00321 * Quick return if possible 00322 * 00323 IF( N.EQ.0 ) 00324 $ RETURN 00325 * 00326 * Get machine constants 00327 * 00328 EPS = DLAMCH( 'E' )*DLAMCH( 'B' ) 00329 SAFMIN = DLAMCH( 'S' ) 00330 SAFMIN = SAFMIN + SAFMIN 00331 SAFMAX = ONE / SAFMIN 00332 ONEPLS = ONE + ( 4*EPS ) 00333 * 00334 * Scale A 00335 * 00336 ANRM = DLANGE( 'M', N, N, A, LDA, WORK ) 00337 ANRM1 = ANRM 00338 ANRM2 = ONE 00339 IF( ANRM.LT.ONE ) THEN 00340 IF( SAFMAX*ANRM.LT.ONE ) THEN 00341 ANRM1 = SAFMIN 00342 ANRM2 = SAFMAX*ANRM 00343 END IF 00344 END IF 00345 * 00346 IF( ANRM.GT.ZERO ) THEN 00347 CALL DLASCL( 'G', -1, -1, ANRM, ONE, N, N, A, LDA, IINFO ) 00348 IF( IINFO.NE.0 ) THEN 00349 INFO = N + 10 00350 RETURN 00351 END IF 00352 END IF 00353 * 00354 * Scale B 00355 * 00356 BNRM = DLANGE( 'M', N, N, B, LDB, WORK ) 00357 BNRM1 = BNRM 00358 BNRM2 = ONE 00359 IF( BNRM.LT.ONE ) THEN 00360 IF( SAFMAX*BNRM.LT.ONE ) THEN 00361 BNRM1 = SAFMIN 00362 BNRM2 = SAFMAX*BNRM 00363 END IF 00364 END IF 00365 * 00366 IF( BNRM.GT.ZERO ) THEN 00367 CALL DLASCL( 'G', -1, -1, BNRM, ONE, N, N, B, LDB, IINFO ) 00368 IF( IINFO.NE.0 ) THEN 00369 INFO = N + 10 00370 RETURN 00371 END IF 00372 END IF 00373 * 00374 * Permute the matrix to make it more nearly triangular 00375 * Workspace layout: (8*N words -- "work" requires 6*N words) 00376 * left_permutation, right_permutation, work... 00377 * 00378 ILEFT = 1 00379 IRIGHT = N + 1 00380 IWORK = IRIGHT + N 00381 CALL DGGBAL( 'P', N, A, LDA, B, LDB, ILO, IHI, WORK( ILEFT ), 00382 $ WORK( IRIGHT ), WORK( IWORK ), IINFO ) 00383 IF( IINFO.NE.0 ) THEN 00384 INFO = N + 1 00385 GO TO 120 00386 END IF 00387 * 00388 * Reduce B to triangular form, and initialize VL and/or VR 00389 * Workspace layout: ("work..." must have at least N words) 00390 * left_permutation, right_permutation, tau, work... 00391 * 00392 IROWS = IHI + 1 - ILO 00393 IF( ILV ) THEN 00394 ICOLS = N + 1 - ILO 00395 ELSE 00396 ICOLS = IROWS 00397 END IF 00398 ITAU = IWORK 00399 IWORK = ITAU + IROWS 00400 CALL DGEQRF( IROWS, ICOLS, B( ILO, ILO ), LDB, WORK( ITAU ), 00401 $ WORK( IWORK ), LWORK+1-IWORK, IINFO ) 00402 IF( IINFO.GE.0 ) 00403 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 ) 00404 IF( IINFO.NE.0 ) THEN 00405 INFO = N + 2 00406 GO TO 120 00407 END IF 00408 * 00409 CALL DORMQR( 'L', 'T', IROWS, ICOLS, IROWS, B( ILO, ILO ), LDB, 00410 $ WORK( ITAU ), A( ILO, ILO ), LDA, WORK( IWORK ), 00411 $ LWORK+1-IWORK, IINFO ) 00412 IF( IINFO.GE.0 ) 00413 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 ) 00414 IF( IINFO.NE.0 ) THEN 00415 INFO = N + 3 00416 GO TO 120 00417 END IF 00418 * 00419 IF( ILVL ) THEN 00420 CALL DLASET( 'Full', N, N, ZERO, ONE, VL, LDVL ) 00421 CALL DLACPY( 'L', IROWS-1, IROWS-1, B( ILO+1, ILO ), LDB, 00422 $ VL( ILO+1, ILO ), LDVL ) 00423 CALL DORGQR( IROWS, IROWS, IROWS, VL( ILO, ILO ), LDVL, 00424 $ WORK( ITAU ), WORK( IWORK ), LWORK+1-IWORK, 00425 $ IINFO ) 00426 IF( IINFO.GE.0 ) 00427 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 ) 00428 IF( IINFO.NE.0 ) THEN 00429 INFO = N + 4 00430 GO TO 120 00431 END IF 00432 END IF 00433 * 00434 IF( ILVR ) 00435 $ CALL DLASET( 'Full', N, N, ZERO, ONE, VR, LDVR ) 00436 * 00437 * Reduce to generalized Hessenberg form 00438 * 00439 IF( ILV ) THEN 00440 * 00441 * Eigenvectors requested -- work on whole matrix. 00442 * 00443 CALL DGGHRD( JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, VL, 00444 $ LDVL, VR, LDVR, IINFO ) 00445 ELSE 00446 CALL DGGHRD( 'N', 'N', IROWS, 1, IROWS, A( ILO, ILO ), LDA, 00447 $ B( ILO, ILO ), LDB, VL, LDVL, VR, LDVR, IINFO ) 00448 END IF 00449 IF( IINFO.NE.0 ) THEN 00450 INFO = N + 5 00451 GO TO 120 00452 END IF 00453 * 00454 * Perform QZ algorithm 00455 * Workspace layout: ("work..." must have at least 1 word) 00456 * left_permutation, right_permutation, work... 00457 * 00458 IWORK = ITAU 00459 IF( ILV ) THEN 00460 CHTEMP = 'S' 00461 ELSE 00462 CHTEMP = 'E' 00463 END IF 00464 CALL DHGEQZ( CHTEMP, JOBVL, JOBVR, N, ILO, IHI, A, LDA, B, LDB, 00465 $ ALPHAR, ALPHAI, BETA, VL, LDVL, VR, LDVR, 00466 $ WORK( IWORK ), LWORK+1-IWORK, IINFO ) 00467 IF( IINFO.GE.0 ) 00468 $ LWKOPT = MAX( LWKOPT, INT( WORK( IWORK ) )+IWORK-1 ) 00469 IF( IINFO.NE.0 ) THEN 00470 IF( IINFO.GT.0 .AND. IINFO.LE.N ) THEN 00471 INFO = IINFO 00472 ELSE IF( IINFO.GT.N .AND. IINFO.LE.2*N ) THEN 00473 INFO = IINFO - N 00474 ELSE 00475 INFO = N + 6 00476 END IF 00477 GO TO 120 00478 END IF 00479 * 00480 IF( ILV ) THEN 00481 * 00482 * Compute Eigenvectors (DTGEVC requires 6*N words of workspace) 00483 * 00484 IF( ILVL ) THEN 00485 IF( ILVR ) THEN 00486 CHTEMP = 'B' 00487 ELSE 00488 CHTEMP = 'L' 00489 END IF 00490 ELSE 00491 CHTEMP = 'R' 00492 END IF 00493 * 00494 CALL DTGEVC( CHTEMP, 'B', LDUMMA, N, A, LDA, B, LDB, VL, LDVL, 00495 $ VR, LDVR, N, IN, WORK( IWORK ), IINFO ) 00496 IF( IINFO.NE.0 ) THEN 00497 INFO = N + 7 00498 GO TO 120 00499 END IF 00500 * 00501 * Undo balancing on VL and VR, rescale 00502 * 00503 IF( ILVL ) THEN 00504 CALL DGGBAK( 'P', 'L', N, ILO, IHI, WORK( ILEFT ), 00505 $ WORK( IRIGHT ), N, VL, LDVL, IINFO ) 00506 IF( IINFO.NE.0 ) THEN 00507 INFO = N + 8 00508 GO TO 120 00509 END IF 00510 DO 50 JC = 1, N 00511 IF( ALPHAI( JC ).LT.ZERO ) 00512 $ GO TO 50 00513 TEMP = ZERO 00514 IF( ALPHAI( JC ).EQ.ZERO ) THEN 00515 DO 10 JR = 1, N 00516 TEMP = MAX( TEMP, ABS( VL( JR, JC ) ) ) 00517 10 CONTINUE 00518 ELSE 00519 DO 20 JR = 1, N 00520 TEMP = MAX( TEMP, ABS( VL( JR, JC ) )+ 00521 $ ABS( VL( JR, JC+1 ) ) ) 00522 20 CONTINUE 00523 END IF 00524 IF( TEMP.LT.SAFMIN ) 00525 $ GO TO 50 00526 TEMP = ONE / TEMP 00527 IF( ALPHAI( JC ).EQ.ZERO ) THEN 00528 DO 30 JR = 1, N 00529 VL( JR, JC ) = VL( JR, JC )*TEMP 00530 30 CONTINUE 00531 ELSE 00532 DO 40 JR = 1, N 00533 VL( JR, JC ) = VL( JR, JC )*TEMP 00534 VL( JR, JC+1 ) = VL( JR, JC+1 )*TEMP 00535 40 CONTINUE 00536 END IF 00537 50 CONTINUE 00538 END IF 00539 IF( ILVR ) THEN 00540 CALL DGGBAK( 'P', 'R', N, ILO, IHI, WORK( ILEFT ), 00541 $ WORK( IRIGHT ), N, VR, LDVR, IINFO ) 00542 IF( IINFO.NE.0 ) THEN 00543 INFO = N + 9 00544 GO TO 120 00545 END IF 00546 DO 100 JC = 1, N 00547 IF( ALPHAI( JC ).LT.ZERO ) 00548 $ GO TO 100 00549 TEMP = ZERO 00550 IF( ALPHAI( JC ).EQ.ZERO ) THEN 00551 DO 60 JR = 1, N 00552 TEMP = MAX( TEMP, ABS( VR( JR, JC ) ) ) 00553 60 CONTINUE 00554 ELSE 00555 DO 70 JR = 1, N 00556 TEMP = MAX( TEMP, ABS( VR( JR, JC ) )+ 00557 $ ABS( VR( JR, JC+1 ) ) ) 00558 70 CONTINUE 00559 END IF 00560 IF( TEMP.LT.SAFMIN ) 00561 $ GO TO 100 00562 TEMP = ONE / TEMP 00563 IF( ALPHAI( JC ).EQ.ZERO ) THEN 00564 DO 80 JR = 1, N 00565 VR( JR, JC ) = VR( JR, JC )*TEMP 00566 80 CONTINUE 00567 ELSE 00568 DO 90 JR = 1, N 00569 VR( JR, JC ) = VR( JR, JC )*TEMP 00570 VR( JR, JC+1 ) = VR( JR, JC+1 )*TEMP 00571 90 CONTINUE 00572 END IF 00573 100 CONTINUE 00574 END IF 00575 * 00576 * End of eigenvector calculation 00577 * 00578 END IF 00579 * 00580 * Undo scaling in alpha, beta 00581 * 00582 * Note: this does not give the alpha and beta for the unscaled 00583 * problem. 00584 * 00585 * Un-scaling is limited to avoid underflow in alpha and beta 00586 * if they are significant. 00587 * 00588 DO 110 JC = 1, N 00589 ABSAR = ABS( ALPHAR( JC ) ) 00590 ABSAI = ABS( ALPHAI( JC ) ) 00591 ABSB = ABS( BETA( JC ) ) 00592 SALFAR = ANRM*ALPHAR( JC ) 00593 SALFAI = ANRM*ALPHAI( JC ) 00594 SBETA = BNRM*BETA( JC ) 00595 ILIMIT = .FALSE. 00596 SCALE = ONE 00597 * 00598 * Check for significant underflow in ALPHAI 00599 * 00600 IF( ABS( SALFAI ).LT.SAFMIN .AND. ABSAI.GE. 00601 $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSB ) ) THEN 00602 ILIMIT = .TRUE. 00603 SCALE = ( ONEPLS*SAFMIN / ANRM1 ) / 00604 $ MAX( ONEPLS*SAFMIN, ANRM2*ABSAI ) 00605 * 00606 ELSE IF( SALFAI.EQ.ZERO ) THEN 00607 * 00608 * If insignificant underflow in ALPHAI, then make the 00609 * conjugate eigenvalue real. 00610 * 00611 IF( ALPHAI( JC ).LT.ZERO .AND. JC.GT.1 ) THEN 00612 ALPHAI( JC-1 ) = ZERO 00613 ELSE IF( ALPHAI( JC ).GT.ZERO .AND. JC.LT.N ) THEN 00614 ALPHAI( JC+1 ) = ZERO 00615 END IF 00616 END IF 00617 * 00618 * Check for significant underflow in ALPHAR 00619 * 00620 IF( ABS( SALFAR ).LT.SAFMIN .AND. ABSAR.GE. 00621 $ MAX( SAFMIN, EPS*ABSAI, EPS*ABSB ) ) THEN 00622 ILIMIT = .TRUE. 00623 SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / ANRM1 ) / 00624 $ MAX( ONEPLS*SAFMIN, ANRM2*ABSAR ) ) 00625 END IF 00626 * 00627 * Check for significant underflow in BETA 00628 * 00629 IF( ABS( SBETA ).LT.SAFMIN .AND. ABSB.GE. 00630 $ MAX( SAFMIN, EPS*ABSAR, EPS*ABSAI ) ) THEN 00631 ILIMIT = .TRUE. 00632 SCALE = MAX( SCALE, ( ONEPLS*SAFMIN / BNRM1 ) / 00633 $ MAX( ONEPLS*SAFMIN, BNRM2*ABSB ) ) 00634 END IF 00635 * 00636 * Check for possible overflow when limiting scaling 00637 * 00638 IF( ILIMIT ) THEN 00639 TEMP = ( SCALE*SAFMIN )*MAX( ABS( SALFAR ), ABS( SALFAI ), 00640 $ ABS( SBETA ) ) 00641 IF( TEMP.GT.ONE ) 00642 $ SCALE = SCALE / TEMP 00643 IF( SCALE.LT.ONE ) 00644 $ ILIMIT = .FALSE. 00645 END IF 00646 * 00647 * Recompute un-scaled ALPHAR, ALPHAI, BETA if necessary. 00648 * 00649 IF( ILIMIT ) THEN 00650 SALFAR = ( SCALE*ALPHAR( JC ) )*ANRM 00651 SALFAI = ( SCALE*ALPHAI( JC ) )*ANRM 00652 SBETA = ( SCALE*BETA( JC ) )*BNRM 00653 END IF 00654 ALPHAR( JC ) = SALFAR 00655 ALPHAI( JC ) = SALFAI 00656 BETA( JC ) = SBETA 00657 110 CONTINUE 00658 * 00659 120 CONTINUE 00660 WORK( 1 ) = LWKOPT 00661 * 00662 RETURN 00663 * 00664 * End of DGEGV 00665 * 00666 END
1.7.3 
