|
LAPACK 3.3.0
|
|
LAPACK 3.3.0
|
00001 REAL FUNCTION SLAMCH( 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 * SLAMCH determines single precision machine parameters. 00015 * 00016 * Arguments 00017 * ========= 00018 * 00019 * CMACH (input) CHARACTER*1 00020 * Specifies the value to be returned by SLAMCH: 00021 * = 'E' or 'e', SLAMCH := eps 00022 * = 'S' or 's , SLAMCH := sfmin 00023 * = 'B' or 'b', SLAMCH := base 00024 * = 'P' or 'p', SLAMCH := eps*base 00025 * = 'N' or 'n', SLAMCH := t 00026 * = 'R' or 'r', SLAMCH := rnd 00027 * = 'M' or 'm', SLAMCH := emin 00028 * = 'U' or 'u', SLAMCH := rmin 00029 * = 'L' or 'l', SLAMCH := emax 00030 * = 'O' or 'o', SLAMCH := 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 REAL ONE, ZERO 00049 PARAMETER ( ONE = 1.0E+0, ZERO = 0.0E+0 ) 00050 * .. 00051 * .. Local Scalars .. 00052 LOGICAL FIRST, LRND 00053 INTEGER BETA, IMAX, IMIN, IT 00054 REAL 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 SLAMC2 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 SLAMC2( 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 SLAMCH = RMACH 00121 FIRST = .FALSE. 00122 RETURN 00123 * 00124 * End of SLAMCH 00125 * 00126 END 00127 * 00128 ************************************************************************ 00129 * 00130 SUBROUTINE SLAMC1( 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 * SLAMC1 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 REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2 00185 * .. 00186 * .. External Functions .. 00187 REAL SLAMC3 00188 EXTERNAL SLAMC3 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 SLAMC3 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 = SLAMC3( A, ONE ) 00221 C = SLAMC3( 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 = SLAMC3( A, B ) 00233 * 00234 *+ WHILE( C.EQ.A )LOOP 00235 20 CONTINUE 00236 IF( C.EQ.A ) THEN 00237 B = 2*B 00238 C = SLAMC3( 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 = SLAMC3( 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 = SLAMC3( B / 2, -B / 100 ) 00258 C = SLAMC3( F, A ) 00259 IF( C.EQ.A ) THEN 00260 LRND = .TRUE. 00261 ELSE 00262 LRND = .FALSE. 00263 END IF 00264 F = SLAMC3( B / 2, B / 100 ) 00265 C = SLAMC3( 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 = SLAMC3( B / 2, A ) 00276 T2 = SLAMC3( 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 = SLAMC3( A, ONE ) 00296 C = SLAMC3( 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 SLAMC1 00311 * 00312 END 00313 * 00314 ************************************************************************ 00315 * 00316 SUBROUTINE SLAMC2( 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 REAL EPS, RMAX, RMIN 00326 * .. 00327 * 00328 * Purpose 00329 * ======= 00330 * 00331 * SLAMC2 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) REAL 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) REAL 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) REAL 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 REAL A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE, 00385 $ SIXTH, SMALL, THIRD, TWO, ZERO 00386 * .. 00387 * .. External Functions .. 00388 REAL SLAMC3 00389 EXTERNAL SLAMC3 00390 * .. 00391 * .. External Subroutines .. 00392 EXTERNAL SLAMC1, SLAMC4, SLAMC5 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 SLAMC3 to ensure 00415 * that relevant values are stored and not held in registers, or 00416 * are not affected by optimizers. 00417 * 00418 * SLAMC1 returns the parameters LBETA, LT, LRND and LIEEE1. 00419 * 00420 CALL SLAMC1( 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 = SLAMC3( B, -HALF ) 00433 THIRD = SLAMC3( SIXTH, SIXTH ) 00434 B = SLAMC3( THIRD, -HALF ) 00435 B = SLAMC3( 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 = SLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) ) 00447 C = SLAMC3( HALF, -C ) 00448 B = SLAMC3( HALF, C ) 00449 C = SLAMC3( HALF, -B ) 00450 B = SLAMC3( 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 = SLAMC3( SMALL*RBASE, ZERO ) 00468 20 CONTINUE 00469 A = SLAMC3( ONE, SMALL ) 00470 CALL SLAMC4( NGPMIN, ONE, LBETA ) 00471 CALL SLAMC4( NGNMIN, -ONE, LBETA ) 00472 CALL SLAMC4( GPMIN, A, LBETA ) 00473 CALL SLAMC4( 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 SLAMC1. 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 = SLAMC3( LRMIN*RBASE, ZERO ) 00543 30 CONTINUE 00544 * 00545 * Finally, call SLAMC5 to compute EMAX and RMAX. 00546 * 00547 CALL SLAMC5( 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 $ ' SLAMC2,', / ' otherwise supply EMIN explicitly.', / ) 00567 * 00568 * End of SLAMC2 00569 * 00570 END 00571 * 00572 ************************************************************************ 00573 * 00574 REAL FUNCTION SLAMC3( 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 REAL A, B 00582 * .. 00583 * 00584 * Purpose 00585 * ======= 00586 * 00587 * SLAMC3 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) REAL 00595 * B (input) REAL 00596 * The values A and B. 00597 * 00598 * ===================================================================== 00599 * 00600 * .. Executable Statements .. 00601 * 00602 SLAMC3 = A + B 00603 * 00604 RETURN 00605 * 00606 * End of SLAMC3 00607 * 00608 END 00609 * 00610 ************************************************************************ 00611 * 00612 SUBROUTINE SLAMC4( 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 00620 INTEGER EMIN 00621 REAL START 00622 * .. 00623 * 00624 * Purpose 00625 * ======= 00626 * 00627 * SLAMC4 is a service routine for SLAMC2. 00628 * 00629 * Arguments 00630 * ========= 00631 * 00632 * EMIN (output) INTEGER 00633 * The minimum exponent before (gradual) underflow, computed by 00634 * setting A = START and dividing by BASE until the previous A 00635 * can not be recovered. 00636 * 00637 * START (input) REAL 00638 * The starting point for determining EMIN. 00639 * 00640 * BASE (input) INTEGER 00641 * The base of the machine. 00642 * 00643 * ===================================================================== 00644 * 00645 * .. Local Scalars .. 00646 INTEGER I 00647 REAL A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO 00648 * .. 00649 * .. External Functions .. 00650 REAL SLAMC3 00651 EXTERNAL SLAMC3 00652 * .. 00653 * .. Executable Statements .. 00654 * 00655 A = START 00656 ONE = 1 00657 RBASE = ONE / BASE 00658 ZERO = 0 00659 EMIN = 1 00660 B1 = SLAMC3( A*RBASE, ZERO ) 00661 C1 = A 00662 C2 = A 00663 D1 = A 00664 D2 = A 00665 *+ WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND. 00666 * $ ( D1.EQ.A ).AND.( D2.EQ.A ) )LOOP 00667 10 CONTINUE 00668 IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND. 00669 $ ( D2.EQ.A ) ) THEN 00670 EMIN = EMIN - 1 00671 A = B1 00672 B1 = SLAMC3( A / BASE, ZERO ) 00673 C1 = SLAMC3( B1*BASE, ZERO ) 00674 D1 = ZERO 00675 DO 20 I = 1, BASE 00676 D1 = D1 + B1 00677 20 CONTINUE 00678 B2 = SLAMC3( A*RBASE, ZERO ) 00679 C2 = SLAMC3( B2 / RBASE, ZERO ) 00680 D2 = ZERO 00681 DO 30 I = 1, BASE 00682 D2 = D2 + B2 00683 30 CONTINUE 00684 GO TO 10 00685 END IF 00686 *+ END WHILE 00687 * 00688 RETURN 00689 * 00690 * End of SLAMC4 00691 * 00692 END 00693 * 00694 ************************************************************************ 00695 * 00696 SUBROUTINE SLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX ) 00697 * 00698 * -- LAPACK auxiliary routine (version 3.3.0) -- 00699 * Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. 00700 * November 2010 00701 * 00702 * .. Scalar Arguments .. 00703 LOGICAL IEEE 00704 INTEGER BETA, EMAX, EMIN, P 00705 REAL RMAX 00706 * .. 00707 * 00708 * Purpose 00709 * ======= 00710 * 00711 * SLAMC5 attempts to compute RMAX, the largest machine floating-point 00712 * number, without overflow. It assumes that EMAX + abs(EMIN) sum 00713 * approximately to a power of 2. It will fail on machines where this 00714 * assumption does not hold, for example, the Cyber 205 (EMIN = -28625, 00715 * EMAX = 28718). It will also fail if the value supplied for EMIN is 00716 * too large (i.e. too close to zero), probably with overflow. 00717 * 00718 * Arguments 00719 * ========= 00720 * 00721 * BETA (input) INTEGER 00722 * The base of floating-point arithmetic. 00723 * 00724 * P (input) INTEGER 00725 * The number of base BETA digits in the mantissa of a 00726 * floating-point value. 00727 * 00728 * EMIN (input) INTEGER 00729 * The minimum exponent before (gradual) underflow. 00730 * 00731 * IEEE (input) LOGICAL 00732 * A logical flag specifying whether or not the arithmetic 00733 * system is thought to comply with the IEEE standard. 00734 * 00735 * EMAX (output) INTEGER 00736 * The largest exponent before overflow 00737 * 00738 * RMAX (output) REAL 00739 * The largest machine floating-point number. 00740 * 00741 * ===================================================================== 00742 * 00743 * .. Parameters .. 00744 REAL ZERO, ONE 00745 PARAMETER ( ZERO = 0.0E0, ONE = 1.0E0 ) 00746 * .. 00747 * .. Local Scalars .. 00748 INTEGER EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP 00749 REAL OLDY, RECBAS, Y, Z 00750 * .. 00751 * .. External Functions .. 00752 REAL SLAMC3 00753 EXTERNAL SLAMC3 00754 * .. 00755 * .. Intrinsic Functions .. 00756 INTRINSIC MOD 00757 * .. 00758 * .. Executable Statements .. 00759 * 00760 * First compute LEXP and UEXP, two powers of 2 that bound 00761 * abs(EMIN). We then assume that EMAX + abs(EMIN) will sum 00762 * approximately to the bound that is closest to abs(EMIN). 00763 * (EMAX is the exponent of the required number RMAX). 00764 * 00765 LEXP = 1 00766 EXBITS = 1 00767 10 CONTINUE 00768 TRY = LEXP*2 00769 IF( TRY.LE.( -EMIN ) ) THEN 00770 LEXP = TRY 00771 EXBITS = EXBITS + 1 00772 GO TO 10 00773 END IF 00774 IF( LEXP.EQ.-EMIN ) THEN 00775 UEXP = LEXP 00776 ELSE 00777 UEXP = TRY 00778 EXBITS = EXBITS + 1 00779 END IF 00780 * 00781 * Now -LEXP is less than or equal to EMIN, and -UEXP is greater 00782 * than or equal to EMIN. EXBITS is the number of bits needed to 00783 * store the exponent. 00784 * 00785 IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN 00786 EXPSUM = 2*LEXP 00787 ELSE 00788 EXPSUM = 2*UEXP 00789 END IF 00790 * 00791 * EXPSUM is the exponent range, approximately equal to 00792 * EMAX - EMIN + 1 . 00793 * 00794 EMAX = EXPSUM + EMIN - 1 00795 NBITS = 1 + EXBITS + P 00796 * 00797 * NBITS is the total number of bits needed to store a 00798 * floating-point number. 00799 * 00800 IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN 00801 * 00802 * Either there are an odd number of bits used to store a 00803 * floating-point number, which is unlikely, or some bits are 00804 * not used in the representation of numbers, which is possible, 00805 * (e.g. Cray machines) or the mantissa has an implicit bit, 00806 * (e.g. IEEE machines, Dec Vax machines), which is perhaps the 00807 * most likely. We have to assume the last alternative. 00808 * If this is true, then we need to reduce EMAX by one because 00809 * there must be some way of representing zero in an implicit-bit 00810 * system. On machines like Cray, we are reducing EMAX by one 00811 * unnecessarily. 00812 * 00813 EMAX = EMAX - 1 00814 END IF 00815 * 00816 IF( IEEE ) THEN 00817 * 00818 * Assume we are on an IEEE machine which reserves one exponent 00819 * for infinity and NaN. 00820 * 00821 EMAX = EMAX - 1 00822 END IF 00823 * 00824 * Now create RMAX, the largest machine number, which should 00825 * be equal to (1.0 - BETA**(-P)) * BETA**EMAX . 00826 * 00827 * First compute 1.0 - BETA**(-P), being careful that the 00828 * result is less than 1.0 . 00829 * 00830 RECBAS = ONE / BETA 00831 Z = BETA - ONE 00832 Y = ZERO 00833 DO 20 I = 1, P 00834 Z = Z*RECBAS 00835 IF( Y.LT.ONE ) 00836 $ OLDY = Y 00837 Y = SLAMC3( Y, Z ) 00838 20 CONTINUE 00839 IF( Y.GE.ONE ) 00840 $ Y = OLDY 00841 * 00842 * Now multiply by BETA**EMAX to get RMAX. 00843 * 00844 DO 30 I = 1, EMAX 00845 Y = SLAMC3( Y*BETA, ZERO ) 00846 30 CONTINUE 00847 * 00848 RMAX = Y 00849 RETURN 00850 * 00851 * End of SLAMC5 00852 * 00853 END
1.7.3 
