# QuadPack -- with OpenMP -- problem/solution
**URL:** https://fortran-lang.discourse.group/t/quadpack-with-openmp-problem-solution/11113
**Category:** Uncategorized
**Created:** [September 22, 2026, 9:50am UTC](https://fortran-lang.discourse.group/t/quadpack-with-openmp-problem-solution/11113 "2026-09-22T09:50:38Z")
**Posts on this page:** 3
**Page:** 1
### Author:  [@sshestov](https://fortran-lang.discourse.group/u/sshestov)
#### Post date: [September 22, 2026, 9:50am UTC](https://fortran-lang.discourse.group/t/quadpack-with-openmp-problem-solution/11113/1 "2026-09-22T09:50:38Z")
Dear all,
This is a continuation of my previous question from [here](https://fortran-lang.discourse.group/t/modern-fortran-quadpack/3714/16) 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.
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
```fortran
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:
```auto
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!
---
### Author:  [@PierU](https://fortran-lang.discourse.group/u/PierU)
#### Post date: [September 22, 2026, 12:52pm UTC](https://fortran-lang.discourse.group/t/quadpack-with-openmp-problem-solution/11113/2 "2026-09-22T12:52:32Z")
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.
---
### Author:  [@PierU](https://fortran-lang.discourse.group/u/PierU)
#### Post date: [September 23, 2026, 7:13am UTC](https://fortran-lang.discourse.group/t/quadpack-with-openmp-problem-solution/11113/4 "2026-09-23T07:13:06Z")
> [@septc](#):
>
> I guess the above approach may be nice (i.e., utilize a thin derived type that contains a few data like rho + a type-bound procedure that wraps/calls an external / QUADPACK routine).
I may be wrong but I don’t think that a wrapper can do without modifying QUADPACK