[BACK]Return to dlamch.f CVS log [TXT] [DIR] Up to [local] / rpl / lapack / lapack

Diff for /rpl/lapack/lapack/dlamch.f between versions 1.6 and 1.7

version 1.6, 2010/08/13 21:03:49 version 1.7, 2010/12/21 13:48:05
Line 1 Line 1
       DOUBLE PRECISION FUNCTION DLAMCH( CMACH )        DOUBLE PRECISION FUNCTION DLAMCH( CMACH )
 *  *
 *  -- LAPACK auxiliary routine (version 3.2) --  *  -- LAPACK auxiliary routine (version 3.3.0) --
 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..  *  -- LAPACK is a software package provided by Univ. of Tennessee,    --
 *     November 2006  *  -- Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd..--
   *     Based on LAPACK DLAMCH but with Fortran 95 query functions
   *     See: http://www.cs.utk.edu/~luszczek/lapack/lamch.html
   *     and  http://www.netlib.org/lapack-dev/lapack-coding/program-style.html#id2537289
   *     July 2010
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       CHARACTER          CMACH        CHARACTER          CMACH
Line 49 Line 53
       PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )        PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
 *     ..  *     ..
 *     .. Local Scalars ..  *     .. Local Scalars ..
       LOGICAL            FIRST, LRND        DOUBLE PRECISION   RND, EPS, SFMIN, SMALL, RMACH
       INTEGER            BETA, IMAX, IMIN, IT  
       DOUBLE PRECISION   BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,  
      $                   RND, SFMIN, SMALL, T  
 *     ..  *     ..
 *     .. External Functions ..  *     .. External Functions ..
       LOGICAL            LSAME        LOGICAL            LSAME
       EXTERNAL           LSAME        EXTERNAL           LSAME
 *     ..  *     ..
 *     .. External Subroutines ..  *     .. Intrinsic Functions ..
       EXTERNAL           DLAMC2        INTRINSIC          DIGITS, EPSILON, HUGE, MAXEXPONENT,
 *     ..       $                   MINEXPONENT, RADIX, TINY
 *     .. Save statement ..  
       SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,  
      $                   EMAX, RMAX, PREC  
 *     ..  
 *     .. Data statements ..  
       DATA               FIRST / .TRUE. /  
 *     ..  *     ..
 *     .. Executable Statements ..  *     .. Executable Statements ..
 *  *
       IF( FIRST ) THEN  
          CALL DLAMC2( BETA, IT, LRND, EPS, IMIN, RMIN, IMAX, RMAX )  
          BASE = BETA  
          T = IT  
          IF( LRND ) THEN  
             RND = ONE  
             EPS = ( BASE**( 1-IT ) ) / 2  
          ELSE  
             RND = ZERO  
             EPS = BASE**( 1-IT )  
          END IF  
          PREC = EPS*BASE  
          EMIN = IMIN  
          EMAX = IMAX  
          SFMIN = RMIN  
          SMALL = ONE / RMAX  
          IF( SMALL.GE.SFMIN ) THEN  
 *  *
 *           Use SMALL plus a bit, to avoid the possibility of rounding  *     Assume rounding, not chopping. Always.
 *           causing overflow when computing  1/sfmin.  
 *  *
             SFMIN = SMALL*( ONE+EPS )        RND = ONE
          END IF  *
         IF( ONE.EQ.RND ) THEN
            EPS = EPSILON(ZERO) * 0.5
         ELSE
            EPS = EPSILON(ZERO)
       END IF        END IF
 *  *
       IF( LSAME( CMACH, 'E' ) ) THEN        IF( LSAME( CMACH, 'E' ) ) THEN
          RMACH = EPS           RMACH = EPS
       ELSE IF( LSAME( CMACH, 'S' ) ) THEN        ELSE IF( LSAME( CMACH, 'S' ) ) THEN
            SFMIN = TINY(ZERO)
            SMALL = ONE / HUGE(ZERO)
            IF( SMALL.GE.SFMIN ) THEN
   *
   *           Use SMALL plus a bit, to avoid the possibility of rounding
   *           causing overflow when computing  1/sfmin.
   *
               SFMIN = SMALL*( ONE+EPS )
            END IF
          RMACH = SFMIN           RMACH = SFMIN
       ELSE IF( LSAME( CMACH, 'B' ) ) THEN        ELSE IF( LSAME( CMACH, 'B' ) ) THEN
          RMACH = BASE           RMACH = RADIX(ZERO)
       ELSE IF( LSAME( CMACH, 'P' ) ) THEN        ELSE IF( LSAME( CMACH, 'P' ) ) THEN
          RMACH = PREC           RMACH = EPS * RADIX(ZERO)
       ELSE IF( LSAME( CMACH, 'N' ) ) THEN        ELSE IF( LSAME( CMACH, 'N' ) ) THEN
          RMACH = T           RMACH = DIGITS(ZERO)
       ELSE IF( LSAME( CMACH, 'R' ) ) THEN        ELSE IF( LSAME( CMACH, 'R' ) ) THEN
          RMACH = RND           RMACH = RND
       ELSE IF( LSAME( CMACH, 'M' ) ) THEN        ELSE IF( LSAME( CMACH, 'M' ) ) THEN
          RMACH = EMIN           RMACH = MINEXPONENT(ZERO)
       ELSE IF( LSAME( CMACH, 'U' ) ) THEN        ELSE IF( LSAME( CMACH, 'U' ) ) THEN
          RMACH = RMIN           RMACH = tiny(zero)
       ELSE IF( LSAME( CMACH, 'L' ) ) THEN        ELSE IF( LSAME( CMACH, 'L' ) ) THEN
          RMACH = EMAX           RMACH = MAXEXPONENT(ZERO)
       ELSE IF( LSAME( CMACH, 'O' ) ) THEN        ELSE IF( LSAME( CMACH, 'O' ) ) THEN
          RMACH = RMAX           RMACH = HUGE(ZERO)
         ELSE
            RMACH = ZERO
       END IF        END IF
 *  *
       DLAMCH = RMACH        DLAMCH = RMACH
       FIRST  = .FALSE.  
       RETURN        RETURN
 *  *
 *     End of DLAMCH  *     End of DLAMCH
 *  *
       END        END
 *  
 ************************************************************************  
 *  
       SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 )  
 *  
 *  -- LAPACK auxiliary routine (version 3.2) --  
 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..  
 *     November 2006  
 *  
 *     .. Scalar Arguments ..  
       LOGICAL            IEEE1, RND  
       INTEGER            BETA, T  
 *     ..  
 *  
 *  Purpose  
 *  =======  
 *  
 *  DLAMC1 determines the machine parameters given by BETA, T, RND, and  
 *  IEEE1.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  BETA    (output) INTEGER  
 *          The base of the machine.  
 *  
 *  T       (output) INTEGER  
 *          The number of ( BETA ) digits in the mantissa.  
 *  
 *  RND     (output) LOGICAL  
 *          Specifies whether proper rounding  ( RND = .TRUE. )  or  
 *          chopping  ( RND = .FALSE. )  occurs in addition. This may not  
 *          be a reliable guide to the way in which the machine performs  
 *          its arithmetic.  
 *  
 *  IEEE1   (output) LOGICAL  
 *          Specifies whether rounding appears to be done in the IEEE  
 *          'round to nearest' style.  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  The routine is based on the routine  ENVRON  by Malcolm and  
 *  incorporates suggestions by Gentleman and Marovich. See  
 *  
 *     Malcolm M. A. (1972) Algorithms to reveal properties of  
 *        floating-point arithmetic. Comms. of the ACM, 15, 949-951.  
 *  
 *     Gentleman W. M. and Marovich S. B. (1974) More on algorithms  
 *        that reveal properties of floating point arithmetic units.  
 *        Comms. of the ACM, 17, 276-277.  
 *  
 * =====================================================================  
 *  
 *     .. Local Scalars ..  
       LOGICAL            FIRST, LIEEE1, LRND  
       INTEGER            LBETA, LT  
       DOUBLE PRECISION   A, B, C, F, ONE, QTR, SAVEC, T1, T2  
 *     ..  
 *     .. External Functions ..  
       DOUBLE PRECISION   DLAMC3  
       EXTERNAL           DLAMC3  
 *     ..  
 *     .. Save statement ..  
       SAVE               FIRST, LIEEE1, LBETA, LRND, LT  
 *     ..  
 *     .. Data statements ..  
       DATA               FIRST / .TRUE. /  
 *     ..  
 *     .. Executable Statements ..  
 *  
       IF( FIRST ) THEN  
          ONE = 1  
 *  
 *        LBETA,  LIEEE1,  LT and  LRND  are the  local values  of  BETA,  
 *        IEEE1, T and RND.  
 *  
 *        Throughout this routine  we use the function  DLAMC3  to ensure  
 *        that relevant values are  stored and not held in registers,  or  
 *        are not affected by optimizers.  
 *  
 *        Compute  a = 2.0**m  with the  smallest positive integer m such  
 *        that  
 *  
 *           fl( a + 1.0 ) = a.  
 *  
          A = 1  
          C = 1  
 *  
 *+       WHILE( C.EQ.ONE )LOOP  
    10    CONTINUE  
          IF( C.EQ.ONE ) THEN  
             A = 2*A  
             C = DLAMC3( A, ONE )  
             C = DLAMC3( C, -A )  
             GO TO 10  
          END IF  
 *+       END WHILE  
 *  
 *        Now compute  b = 2.0**m  with the smallest positive integer m  
 *        such that  
 *  
 *           fl( a + b ) .gt. a.  
 *  
          B = 1  
          C = DLAMC3( A, B )  
 *  
 *+       WHILE( C.EQ.A )LOOP  
    20    CONTINUE  
          IF( C.EQ.A ) THEN  
             B = 2*B  
             C = DLAMC3( A, B )  
             GO TO 20  
          END IF  
 *+       END WHILE  
 *  
 *        Now compute the base.  a and c  are neighbouring floating point  
 *        numbers  in the  interval  ( beta**t, beta**( t + 1 ) )  and so  
 *        their difference is beta. Adding 0.25 to c is to ensure that it  
 *        is truncated to beta and not ( beta - 1 ).  
 *  
          QTR = ONE / 4  
          SAVEC = C  
          C = DLAMC3( C, -A )  
          LBETA = C + QTR  
 *  
 *        Now determine whether rounding or chopping occurs,  by adding a  
 *        bit  less  than  beta/2  and a  bit  more  than  beta/2  to  a.  
 *  
          B = LBETA  
          F = DLAMC3( B / 2, -B / 100 )  
          C = DLAMC3( F, A )  
          IF( C.EQ.A ) THEN  
             LRND = .TRUE.  
          ELSE  
             LRND = .FALSE.  
          END IF  
          F = DLAMC3( B / 2, B / 100 )  
          C = DLAMC3( F, A )  
          IF( ( LRND ) .AND. ( C.EQ.A ) )  
      $      LRND = .FALSE.  
 *  
 *        Try and decide whether rounding is done in the  IEEE  'round to  
 *        nearest' style. B/2 is half a unit in the last place of the two  
 *        numbers A and SAVEC. Furthermore, A is even, i.e. has last  bit  
 *        zero, and SAVEC is odd. Thus adding B/2 to A should not  change  
 *        A, but adding B/2 to SAVEC should change SAVEC.  
 *  
          T1 = DLAMC3( B / 2, A )  
          T2 = DLAMC3( B / 2, SAVEC )  
          LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND  
 *  
 *        Now find  the  mantissa, t.  It should  be the  integer part of  
 *        log to the base beta of a,  however it is safer to determine  t  
 *        by powering.  So we find t as the smallest positive integer for  
 *        which  
 *  
 *           fl( beta**t + 1.0 ) = 1.0.  
 *  
          LT = 0  
          A = 1  
          C = 1  
 *  
 *+       WHILE( C.EQ.ONE )LOOP  
    30    CONTINUE  
          IF( C.EQ.ONE ) THEN  
             LT = LT + 1  
             A = A*LBETA  
             C = DLAMC3( A, ONE )  
             C = DLAMC3( C, -A )  
             GO TO 30  
          END IF  
 *+       END WHILE  
 *  
       END IF  
 *  
       BETA = LBETA  
       T = LT  
       RND = LRND  
       IEEE1 = LIEEE1  
       FIRST = .FALSE.  
       RETURN  
 *  
 *     End of DLAMC1  
 *  
       END  
 *  
 ************************************************************************  
 *  
       SUBROUTINE DLAMC2( BETA, T, RND, EPS, EMIN, RMIN, EMAX, RMAX )  
 *  
 *  -- LAPACK auxiliary routine (version 3.2) --  
 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..  
 *     November 2006  
 *  
 *     .. Scalar Arguments ..  
       LOGICAL            RND  
       INTEGER            BETA, EMAX, EMIN, T  
       DOUBLE PRECISION   EPS, RMAX, RMIN  
 *     ..  
 *  
 *  Purpose  
 *  =======  
 *  
 *  DLAMC2 determines the machine parameters specified in its argument  
 *  list.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  BETA    (output) INTEGER  
 *          The base of the machine.  
 *  
 *  T       (output) INTEGER  
 *          The number of ( BETA ) digits in the mantissa.  
 *  
 *  RND     (output) LOGICAL  
 *          Specifies whether proper rounding  ( RND = .TRUE. )  or  
 *          chopping  ( RND = .FALSE. )  occurs in addition. This may not  
 *          be a reliable guide to the way in which the machine performs  
 *          its arithmetic.  
 *  
 *  EPS     (output) DOUBLE PRECISION  
 *          The smallest positive number such that  
 *  
 *             fl( 1.0 - EPS ) .LT. 1.0,  
 *  
 *          where fl denotes the computed value.  
 *  
 *  EMIN    (output) INTEGER  
 *          The minimum exponent before (gradual) underflow occurs.  
 *  
 *  RMIN    (output) DOUBLE PRECISION  
 *          The smallest normalized number for the machine, given by  
 *          BASE**( EMIN - 1 ), where  BASE  is the floating point value  
 *          of BETA.  
 *  
 *  EMAX    (output) INTEGER  
 *          The maximum exponent before overflow occurs.  
 *  
 *  RMAX    (output) DOUBLE PRECISION  
 *          The largest positive number for the machine, given by  
 *          BASE**EMAX * ( 1 - EPS ), where  BASE  is the floating point  
 *          value of BETA.  
 *  
 *  Further Details  
 *  ===============  
 *  
 *  The computation of  EPS  is based on a routine PARANOIA by  
 *  W. Kahan of the University of California at Berkeley.  
 *  
 * =====================================================================  
 *  
 *     .. Local Scalars ..  
       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  
 *     ..  
 *     .. External Functions ..  
       DOUBLE PRECISION   DLAMC3  
       EXTERNAL           DLAMC3  
 *     ..  
 *     .. External Subroutines ..  
       EXTERNAL           DLAMC1, DLAMC4, DLAMC5  
 *     ..  
 *     .. Intrinsic Functions ..  
       INTRINSIC          ABS, MAX, MIN  
 *     ..  
 *     .. Save statement ..  
       SAVE               FIRST, IWARN, LBETA, LEMAX, LEMIN, LEPS, LRMAX,  
      $                   LRMIN, LT  
 *     ..  
 *     .. Data statements ..  
       DATA               FIRST / .TRUE. / , IWARN / .FALSE. /  
 *     ..  
 *     .. Executable Statements ..  
 *  
       IF( FIRST ) THEN  
          ZERO = 0  
          ONE = 1  
          TWO = 2  
 *  
 *        LBETA, LT, LRND, LEPS, LEMIN and LRMIN  are the local values of  
 *        BETA, T, RND, EPS, EMIN and RMIN.  
 *  
 *        Throughout this routine  we use the function  DLAMC3  to ensure  
 *        that relevant values are stored  and not held in registers,  or  
 *        are not affected by optimizers.  
 *  
 *        DLAMC1 returns the parameters  LBETA, LT, LRND and LIEEE1.  
 *  
          CALL DLAMC1( LBETA, LT, LRND, LIEEE1 )  
 *  
 *        Start to find EPS.  
 *  
          B = LBETA  
          A = B**( -LT )  
          LEPS = A  
 *  
 *        Try some tricks to see whether or not this is the correct  EPS.  
 *  
          B = TWO / 3  
          HALF = ONE / 2  
          SIXTH = DLAMC3( B, -HALF )  
          THIRD = DLAMC3( SIXTH, SIXTH )  
          B = DLAMC3( THIRD, -HALF )  
          B = DLAMC3( B, SIXTH )  
          B = ABS( B )  
          IF( B.LT.LEPS )  
      $      B = LEPS  
 *  
          LEPS = 1  
 *  
 *+       WHILE( ( LEPS.GT.B ).AND.( B.GT.ZERO ) )LOOP  
    10    CONTINUE  
          IF( ( LEPS.GT.B ) .AND. ( B.GT.ZERO ) ) THEN  
             LEPS = B  
             C = DLAMC3( HALF*LEPS, ( TWO**5 )*( LEPS**2 ) )  
             C = DLAMC3( HALF, -C )  
             B = DLAMC3( HALF, C )  
             C = DLAMC3( HALF, -B )  
             B = DLAMC3( HALF, C )  
             GO TO 10  
          END IF  
 *+       END WHILE  
 *  
          IF( A.LT.LEPS )  
      $      LEPS = A  
 *  
 *        Computation of EPS complete.  
 *  
 *        Now find  EMIN.  Let A = + or - 1, and + or - (1 + BASE**(-3)).  
 *        Keep dividing  A by BETA until (gradual) underflow occurs. This  
 *        is detected when we cannot recover the previous A.  
 *  
          RBASE = ONE / LBETA  
          SMALL = ONE  
          DO 20 I = 1, 3  
             SMALL = DLAMC3( SMALL*RBASE, ZERO )  
    20    CONTINUE  
          A = DLAMC3( ONE, SMALL )  
          CALL DLAMC4( NGPMIN, ONE, LBETA )  
          CALL DLAMC4( NGNMIN, -ONE, LBETA )  
          CALL DLAMC4( GPMIN, A, LBETA )  
          CALL DLAMC4( GNMIN, -A, LBETA )  
          IEEE = .FALSE.  
 *  
          IF( ( NGPMIN.EQ.NGNMIN ) .AND. ( GPMIN.EQ.GNMIN ) ) THEN  
             IF( NGPMIN.EQ.GPMIN ) THEN  
                LEMIN = NGPMIN  
 *            ( Non twos-complement machines, no gradual underflow;  
 *              e.g.,  VAX )  
             ELSE IF( ( GPMIN-NGPMIN ).EQ.3 ) THEN  
                LEMIN = NGPMIN - 1 + LT  
                IEEE = .TRUE.  
 *            ( Non twos-complement machines, with gradual underflow;  
 *              e.g., IEEE standard followers )  
             ELSE  
                LEMIN = MIN( NGPMIN, GPMIN )  
 *            ( A guess; no known machine )  
                IWARN = .TRUE.  
             END IF  
 *  
          ELSE IF( ( NGPMIN.EQ.GPMIN ) .AND. ( NGNMIN.EQ.GNMIN ) ) THEN  
             IF( ABS( NGPMIN-NGNMIN ).EQ.1 ) THEN  
                LEMIN = MAX( NGPMIN, NGNMIN )  
 *            ( Twos-complement machines, no gradual underflow;  
 *              e.g., CYBER 205 )  
             ELSE  
                LEMIN = MIN( NGPMIN, NGNMIN )  
 *            ( A guess; no known machine )  
                IWARN = .TRUE.  
             END IF  
 *  
          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  
 *            ( Twos-complement machines with gradual underflow;  
 *              no known machine )  
             ELSE  
                LEMIN = MIN( NGPMIN, NGNMIN )  
 *            ( A guess; no known machine )  
                IWARN = .TRUE.  
             END IF  
 *  
          ELSE  
             LEMIN = MIN( NGPMIN, NGNMIN, GPMIN, GNMIN )  
 *         ( A guess; no known machine )  
             IWARN = .TRUE.  
          END IF  
          FIRST = .FALSE.  
 ***  
 * Comment out this if block if EMIN is ok  
          IF( IWARN ) THEN  
             FIRST = .TRUE.  
             WRITE( 6, FMT = 9999 )LEMIN  
          END IF  
 ***  
 *  
 *        Assume IEEE arithmetic if we found denormalised  numbers above,  
 *        or if arithmetic seems to round in the  IEEE style,  determined  
 *        in routine DLAMC1. A true IEEE machine should have both  things  
 *        true; however, faulty machines may have one or the other.  
 *  
          IEEE = IEEE .OR. LIEEE1  
 *  
 *        Compute  RMIN by successive division by  BETA. We could compute  
 *        RMIN as BASE**( EMIN - 1 ),  but some machines underflow during  
 *        this computation.  
 *  
          LRMIN = 1  
          DO 30 I = 1, 1 - LEMIN  
             LRMIN = DLAMC3( LRMIN*RBASE, ZERO )  
    30    CONTINUE  
 *  
 *        Finally, call DLAMC5 to compute EMAX and RMAX.  
 *  
          CALL DLAMC5( LBETA, LT, LEMIN, IEEE, LEMAX, LRMAX )  
       END IF  
 *  
       BETA = LBETA  
       T = LT  
       RND = LRND  
       EPS = LEPS  
       EMIN = LEMIN  
       RMIN = LRMIN  
       EMAX = LEMAX  
       RMAX = LRMAX  
 *  
       RETURN  
 *  
  9999 FORMAT( / / ' WARNING. The value EMIN may be incorrect:-',  
      $      '  EMIN = ', I8, /  
      $      ' If, after inspection, the value EMIN looks',  
      $      ' acceptable please comment out ',  
      $      / ' the IF block as marked within the code of routine',  
      $      ' DLAMC2,', / ' otherwise supply EMIN explicitly.', / )  
 *  
 *     End of DLAMC2  
 *  
       END  
 *  
 ************************************************************************  ************************************************************************
 *  *
       DOUBLE PRECISION FUNCTION DLAMC3( A, B )        DOUBLE PRECISION FUNCTION DLAMC3( A, B )
 *  *
 *  -- LAPACK auxiliary routine (version 3.2) --  *  -- LAPACK auxiliary routine (version 3.3.0) --
 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..  *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
 *     November 2006  *     November 2010
 *  *
 *     .. Scalar Arguments ..  *     .. Scalar Arguments ..
       DOUBLE PRECISION   A, B        DOUBLE PRECISION   A, B
