dlamch.f.html

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

</span><span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DLAMC4.689"></a><a href="dlamch.f.html#DLAMC4.612">DLAMC4</a>
</span><span class="comment">*</span><span class="comment">
</span>      END
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">***********************************************************************
</span><span class="comment">*</span><span class="comment">
</span><a name="DLAMC5.695"></a>      SUBROUTINE <a name="DLAMC5.695"></a><a href="dlamch.f.html#DLAMC5.695">DLAMC5</a>( BETA, P, EMIN, IEEE, EMAX, RMAX )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK auxiliary routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      LOGICAL            IEEE
      INTEGER            BETA, EMAX, EMIN, P
      DOUBLE PRECISION   RMAX
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="DLAMC5.710"></a><a href="dlamch.f.html#DLAMC5.695">DLAMC5</a> attempts to compute RMAX, the largest machine floating-point
</span><span class="comment">*</span><span class="comment">  number, without overflow.  It assumes that EMAX + abs(EMIN) sum
</span><span class="comment">*</span><span class="comment">  approximately to a power of 2.  It will fail on machines where this
</span><span class="comment">*</span><span class="comment">  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,
</span><span class="comment">*</span><span class="comment">  EMAX = 28718).  It will also fail if the value supplied for EMIN is
</span><span class="comment">*</span><span class="comment">  too large (i.e. too close to zero), probably with overflow.
</span><span class="comment">*</span><span class="comment">
</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    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The base of floating-point arithmetic.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  P       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of base BETA digits in the mantissa of a
</span><span class="comment">*</span><span class="comment">          floating-point value.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  EMIN    (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The minimum exponent before (gradual) underflow.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  IEEE    (input) LOGICAL
</span><span class="comment">*</span><span class="comment">          A logical flag specifying whether or not the arithmetic
</span><span class="comment">*</span><span class="comment">          system is thought to comply with the IEEE standard.
</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 largest exponent before overflow
</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 machine floating-point number.
</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">     .. Parameters ..
</span>      DOUBLE PRECISION   ZERO, ONE
      PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP
      DOUBLE PRECISION   OLDY, RECBAS, Y, Z
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMC3.751"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>
      EXTERNAL           <a name="DLAMC3.752"></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">     .. Intrinsic Functions ..
</span>      INTRINSIC          MOD
<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><span class="comment">*</span><span class="comment">     First compute LEXP and UEXP, two powers of 2 that bound
</span><span class="comment">*</span><span class="comment">     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum
</span><span class="comment">*</span><span class="comment">     approximately to the bound that is closest to abs(EMIN).
</span><span class="comment">*</span><span class="comment">     (EMAX is the exponent of the required number RMAX).
</span><span class="comment">*</span><span class="comment">
</span>      LEXP = 1
      EXBITS = 1
   10 CONTINUE
      TRY = LEXP*2
      IF( TRY.LE.( -EMIN ) ) THEN
         LEXP = TRY
         EXBITS = EXBITS + 1
         GO TO 10
      END IF
      IF( LEXP.EQ.-EMIN ) THEN
         UEXP = LEXP
      ELSE
         UEXP = TRY
         EXBITS = EXBITS + 1
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Now -LEXP is less than or equal to EMIN, and -UEXP is greater
</span><span class="comment">*</span><span class="comment">     than or equal to EMIN. EXBITS is the number of bits needed to
</span><span class="comment">*</span><span class="comment">     store the exponent.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN
         EXPSUM = 2*LEXP
      ELSE
         EXPSUM = 2*UEXP
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     EXPSUM is the exponent range, approximately equal to
</span><span class="comment">*</span><span class="comment">     EMAX - EMIN + 1 .
</span><span class="comment">*</span><span class="comment">
</span>      EMAX = EXPSUM + EMIN - 1
      NBITS = 1 + EXBITS + P
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     NBITS is the total number of bits needed to store a
</span><span class="comment">*</span><span class="comment">     floating-point number.
</span><span class="comment">*</span><span class="comment">
</span>      IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Either there are an odd number of bits used to store a
</span><span class="comment">*</span><span class="comment">        floating-point number, which is unlikely, or some bits are
</span><span class="comment">*</span><span class="comment">        not used in the representation of numbers, which is possible,
</span><span class="comment">*</span><span class="comment">        (e.g. Cray machines) or the mantissa has an implicit bit,
</span><span class="comment">*</span><span class="comment">        (e.g. IEEE machines, Dec Vax machines), which is perhaps the
</span><span class="comment">*</span><span class="comment">        most likely. We have to assume the last alternative.
</span><span class="comment">*</span><span class="comment">        If this is true, then we need to reduce EMAX by one because
</span><span class="comment">*</span><span class="comment">        there must be some way of representing zero in an implicit-bit
</span><span class="comment">*</span><span class="comment">        system. On machines like Cray, we are reducing EMAX by one
</span><span class="comment">*</span><span class="comment">        unnecessarily.
</span><span class="comment">*</span><span class="comment">
</span>         EMAX = EMAX - 1
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( IEEE ) THEN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Assume we are on an IEEE machine which reserves one exponent
</span><span class="comment">*</span><span class="comment">        for infinity and NaN.
</span><span class="comment">*</span><span class="comment">
</span>         EMAX = EMAX - 1
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Now create RMAX, the largest machine number, which should
</span><span class="comment">*</span><span class="comment">     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     First compute 1.0 - BETA**(-P), being careful that the
</span><span class="comment">*</span><span class="comment">     result is less than 1.0 .
</span><span class="comment">*</span><span class="comment">
</span>      RECBAS = ONE / BETA
      Z = BETA - ONE
      Y = ZERO
      DO 20 I = 1, P
         Z = Z*RECBAS
         IF( Y.LT.ONE )
     $      OLDY = Y
         Y = <a name="DLAMC3.836"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( Y, Z )
   20 CONTINUE
      IF( Y.GE.ONE )
     $   Y = OLDY
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Now multiply by BETA**EMAX to get RMAX.
</span><span class="comment">*</span><span class="comment">
</span>      DO 30 I = 1, EMAX
         Y = <a name="DLAMC3.844"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( Y*BETA, ZERO )
   30 CONTINUE
<span class="comment">*</span><span class="comment">
</span>      RMAX = Y
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DLAMC5.850"></a><a href="dlamch.f.html#DLAMC5.695">DLAMC5</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>