QuadPack -- with OpenMP -- problem/solution

Uncategorized
1

Dear all,

This is a continuation of my previous question from here about QuadPack’s QAWO/OpenMP troubles. AI provided some hints/solutions, but may be a) it is hallucinating; b) there are better solutions?

So before description of the problem:

Diffraction behind a circular occulter can be calculated using Hankel integral, which has a form of (1), but let’s stay with trigonometric functions for a while (2). I need to calculate the integral for a range of rho-s, thus some speed is desirable. Since it has the oscillating terms, some care should be taken. With variable substitution we can obtain (3) for which there are standard approaches.

description
description768×368 35.2 KB

There exist standard solutions-routines, like QAWO from QuadPack (or Levin method). Now I’m using QAWO, but I need to provide the second function (Bessel J0 in my case). I’ve implemented it with the following code

program main_prog
  ...
  !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(rho, val,other_params)
  do i=10,100
    rho=rho_axis(i)
    !!! a,b are iterated in a loop to avoid J0 oscillations
    call dqawo(qawo_f,a,b,dpi*lamz0_1,1,val,other_params)
  end do
  !$OMP END PARALLEL DO
  ...
contains

  real(wp) function qawo_f(x)
    real(wp), intent(in) :: x
    qawo_f = BESSEL_J0(2.*dpi*rho*sqrt(x)*lamz0_1)
  end function qawo_f

end main_prog

The code was working without OpenMP or with a single thread. But was breaking (wrong results) with many threads. After some investigations, AI suggested the following, I’m copy-pasting his text:

*******************************************************************************************************

Good, that set of checks actually narrows things down a lot. Since NUM_THREADS(1) fixes it, and you’ve ruled out SAVEd state in the library itself, the strongest remaining candidate is a well-known Fortran+OpenMP gotcha:

Host association through an internal (CONTAINS) procedure does not reliably see per-thread OpenMP-private copies when that procedure is called back indirectly through an external library.

Here’s the mechanism: qags_fc is internal to your subroutine, so it accesses rho via a compiler-generated “static link”/closure (often implemented as a runtime trampoline) pointing back to the host subroutine’s stack frame. But OMP PARALLEL DO PRIVATE(rho) is a different mechanism — it usually creates each thread’s private copy inside the outlined parallel-region microtask, which is technically a separate function frame from the “original” host subroutine frame that the internal procedure’s closure was built to point at. These two mechanisms (Fortran host association vs. OpenMP privatization) don’t always compose correctly, especially once you go through an indirect call (external dqags calling back into qags_fc as a passed dummy procedure). Depending on compiler and version, the callback can end up seeing a stale/shared value instead of the thread’s private one — which matches your exact symptom (works serially, breaks only under threading, no crash, just wrong numbers).

This is a documented pitfall with gfortran in particular (there are old gcc Bugzilla threads about internal procedures + trampolines + OpenMP), though it’s not strictly guaranteed-safe behavior in any implementation.

Fix: get rho out of host-association entirely. Move the integrand functions out of CONTAINS into a module, and make the shared parameter THREADPRIVATE:

module integrand_mod
  use, intrinsic :: iso_fortran_env, only: wp => real64
  implicit none
  real(wp) :: rho
  !$omp threadprivate(rho)

contains

  real(wp) function qawo_f(x)
    real(wp), intent(in) :: x
    qawo_f =  BESSEL_J0(2._wp*dpi*rho*sqrt(x)*lamz0_1)
  end function qawo_f

end module integrand_mod

Then in your driver, use integrand_mod, and inside the parallel loop just assign rho before each dqags call as before — but now it’s writing to the thread’s own THREADPRIVATE copy, and qags_fc/qags_fs are plain module procedures passed by address with no closure/trampoline involved at all. This removes the ambiguous interaction entirely rather than relying on it working out.

A couple of things worth telling me to pin this down further:

Which compiler/version are you using (gfortran, ifort/ifx, other)?
Does jacobwilliams/quadpack in that repo also expose an OOP-style interface (a derived type you extend/instantiate with the integrand as a type-bound procedure)? If so, that’s actually the cleanest fix of all — you’d hold rho as a component of a per-thread-local object instance instead of module state, avoiding THREADPRIVATE altogether. Worth checking their README if you want the more modern-Fortran-idiomatic route rather than the THREADPRIVATE patch above.

*******************************************************************************************************

(originally I was using QAGS and r_Aplane instead of rho, thus typos can remain)

And surprise-surprise, with these suggestions it works.

So my questions are: 1) is this explanation reasonable, does it invent things? 2) how do I pass rho and lamz0_1 into the function in a good way? I’m using another module-wide variable right now; 3) is there a better way to repair the original problem?

Thanks!

2

The AI is essentially right: in your initial version, when compiling qawo_f the compiler has no way to know that rho should be a private variable of an OPENMP parallel region.

The THREADPRIVATE solution is a possible fix, but it’s not elegant. And I’m not sure it’s very robust. But I don’t think it’s possible to do much better without modifying QUADPACK.

2 Likes