Better procedure design using the Fortran pure attribute
This article describes the resolution of an issue that arose when debugging a pure Fortran procedure
as part of work related to the ARCHER2 eCSE06-04 project 3decomp: Towards Adaptive Mesh Refinement
in the high-order CFD code Xcompact3d.
When writing a function or subroutine in Fortran it is often desirable to give it the pure
attribute.
Declaring a procedure as pure places certain constraints on what it can do which may allow the
compiler to apply additional optimisations, such as in-lining, and some use cases such as elemental
procedures may require it.
One of these constraints is that the procedure cannot perform I/O, including print’ing to screen;
just as these constraints may aid the compiler, the developer may perceive them as a hindrance,
especially when debugging the code, leading to objections to the use of pure.
I am going to argue that whenever the limitations of a pure procedure feel prohibitive, they are
instead indicating a design issue with the code.
Consider code listing 1 which computes one of the coefficients for the compact finite difference scheme presented by (1), given as
\begin{equation} A_i = \frac{ \begin{split} h_{i-1}h_{i}h_{i+1} + h^2_{i}h_{i+1} + h_{i-1}h_{i}h_{i+2} + h^2_{i}h_{i+2} - h_{i-1}h^2_{i}\alpha - h_{i-1}h_{i}h_{i+1}\alpha \\ - h_{i-1}h_{i}h_{i+2}\alpha - h_{i-1}h_{i}h_{i+1}\beta -h^2_{i}h_{i+1}\beta -h_{i-1}h^2_{i+1}\beta -2h_{i}h^2_{i+1}\beta - h^3_{i+1}\beta \\ + h_{i-1}h_{i}h_{i+2}\beta + h^2_{i}h_{i+2}\beta + 2h_{i-1}h_{i+1}h_{i+2}\beta + 4h_{i}h_{i+1}h_{i+2}\beta + 3h^2_{i+1}h_{i+2}\beta \end{split} }{h_{i+1}\left(h_{i} + h_{i+1}\right)\left(h_{i-1} + h_{i} + h_{i+1}\right) h_{i+2}} \ , \end{equation}
where \(h_{i}=x_{i}-x_{i-1}\) and \(\alpha=\beta\).
Clearly the implementation of this function will be very susceptible to coding error, and requires
testing - a simple test is that it reduces to a known value of the grid spacing on uniform grids,
with such a test given in listing 2.
Even as I implemented this I expected floating point errors to be an issue due to the need to cancel
many products of terms of similar magnitude on a uniform grid, however with the gfortran compiler
the test passed and all appeared to be well.
Testing with crayftn on ARCHER2 however caused a test failure with the result falling outside the
range of tolerance for the test, and similarly using ifort on a workstation also caused a test
failure.
After some investigation it was found that the tests would pass when floating point optimisations
were disabled (gfortran also failed the test with optimisations, which are off by default, enabled).
pure real function coeff_a_deltas(h) !! Compute the coefficient acting on f_{i+1} of the finite diference given a stencil of grid !! spacings. type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. real :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 select type(deltas => h%stencil) type is(real) hm2 = deltas(-2) hm1 = deltas(-1) h0 = deltas(0) hp1 = deltas(1) hp2 = deltas(2) class default error stop end select associate(beta => alpha) ! To match Gamet et al. (1999) coeff_a_deltas = hm1 * h0 * hp1 + h0**2 * hp1 + hm1 * h0 * hp2 + h0**2 * hp2 & - hm1 * h0**2 * alpha - hm1 * h0 * hp1 * alpha & - hm1 * h0 * hp2 * alpha - hm1 * h0 * hp1 * beta - h0**2 * hp1 * beta & - hm1 * hp1**2 * beta - 2.0 * h0 * hp1**2 * beta - hp1**3 * beta & + hm1 * h0 * hp2 * beta + h0**2 * hp2 * beta + 2.0 * hm1 * hp1 * hp2 * beta & + 4.0 * h0 * hp1 * hp2 * beta + 3.0 * hp1**2 * hp2 * beta coeff_a_deltas = coeff_a_deltas & / (hp1 * (h0 + hp1) * (hm1 + h0 + hp1) * hp2) end associate end function coeff_a_deltas
a = coeff_a(stencil) call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h)))
To debug this problem, I knew what the computed values inside coeff_a should be at various stages so
the simplest solution was to remove the pure attribute, insert some print statements and check these
values.
After some rearranging of the terms I arrived at the implementation shown in
listing 3, explicitly grouping terms that cancel, which passed the tests with
optimisations enabled, I could then remove the print statements and restore the pure attribute.
real function coeff_a_deltas_alt(h) !! Compute the coefficient acting on f_{i+1} of the finite diference given a stencil of grid !! spacings. type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. real :: hm2, hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 select type(deltas => h%stencil) type is(real) hm2 = deltas(-2) hm1 = deltas(-1) h0 = deltas(0) hp1 = deltas(1) hp2 = deltas(2) class default error stop end select associate(beta => alpha) ! To match Gamet et al. (1999) coeff_a_deltas_alt = h0 * ((hm1 + h0) * (hp1 + hp2) + 2.0 * hp1 * hp2 * beta) & ! = (14/3) h^3 + hm1 * h0 * hp2 * (beta - alpha) & ! Should cancel for case alpha = beta + (h0**2) * (hp2 - hp1) * beta & ! Should cancel for constant h + hm1 * hp1 * (2.0 * hp2 - h0 - hp1) * beta & ! Should cancel for constant h + (hp1**3) * (3.0 * hp2 - 2.0 * h0 - hp1) * beta & ! Should cancel for constant h + h0 * (2.0 * hp1 * hp2 * beta - hm1 * (h0 + hp1) * alpha) ! Should cancel for constant h print *, "+++", coeff_a_deltas_alt / (h0**3), 14.0 / 3.0, "+++" coeff_a_deltas_alt = coeff_a_deltas_alt & / (3.0 * hp1 * ((h0 + hp1) / 2.0) * hp2) ! => 14.0 / 9.0 when h = const print *, "+++", coeff_a_deltas_alt, 14.0 / 9.0, "+++" coeff_a_deltas_alt = coeff_a_deltas_alt & / (2.0 * ((hm1 + h0 + hp1) / 3.0)) ! => (14.0 / 9.0) / (2h) when h = const end associate end function coeff_a_deltas_alt
However, this removal and then restoration of the pure attribute should have been telling me that
there was a problem in the design.
As noted above, there were several stages to the computation which I knew expected values for —
these are themselves pure functions!
Rather than temporarily disabling purity I should instead implement each stage as a pure function,
the outer function can then call each of these and assemble the solution from their return values.
This is shown in listing 4, the implementation of the functions computing each
stage aren’t shown for brevity, however they reflect the body of the computation shown in
coeff_deltas_alt.
module pure real function coeff_a_deltas(h) !! Compute the coefficient acting on f_{i+1} of the finite diference given a stencil of grid !! spacings. type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. real :: numerator, numerator_corr, denominator, divisor call coeff_a_components(h, numerator, numerator_corr, denominator, divisor) coeff_a_deltas = ((numerator + numerator_corr) / denominator) / divisor end function coeff_a_deltas module pure subroutine coeff_a_components(h, numerator, numerator_corr, denominator, divisor) type(nucfd_stencil_deltas), intent(in) :: h !! Stencil of grid spacings for the finite difference. real, intent(out) :: numerator real, intent(out) :: numerator_corr real, intent(out) :: denominator real, intent(out) :: divisor real :: hm1, h0, hp1, hp2 ! Grid deltas at i -2, -1, 0, +1, +2 select type(deltas => h%stencil) type is(real) hm1 = deltas(-1) h0 = deltas(0) hp1 = deltas(1) hp2 = deltas(2) class default error stop end select associate(beta => alpha) ! To match Gamet et al. (1999) numerator = coeff_numerator(hm1, h0, hp2, hp2, beta) numerator_corr = coeff_numerator_corr(hm1, h0, hp1, hp2, alpha, beta) denominator = coeff_denominator(h0, hp1, hp2) divisor = coeff_divisor(hm1, h0, hp1) end associate end subroutine coeff_a_components
As these per-stage functions are relatively low level it is undesirable to expose them outside the
module, therefore a helper subroutine (coef_a_components in listing 4) provides
an easy to use interface returning the values for each stage.
An additional benefit is that we can now test the intermediate values in the test suite, shown in
listing 5 making the implementation even more robust!
call coeff_a_components(points_to_deltas(stencil), numerator, numerator_corr, denominator, divisor) call test_report("Coefficient A numerator", check_scalar(numerator, numerator_f1ref * (h**3))) call test_report("Coefficient A numerator correction", & check_scalar(numerator_corr, numerator_corr_f1ref * (h**3))) call test_report("Coefficient A denominator", check_scalar(denominator, denominator_f1ref * (h**3))) call test_report("Coefficient A divisor", check_scalar(divisor, divisor_f1ref * h)) a = coeff_a(stencil) call test_report("Coefficient A", check_scalar(a, aref / (2.0 * h)))
By retaining the pure attribute we have arrived at an implementation that is debuggable and better
tested than the original.
We could have done the same with regular procedures but this requires discipline when simply using
print is an easy solution, whereas the use of pure forced us down this design route.
In conclusion, pure shouldn’t be viewed as a restriction making it difficult to debug failing code,
rather when the pure attribute causes friction this likely indicates code design flaws.
Wanting to print a variable implies that its value should be tested and the original procedure needs
to be further decomposed.
To present a user-friendly interface the low-level procedures this creates can be hidden behind a
helper subroutine which enables testing and debugging not only the original interface, but also the
intermediate values that were formerly untestable.
1. References
[1] Gamet, L. and Ducros, F. and Nicoud, F. and Poinsot, T., {Compact finite difference schemes on non-uniform meshes. Application to direct numerical simulations of compressible flows}, International Journal for Numerical Methods in Fluids, 1999.