Diffusion lammps (efficient_neuronet) Browse

dlamch.f
753367408642efficient_neuronet

dlamch.f
View Options

			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 ( betat, 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( HALFLEPS, ( TWO5 )( 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

lammps/lib/linalg/dlamch.f753367408642efficient_neuronet

dlamch.fView Options

lammps/lib/linalg/dlamch.f
753367408642efficient_neuronet

dlamch.f
View Options