[430] | 1 | !*********************************************************************** |
---|
| 2 | ! |
---|
| 3 | SUBROUTINE DLAMC1( BETA, T, RND, IEEE1 ) |
---|
| 4 | implicit none |
---|
| 5 | ! |
---|
| 6 | ! -- LAPACK auxiliary routine (version 3.1) -- |
---|
| 7 | ! Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd.. |
---|
| 8 | ! November 2006 |
---|
| 9 | ! |
---|
| 10 | ! .. Scalar Arguments .. |
---|
| 11 | LOGICAL IEEE1, RND |
---|
| 12 | INTEGER BETA, T |
---|
| 13 | ! .. |
---|
| 14 | ! |
---|
| 15 | ! Purpose |
---|
| 16 | ! ======= |
---|
| 17 | ! |
---|
| 18 | ! DLAMC1 determines the machine parameters given by BETA, T, RND, and |
---|
| 19 | ! IEEE1. |
---|
| 20 | ! |
---|
| 21 | ! Arguments |
---|
| 22 | ! ========= |
---|
| 23 | ! |
---|
| 24 | ! BETA (output) INTEGER |
---|
| 25 | ! The base of the machine. |
---|
| 26 | ! |
---|
| 27 | ! T (output) INTEGER |
---|
| 28 | ! The number of ( BETA ) digits in the mantissa. |
---|
| 29 | ! |
---|
| 30 | ! RND (output) LOGICAL |
---|
| 31 | ! Specifies whether proper rounding ( RND = .TRUE. ) or |
---|
| 32 | ! chopping ( RND = .FALSE. ) occurs in addition. This may not |
---|
| 33 | ! be a reliable guide to the way in which the machine performs |
---|
| 34 | ! its arithmetic. |
---|
| 35 | ! |
---|
| 36 | ! IEEE1 (output) LOGICAL |
---|
| 37 | ! Specifies whether rounding appears to be done in the IEEE |
---|
| 38 | ! 'round to nearest' style. |
---|
| 39 | ! |
---|
| 40 | ! Further Details |
---|
| 41 | ! =============== |
---|
| 42 | ! |
---|
| 43 | ! The routine is based on the routine ENVRON by Malcolm and |
---|
| 44 | ! incorporates suggestions by Gentleman and Marovich. See |
---|
| 45 | ! |
---|
| 46 | ! Malcolm M. A. (1972) Algorithms to reveal properties of |
---|
| 47 | ! floating-point arithmetic. Comms. of the ACM, 15, 949-951. |
---|
| 48 | ! |
---|
| 49 | ! Gentleman W. M. and Marovich S. B. (1974) More on algorithms |
---|
| 50 | ! that reveal properties of floating point arithmetic units. |
---|
| 51 | ! Comms. of the ACM, 17, 276-277. |
---|
| 52 | ! |
---|
| 53 | ! ===================================================================== |
---|
| 54 | ! |
---|
| 55 | ! .. Local Scalars .. |
---|
| 56 | LOGICAL FIRST, LIEEE1, LRND |
---|
| 57 | INTEGER LBETA, LT |
---|
| 58 | DOUBLE PRECISION A, B, C, F, ONE, QTR, SAVEC, T1, T2 |
---|
| 59 | ! .. |
---|
| 60 | ! .. External Functions .. |
---|
| 61 | DOUBLE PRECISION DLAMC3 |
---|
| 62 | EXTERNAL DLAMC3 |
---|
| 63 | ! .. |
---|
| 64 | ! .. Save statement .. |
---|
| 65 | SAVE FIRST, LIEEE1, LBETA, LRND, LT |
---|
| 66 | ! .. |
---|
| 67 | ! .. Data statements .. |
---|
| 68 | DATA FIRST / .TRUE. / |
---|
| 69 | ! .. |
---|
| 70 | ! .. Executable Statements .. |
---|
| 71 | ! |
---|
| 72 | IF( FIRST ) THEN |
---|
| 73 | ONE = 1 |
---|
| 74 | ! |
---|
| 75 | ! LBETA, LIEEE1, LT and LRND are the local values of BETA, |
---|
| 76 | ! IEEE1, T and RND. |
---|
| 77 | ! |
---|
| 78 | ! Throughout this routine we use the function DLAMC3 to ensure |
---|
| 79 | ! that relevant values are stored and not held in registers, or |
---|
| 80 | ! are not affected by optimizers. |
---|
| 81 | ! |
---|
| 82 | ! Compute a = 2.0**m with the smallest positive integer m such |
---|
| 83 | ! that |
---|
| 84 | ! |
---|
| 85 | ! fl( a + 1.0 ) = a. |
---|
| 86 | ! |
---|
| 87 | A = 1 |
---|
| 88 | C = 1 |
---|
| 89 | ! |
---|
| 90 | !+ WHILE( C.EQ.ONE )LOOP |
---|
| 91 | 10 CONTINUE |
---|
| 92 | IF( C.EQ.ONE ) THEN |
---|
| 93 | A = 2*A |
---|
| 94 | C = DLAMC3( A, ONE ) |
---|
| 95 | C = DLAMC3( C, -A ) |
---|
| 96 | GO TO 10 |
---|
| 97 | END IF |
---|
| 98 | !+ END WHILE |
---|
| 99 | ! |
---|
| 100 | ! Now compute b = 2.0**m with the smallest positive integer m |
---|
| 101 | ! such that |
---|
| 102 | ! |
---|
| 103 | ! fl( a + b ) .gt. a. |
---|
| 104 | ! |
---|
| 105 | B = 1 |
---|
| 106 | C = DLAMC3( A, B ) |
---|
| 107 | ! |
---|
| 108 | !+ WHILE( C.EQ.A )LOOP |
---|
| 109 | 20 CONTINUE |
---|
| 110 | IF( C.EQ.A ) THEN |
---|
| 111 | B = 2*B |
---|
| 112 | C = DLAMC3( A, B ) |
---|
| 113 | GO TO 20 |
---|
| 114 | END IF |
---|
| 115 | !+ END WHILE |
---|
| 116 | ! |
---|
| 117 | ! Now compute the base. a and c are neighbouring floating point |
---|
| 118 | ! numbers in the interval ( beta**t, beta**( t + 1 ) ) and so |
---|
| 119 | ! their difference is beta. Adding 0.25 to c is to ensure that it |
---|
| 120 | ! is truncated to beta and not ( beta - 1 ). |
---|
| 121 | ! |
---|
| 122 | QTR = ONE / 4 |
---|
| 123 | SAVEC = C |
---|
| 124 | C = DLAMC3( C, -A ) |
---|
| 125 | LBETA = C + QTR |
---|
| 126 | ! |
---|
| 127 | ! Now determine whether rounding or chopping occurs, by adding a |
---|
| 128 | ! bit less than beta/2 and a bit more than beta/2 to a. |
---|
| 129 | ! |
---|
| 130 | B = LBETA |
---|
| 131 | F = DLAMC3( B / 2, -B / 100 ) |
---|
| 132 | C = DLAMC3( F, A ) |
---|
| 133 | IF( C.EQ.A ) THEN |
---|
| 134 | LRND = .TRUE. |
---|
| 135 | ELSE |
---|
| 136 | LRND = .FALSE. |
---|
| 137 | END IF |
---|
| 138 | F = DLAMC3( B / 2, B / 100 ) |
---|
| 139 | C = DLAMC3( F, A ) |
---|
| 140 | IF( ( LRND ) .AND. ( C.EQ.A ) ) & |
---|
| 141 | &LRND = .FALSE. |
---|
| 142 | ! |
---|
| 143 | ! Try and decide whether rounding is done in the IEEE 'round to |
---|
| 144 | ! nearest' style. B/2 is half a unit in the last place of the two |
---|
| 145 | ! numbers A and SAVEC. Furthermore, A is even, i.e. has last bit |
---|
| 146 | ! zero, and SAVEC is odd. Thus adding B/2 to A should not change |
---|
| 147 | ! A, but adding B/2 to SAVEC should change SAVEC. |
---|
| 148 | ! |
---|
| 149 | T1 = DLAMC3( B / 2, A ) |
---|
| 150 | T2 = DLAMC3( B / 2, SAVEC ) |
---|
| 151 | LIEEE1 = ( T1.EQ.A ) .AND. ( T2.GT.SAVEC ) .AND. LRND |
---|
| 152 | ! |
---|
| 153 | ! Now find the mantissa, t. It should be the integer part of |
---|
| 154 | ! log to the base beta of a, however it is safer to determine t |
---|
| 155 | ! by powering. So we find t as the smallest positive integer for |
---|
| 156 | ! which |
---|
| 157 | ! |
---|
| 158 | ! fl( beta**t + 1.0 ) = 1.0. |
---|
| 159 | ! |
---|
| 160 | LT = 0 |
---|
| 161 | A = 1 |
---|
| 162 | C = 1 |
---|
| 163 | ! |
---|
| 164 | !+ WHILE( C.EQ.ONE )LOOP |
---|
| 165 | 30 CONTINUE |
---|
| 166 | IF( C.EQ.ONE ) THEN |
---|
| 167 | LT = LT + 1 |
---|
| 168 | A = A*LBETA |
---|
| 169 | C = DLAMC3( A, ONE ) |
---|
| 170 | C = DLAMC3( C, -A ) |
---|
| 171 | GO TO 30 |
---|
| 172 | END IF |
---|
| 173 | !+ END WHILE |
---|
| 174 | ! |
---|
| 175 | END IF |
---|
| 176 | ! |
---|
| 177 | BETA = LBETA |
---|
| 178 | T = LT |
---|
| 179 | RND = LRND |
---|
| 180 | IEEE1 = LIEEE1 |
---|
| 181 | FIRST = .FALSE. |
---|
| 182 | RETURN |
---|
| 183 | ! |
---|
| 184 | ! End of DLAMC1 |
---|
| 185 | ! |
---|
| 186 | END |
---|
| 187 | ! |
---|