Chimie Théorique » Opt'n Path

      DOUBLE PRECISION FUNCTION DLAMCH( CMACH )

*

*  -- LAPACK auxiliary routine (version 3.2) --

*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..

*     November 2006

*

*     .. Scalar Arguments ..

      CHARACTER          CMACH

*     ..

*

*  Purpose

*  =======

*

*  DLAMCH determines double precision machine parameters.

*

*  Arguments

*  =========

*

*  CMACH   (input) CHARACTER*1

*          Specifies the value to be returned by DLAMCH:

*          = 'E' or 'e',   DLAMCH := eps

*          = 'S' or 's ,   DLAMCH := sfmin

*          = 'B' or 'b',   DLAMCH := base

*          = 'P' or 'p',   DLAMCH := eps*base

*          = 'N' or 'n',   DLAMCH := t

*          = 'R' or 'r',   DLAMCH := rnd

*          = 'M' or 'm',   DLAMCH := emin

*          = 'U' or 'u',   DLAMCH := rmin

*          = 'L' or 'l',   DLAMCH := emax

*          = 'O' or 'o',   DLAMCH := rmax

*

*          where

*

*          eps   = relative machine precision

*          sfmin = safe minimum, such that 1/sfmin does not overflow

*          base  = base of the machine

*          prec  = eps*base

*          t     = number of (base) digits in the mantissa

*          rnd   = 1.0 when rounding occurs in addition, 0.0 otherwise

*          emin  = minimum exponent before (gradual) underflow

*          rmin  = underflow threshold - base**(emin-1)

*          emax  = largest exponent before overflow

*          rmax  = overflow threshold  - (base**emax)*(1-eps)

*

* =====================================================================

*

*     .. Parameters ..

      DOUBLE PRECISION   ONE, ZERO

      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )

*     ..

*     .. Local Scalars ..

      LOGICAL            FIRST, LRND

      INTEGER            BETA, IMAX, IMIN, IT

      DOUBLE PRECISION   BASE, EMAX, EMIN, EPS, PREC, RMACH, RMAX, RMIN,

     $                   RND, SFMIN, SMALL, T

*     ..

*     .. External Functions ..

      LOGICAL            LSAME

      EXTERNAL           LSAME

*     ..

*     .. External Subroutines ..

      EXTERNAL           DLAMC2

*     ..

*     .. Save statement ..

      SAVE               FIRST, EPS, SFMIN, BASE, T, RND, EMIN, RMIN,

     $                   EMAX, RMAX, PREC

*     ..

*     .. Data statements ..

      DATA               FIRST / .TRUE. /

*     ..

*     .. 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

*           causing overflow when computing  1/sfmin.

*

            SFMIN = SMALL*( ONE+EPS )

         END IF

      END IF

*

      IF( LSAME( CMACH, 'E' ) ) THEN

         RMACH = EPS

      ELSE IF( LSAME( CMACH, 'S' ) ) THEN

         RMACH = SFMIN

      ELSE IF( LSAME( CMACH, 'B' ) ) THEN

         RMACH = BASE

      ELSE IF( LSAME( CMACH, 'P' ) ) THEN

         RMACH = PREC

      ELSE IF( LSAME( CMACH, 'N' ) ) THEN

         RMACH = T

      ELSE IF( LSAME( CMACH, 'R' ) ) THEN

         RMACH = RND

      ELSE IF( LSAME( CMACH, 'M' ) ) THEN

         RMACH = EMIN

      ELSE IF( LSAME( CMACH, 'U' ) ) THEN

         RMACH = RMIN

      ELSE IF( LSAME( CMACH, 'L' ) ) THEN

         RMACH = EMAX

      ELSE IF( LSAME( CMACH, 'O' ) ) THEN

         RMACH = RMAX

      END IF

*

      DLAMCH = RMACH

      FIRST  = .FALSE.

      RETURN

*

*     End of DLAMCH

*

      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 )

*

*  -- LAPACK auxiliary routine (version 3.2) --

*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..

*     November 2006

*

*     .. Scalar Arguments ..

      DOUBLE PRECISION   A, B

*     ..

*

*  Purpose

*  =======

*

*  DLAMC3  is intended to force  A  and  B  to be stored prior to doing

*  the addition of  A  and  B ,  for use in situations where optimizers

*  might hold one of these in a register.

*

*  Arguments

*  =========

*

*  A       (input) DOUBLE PRECISION

*  B       (input) DOUBLE PRECISION

*          The values A and B.

*

* =====================================================================

*

*     .. Executable Statements ..

*

      DLAMC3 = A + B

*

      RETURN

*

*     End of DLAMC3

*

      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