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