Dear all, I’m having troubles with OpenMP-parallelized do-loops with QAWO or QAGS. My minimal (not working) example is below
!$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(r, r_Aplane, int_cos, int_sin, f_cos, f_sin, val, eps, neval, ierr, last, iwork, work)
do j=1, 500
r=A_axis(j) !! A-plane coord.
!! ************ my straightforward sampling of the 1st plane works well with OpenMP
f_cos = 2. * dpi * r_axis * cos(dpi*r_axis**2 *lamz0_1) * BESSEL_J0(2.*dpi*r*r_axis*lamz0_1)
f_sin = 2. * dpi * r_axis * sin(dpi*r_axis**2 *lamz0_1) * BESSEL_J0(2.*dpi*r*r_axis*lamz0_1)
int_cos = integrate_trap(real(r_axis,kind=8), f_cos)
int_sin = integrate_trap(real(r_axis,kind=8), f_sin)
!! ******** but I want to implement it via QUAD_PACK routines **************************
!! ************** variable substitution in the integral ksi <=> r_axis *****************
!! Int ( ksi*cos( .. ksi^2) * J0(2pi r ksi/lamz0) dksi ) =>
!! 1/2 Int ( cos(ksi^2) * J0(2pi r sqrt(ksi^2) / lamz0) dksi^2
r_Aplane=r
call dqags(qags_fc,0._wp,710.0_wp,1d-2, 1d-4 ,val, eps, neval, ierr, limit, lenw, last, iwork, work)
int_cos = val*0.5*2. ! 2 to conform the definition in my original code
r_Aplane=r
call dqags(qags_fs,0._wp,710.0_wp,1d-2, 1d-4 ,val, eps, neval, ierr, limit, lenw, last, iwork, work)
int_sin = val*0.5*2.
!!!! temporarily put newly computed real part to imaginary or vise versa
!!psi_imag(j) = 1. - ( cos(dpi*r**2*lamz0_1)*int_sin + sin(dpi*r**2*lamz0_1)*int_cos ) * lamz0_1
psi_real(j) = ( cos(dpi*r**2*lamz0_1)*int_cos - sin(dpi*r**2*lamz0_1)*int_sin ) * lamz0_1
if (mod(j,100) == 0) write(*,'(I4,A)', advance='no') int(j/100)," "
end do
!$OMP END PARALLEL DO
...
contains
real(wp) function qags_fc(x)
real(wp), intent(in) :: x
qags_fc = 2. * dpi * x * cos(dpi * x**2 * lamz0_1) * BESSEL_J0(2.*dpi*r_Aplane*x*lamz0_1)
end function qags_fc
end
When I use self-implemented trapezoid implementation and explicit sampling (in r_axis) it works flawlessly. But I wanted to switch to QUAD_PACK routines or a separate Levin method, with a hope that presence of cos/sin terms can speed up.
I’ve started from Levin method - it was working but giving spurious outliers; then QAWO was working only until some j. I guess this is because Bessel J0 is also oscillating. Finally I switched to QAGS, which seems to work, but I struggling with OpenMP.
It works, but gives completely wrong results.
Do you have any ideas/suggestions?
PS: this is to calculate diffraction behind a circular occulter.