5. Is a refactoring harmless?
Simpson's rule is rewritten to sum the odd and the even nodes apart, without the merge in the loop:
fortran
!< Composite quadrature rules on a uniform grid.
module quad_rules
!< Composite quadrature rules on a uniform grid.
use, intrinsic :: iso_fortran_env, only : R8P => real64
implicit none
private
public :: R8P, integrand, midpoint, trapezoid, simpson
abstract interface
pure function integrand(x) result(y)
!< A function of one real variable.
import :: R8P
implicit none
real(R8P), intent(in) :: x !< Abscissa.
real(R8P) :: y !< Value at x.
endfunction integrand
endinterface
contains
pure function midpoint(f, a, b, n) result(s)
!< Composite midpoint rule with n intervals.
procedure(integrand) :: f !< Integrand.
real(R8P), intent(in) :: a !< Lower bound.
real(R8P), intent(in) :: b !< Upper bound.
integer, intent(in) :: n !< Number of intervals.
real(R8P) :: s !< Integral estimate.
real(R8P) :: h !< Interval width.
integer :: i !< Counter.
h = (b - a) / n
s = 0._R8P
do i = 1, n
s = s + f(a + (i - 0.5_R8P) * h)
enddo
s = s * h
endfunction midpoint
pure function trapezoid(f, a, b, n) result(s)
!< Composite trapezoidal rule with n intervals.
!<```fortran
!< use integrands, only : sine
!< print '(F8.5)', trapezoid(sine, 0._R8P, acos(-1._R8P), 64)
!<```
!=> 1.99960 <<<
procedure(integrand) :: f !< Integrand.
real(R8P), intent(in) :: a !< Lower bound.
real(R8P), intent(in) :: b !< Upper bound.
integer, intent(in) :: n !< Number of intervals.
real(R8P) :: s !< Integral estimate.
real(R8P) :: h !< Interval width.
integer :: i !< Counter.
h = (b - a) / n
s = 0.5_R8P * (f(a) + f(b))
do i = 1, n - 1
s = s + f(a + i * h)
enddo
s = s * h
endfunction trapezoid
pure function simpson(f, a, b, n) result(s)
!< Composite Simpson rule with n intervals (n even).
!<```fortran
!< use integrands, only : sine
!< print '(F8.5)', simpson(sine, 0._R8P, acos(-1._R8P), 64)
!<```
!=> 2.00000 <<<
procedure(integrand) :: f !< Integrand.
real(R8P), intent(in) :: a !< Lower bound.
real(R8P), intent(in) :: b !< Upper bound.
integer, intent(in) :: n !< Number of intervals, even.
real(R8P) :: s !< Integral estimate.
real(R8P) :: h !< Interval width.
real(R8P) :: odd !< Sum at the odd nodes.
real(R8P) :: even !< Sum at the even inner nodes.
integer :: i !< Counter.
#ifdef QUAD_CHECKS
if (mod(n, 2) /= 0) error stop 'simpson: the number of intervals must be even'
#endif
h = (b - a) / n
odd = 0._R8P
do i = 1, n - 1, 2
odd = odd + f(a + i * h)
enddo
even = 0._R8P
do i = 2, n - 2, 2
even = even + f(a + i * h)
enddo
s = (f(a) + f(b) + 4._R8P * odd + 2._R8P * even) * h / 3._R8P
endfunction simpson
endmodule quad_rulesThe tests pass:
console
$ fobis test --quiet
Building test suite (4 target(s))...
PASS test_exactness (0.00 s)
PASS test_order (0.00 s)
PASS test_romberg (0.00 s)
PASS test_symmetry (0.00 s)
Results: 4 passed, 0 failed, 0 skipped
Total: 0.00 sThe reference does not:
console
$ fobis verify
Verifying case area...
area: FAIL
n: equal (1 value(s), exact)
trapezoid: equal (1 value(s), exact)
simpson: 1/1 out of tolerance (exact): max abs 4.44089e-16, max rel 2.22043e-16, L2 4.44089e-16
got 2.0000165910479355, reference 2.000016591047936
changed since the reference:
git.commit: 'b15eb6be644457b2e64c9cc90d3fb60599820baa' -> 'd5caf8a979673ea1a976a138c43802b9df7b3ad3'
git.diff_sha256: None -> '621b6a31a343078c1ce8b693cde856942a6dfe2ca987bbf6de0d15cf494a42de'
git.dirty: False -> True
0 of 1 case(s) verified.
[exit status 1]The result moved by one unit in the last place (2.2e-16 relative): the sums are added in another order. Tests with a tolerance of 1e-14 cannot see this; an exact reference can.
Against a revision
A reference is one point in history. fobis diff-run REF compares the working tree with any revision directly: both are copied to a scratch directory (your tree and repository are never touched), built, and run on every case:
console
$ fobis diff-run HEAD
diff-run: extracting HEAD (d5caf8a97967) and copying the working tree to ~/tmp/fobis-diff-run-u_zmeskf
diff-run: case area: build and run on the revision
diff-run: case area: build and run on the working tree
diff-run: HEAD (d5caf8a97967) vs the working tree, 1 file(s) changed
M src/quad_rules.F90: quad_rules::simpson
note: both sides use the working copy's fetched dependencies (.fobis_deps)
area: DIFFERS
n: equal (1 value(s), exact)
trapezoid: equal (1 value(s), exact)
simpson: 1/1 out of tolerance (exact): max abs 4.44089e-16, max rel 2.22043e-16, L2 4.44089e-16
got 2.0000165910479355, HEAD 2.000016591047936
0 of 1 case(s) equivalent.
[exit status 1]It also names the procedures whose lines changed (quad_rules::simpson), so a difference points to its cause. --assess weighs the difference against the working tree's own spread under compiler perturbations (the rounding and precision families of the case's assessment):
console
$ fobis diff-run HEAD --assess
diff-run: extracting HEAD (d5caf8a97967) and copying the working tree to ~/tmp/fobis-diff-run-oj3z2pqy
diff-run: case area: build and run on the revision
diff-run: case area: build and run on the working tree
diff-run: case area: assessing the working tree's own spread
diff-run: HEAD (d5caf8a97967) vs the working tree, 1 file(s) changed
M src/quad_rules.F90: quad_rules::simpson
note: both sides use the working copy's fetched dependencies (.fobis_deps)
area: DIFFERS (simpson: beyond the working tree's own spread (every digit shared by its variants))
n: equal (1 value(s), exact)
trapezoid: equal (1 value(s), exact)
simpson: 1/1 out of tolerance (exact): max abs 4.44089e-16, max rel 2.22043e-16, L2 4.44089e-16
got 2.0000165910479355, HEAD 2.000016591047936
0 of 1 case(s) equivalent.
[exit status 1]Here every variant of the new code gives the same simpson value, so even a one-ulp difference is beyond the spread: the change is real, not compiler noise. Whether it is acceptable is a decision the case has to state.
Saying what equivalent means
A tolerance in the case: four ulps for every quantity (*). Absolute (abs A), relative (rel R) and ulp tolerances can be combined, per quantity:
ini
[case:area]
description = sin(x) over [0, pi] with 16 intervals, both rules
target = src/area.F90
mode = release
args = 16
quantities = stdout
tolerance = *: ulp 4console
$ fobis diff-run HEAD
diff-run: extracting HEAD (d5caf8a97967) and copying the working tree to ~/tmp/fobis-diff-run-wf65wg4_
diff-run: case area: build and run on the revision
diff-run: case area: build and run on the working tree
diff-run: HEAD (d5caf8a97967) vs the working tree, 2 file(s) changed
M fobos
M src/quad_rules.F90: quad_rules::simpson
note: both sides use the working copy's fetched dependencies (.fobis_deps)
area: EQUIVALENT
n: equal (1 value(s), ulp 4)
trapezoid: equal (1 value(s), ulp 4)
simpson: within tolerance (ulp 4): max abs 4.44089e-16, max rel 2.22043e-16, L2 4.44089e-16
got 2.0000165910479355, HEAD 2.000016591047936
1 of 1 case(s) equivalent.console
$ fobis verify
Verifying case area...
area: PASS
n: equal (1 value(s), ulp 4)
trapezoid: equal (1 value(s), ulp 4)
simpson: within tolerance (ulp 4): max abs 4.44089e-16, max rel 2.22043e-16, L2 4.44089e-16
got 2.0000165910479355, reference 2.000016591047936
changed since the reference (results unaffected):
git.commit: 'b15eb6be644457b2e64c9cc90d3fb60599820baa' -> 'd5caf8a979673ea1a976a138c43802b9df7b3ad3'
git.diff_sha256: None -> 'b6fd3581f1a4936c38d8063b242c29cbec3da545734ca0a61046eab76ee31f64'
git.dirty: False -> True
1 of 1 case(s) verified.The reference stays as it was: the new result is within the declared tolerance of it. "Results unaffected" is the verdict on the listed changes.
What you learned
fobis diff-run REF, --assess, tolerance = q: spec; *: spec (exact, abs, rel, ulp). Reference: differential runs, tolerances.
Next: 6. Is it faster?.