home *** CD-ROM | disk | FTP | other *** search
/ Geek Gadgets 1 / ADE-1.bin / ade-dist / octave-1.1.1p1-src.tgz / tar.out / fsf / octave / libcruft / lapack / dlamc2.f < prev    next >
Text File  |  1996-09-28  |  8KB  |  260 lines
  1. *
  2. ************************************************************************
  3. *
  4.       SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )
  5. *
  6. *  -- LAPACK auxiliary routine (version 2.0) --
  7. *     Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
  8. *     Courant Institute, Argonne National Lab, and Rice University
  9. *     October 31, 1992
  10. *
  11. *     .. Scalar Arguments ..
  12.       LOGICAL            RND
  13.       INTEGER            BETA, EMAX, EMIN, T
  14.       DOUBLE PRECISION   EPS, RMAX, RMIN
  15. *     ..
  16. *
  17. *  Purpose
  18. *  =======
  19. *
  20. *  DLAMC2 determines the machine parameters specified in its argument
  21. *  list.
  22. *
  23. *  Arguments
  24. *  =========
  25. *
  26. *  BETA    (output) INTEGER
  27. *          The base of the machine.
  28. *
  29. *  T       (output) INTEGER
  30. *          The number of ( BETA ) digits in the mantissa.
  31. *
  32. *  RND     (output) LOGICAL
  33. *          Specifies whether proper rounding  ( RND = .TRUE. )  or
  34. *          chopping  ( RND = .FALSE. )  occurs in addition. This may not
  35. *          be a reliable guide to the way in which the machine performs
  36. *          its arithmetic.
  37. *
  38. *  EPS     (output) DOUBLE PRECISION
  39. *          The smallest positive number such that
  40. *
  41. *             fl( 1.0 - EPS ) .LT. 1.0,
  42. *
  43. *          where fl denotes the computed value.
  44. *
  45. *  EMIN    (output) INTEGER
  46. *          The minimum exponent before (gradual) underflow occurs.
  47. *
  48. *  RMIN    (output) DOUBLE PRECISION
  49. *          The smallest normalized number for the machine, given by
  50. *          BASE**( EMIN - 1 ), where  BASE  is the floating point value
  51. *          of BETA.
  52. *
  53. *  EMAX    (output) INTEGER
  54. *          The maximum exponent before overflow occurs.
  55. *
  56. *  RMAX    (output) DOUBLE PRECISION
  57. *          The largest positive number for the machine, given by
  58. *          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point
  59. *          value of BETA.
  60. *
  61. *  Further Details
  62. *  ===============
  63. *
  64. *  The computation of  EPS  is based on a routine PARANOIA by
  65. *  W. Kahan of the University of California at Berkeley.
  66. *
  67. * =====================================================================
  68. *
  69. *     .. Local Scalars ..
  70.       LOGICAL            FIRST, IEEE, IWARN, LIEEE1, LRND
  71.       INTEGER            GNMIN, GPMIN, I, LBETA, LEMAX, LEMIN, LT,
  72.      $                   NGNMIN, NGPMIN
  73.       DOUBLE PRECISION   A, B, C, HALF, LEPS, LRMAX, LRMIN, ONE, RBASE,
  74.      $                   SIXTH, SMALL, THIRD, TWO, ZERO
  75. *     ..
  76. *     .. External Functions ..
  77.       DOUBLE PRECISION   DLAMC3
  78.       EXTERNAL           DLAMC3
  79. *     ..
  80. *     .. External Subroutines ..
  81.       EXTERNAL           DLAMC1, DLAMC4, DLAMC5
  82. *     ..
  83. *     .. Intrinsic Functions ..
  84.       INTRINSIC          ABS, MAX, MIN
  85. *     ..
  86. *     .. Save statement ..
  87.       SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,
  88.      $                   LRMIN, LT
  89. *     ..
  90. *     .. Data statements ..
  91.       DATA               FIRST / .TRUE. / , IWARN / .FALSE. /
  92. *     ..
  93. *     .. Executable Statements ..
  94. *
  95.       IF( FIRST ) THEN
  96.          FIRST = .FALSE.
  97.          ZERO = 0
  98.          ONE = 1
  99.          TWO = 2
  100. *
  101. *        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of
  102. *        BETA, T, RND, EPS, EMIN and RMIN.
  103. *
  104. *        Throughout this routine  we use the function  DLAMC3  to ensure
  105. *        that relevant values are stored  and not held in registers,  or
  106. *        are not affected by optimizers.
  107. *
  108. *        DLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.
  109. *
  110.          CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )
  111. *
  112. *        Start to find EPS.
  113. *
  114.          B = LBETA
  115.          A = B**( -LT )
  116.          LEPS = A
  117. *
  118. *        Try some tricks to see whether or not this is the correct  EPS.
  119. *
  120.          B = TWO / 3
  121.          HALF = ONE / 2
  122.          SIXTH = DLAMC3( B, -HALF )
  123.          THIRD = DLAMC3( SIXTH, SIXTH )
  124.          B = DLAMC3( THIRD, -HALF )
  125.          B = DLAMC3( B, SIXTH )
  126.          B = ABS( B )
  127.          IF( B.LT.LEPS )
  128.      $      B = LEPS
  129. *
  130.          LEPS = 1
  131. *
  132. *+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP
  133.    10    CONTINUE
  134.          IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN
  135.             LEPS = B
  136.             C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )
  137.             C = DLAMC3( HALF, -C )
  138.             B = DLAMC3( HALF, C )
  139.             C = DLAMC3( HALF, -B )
  140.             B = DLAMC3( HALF, C )
  141.             GO TO 10
  142.          END IF
  143. *+       END WHILE
  144. *
  145.          IF( A.LT.LEPS )
  146.      $      LEPS = A
  147. *
  148. *        Computation of EPS complete.
  149. *
  150. *        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).
  151. *        Keep dividing  A by BETA until (gradual) underflow occurs. This
  152. *        is detected when we cannot recover the previous A.
  153. *
  154.          RBASE = ONE / LBETA
  155.          SMALL = ONE
  156.          DO 20 I = 1, 3
  157.             SMALL = DLAMC3( SMALL*RBASE, ZERO )
  158.    20    CONTINUE
  159.          A = DLAMC3( ONE, SMALL )
  160.          CALL DLAMC4( NGPMIN, ONE, LBETA )
  161.          CALL DLAMC4( NGNMIN, -ONE, LBETA )
  162.          CALL DLAMC4( GPMIN, A, LBETA )
  163.          CALL DLAMC4( GNMIN, -A, LBETA )
  164.          IEEE = .FALSE.
  165. *
  166.          IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN
  167.             IF( NGPMIN.EQ.GPMIN ) THEN
  168.                LEMIN = NGPMIN
  169. *            ( Non twos-complement machines, no gradual underflow;
  170. *              e.g.,  VAX )
  171.             ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN
  172.                LEMIN = NGPMIN - 1 + LT
  173.                IEEE = .TRUE.
  174. *            ( Non twos-complement machines, with gradual underflow;
  175. *              e.g., IEEE standard followers )
  176.             ELSE
  177.                LEMIN = MIN( NGPMIN, GPMIN )
  178. *            ( A guess; no known machine )
  179.                IWARN = .TRUE.
  180.             END IF
  181. *
  182.          ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN
  183.             IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN
  184.                LEMIN = MAX( NGPMIN, NGNMIN )
  185. *            ( Twos-complement machines, no gradual underflow;
  186. *              e.g., CYBER 205 )
  187.             ELSE
  188.                LEMIN = MIN( NGPMIN, NGNMIN )
  189. *            ( A guess; no known machine )
  190.                IWARN = .TRUE.
  191.             END IF
  192. *
  193.          ELSE IF( ( ABS( NGPMIN-NGNMIN ).EQ.1 ) .AND.
  194.      $            ( GPMIN.EQ.GNMIN ) ) THEN
  195.             IF( ( GPMIN-MIN( NGPMIN, NGNMIN ) ).EQ.3 ) THEN
  196.                LEMIN = MAX( NGPMIN, NGNMIN ) - 1 + LT
  197. *            ( Twos-complement machines with gradual underflow;
  198. *              no known machine )
  199.             ELSE
  200.                LEMIN = MIN( NGPMIN, NGNMIN )
  201. *            ( A guess; no known machine )
  202.                IWARN = .TRUE.
  203.             END IF
  204. *
  205.          ELSE
  206.             LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )
  207. *         ( A guess; no known machine )
  208.             IWARN = .TRUE.
  209.          END IF
  210. ***
  211. * Comment out this if block if EMIN is ok
  212.          IF( IWARN ) THEN
  213.             FIRST = .TRUE.
  214.             WRITE( 6, FMT = 9999 )LEMIN
  215.          END IF
  216. ***
  217. *
  218. *        Assume IEEE arithmetic if we found denormalised  numbers above,
  219. *        or if arithmetic seems to round in the  IEEE style,  determined
  220. *        in routine DLAMC1. A true IEEE machine should have both  things
  221. *        true; however, faulty machines may have one or the other.
  222. *
  223.          IEEE = IEEE .OR. LIEEE1
  224. *
  225. *        Compute  RMIN by successive division by  BETA. We could compute
  226. *        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during
  227. *        this computation.
  228. *
  229.          LRMIN = 1
  230.          DO 30 I = 1, 1 - LEMIN
  231.             LRMIN = DLAMC3( LRMIN*RBASE, ZERO )
  232.    30    CONTINUE
  233. *
  234. *        Finally, call DLAMC5 to compute EMAX and RMAX.
  235. *
  236.          CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )
  237.       END IF
  238. *
  239.       BETA = LBETA
  240.       T = LT
  241.       RND = LRND
  242.       EPS = LEPS
  243.       EMIN = LEMIN
  244.       RMIN = LRMIN
  245.       EMAX = LEMAX
  246.       RMAX = LRMAX
  247. *
  248.       RETURN
  249. *
  250.  9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',
  251.      $      '  EMIN = ', I8, /
  252.      $      ' If, after inspection, the value EMIN looks',
  253.      $      ' acceptable please comment out ',
  254.      $      / ' the IF block as marked within the code of routine',
  255.      $      ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )
  256. *
  257. *     End of DLAMC2
  258. *
  259.       END
  260.