Line 608 Line 154
       END        END
 *  *
 ************************************************************************  ************************************************************************
 *  
       SUBROUTINE DLAMC4( EMIN, START, BASE )  
 *  
 *  -- LAPACK auxiliary routine (version 3.2) --  
 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..  
 *     November 2006  
 *  
 *     .. Scalar Arguments ..  
       INTEGER            BASE, EMIN  
       DOUBLE PRECISION   START  
 *     ..  
 *  
 *  Purpose  
 *  =======  
 *  
 *  DLAMC4 is a service routine for DLAMC2.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  EMIN    (output) INTEGER   
 *          The minimum exponent before (gradual) underflow, computed by  
 *          setting A = START and dividing by BASE until the previous A  
 *          can not be recovered.  
 *  
 *  START   (input) DOUBLE PRECISION  
 *          The starting point for determining EMIN.  
 *  
 *  BASE    (input) INTEGER  
 *          The base of the machine.  
 *  
 * =====================================================================  
 *  
 *     .. Local Scalars ..  
       INTEGER            I  
       DOUBLE PRECISION   A, B1, B2, C1, C2, D1, D2, ONE, RBASE, ZERO  
 *     ..  
 *     .. External Functions ..  
       DOUBLE PRECISION   DLAMC3  
       EXTERNAL           DLAMC3  
 *     ..  
 *     .. Executable Statements ..  
 *  
       A = START  
       ONE = 1  
       RBASE = ONE / BASE  
       ZERO = 0  
       EMIN = 1  
       B1 = DLAMC3( A*RBASE, ZERO )  
       C1 = A  
       C2 = A  
       D1 = A  
       D2 = A  
 *+    WHILE( ( C1.EQ.A ).AND.( C2.EQ.A ).AND.  
 *    $       ( D1.EQ.A ).AND.( D2.EQ.A )      )LOOP  
    10 CONTINUE  
       IF( ( C1.EQ.A ) .AND. ( C2.EQ.A ) .AND. ( D1.EQ.A ) .AND.  
      $    ( D2.EQ.A ) ) THEN  
          EMIN = EMIN - 1  
          A = B1  
          B1 = DLAMC3( A / BASE, ZERO )  
          C1 = DLAMC3( B1*BASE, ZERO )  
          D1 = ZERO  
          DO 20 I = 1, BASE  
             D1 = D1 + B1  
    20    CONTINUE  
          B2 = DLAMC3( A*RBASE, ZERO )  
          C2 = DLAMC3( B2 / RBASE, ZERO )  
          D2 = ZERO  
          DO 30 I = 1, BASE  
             D2 = D2 + B2  
    30    CONTINUE  
          GO TO 10  
       END IF  
 *+    END WHILE  
 *  
       RETURN  
 *  
 *     End of DLAMC4  
 *  
       END  
 *  
 ************************************************************************  
 *  
       SUBROUTINE DLAMC5( BETA, P, EMIN, IEEE, EMAX, RMAX )  
 *  
 *  -- LAPACK auxiliary routine (version 3.2) --  
 *     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..  
 *     November 2006  
 *  
 *     .. Scalar Arguments ..  
       LOGICAL            IEEE  
       INTEGER            BETA, EMAX, EMIN, P  
       DOUBLE PRECISION   RMAX  
 *     ..  
 *  
 *  Purpose  
 *  =======  
 *  
 *  DLAMC5 attempts to compute RMAX, the largest machine floating-point  
 *  number, without overflow.  It assumes that EMAX + abs(EMIN) sum  
 *  approximately to a power of 2.  It will fail on machines where this  
 *  assumption does not hold, for example, the Cyber 205 (EMIN = -28625,  
 *  EMAX = 28718).  It will also fail if the value supplied for EMIN is  
 *  too large (i.e. too close to zero), probably with overflow.  
 *  
 *  Arguments  
 *  =========  
 *  
 *  BETA    (input) INTEGER  
 *          The base of floating-point arithmetic.  
 *  
 *  P       (input) INTEGER  
 *          The number of base BETA digits in the mantissa of a  
 *          floating-point value.  
 *  
 *  EMIN    (input) INTEGER  
 *          The minimum exponent before (gradual) underflow.  
 *  
 *  IEEE    (input) LOGICAL  
 *          A logical flag specifying whether or not the arithmetic  
 *          system is thought to comply with the IEEE standard.  
 *  
 *  EMAX    (output) INTEGER  
 *          The largest exponent before overflow  
 *  
 *  RMAX    (output) DOUBLE PRECISION  
 *          The largest machine floating-point number.  
 *  
 * =====================================================================  
 *  
 *     .. Parameters ..  
       DOUBLE PRECISION   ZERO, ONE  
       PARAMETER          ( ZERO = 0.0D0, ONE = 1.0D0 )  
 *     ..  
 *     .. Local Scalars ..  
       INTEGER            EXBITS, EXPSUM, I, LEXP, NBITS, TRY, UEXP  
       DOUBLE PRECISION   OLDY, RECBAS, Y, Z  
 *     ..  
 *     .. External Functions ..  
       DOUBLE PRECISION   DLAMC3  
       EXTERNAL           DLAMC3  
 *     ..  
 *     .. Intrinsic Functions ..  
       INTRINSIC          MOD  
 *     ..  
 *     .. Executable Statements ..  
 *  
 *     First compute LEXP and UEXP, two powers of 2 that bound  
 *     abs(EMIN). We then assume that EMAX + abs(EMIN) will sum  
 *     approximately to the bound that is closest to abs(EMIN).  
 *     (EMAX is the exponent of the required number RMAX).  
 *  
       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  
 *  
 *     Now -LEXP is less than or equal to EMIN, and -UEXP is greater  
 *     than or equal to EMIN. EXBITS is the number of bits needed to  
 *     store the exponent.  
 *  
       IF( ( UEXP+EMIN ).GT.( -LEXP-EMIN ) ) THEN  
          EXPSUM = 2*LEXP  
       ELSE  
          EXPSUM = 2*UEXP  
       END IF  
 *  
 *     EXPSUM is the exponent range, approximately equal to  
 *     EMAX - EMIN + 1 .  
 *  
       EMAX = EXPSUM + EMIN - 1  
       NBITS = 1 + EXBITS + P  
 *  
 *     NBITS is the total number of bits needed to store a  
 *     floating-point number.  
 *  
       IF( ( MOD( NBITS, 2 ).EQ.1 ) .AND. ( BETA.EQ.2 ) ) THEN  
 *  
 *        Either there are an odd number of bits used to store a  
 *        floating-point number, which is unlikely, or some bits are  
 *        not used in the representation of numbers, which is possible,  
 *        (e.g. Cray machines) or the mantissa has an implicit bit,  
 *        (e.g. IEEE machines, Dec Vax machines), which is perhaps the  
 *        most likely. We have to assume the last alternative.  
 *        If this is true, then we need to reduce EMAX by one because  
 *        there must be some way of representing zero in an implicit-bit  
 *        system. On machines like Cray, we are reducing EMAX by one  
 *        unnecessarily.  
 *  
          EMAX = EMAX - 1  
       END IF  
 *  
       IF( IEEE ) THEN  
 *  
 *        Assume we are on an IEEE machine which reserves one exponent  
 *        for infinity and NaN.  
 *  
          EMAX = EMAX - 1  
       END IF  
 *  
 *     Now create RMAX, the largest machine number, which should  
 *     be equal to (1.0 - BETA**(-P)) * BETA**EMAX .  
 *  
 *     First compute 1.0 - BETA**(-P), being careful that the  
 *     result is less than 1.0 .  
 *  
       RECBAS = ONE / BETA  
       Z = BETA - ONE  
       Y = ZERO  
       DO 20 I = 1, P  
          Z = Z*RECBAS  
          IF( Y.LT.ONE )  
      $      OLDY = Y  
          Y = DLAMC3( Y, Z )  
    20 CONTINUE  
       IF( Y.GE.ONE )  
      $   Y = OLDY  
 *  
 *     Now multiply by BETA**EMAX to get RMAX.  
 *  
       DO 30 I = 1, EMAX  
          Y = DLAMC3( Y*BETA, ZERO )  
    30 CONTINUE  
 *  
       RMAX = Y  
       RETURN  
 *  
 *     End of DLAMC5  
 *  
       END  
Removed from v.1.6  
changed lines
  Added in v.1.7

CVSweb interface <joel.bertrand@systella.fr>