|
| 1 | +module assert |
| 2 | + |
| 3 | +use, intrinsic:: iso_fortran_env, only: stderr=>error_unit, real64, real32 |
| 4 | +use, intrinsic:: ieee_arithmetic, only: ieee_is_finite, ieee_is_nan |
| 5 | + |
| 6 | +implicit none (type, external) |
| 7 | + |
| 8 | +private |
| 9 | + |
| 10 | +public :: isclose, assert_allclose |
| 11 | + |
| 12 | +contains |
| 13 | + |
| 14 | +elemental logical function isclose(actual, desired, rtol, atol, equal_nan) |
| 15 | +!! ## inputs |
| 16 | +!! |
| 17 | +!! * actual: value "measured" |
| 18 | +!! * desired: value "wanted" |
| 19 | +!! * rtol: relative tolerance |
| 20 | +!! * atol: absolute tolerance |
| 21 | +!! * equal_nan: consider NaN to be equal |
| 22 | +!! |
| 23 | +!! https://www.python.org/dev/peps/pep-0485/#proposed-implementation |
| 24 | +!! https://github.com/PythonCHB/close_pep/blob/master/is_close.py |
| 25 | + |
| 26 | +class(*), intent(in) :: actual, desired |
| 27 | +class(*), intent(in), optional :: rtol, atol |
| 28 | +logical, intent(in), optional :: equal_nan |
| 29 | + |
| 30 | +real(real64) :: r,a, act, des |
| 31 | +logical :: n |
| 32 | + |
| 33 | +isclose = .false. !< ensure it's defined |
| 34 | + |
| 35 | +!> INSTEAD OF merge(), since non present values aren't defined. |
| 36 | +r = 1e-6 |
| 37 | +a = 1e-12 |
| 38 | +n = .false. |
| 39 | + |
| 40 | +select type (actual) |
| 41 | +type is (real(real32)) |
| 42 | + act = real(actual, real32) |
| 43 | +type is (real(real64)) |
| 44 | + act = actual |
| 45 | +class default |
| 46 | + error stop "assert: actual must be real32 or real64" |
| 47 | +end select |
| 48 | + |
| 49 | +select type (desired) |
| 50 | +type is (real(real32)) |
| 51 | + des = real(desired, real32) |
| 52 | +type is (real(real64)) |
| 53 | + des = desired |
| 54 | +class default |
| 55 | + error stop "assert: desired must be real32 or real64" |
| 56 | +end select |
| 57 | + |
| 58 | +if (present(rtol)) then |
| 59 | + select type (rtol) |
| 60 | + type is (real(real64)) |
| 61 | + r = rtol |
| 62 | + type is (real(real32)) |
| 63 | + r = real(rtol, real64) |
| 64 | + class default |
| 65 | + error stop "assert: rtol needs real32 or real64" |
| 66 | + end select |
| 67 | +endif |
| 68 | + |
| 69 | +if (present(atol)) then |
| 70 | + select type (atol) |
| 71 | + type is (real(real64)) |
| 72 | + a = atol |
| 73 | + type is (real(real32)) |
| 74 | + a = real(atol, real64) |
| 75 | + class default |
| 76 | + error stop "assert: atol needs real32 or real64" |
| 77 | + end select |
| 78 | +endif |
| 79 | + |
| 80 | +if (present(equal_nan)) n = equal_nan |
| 81 | + |
| 82 | +!print*,r,a,n,act,des |
| 83 | + |
| 84 | +!> sanity check |
| 85 | +if ((r < 0).or.(a < 0)) error stop 'improper tolerances specified' |
| 86 | +!> simplest case -- too unlikely, especially for arrays? |
| 87 | +!isclose = (act == des) |
| 88 | +!if (isclose) return |
| 89 | + |
| 90 | +!> equal nan |
| 91 | +if (n) then ! fortran is NOT short circuit logic in general |
| 92 | + isclose = (ieee_is_nan(act) .and. ieee_is_nan(des)) |
| 93 | + if (isclose) return |
| 94 | +endif |
| 95 | + |
| 96 | +!> Inf /= -Inf, unequal NaN |
| 97 | +if (.not.ieee_is_finite(act) .or. .not.ieee_is_finite(des)) return |
| 98 | + |
| 99 | +!> floating point closeness check |
| 100 | +isclose = abs(act-des) <= max(r * max(abs(act), abs(des)), a) |
| 101 | + |
| 102 | +end function isclose |
| 103 | + |
| 104 | + |
| 105 | +impure elemental subroutine assert_allclose(actual, desired, rtol, atol, equal_nan, err_msg) |
| 106 | + |
| 107 | +!! ## inputs |
| 108 | +!! |
| 109 | +!! * actual: value "measured" |
| 110 | +!! * desired: value "wanted" |
| 111 | +!! * rtol: relative tolerance |
| 112 | +!! * atol: absolute tolerance |
| 113 | +!! * equal_nan: consider NaN to be equal |
| 114 | +!! * err_msg: message to print on mismatch |
| 115 | + |
| 116 | +class(*), intent(in) :: actual, desired |
| 117 | +class(*), intent(in), optional :: rtol, atol |
| 118 | +logical, intent(in), optional :: equal_nan |
| 119 | +character(*), intent(in), optional :: err_msg |
| 120 | +character(:), allocatable :: emsg |
| 121 | + |
| 122 | +real(real64) :: act, des |
| 123 | + |
| 124 | +select type (actual) |
| 125 | +type is (real(real32)) |
| 126 | + act = real(actual, real32) |
| 127 | +type is (real(real64)) |
| 128 | + act = actual |
| 129 | +class default |
| 130 | + error stop "assert: actual must be real32 or real64" |
| 131 | +end select |
| 132 | + |
| 133 | +select type (desired) |
| 134 | +type is (real(real32)) |
| 135 | + des = real(desired, real32) |
| 136 | +type is (real(real64)) |
| 137 | + des = desired |
| 138 | +class default |
| 139 | + error stop "assert: desired must be real32 or real64" |
| 140 | +end select |
| 141 | + |
| 142 | +if (present(err_msg)) then |
| 143 | + emsg = err_msg |
| 144 | +else |
| 145 | + emsg = 'assert: MISMATCH' |
| 146 | +endif |
| 147 | + |
| 148 | +if (.not.isclose(act, des, rtol,atol,equal_nan)) then |
| 149 | + write(stderr,*) emsg // ': actual',act,'desired',des |
| 150 | + error stop |
| 151 | +endif |
| 152 | + |
| 153 | +end subroutine assert_allclose |
| 154 | + |
| 155 | +end module assert |
0 commit comments