00001 DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) 00002 * 00003 * -- LAPACK auxiliary routine (version 3.3.0) -- 00004 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00005 * November 2010 00006 * 00007 * .. Scalar Arguments .. 00008 CHARACTER CMACH 00009 * .. 00010 * 00011 * Purpose 00012 * ======= 00013 * 00014 * DLAMCH determines double precision machine parameters. 00015 * 00016 * Arguments 00017 * ========= 00018 * 00019 * CMACH (input) CHARACTER*1 00020 * Specifies the value to be returned by DLAMCH: 00021 * = 'E' or 'e', DLAMCH := eps 00022 * = 'S' or 's , DLAMCH := sfmin 00023 * = 'B' or 'b', DLAMCH := base 00024 * = 'P' or 'p', DLAMCH := eps*base 00025 * = 'N' or 'n', DLAMCH := t 00026 * = 'R' or 'r', DLAMCH := rnd 00027 * = 'M' or 'm', DLAMCH := emin 00028 * = 'U' or 'u', DLAMCH := rmin 00029 * = 'L' or 'l', DLAMCH := emax 00030 * = 'O' or 'o', DLAMCH := rmax 00031 * 00032 * where 00033 * 00034 * eps = relative machine precision 00035 * sfmin = safe minimum, such that 1/sfmin does not overflow 00036 * base = base of the machine 00037 * prec = eps*base 00038 * t = number of (base) digits in the mantissa 00039 * rnd = 1.0 when rounding occurs in addition, 0.0 otherwise 00040 * emin = minimum exponent before (gradual) underflow 00041 * rmin = underflow threshold - base**(emin-1) 00042 * emax = largest exponent before overflow 00043 * rmax = overflow threshold - (base**emax)*(1-eps) 00044 * 00045 * ===================================================================== 00046 * 00047 * .. Parameters .. 00048 DOUBLE PRECISION ONE, ZERO 00049 PARAMETER ( ONE = 1.0D+0, ZERO = 0.0D+0 ) 00050 * .. 00051 * .. Local Scalars .. 00052 LOGICAL FIRST, LRND 00053 INTEGER BETA, IMAX, IMIN, IT 00054 DOUBLE PRECISION BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN, 00055 $ RND, SFMIN, SMALL, T 00056 * .. 00057 * .. External Functions .. 00058 LOGICAL LSAME 00059 EXTERNAL LSAME 00060 * .. 00061 * .. External Subroutines .. 00062 EXTERNAL DLAMC2 00063 * .. 00064 * .. Save statement .. 00065 SAVE FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN, 00066 $ EMAX, RMAX, PREC 00067 * .. 00068 * .. Data statements .. 00069 DATA FIRST / .TRUE. / 00070 * .. 00071 * .. Executable Statements .. 00072 * 00073 IF( FIRST ) THEN 00074 CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX ) 00075 BASE = BETA 00076 T = IT 00077 IF( LRND ) THEN 00078 RND = ONE 00079 EPS = ( BASE**( 1-IT ) ) / 2 00080 ELSE 00081 RND = ZERO 00082 EPS = BASE**( 1-IT ) 00083 END IF 00084 PREC = EPS*BASE 00085 EMIN = IMIN 00086 EMAX = IMAX 00087 SFMIN = RMIN 00088 SMALL = ONE / RMAX 00089 IF( SMALL.GE.SFMIN ) THEN 00090 * 00091 * Use SMALL plus a bit, to avoid the possibility of rounding 00092 * causing overflow when computing 1/sfmin. 00093 * 00094 SFMIN = SMALL*( ONE+EPS ) 00095 END IF 00096 END IF 00097 * 00098 IF( LSAME( CMACH, 'E' ) ) THEN 00099 RMACH = EPS 00100 ELSE IF( LSAME( CMACH, 'S' ) ) THEN 00101 RMACH = SFMIN 00102 ELSE IF( LSAME( CMACH, 'B' ) ) THEN 00103 RMACH = BASE 00104 ELSE IF( LSAME( CMACH, 'P' ) ) THEN 00105 RMACH = PREC 00106 ELSE IF( LSAME( CMACH, 'N' ) ) THEN 00107 RMACH = T 00108 ELSE IF( LSAME( CMACH, 'R' ) ) THEN 00109 RMACH = RND 00110 ELSE IF( LSAME( CMACH, 'M' ) ) THEN 00111 RMACH = EMIN 00112 ELSE IF( LSAME( CMACH, 'U' ) ) THEN 00113 RMACH = RMIN 00114 ELSE IF( LSAME( CMACH, 'L' ) ) THEN 00115 RMACH = EMAX 00116 ELSE IF( LSAME( CMACH, 'O' ) ) THEN 00117 RMACH = RMAX 00118 END IF 00119 * 00120 DLAMCH = RMACH 00121 FIRST = .FALSE. 00122 RETURN 00123 * 00124 * End of DLAMCH 00125 * 00126 END 00127 * 00128 ************************************************************************ 00129 * 00130 SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) 00131 * 00132 * -- LAPACK auxiliary routine (version 3.3.0) -- 00133 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00134 * November 2010 00135 * 00136 * .. Scalar Arguments .. 00137 LOGICAL IEEE1, RND 00138 INTEGER BETA, T 00139 * .. 00140 * 00141 * Purpose 00142 * ======= 00143 * 00144 * DLAMC1 determines the machine parameters given by BETA, T, RND, and 00145 * IEEE1. 00146 * 00147 * Arguments 00148 * ========= 00149 * 00150 * BETA (output) INTEGER 00151 * The base of the machine. 00152 * 00153 * T (output) INTEGER 00154 * The number of ( BETA ) digits in the mantissa. 00155 * 00156 * RND (output) LOGICAL 00157 * Specifies whether proper rounding ( RND = .TRUE. ) or 00158 * chopping ( RND = .FALSE. ) occurs in addition. This may not 00159 * be a reliable guide to the way in which the machine performs 00160 * its arithmetic. 00161 * 00162 * IEEE1 (output) LOGICAL 00163 * Specifies whether rounding appears to be done in the IEEE 00164 * 'round to nearest' style. 00165 * 00166 * Further Details 00167 * =============== 00168 * 00169 * The routine is based on the routine ENVRON by Malcolm and 00170 * incorporates suggestions by Gentleman and Marovich. See 00171 * 00172 * Malcolm M. A. (1972) Algorithms to reveal properties of 00173 * floating-point arithmetic. Comms. of the ACM, 15, 949-951. 00174 * 00175 * Gentleman W. M. and Marovich S. B. (1974) More on algorithms 00176 * that reveal properties of floating point arithmetic units. 00177 * Comms. of the ACM, 17, 276-277. 00178 * 00179 * ===================================================================== 00180 * 00181 * .. Local Scalars .. 00182 LOGICAL FIRST, LIEEE1, LRND 00183 INTEGER LBETA, LT 00184 DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 00185 * .. 00186 * .. External Functions .. 00187 DOUBLE PRECISION DLAMC3 00188 EXTERNAL DLAMC3 00189 * .. 00190 * .. Save statement .. 00191 SAVE FIRST, LIEEE1, LBETA, LRND, LT 00192 * .. 00193 * .. Data statements .. 00194 DATA FIRST / .TRUE. / 00195 * .. 00196 * .. Executable Statements .. 00197 * 00198 IF( FIRST ) THEN 00199 ONE = 1 00200 * 00201 * LBETA, LIEEE1, LT and LRND are the local values of BETA, 00202 * IEEE1, T and RND. 00203 * 00204 * Throughout this routine we use the function DLAMC3 to ensure 00205 * that relevant values are stored and not held in registers, or 00206 * are not affected by optimizers. 00207 * 00208 * Compute a = 2.0**m with the smallest positive integer m such 00209 * that 00210 * 00211 * fl( a + 1.0 ) = a. 00212 * 00213 A = 1 00214 C = 1 00215 * 00216 *+ WHILE( C.EQ.ONE )LOOP 00217 10 CONTINUE 00218 IF( C.EQ.ONE ) THEN 00219 A = 2*A 00220 C = DLAMC3( A, ONE ) 00221 C = DLAMC3( C, -A ) 00222 GO TO 10 00223 END IF 00224 *+ END WHILE 00225 * 00226 * Now compute b = 2.0**m with the smallest positive integer m 00227 * such that 00228 * 00229 * fl( a + b ) .gt. a. 00230 * 00231 B = 1 00232 C = DLAMC3( A, B ) 00233 * 00234 *+ WHILE( C.EQ.A )LOOP 00235 20 CONTINUE 00236 IF( C.EQ.A ) THEN 00237 B = 2*B 00238 C = DLAMC3( A, B ) 00239 GO TO 20 00240 END IF 00241 *+ END WHILE 00242 * 00243 * Now compute the base. a and c are neighbouring floating point 00244 * numbers in the interval ( beta**t, beta**( t + 1 ) ) and so 00245 * their difference is beta. Adding 0.25 to c is to ensure that it 00246 * is truncated to beta and not ( beta - 1 ). 00247 * 00248 QTR = ONE / 4 00249 SAVEC = C 00250 C = DLAMC3( C, -A ) 00251 LBETA = C + QTR 00252 * 00253 * Now determine whether rounding or chopping occurs, by adding a 00254 * bit less than beta/2 and a bit more than beta/2 to a. 00255 * 00256 B = LBETA 00257 F = DLAMC3( B / 2, -B / 100 ) 00258 C = DLAMC3( F, A ) 00259 IF( C.EQ.A ) THEN 00260 LRND = .TRUE. 00261 ELSE 00262 LRND = .FALSE. 00263 END IF 00264 F = DLAMC3( B / 2, B / 100 ) 00265 C = DLAMC3( F, A ) 00266 IF( ( LRND ) .AND. ( C.EQ.A ) ) 00267 $ LRND = .FALSE. 00268 * 00269 * Try and decide whether rounding is done in the IEEE 'round to 00270 * nearest' style. B/2 is half a unit in the last place of the two 00271 * numbers A and SAVEC. Furthermore, A is even, i.e. has last bit 00272 * zero, and SAVEC is odd. Thus adding B/2 to A should not change 00273 * A, but adding B/2 to SAVEC should change SAVEC. 00274 * 00275 T1 = DLAMC3( B / 2, A ) 00276 T2 = DLAMC3( B / 2, SAVEC ) 00277 LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND 00278 * 00279 * Now find the mantissa, t. It should be the integer part of 00280 * log to the base beta of a, however it is safer to determine t 00281 * by powering. So we find t as the smallest positive integer for 00282 * which 00283 * 00284 * fl( beta**t + 1.0 ) = 1.0. 00285 * 00286 LT = 0 00287 A = 1 00288 C = 1 00289 * 00290 *+ WHILE( C.EQ.ONE )LOOP 00291 30 CONTINUE 00292 IF( C.EQ.ONE ) THEN 00293 LT = LT + 1 00294 A = A*LBETA 00295 C = DLAMC3( A, ONE ) 00296 C = DLAMC3( C, -A ) 00297 GO TO 30 00298 END IF 00299 *+ END WHILE 00300 * 00301 END IF 00302 * 00303 BETA = LBETA 00304 T = LT 00305 RND = LRND 00306 IEEE1 = LIEEE1 00307 FIRST = .FALSE. 00308 RETURN 00309 * 00310 * End of DLAMC1 00311 * 00312 END 00313 * 00314 ************************************************************************ 00315 * 00316 SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX ) 00317 * 00318 * -- LAPACK auxiliary routine (version 3.3.0) -- 00319 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00320 * November 2010 00321 * 00322 * .. Scalar Arguments .. 00323 LOGICAL RND 00324 INTEGER BETA, EMAX, EMIN, T 00325 DOUBLE PRECISION EPS, RMAX, RMIN 00326 * .. 00327 * 00328 * Purpose 00329 * ======= 00330 * 00331 * DLAMC2 determines the machine parameters specified in its argument 00332 * list. 00333 * 00334 * Arguments 00335 * ========= 00336 * 00337 * BETA (output) INTEGER 00338 * The base of the machine. 00339 * 00340 * T (output) INTEGER 00341 * The number of ( BETA ) digits in the mantissa. 00342 * 00343 * RND (output) LOGICAL 00344 * Specifies whether proper rounding ( RND = .TRUE. ) or 00345 * chopping ( RND = .FALSE. ) occurs in addition. This may not 00346 * be a reliable guide to the way in which the machine performs 00347 * its arithmetic. 00348 * 00349 * EPS (output) DOUBLE PRECISION 00350 * The smallest positive number such that 00351 * 00352 * fl( 1.0 - EPS ) .LT. 1.0, 00353 * 00354 * where fl denotes the computed value. 00355 * 00356 * EMIN (output) INTEGER 00357 * The minimum exponent before (gradual) underflow occurs. 00358 * 00359 * RMIN (output) DOUBLE PRECISION 00360 * The smallest normalized number for the machine, given by 00361 * BASE**( EMIN - 1 ), where BASE is the floating point value 00362 * of BETA. 00363 * 00364 * EMAX (output) INTEGER 00365 * The maximum exponent before overflow occurs. 00366 * 00367 * RMAX (output) DOUBLE PRECISION 00368 * The largest positive number for the machine, given by 00369 * BASE**EMAX * ( 1 - EPS ), where BASE is the floating point 00370 * value of BETA. 00371 * 00372 * Further Details 00373 * =============== 00374 * 00375 * The computation of EPS is based on a routine PARANOIA by 00376 * W. Kahan of the University of California at Berkeley. 00377 * 00378 * ===================================================================== 00379 * 00380 * .. Local Scalars .. 00381 LOGICAL FIRST, IEEE, IWARN, LIEEE1, LRND 00382 INTEGER GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT, 00383 $ NGNMIN, NGPMIN 00384 DOUBLE PRECISION A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, 00385 $ SIXTH, SMALL, THIRD, TWO, ZERO 00386 * .. 00387 * .. External Functions .. 00388 DOUBLE PRECISION DLAMC3 00389 EXTERNAL DLAMC3 00390 * .. 00391 * .. External Subroutines .. 00392 EXTERNAL DLAMC1, DLAMC4, DLAMC5 00393 * .. 00394 * .. Intrinsic Functions .. 00395 INTRINSIC ABS, MAX, MIN 00396 * .. 00397 * .. Save statement .. 00398 SAVE FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX, 00399 $ LRMIN, LT 00400 * .. 00401 * .. Data statements .. 00402 DATA FIRST / .TRUE. / , IWARN / .FALSE. / 00403 * .. 00404 * .. Executable Statements .. 00405 * 00406 IF( FIRST ) THEN 00407 ZERO = 0 00408 ONE = 1 00409 TWO = 2 00410 * 00411 * LBETA, LT, LRND, LEPS, LEMIN and LRMIN are the local values of 00412 * BETA, T, RND, EPS, EMIN and RMIN. 00413 * 00414 * Throughout this routine we use the function DLAMC3 to ensure 00415 * that relevant values are stored and not held in registers, or 00416 * are not affected by optimizers. 00417 * 00418 * DLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. 00419 * 00420 CALL DLAMC1( LBETA, LT, LRND, LIEEE1 ) 00421 * 00422 * Start to find EPS. 00423 * 00424 B = LBETA 00425 A = B**( -LT ) 00426 LEPS = A 00427 * 00428 * Try some tricks to see whether or not this is the correct EPS. 00429 * 00430 B = TWO / 3 00431 HALF = ONE / 2 00432 SIXTH = DLAMC3( B, -HALF ) 00433 THIRD = DLAMC3( SIXTH, SIXTH ) 00434 B = DLAMC3( THIRD, -HALF ) 00435 B = DLAMC3( B, SIXTH ) 00436 B = ABS( B ) 00437 IF( B.LT.LEPS ) 00438 $ B = LEPS 00439 * 00440 LEPS = 1 00441 * 00442 *+ WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP 00443 10 CONTINUE 00444 IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN 00445 LEPS = B 00446 C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) 00447 C = DLAMC3( HALF, -C ) 00448 B = DLAMC3( HALF, C ) 00449 C = DLAMC3( HALF, -B ) 00450 B = DLAMC3( HALF, C ) 00451 GO TO 10 00452 END IF 00453 *+ END WHILE 00454 * 00455 IF( A.LT.LEPS ) 00456 $ LEPS = A 00457 * 00458 * Computation of EPS complete. 00459 * 00460 * Now find EMIN. Let A = + or - 1, and + or - (1 + BASE**(-3)). 00461 * Keep dividing A by BETA until (gradual) underflow occurs. This 00462 * is detected when we cannot recover the previous A. 00463 * 00464 RBASE = ONE / LBETA 00465 SMALL = ONE 00466 DO 20 I = 1, 3 00467 SMALL = DLAMC3( SMALL*RBASE, ZERO ) 00468 20 CONTINUE 00469 A = DLAMC3( ONE, SMALL ) 00470 CALL DLAMC4( NGPMIN, ONE, LBETA ) 00471 CALL DLAMC4( NGNMIN, -ONE, LBETA ) 00472 CALL DLAMC4( GPMIN, A, LBETA ) 00473 CALL DLAMC4( GNMIN, -A, LBETA ) 00474 IEEE = .FALSE. 00475 * 00476 IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN 00477 IF( NGPMIN.EQ.GPMIN ) THEN 00478 LEMIN = NGPMIN 00479 * ( Non twos-complement machines, no gradual underflow; 00480 * e.g., VAX ) 00481 ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN 00482 LEMIN = NGPMIN - 1 + LT 00483 IEEE = .TRUE. 00484 * ( Non twos-complement machines, with gradual underflow; 00485 * e.g., IEEE standard followers ) 00486 ELSE 00487 LEMIN = MIN( NGPMIN, GPMIN ) 00488 * ( A guess; no known machine ) 00489 IWARN = .TRUE. 00490 END IF 00491 * 00492 ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN 00493 IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN 00494 LEMIN = MAX( NGPMIN, NGNMIN ) 00495 * ( Twos-complement machines, no gradual underflow; 00496 * e.g., CYBER 205 ) 00497 ELSE 00498 LEMIN = MIN( NGPMIN, NGNMIN ) 00499 * ( A guess; no known machine ) 00500 IWARN = .TRUE. 00501 END IF 00502 * 00503 ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND. 00504 $ ( GPMIN.EQ.GNMIN ) ) THEN 00505 IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN 00506 LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT 00507 * ( Twos-complement machines with gradual underflow; 00508 * no known machine ) 00509 ELSE 00510 LEMIN = MIN( NGPMIN, NGNMIN ) 00511 * ( A guess; no known machine ) 00512 IWARN = .TRUE. 00513 END IF 00514 * 00515 ELSE 00516 LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN ) 00517 * ( A guess; no known machine ) 00518 IWARN = .TRUE. 00519 END IF 00520 FIRST = .FALSE. 00521 *** 00522 * Comment out this if block if EMIN is ok 00523 IF( IWARN ) THEN 00524 FIRST = .TRUE. 00525 WRITE( 6, FMT = 9999 )LEMIN 00526 END IF 00527 *** 00528 * 00529 * Assume IEEE arithmetic if we found denormalised numbers above, 00530 * or if arithmetic seems to round in the IEEE style, determined 00531 * in routine DLAMC1. A true IEEE machine should have both things 00532 * true; however, faulty machines may have one or the other. 00533 * 00534 IEEE = IEEE .OR. LIEEE1 00535 * 00536 * Compute RMIN by successive division by BETA. We could compute 00537 * RMIN as BASE**( EMIN - 1 ), but some machines underflow during 00538 * this computation. 00539 * 00540 LRMIN = 1 00541 DO 30 I = 1, 1 - LEMIN 00542 LRMIN = DLAMC3( LRMIN*RBASE, ZERO ) 00543 30 CONTINUE 00544 * 00545 * Finally, call DLAMC5 to compute EMAX and RMAX. 00546 * 00547 CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX ) 00548 END IF 00549 * 00550 BETA = LBETA 00551 T = LT 00552 RND = LRND 00553 EPS = LEPS 00554 EMIN = LEMIN 00555 RMIN = LRMIN 00556 EMAX = LEMAX 00557 RMAX = LRMAX 00558 * 00559 RETURN 00560 * 00561 9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-', 00562 $ ' EMIN = ', I8, / 00563 $ ' If, after inspection, the value EMIN looks', 00564 $ ' acceptable please comment out ', 00565 $ / ' the IF block as marked within the code of routine', 00566 $ ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / ) 00567 * 00568 * End of DLAMC2 00569 * 00570 END 00571 * 00572 ************************************************************************ 00573 * 00574 DOUBLE PRECISION FUNCTION DLAMC3( A, B ) 00575 * 00576 * -- LAPACK auxiliary routine (version 3.3.0) -- 00577 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00578 * November 2010 00579 * 00580 * .. Scalar Arguments .. 00581 DOUBLE PRECISION A, B 00582 * .. 00583 * 00584 * Purpose 00585 * ======= 00586 * 00587 * DLAMC3 is intended to force A and B to be stored prior to doing 00588 * the addition of A and B , for use in situations where optimizers 00589 * might hold one of these in a register. 00590 * 00591 * Arguments 00592 * ========= 00593 * 00594 * A (input) DOUBLE PRECISION 00595 * B (input) DOUBLE PRECISION 00596 * The values A and B. 00597 * 00598 * ===================================================================== 00599 * 00600 * .. Executable Statements .. 00601 * 00602 DLAMC3 = A + B 00603 * 00604 RETURN 00605 * 00606 * End of DLAMC3 00607 * 00608 END 00609 * 00610 ************************************************************************ 00611 * 00612 SUBROUTINE DLAMC4( EMIN, START, BASE ) 00613 * 00614 * -- LAPACK auxiliary routine (version 3.3.0) -- 00615 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00616 * November 2010 00617 * 00618 * .. Scalar Arguments .. 00619 INTEGER BASE, EMIN 00620 DOUBLE PRECISION START 00621 * .. 00622 * 00623 * Purpose 00624 * ======= 00625 * 00626 * DLAMC4 is a service routine for DLAMC2. 00627 * 00628 * Arguments 00629 * ========= 00630 * 00631 * EMIN (output) INTEGER 00632 * The minimum exponent before (gradual) underflow, computed by 00633 * setting A = START and dividing by BASE until the previous A 00634 * can not be recovered. 00635 * 00636 * START (input) DOUBLE PRECISION 00637 * The starting point for determining EMIN. 00638 * 00639 * BASE (input) INTEGER 00640 * The base of the machine. 00641 * 00642 * ===================================================================== 00643 * 00644 * .. Local Scalars .. 00645 INTEGER I 00646 DOUBLE PRECISION A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO 00647 * .. 00648 * .. External Functions .. 00649 DOUBLE PRECISION DLAMC3 00650 EXTERNAL DLAMC3 00651 * .. 00652 * .. Executable Statements .. 00653 * 00654 A = START 00655 ONE = 1 00656 RBASE = ONE / BASE 00657 ZERO = 0 00658 EMIN = 1 00659 B1 = DLAMC3( A*RBASE, ZERO ) 00660 C1 = A 00661 C2 = A 00662 D1 = A 00663 D2 = A 00664 *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. 00665 * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP 00666 10 CONTINUE 00667 IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. 00668 $ ( D2.EQ.A ) ) THEN 00669 EMIN = EMIN - 1 00670 A = B1 00671 B1 = DLAMC3( A / BASE, ZERO ) 00672 C1 = DLAMC3( B1*BASE, ZERO ) 00673 D1 = ZERO 00674 DO 20 I = 1, BASE 00675 D1 = D1 + B1 00676 20 CONTINUE 00677 B2 = DLAMC3( A*RBASE, ZERO ) 00678 C2 = DLAMC3( B2 / RBASE, ZERO ) 00679 D2 = ZERO 00680 DO 30 I = 1, BASE 00681 D2 = D2 + B2 00682 30 CONTINUE 00683 GO TO 10 00684 END IF 00685 *+ END WHILE 00686 * 00687 RETURN 00688 * 00689 * End of DLAMC4 00690 * 00691 END 00692 * 00693 ************************************************************************ 00694 * 00695 SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) 00696 * 00697 * -- LAPACK auxiliary routine (version 3.3.0) -- 00698 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00699 * November 2010 00700 * 00701 * .. Scalar Arguments .. 00702 LOGICAL IEEE 00703 INTEGER BETA, EMAX, EMIN, P 00704 DOUBLE PRECISION RMAX 00705 * .. 00706 * 00707 * Purpose 00708 * ======= 00709 * 00710 * DLAMC5 attempts to compute RMAX, the largest machine floating-point 00711 * number, without overflow. It assumes that EMAX + abs(EMIN) sum 00712 * approximately to a power of 2. It will fail on machines where this 00713 * assumption does not hold, for example, the Cyber 205 (EMIN = -28625, 00714 * EMAX = 28718). It will also fail if the value supplied for EMIN is 00715 * too large (i.e. too close to zero), probably with overflow. 00716 * 00717 * Arguments 00718 * ========= 00719 * 00720 * BETA (input) INTEGER 00721 * The base of floating-point arithmetic. 00722 * 00723 * P (input) INTEGER 00724 * The number of base BETA digits in the mantissa of a 00725 * floating-point value. 00726 * 00727 * EMIN (input) INTEGER 00728 * The minimum exponent before (gradual) underflow. 00729 * 00730 * IEEE (input) LOGICAL 00731 * A logical flag specifying whether or not the arithmetic 00732 * system is thought to comply with the IEEE standard. 00733 * 00734 * EMAX (output) INTEGER 00735 * The largest exponent before overflow 00736 * 00737 * RMAX (output) DOUBLE PRECISION 00738 * The largest machine floating-point number. 00739 * 00740 * ===================================================================== 00741 * 00742 * .. Parameters .. 00743 DOUBLE PRECISION ZERO, ONE 00744 PARAMETER ( ZERO = 0.0D0, ONE = 1.0D0 ) 00745 * .. 00746 * .. Local Scalars .. 00747 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP 00748 DOUBLE PRECISION OLDY, RECBAS, Y, Z 00749 * .. 00750 * .. External Functions .. 00751 DOUBLE PRECISION DLAMC3 00752 EXTERNAL DLAMC3 00753 * .. 00754 * .. Intrinsic Functions .. 00755 INTRINSIC MOD 00756 * .. 00757 * .. Executable Statements .. 00758 * 00759 * First compute LEXP and UEXP, two powers of 2 that bound 00760 * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum 00761 * approximately to the bound that is closest to abs(EMIN). 00762 * (EMAX is the exponent of the required number RMAX). 00763 * 00764 LEXP = 1 00765 EXBITS = 1 00766 10 CONTINUE 00767 TRY = LEXP*2 00768 IF( TRY.LE.( -EMIN ) ) THEN 00769 LEXP = TRY 00770 EXBITS = EXBITS + 1 00771 GO TO 10 00772 END IF 00773 IF( LEXP.EQ.-EMIN ) THEN 00774 UEXP = LEXP 00775 ELSE 00776 UEXP = TRY 00777 EXBITS = EXBITS + 1 00778 END IF 00779 * 00780 * Now -LEXP is less than or equal to EMIN, and -UEXP is greater 00781 * than or equal to EMIN. EXBITS is the number of bits needed to 00782 * store the exponent. 00783 * 00784 IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN 00785 EXPSUM = 2*LEXP 00786 ELSE 00787 EXPSUM = 2*UEXP 00788 END IF 00789 * 00790 * EXPSUM is the exponent range, approximately equal to 00791 * EMAX - EMIN + 1 . 00792 * 00793 EMAX = EXPSUM + EMIN - 1 00794 NBITS = 1 + EXBITS + P 00795 * 00796 * NBITS is the total number of bits needed to store a 00797 * floating-point number. 00798 * 00799 IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN 00800 * 00801 * Either there are an odd number of bits used to store a 00802 * floating-point number, which is unlikely, or some bits are 00803 * not used in the representation of numbers, which is possible, 00804 * (e.g. Cray machines) or the mantissa has an implicit bit, 00805 * (e.g. IEEE machines, Dec Vax machines), which is perhaps the 00806 * most likely. We have to assume the last alternative. 00807 * If this is true, then we need to reduce EMAX by one because 00808 * there must be some way of representing zero in an implicit-bit 00809 * system. On machines like Cray, we are reducing EMAX by one 00810 * unnecessarily. 00811 * 00812 EMAX = EMAX - 1 00813 END IF 00814 * 00815 IF( IEEE ) THEN 00816 * 00817 * Assume we are on an IEEE machine which reserves one exponent 00818 * for infinity and NaN. 00819 * 00820 EMAX = EMAX - 1 00821 END IF 00822 * 00823 * Now create RMAX, the largest machine number, which should 00824 * be equal to (1.0 - BETA**(-P)) * BETA**EMAX . 00825 * 00826 * First compute 1.0 - BETA**(-P), being careful that the 00827 * result is less than 1.0 . 00828 * 00829 RECBAS = ONE / BETA 00830 Z = BETA - ONE 00831 Y = ZERO 00832 DO 20 I = 1, P 00833 Z = Z*RECBAS 00834 IF( Y.LT.ONE ) 00835 $ OLDY = Y 00836 Y = DLAMC3( Y, Z ) 00837 20 CONTINUE 00838 IF( Y.GE.ONE ) 00839 $ Y = OLDY 00840 * 00841 * Now multiply by BETA**EMAX to get RMAX. 00842 * 00843 DO 30 I = 1, EMAX 00844 Y = DLAMC3( Y*BETA, ZERO ) 00845 30 CONTINUE 00846 * 00847 RMAX = Y 00848 RETURN 00849 * 00850 * End of DLAMC5 00851 * 00852 END
1.7.2