dlamch.f.html

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

HTML
877 行
字号
12345
</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">  IEEE1   (output) LOGICAL
</span><span class="comment">*</span><span class="comment">          Specifies whether rounding appears to be done in the IEEE
</span><span class="comment">*</span><span class="comment">          'round to nearest' style.
</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 routine is based on the routine  ENVRON  by Malcolm and
</span><span class="comment">*</span><span class="comment">  incorporates suggestions by Gentleman and Marovich. See
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Malcolm M. A. (1972) Algorithms to reveal properties of
</span><span class="comment">*</span><span class="comment">        floating-point arithmetic. Comms. of the ACM, 15, 949-951.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Gentleman W. M. and Marovich S. B. (1974) More on algorithms
</span><span class="comment">*</span><span class="comment">        that reveal properties of floating point arithmetic units.
</span><span class="comment">*</span><span class="comment">        Comms. of the ACM, 17, 276-277.
</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, LIEEE1, LRND
      INTEGER            LBETA, LT
      DOUBLE PRECISION   A, B, C, F, ONE, QTR, SAVEC, T1, T2
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      DOUBLE PRECISION   <a name="DLAMC3.187"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>
      EXTERNAL           <a name="DLAMC3.188"></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">     .. Save statement ..
</span>      SAVE               FIRST, LIEEE1, LBETA, LRND, LT
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Data statements ..
</span>      DATA               FIRST / .TRUE. /
<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
         ONE = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,
</span><span class="comment">*</span><span class="comment">        IEEE1, T and RND.
</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.204"></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">        Compute  a = 2.0**m  with the  smallest positive integer m such
</span><span class="comment">*</span><span class="comment">        that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           fl( a + 1.0 ) = a.
</span><span class="comment">*</span><span class="comment">
</span>         A = 1
         C = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">+       WHILE( C.EQ.ONE )LOOP
</span>   10    CONTINUE
         IF( C.EQ.ONE ) THEN
            A = 2*A
            C = <a name="DLAMC3.220"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( A, ONE )
            C = <a name="DLAMC3.221"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( C, -A )
            GO TO 10
         END IF
<span class="comment">*</span><span class="comment">+       END WHILE
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now compute  b = 2.0**m  with the smallest positive integer m
</span><span class="comment">*</span><span class="comment">        such that
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           fl( a + b ) .gt. a.
</span><span class="comment">*</span><span class="comment">
</span>         B = 1
         C = <a name="DLAMC3.232"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( A, B )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">+       WHILE( C.EQ.A )LOOP
</span>   20    CONTINUE
         IF( C.EQ.A ) THEN
            B = 2*B
            C = <a name="DLAMC3.238"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( A, B )
            GO TO 20
         END IF
<span class="comment">*</span><span class="comment">+       END WHILE
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now compute the base.  a and c  are neighbouring floating point
</span><span class="comment">*</span><span class="comment">        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so
</span><span class="comment">*</span><span class="comment">        their difference is beta. Adding 0.25 to c is to ensure that it
</span><span class="comment">*</span><span class="comment">        is truncated to beta and not ( beta - 1 ).
</span><span class="comment">*</span><span class="comment">
</span>         QTR = ONE / 4
         SAVEC = C
         C = <a name="DLAMC3.250"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( C, -A )
         LBETA = C + QTR
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now determine whether rounding or chopping occurs,  by adding a
</span><span class="comment">*</span><span class="comment">        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.
</span><span class="comment">*</span><span class="comment">
</span>         B = LBETA
         F = <a name="DLAMC3.257"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( B / 2, -B / 100 )
         C = <a name="DLAMC3.258"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( F, A )
         IF( C.EQ.A ) THEN
            LRND = .TRUE.
         ELSE
            LRND = .FALSE.
         END IF
         F = <a name="DLAMC3.264"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( B / 2, B / 100 )
         C = <a name="DLAMC3.265"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( F, A )
         IF( ( LRND ) .AND. ( C.EQ.A ) )
     $      LRND = .FALSE.
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Try and decide whether rounding is done in the  IEEE  'round to
</span><span class="comment">*</span><span class="comment">        nearest' style. B/2 is half a unit in the last place of the two
</span><span class="comment">*</span><span class="comment">        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit
</span><span class="comment">*</span><span class="comment">        zero, and SAVEC is odd. Thus adding B/2 to A should not  change
</span><span class="comment">*</span><span class="comment">        A, but adding B/2 to SAVEC should change SAVEC.
</span><span class="comment">*</span><span class="comment">
</span>         T1 = <a name="DLAMC3.275"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( B / 2, A )
         T2 = <a name="DLAMC3.276"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( B / 2, SAVEC )
         LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Now find  the  mantissa, t.  It should  be the  integer part of
</span><span class="comment">*</span><span class="comment">        log to the base beta of a,  however it is safer to determine  t
</span><span class="comment">*</span><span class="comment">        by powering.  So we find t as the smallest positive integer for
</span><span class="comment">*</span><span class="comment">        which
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">           fl( beta**t + 1.0 ) = 1.0.
</span><span class="comment">*</span><span class="comment">
</span>         LT = 0
         A = 1
         C = 1
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">+       WHILE( C.EQ.ONE )LOOP
</span>   30    CONTINUE
         IF( C.EQ.ONE ) THEN
            LT = LT + 1
            A = A*LBETA
            C = <a name="DLAMC3.295"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( A, ONE )
            C = <a name="DLAMC3.296"></a><a href="dlamch.f.html#DLAMC3.574">DLAMC3</a>( C, -A )
            GO TO 30
         END IF
<span class="comment">*</span><span class="comment">+       END WHILE
</span><span class="comment">*</span><span class="comment">
</span>      END IF
<span class="comment">*</span><span class="comment">
</span>      BETA = LBETA
      T = LT
      RND = LRND
      IEEE1 = LIEEE1
      FIRST = .FALSE.
      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="DLAMC1.310"></a><a href="dlamch.f.html#DLAMC1.130">DLAMC1</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="DLAMC2.316"></a>      SUBROUTINE <a name="DLAMC2.316"></a><a href="dlamch.f.html#DLAMC2.316">DLAMC2</a>( BETA, T, RND, EPS, EMIN, RMIN, 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            RND
      INTEGER            BETA, EMAX, EMIN, T
      DOUBLE PRECISION   EPS, RMAX, RMIN
<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="DLAMC2.331"></a><a href="dlamch.f.html#DLAMC2.316">DLAMC2</a> determines the machine parameters specified in its argument
</span><span class="comment">*</span><span class="comment">  list.
</span><span class="comment">*</span><span class="comment">
12345

dlamch.f.html - 源码说明

本页面展示了「famous linear algebra library (LAPACK) ports to windows」中的 dlamch.f.html 源码文件，采用 HTML 编程语言编写，共 877 行代码。您可以在线阅读完整代码内容，也可以返回资源详情页下载完整源码包进行本地学习和开发。

虫虫下载站收录了大量与algebra相关的技术资源，包括源代码、技术文档、电路图等，是电子工程师和嵌入式开发者的专业学习平台。

⌨️ 快捷键说明

复制代码Ctrl + C
搜索代码Ctrl + F
全屏模式F11
增大字号Ctrl + =
减小字号Ctrl + -
显示快捷键?