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