version 1.4, 2010/08/06 15:32:27
|
version 1.10, 2011/11/21 20:42:55
|
Line 1
|
Line 1
|
|
*> \brief \b DLAMCH |
|
* |
|
* =========== DOCUMENTATION =========== |
|
* |
|
* Online html documentation available at |
|
* http://www.netlib.org/lapack/explore-html/ |
|
* |
|
* Definition: |
|
* =========== |
|
* |
|
* DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) |
|
* |
|
* |
|
*> \par Purpose: |
|
* ============= |
|
*> |
|
*> \verbatim |
|
*> |
|
*> DLAMCH determines double precision machine parameters. |
|
*> \endverbatim |
|
* |
|
* Arguments: |
|
* ========== |
|
* |
|
*> \param[in] CMACH |
|
*> \verbatim |
|
*> 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) |
|
*> \endverbatim |
|
* |
|
* Authors: |
|
* ======== |
|
* |
|
*> \author Univ. of Tennessee |
|
*> \author Univ. of California Berkeley |
|
*> \author Univ. of Colorado Denver |
|
*> \author NAG Ltd. |
|
* |
|
*> \date November 2011 |
|
* |
|
*> \ingroup auxOTHERauxiliary |
|
* |
|
* ===================================================================== |
DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) |
DOUBLE PRECISION FUNCTION DLAMCH( CMACH ) |
* |
* |
* -- LAPACK auxiliary routine (version 3.2) -- |
* -- LAPACK auxiliary routine (version 3.4.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..-- |
|
* November 2011 |
* |
* |
* .. Scalar Arguments .. |
* .. Scalar Arguments .. |
CHARACTER CMACH |
CHARACTER CMACH |
* .. |
* .. |
* |
* |
* Purpose |
* .. Scalar Arguments .. |
* ======= |
DOUBLE PRECISION A, B |
* |
* .. |
* 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) |
|
* |
* |
* ===================================================================== |
* ===================================================================== |
* |
* |
Line 49
|
Line 83
|
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 |
|
* |
|
************************************************************************ |
************************************************************************ |
* |
*> \brief \b DLAMC3 |
|
*> \details |
|
*> \b Purpose: |
|
*> \verbatim |
|
*> 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. |
|
*> \endverbatim |
|
*> \author LAPACK is a software package provided by Univ. of Tennessee, Univ. of California Berkeley, Univ. of Colorado Denver and NAG Ltd.. |
|
*> \date November 2011 |
|
*> \ingroup auxOTHERauxiliary |
|
*> |
|
*> \param[in] A |
|
*> \verbatim |
|
*> A is a DOUBLE PRECISION |
|
*> \endverbatim |
|
*> |
|
*> \param[in] B |
|
*> \verbatim |
|
*> B is a DOUBLE PRECISION |
|
*> The values A and B. |
|
*> \endverbatim |
|
*> |
DOUBLE PRECISION FUNCTION DLAMC3( A, B ) |
DOUBLE PRECISION FUNCTION DLAMC3( A, B ) |
* |
* |
* -- LAPACK auxiliary routine (version 3.2) -- |
* -- LAPACK auxiliary routine (version 3.4.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 |
* .. |
* .. |
* |
|
* 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 .. |
* .. Executable Statements .. |
Line 608
|
Line 191
|
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 |
|