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