SCALAPACK 2.2.2
LAPACK: Linear Algebra PACKage
Loading...
Searching...
No Matches

◆ slamc1()

subroutine slamc1 ( integer  beta,
integer  t,
logical  rnd,
logical  ieee1 
)

Definition at line 996 of file tools.f.

997*
998* -- LAPACK auxiliary routine (version 2.0) --
999* Univ. of Tennessee, Univ. of California Berkeley, NAG Ltd.,
1000* Courant Institute, Argonne National Lab, and Rice University
1001* October 31, 1992
1002*
1003* .. Scalar Arguments ..
1004 LOGICAL IEEE1, RND
1005 INTEGER BETA, T
1006* ..
1007*
1008* Purpose
1009* =======
1010*
1011* SLAMC1 determines the machine parameters given by BETA, T, RND, and
1012* IEEE1.
1013*
1014* Arguments
1015* =========
1016*
1017* BETA (output) INTEGER
1018* The base of the machine.
1019*
1020* T (output) INTEGER
1021* The number of ( BETA ) digits in the mantissa.
1022*
1023* RND (output) LOGICAL
1024* Specifies whether proper rounding ( RND = .TRUE. ) or
1025* chopping ( RND = .FALSE. ) occurs in addition. This may not
1026* be a reliable guide to the way in which the machine performs
1027* its arithmetic.
1028*
1029* IEEE1 (output) LOGICAL
1030* Specifies whether rounding appears to be done in the IEEE
1031* 'round to nearest' style.
1032*
1033* Further Details
1034* ===============
1035*
1036* The routine is based on the routine ENVRON by Malcolm and
1037* incorporates suggestions by Gentleman and Marovich. See
1038*
1039* Malcolm M. A. (1972) Algorithms to reveal properties of
1040* floating-point arithmetic. Comms. of the ACM, 15, 949-951.
1041*
1042* Gentleman W. M. and Marovich S. B. (1974) More on algorithms
1043* that reveal properties of floating point arithmetic units.
1044* Comms. of the ACM, 17, 276-277.
1045*
1046* =====================================================================
1047*
1048* .. Local Scalars ..
1049 LOGICAL FIRST, LIEEE1, LRND
1050 INTEGER LBETA, LT
1051 REAL A, B, C, F, ONE, QTR, SAVEC, T1, T2
1052* ..
1053* .. External Functions ..
1054 REAL SLAMC3
1055 EXTERNAL slamc3
1056* ..
1057* .. Save statement ..
1058 SAVE first, lieee1, lbeta, lrnd, lt
1059* ..
1060* .. Data statements ..
1061 DATA first / .true. /
1062* ..
1063* .. Executable Statements ..
1064*
1065 IF( first ) THEN
1066 first = .false.
1067 one = 1
1068*
1069* LBETA, LIEEE1, LT and LRND are the local values of BETA,
1070* IEEE1, T and RND.
1071*
1072* Throughout this routine we use the function SLAMC3 to ensure
1073* that relevant values are stored and not held in registers, or
1074* are not affected by optimizers.
1075*
1076* Compute a = 2.0**m with the smallest positive integer m such
1077* that
1078*
1079* fl( a + 1.0 ) = a.
1080*
1081 a = 1
1082 c = 1
1083*
1084*+ WHILE( C.EQ.ONE )LOOP
1085 10 CONTINUE
1086 IF( c.EQ.one ) THEN
1087 a = 2*a
1088 c = slamc3( a, one )
1089 c = slamc3( c, -a )
1090 GO TO 10
1091 END IF
1092*+ END WHILE
1093*
1094* Now compute b = 2.0**m with the smallest positive integer m
1095* such that
1096*
1097* fl( a + b ) .gt. a.
1098*
1099 b = 1
1100 c = slamc3( a, b )
1101*
1102*+ WHILE( C.EQ.A )LOOP
1103 20 CONTINUE
1104 IF( c.EQ.a ) THEN
1105 b = 2*b
1106 c = slamc3( a, b )
1107 GO TO 20
1108 END IF
1109*+ END WHILE
1110*
1111* Now compute the base. a and c are neighbouring floating point
1112* numbers in the interval ( beta**t, beta**( t + 1 ) ) and so
1113* their difference is beta. Adding 0.25 to c is to ensure that it
1114* is truncated to beta and not ( beta - 1 ).
1115*
1116 qtr = one / 4
1117 savec = c
1118 c = slamc3( c, -a )
1119 lbeta = c + qtr
1120*
1121* Now determine whether rounding or chopping occurs, by adding a
1122* bit less than beta/2 and a bit more than beta/2 to a.
1123*
1124 b = lbeta
1125 f = slamc3( b / 2, -b / 100 )
1126 c = slamc3( f, a )
1127 IF( c.EQ.a ) THEN
1128 lrnd = .true.
1129 ELSE
1130 lrnd = .false.
1131 END IF
1132 f = slamc3( b / 2, b / 100 )
1133 c = slamc3( f, a )
1134 IF( ( lrnd ) .AND. ( c.EQ.a ) )
1135 $ lrnd = .false.
1136*
1137* Try and decide whether rounding is done in the IEEE 'round to
1138* nearest' style. B/2 is half a unit in the last place of the two
1139* numbers A and SAVEC. Furthermore, A is even, i.e. has last bit
1140* zero, and SAVEC is odd. Thus adding B/2 to A should not change
1141* A, but adding B/2 to SAVEC should change SAVEC.
1142*
1143 t1 = slamc3( b / 2, a )
1144 t2 = slamc3( b / 2, savec )
1145 lieee1 = ( t1.EQ.a ) .AND. ( t2.GT.savec ) .AND. lrnd
1146*
1147* Now find the mantissa, t. It should be the integer part of
1148* log to the base beta of a, however it is safer to determine t
1149* by powering. So we find t as the smallest positive integer for
1150* which
1151*
1152* fl( beta**t + 1.0 ) = 1.0.
1153*
1154 lt = 0
1155 a = 1
1156 c = 1
1157*
1158*+ WHILE( C.EQ.ONE )LOOP
1159 30 CONTINUE
1160 IF( c.EQ.one ) THEN
1161 lt = lt + 1
1162 a = a*lbeta
1163 c = slamc3( a, one )
1164 c = slamc3( c, -a )
1165 GO TO 30
1166 END IF
1167*+ END WHILE
1168*
1169 END IF
1170*
1171 beta = lbeta
1172 t = lt
1173 rnd = lrnd
1174 ieee1 = lieee1
1175 RETURN
1176*
1177* End of SLAMC1
1178*
real function slamc3(a, b)
Definition tools.f:1443
Here is the caller graph for this function: