LAPACK  3.4.1
LAPACK: Linear Algebra PACKage
dlamchf77.f
Go to the documentation of this file.
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
 All Files Functions