|
| 1 | +module assert |
| 2 | + |
| 3 | + use, intrinsic:: iso_c_binding, only: sp=>c_float, dp=>c_double |
| 4 | + use, intrinsic:: iso_fortran_env, only: stderr=>error_unit |
| 5 | + use, intrinsic:: ieee_arithmetic |
| 6 | + implicit none |
| 7 | + private |
| 8 | + |
| 9 | + integer,parameter :: wp = sp |
| 10 | + |
| 11 | + public :: wp,isclose, assert_isclose, err |
| 12 | + |
| 13 | +contains |
| 14 | + |
| 15 | +elemental logical function isclose(actual, desired, rtol, atol, equal_nan) |
| 16 | +! inputs |
| 17 | +! ------ |
| 18 | +! actual: value "measured" |
| 19 | +! desired: value "wanted" |
| 20 | +! rtol: relative tolerance |
| 21 | +! atol: absolute tolerance |
| 22 | +! equal_nan: consider NaN to be equal? |
| 23 | +! |
| 24 | +! rtol overrides atol when both are specified |
| 25 | +! |
| 26 | +! https://www.python.org/dev/peps/pep-0485/#proposed-implementation |
| 27 | +! https://github.com/PythonCHB/close_pep/blob/master/is_close.py |
| 28 | + |
| 29 | + real(wp), intent(in) :: actual, desired |
| 30 | + real(wp), intent(in), optional :: rtol, atol |
| 31 | + logical, intent(in), optional :: equal_nan |
| 32 | + |
| 33 | + real(wp) :: r,a |
| 34 | + logical :: n |
| 35 | + ! this is appropriate INSTEAD OF merge(), since non present values aren't defined. |
| 36 | + r = 1e-5_wp |
| 37 | + a = 0._wp |
| 38 | + n = .false. |
| 39 | + if (present(rtol)) r = rtol |
| 40 | + if (present(atol)) a = atol |
| 41 | + if (present(equal_nan)) n = equal_nan |
| 42 | + |
| 43 | + !print*,r,a,n,actual,desired |
| 44 | + |
| 45 | +!--- sanity check |
| 46 | + if ((r < 0._wp).or.(a < 0._wp)) stop 'tolerances must be non-negative' |
| 47 | +!--- simplest case |
| 48 | + isclose = (actual == desired) |
| 49 | + if (isclose) return |
| 50 | +!--- equal nan |
| 51 | + isclose = n.and.(ieee_is_nan(actual).and.ieee_is_nan(desired)) |
| 52 | + if (isclose) return |
| 53 | +!--- Inf /= -Inf, unequal NaN |
| 54 | + if (.not.ieee_is_finite(actual) .or. .not.ieee_is_finite(desired)) return |
| 55 | +!--- floating point closeness check |
| 56 | + isclose = abs(actual-desired) <= max(r * max(abs(actual), abs(desired)), a) |
| 57 | + |
| 58 | +end function isclose |
| 59 | + |
| 60 | + |
| 61 | +impure elemental subroutine assert_isclose(actual, desired, rtol, atol, equal_nan, err_msg) |
| 62 | +! inputs |
| 63 | +! ------ |
| 64 | +! actual: value "measured" |
| 65 | +! desired: value "wanted" |
| 66 | +! rtol: relative tolerance |
| 67 | +! atol: absolute tolerance |
| 68 | +! equal_nan: consider NaN to be equal? |
| 69 | +! err_msg: message to print on mismatch |
| 70 | +! |
| 71 | +! rtol overrides atol when both are specified |
| 72 | + |
| 73 | + real(wp), intent(in) :: actual, desired |
| 74 | + real(wp), intent(in), optional :: rtol, atol |
| 75 | + logical, intent(in), optional :: equal_nan |
| 76 | + character(*), intent(in), optional :: err_msg |
| 77 | + |
| 78 | + if (.not.isclose(actual,desired,rtol,atol,equal_nan)) then |
| 79 | + write(stderr,*) merge(err_msg,'',present(err_msg)),': actual',actual,'desired',desired |
| 80 | + stop |
| 81 | + endif |
| 82 | + |
| 83 | +end subroutine assert_isclose |
| 84 | + |
| 85 | + |
| 86 | +pure subroutine err(msg) |
| 87 | + character, intent(in) :: msg |
| 88 | + stop msg |
| 89 | +end subroutine err |
| 90 | + |
| 91 | +end module assert |
0 commit comments