Path: utzoo!utgpu!news-server.csri.toronto.edu!rpi!zaphod.mps.ohio-state.edu!casbah.acns.nwu.edu!ftpbox!mothost!motcid!reilly From: reilly@motcid.UUCP (Patrick L. Reilly) Newsgroups: comp.lang.fortran Subject: Re: validating fortran compiler's numerical results Keywords: fortran, validating Message-ID: <6195@gypsum.UUCP> Date: 12 Apr 91 19:24:30 GMT References: <1991Apr9.023737.18140@mlb.semi.harris.com> Reply-To: motcid!reillyp@uunet.uu.net Organization: Motorola Inc., Cellular Infrastructure Div., Arlington Heights, IL Lines: 176 jdr@sloth.mlb.semi.harris.com (Jim Ray) writes: :As the topic states, I am looking for some fortran program that one :could use to compare the numeric (floating point) results with some :"known" values? Specifically, I am looking for such a program that :would exercise (extensively and repeditively ) most of the :transcendental functions plus the normal run of the mill ones. This :"program" would checking for roundoff and complience with the IEEE :standard floating point specs. You need to get a copy of "Testing Unconstrained Optimization Software" by More, J. et. al, ACM Trans. on Mathematical Software, v. 7,n. 1, March 1981, pages 17-41. Although what you specifically seek is a bit different from the referenced paper, you will find much insight and routines to test with in the paper. Being up to date with Rosenbrock's, Freudenstein and Roth's, Beale, Helical valley, etc. functions can go a long way towards your problem. As for roundoff, you might try the i1mach, r1mach, and d1mach routines from the SLATEC library to get important parameters about your target machine. Alternatively, give some of the code below a whirl. Good luck, ========================== `` '' =================================== Patrick Reilly, Ph.D. < @ @ > tel: 1+708+632-3109 Motorola CIG/Switch Dev. ( > ) fax: 1+708+632-2413 1501 W. Shure Dr. -BC569 \~/ UUNET: uunet!motcid!reilly Arlington Hts., IL 60004 Any ill-formed premise may be REilly's CAutioNary Truth : modeled and validated with -- -- - - simulation. ==================================================================== Roundoff.f: integer i real h double x i = 0 x = 0.0 h = 0.1 10 continue i = i + 1 x = x + h print *, i, x if(x.le.1.0)go to 10 end Epsilon.f: implicit double precision(a-h,o-z) eps= 1.0 10 continue eps= eps/2.0 t= 1.0 + eps if(t.gt.1.0) go to 10 eps = 2.0*eps print *,eps,'.le.machine epsilon.le.',2.0 * eps end **************************************************************** First results from the CrayX.... Output from epsilon.f: 2.5243548967072377773175314089E-29.le.machine epsilon.le.5.0487097934144755546350628178E-29 (Epsilon.f output can be considered, roughly speaking, to be the fractional accuracy to which floating point numbers are represented and corresponds to a change of one in the least significant bit of the mantissa. That is, machine epsilon is a positive floating-point number for which fl(1+epsilon)>= 1, where fl(x) is the floating-point representation of a real number x. This code yields this value, epsilon, within a factor of two.) Output from roundoff.f: 1, 1.0000000000000008881784197001E-1 2, 2.0000000000000017763568394002E-1 3, 3.0000000000000026645352591003E-1 4, 4.0000000000000035527136788005E-1 5, 5.0000000000000044408920985006E-1 6, 6.0000000000000053290705182007E-1 7, 7.0000000000000062172489379008E-1 8, 8.000000000000007105427357601E-1 9, 9.0000000000000079936057773011E-1 10, 1.0000000000000008881784197001E+0 Compare the above with the following after changing x = x + h to read x = i * h in the roundoff.f program: 1, 1.0000000000000008881784197001E-1 2, 2.0000000000000017763568394002E-1 3, 3.000000000000007105427357601E-1 4, 4.0000000000000035527136788005E-1 5, 5.E-1 6, 6.000000000000014210854715202E-1 7, 6.9999999999999928945726423989E-1 8, 8.000000000000007105427357601E-1 9, 9.000000000000021316282072803E-1 10, 1.E+0 11, 1.1000000000000014210854715202E+0 If exact math were possible on the machine, both outputs should be 11 lines! The problem with the 1st run is that the fraction 1/10 cannot be exactly represented by the bits used in the machine word length. Note the direction of the rounding errors...when i=10, x is greater than 1 and the loop stops. Now for the Cray2: Output from epsilon.f: 2.5243548967072377773175314D-29.le.machine epsilon.le.5.0487097934144755546350628D-29 Output from roundoff.f: 1, 1.0000000000000008881784197D-1 2, 2.0000000000000017763568394D-1 3, 3.0000000000000026645352591D-1 4, 4.0000000000000035527136788D-1 5, 5.0000000000000044408920985D-1 6, 6.0000000000000053290705182D-1 7, 7.0000000000000062172489379D-1 8, 8.0000000000000071054273576D-1 9, 9.0000000000000079936057773D-1 10, 1.0000000000000008881784197D+0 Output from roundoff.f with same modification as above: 1, 1.0000000000000008881784197D-1 2, 2.0000000000000017763568394D-1 3, 3.0000000000000071054273576D-1 4, 4.0000000000000035527136788D-1 5, 5.0D-1 6, 6.0000000000000142108547152D-1 7, 6.9999999999999928945726423D-1 8, 8.0000000000000071054273576D-1 9, 9.0000000000000213162820728D-1 10, 1.0D+0 11, 1.1000000000000014210854715D+0 On a Sun 3/80, epsilon.f produces: 2.2204460492503D-16.le.machine epsilon.le. 4.4408920985006D-16 roundoff.f output: 1 1.0000000149012D-01 2 0.20000000298023 3 0.30000000447035 4 0.40000000596046 5 0.50000000745058 6 0.60000000894070 7 0.70000001043081 8 0.80000001192093 9 0.90000001341105 10 1.0000000149012 Now, making the same modifications to the code as above: 1 1.0000000149012D-01 2 0.20000000298023 3 0.30000001192093 4 0.40000000596046 5 0.50000000000000 6 0.60000002384186 7 0.69999998807907 8 0.80000001192093 9 0.90000003576279 10 1.0000000000000 11 1.1000000238419 The epsilon value is useful for terminating iterations instead of waiting for an output that will never come while the computer churns along seeking an exact answer. Epsilon could be the stopping criteria. The programs used were taken from Kalbasi, "Can You Trust Your Computer?", IEEE Potentials, April 1990, pp15-18.