dlamch.f.html

来自「famous linear algebra library (LAPACK) p」· HTML 代码 · 共 877 行 · 第 1/5 页

</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  BETA    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The base of the machine.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  T       (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of ( BETA ) digits in the mantissa.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RND     (output) LOGICAL
</span><span class="comment">*</span><span class="comment">          Specifies whether proper rounding  ( RND = .TRUE. )  or
</span><span class="comment">*</span><span class="comment">          chopping  ( RND = .FALSE. )  occurs in addition. This may not
</span><span class="comment">*</span><span class="comment">          be a reliable guide to the way in which the machine performs
</span><span class="comment">*</span><span class="comment">          its arithmetic.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  EPS     (output) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          The smallest positive number such that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">             fl( 1.0 - EPS ) .LT. 1.0,
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          where fl denotes the computed value.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  EMIN    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The minimum exponent before (gradual) underflow occurs.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RMIN    (output) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          The smallest normalized number for the machine, given by
</span><span class="comment">*</span><span class="comment">          BASE**( EMIN - 1 ), where  BASE  is the floating point value
</span><span class="comment">*</span><span class="comment">          of BETA.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  EMAX    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          The maximum exponent before overflow occurs.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  RMAX    (output) DOUBLE PRECISION
</span><span class="comment">*</span><span class="comment">          The largest positive number for the machine, given by
</span><span class="comment">*</span><span class="comment">          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
</span><span class="comment">*</span><span class="comment">          value of BETA.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Further Details
</span><span class="comment">*</span><span class="comment">  ===============
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  The computation of  EPS  is based on a routine PARANOIA by
</span><span class="comment">*</span><span class="comment">  W. Kahan of the University of California at Berkeley.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment"> =====================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
      INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
     $                   NGNMIN, NGPMIN
      DOUBLE PRECISION   A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
     $                   SIXTH, SMALL, THIRD, TWO, ZERO
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMC3.388"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>
      EXTERNAL           <a name="DLAMC3.389"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           <a name="DLAMC1.392"></a><a href="dlamch.f.html#DLAMC1.130">DLAMC1</a>, <a name="DLAMC4.392"></a><a href="dlamch.f.html#DLAMC4.612">DLAMC4</a>, <a name="DLAMC5.392"></a><a href="dlamch.f.html#DLAMC5.695">DLAMC5</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          ABS, MAX, MIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Save statement ..
</span>      SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
     $                   LRMIN, LT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Data statements ..
</span>      DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span>      IF( FIRST ) THEN
         ZERO = 0
         ONE = 1
         TWO = 2
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
</span><span class="comment">*</span><span class="comment">        BETA, T, RND, EPS, EMIN and RMIN.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Throughout this routine  we use the function  <a name="DLAMC3.414"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>  to ensure
</span><span class="comment">*</span><span class="comment">        that relevant values are stored  and not held in registers,  or
</span><span class="comment">*</span><span class="comment">        are not affected by optimizers.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        <a name="DLAMC1.418"></a><a href="dlamch.f.html#DLAMC1.130">DLAMC1</a> returns the parameters  LBETA, LT, LRND and LIEEE1.
</span><span class="comment">*</span><span class="comment">
</span>         CALL <a name="DLAMC1.420"></a><a href="dlamch.f.html#DLAMC1.130">DLAMC1</a>( LBETA, LT, LRND, LIEEE1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Start to find EPS.
</span><span class="comment">*</span><span class="comment">
</span>         B = LBETA
         A = B**( -LT )
         LEPS = A
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Try some tricks to see whether or not this is the correct  EPS.
</span><span class="comment">*</span><span class="comment">
</span>         B = TWO / 3
         HALF = ONE / 2
         SIXTH = <a name="DLAMC3.432"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( B, -HALF )
         THIRD = <a name="DLAMC3.433"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( SIXTH, SIXTH )
         B = <a name="DLAMC3.434"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( THIRD, -HALF )
         B = <a name="DLAMC3.435"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( B, SIXTH )
         B = ABS( B )
         IF( B.LT.LEPS )
     $      B = LEPS
<span class="comment">*</span><span class="comment">
</span>         LEPS = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
</span>   10    CONTINUE
         IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
            LEPS = B
            C = <a name="DLAMC3.446"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
            C = <a name="DLAMC3.447"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( HALF, -C )
            B = <a name="DLAMC3.448"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( HALF, C )
            C = <a name="DLAMC3.449"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( HALF, -B )
            B = <a name="DLAMC3.450"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( HALF, C )
            GO TO 10
         END IF
<span class="comment">*</span><span class="comment">+       END WHILE
</span><span class="comment">*</span><span class="comment">
</span>         IF( A.LT.LEPS )
     $      LEPS = A
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Computation of EPS complete.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
</span><span class="comment">*</span><span class="comment">        Keep dividing  A by BETA until (gradual) underflow occurs. This
</span><span class="comment">*</span><span class="comment">        is detected when we cannot recover the previous A.
</span><span class="comment">*</span><span class="comment">
</span>         RBASE = ONE / LBETA
         SMALL = ONE
         DO 20 I = 1, 3
            SMALL = <a name="DLAMC3.467"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( SMALL*RBASE, ZERO )
   20    CONTINUE
         A = <a name="DLAMC3.469"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( ONE, SMALL )
         CALL <a name="DLAMC4.470"></a><a href="dlamch.f.html#DLAMC4.612">DLAMC4</a>( NGPMIN, ONE, LBETA )
         CALL <a name="DLAMC4.471"></a><a href="dlamch.f.html#DLAMC4.612">DLAMC4</a>( NGNMIN, -ONE, LBETA )
         CALL <a name="DLAMC4.472"></a><a href="dlamch.f.html#DLAMC4.612">DLAMC4</a>( GPMIN, A, LBETA )
         CALL <a name="DLAMC4.473"></a><a href="dlamch.f.html#DLAMC4.612">DLAMC4</a>( GNMIN, -A, LBETA )
         IEEE = .FALSE.
<span class="comment">*</span><span class="comment">
</span>         IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
            IF( NGPMIN.EQ.GPMIN ) THEN
               LEMIN = NGPMIN
<span class="comment">*</span><span class="comment">            ( Non twos-complement machines, no gradual underflow;
</span><span class="comment">*</span><span class="comment">              e.g.,  VAX )
</span>            ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
               LEMIN = NGPMIN - 1 + LT
               IEEE = .TRUE.
<span class="comment">*</span><span class="comment">            ( Non twos-complement machines, with gradual underflow;
</span><span class="comment">*</span><span class="comment">              e.g., IEEE standard followers )
</span>            ELSE
               LEMIN = MIN( NGPMIN, GPMIN )
<span class="comment">*</span><span class="comment">            ( A guess; no known machine )
</span>               IWARN = .TRUE.
            END IF
<span class="comment">*</span><span class="comment">
</span>         ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
            IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
               LEMIN = MAX( NGPMIN, NGNMIN )
<span class="comment">*</span><span class="comment">            ( Twos-complement machines, no gradual underflow;
</span><span class="comment">*</span><span class="comment">              e.g., CYBER 205 )
</span>            ELSE
               LEMIN = MIN( NGPMIN, NGNMIN )
<span class="comment">*</span><span class="comment">            ( A guess; no known machine )
</span>               IWARN = .TRUE.
            END IF
<span class="comment">*</span><span class="comment">
</span>         ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
     $            ( GPMIN.EQ.GNMIN ) ) THEN
            IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
               LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
<span class="comment">*</span><span class="comment">            ( Twos-complement machines with gradual underflow;
</span><span class="comment">*</span><span class="comment">              no known machine )
</span>            ELSE