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