I have a numerical function
f(x, y) returning a double floating point number that implements some formula and I want to check that it is correct against analytic expressions for all combination of the parameters
y that I am interested in. What is the proper way to compare the computed and analytical floating point numbers?
Let's say the two numbers are
b. So far I've been making sure that both absolute (
abs(a-b) < eps) as well as relative (
abs(a-b)/max(abs(a), abs(b)) < eps) errors are less than eps. That way it will catch numerical inaccuracies even if the numbers are let's say around 1e-20.
However, today I discovered a problem, the numerical value
a and analytic value
In : a Out: 5.9781943146790832e-322 In : b Out: 6.0276008792632078e-322 In : abs(a-b) Out: 4.9406564584124654e-324 In : abs(a-b) / max(a, b) Out: 0.0081967213114754103
So the absolute error  is (obviously) small, but the relative error  is large. So I thought that I have a bug in my program. By debugging, I realized, that these numbers are denormal. As such, I wrote the following routine to do the proper relative comparison:
real(dp) elemental function rel_error(a, b) result(r) real(dp), intent(in) :: a, b real(dp) :: m, d d = abs(a-b) m = max(abs(a), abs(b)) if (d < tiny(1._dp)) then r = 0 else r = d / m end if end function
tiny(1._dp) returns 2.22507385850720138E-308 on my computer. Now everything works and I simply get 0 as the relative error and all is ok.
In particular, the above relative error  is wrong, it's simply caused by insufficient accuracy of the denormal numbers. Is my implementation of the
rel_error function correct? Should I just check that
abs(a-b) is less than tiny (=denormal), and return 0? Or should I check some other combination, like
I would just like to know what the "proper" way is.