OpenMP多线程（>1）运行结果异常，寻求排查建议

我了解有许多类似标题的问题，已查阅大部分，但本次问题成因可能不同，故发帖求助，期待您的判断与建议。

重要说明：我本想提供最小复现示例，但涉及代码规模大且包含不可共享的依赖，故将尽可能详细说明代码结构，并提供一个无错误的示例以展示整体结构。

初始代码 #

我尝试并行化的核心代码是多层嵌套循环，仅显式并行化外层1到4的循环：

!$omp parallel do &
!$omp   default(firstprivate), &
!$omp   private(idir, timer, zone_title), &
!$omp   shared(DIRS_CH, ROTATIONS, NLimsP1, COORDS_DIR_CH, LIM_SIGN_DIRS &
!$omp          , policies, limits, df_I_ref, df_J_ref, msh_ZoneLimsInterestModes &
!$omp          , bases_ch, refmts, deltas, basePts, base_i &
!$omp          , struct_data, wd, settings, logger_debug &
!$omp          , id_im_last, maxF, NLims, getBFM_msh &
!$omp          , NFREQS, NNODES, NNODESL, NLIBS, NLIBSL &
!$omp          , NMODES, NMODES_EFF, MODES &
!$omp          , NPSDEL, NTCOMPS, NDIRS, TCOMPS, DIRS &
!$omp          , msh_izone, msh_totNPts, m3mf_msh_ptr_), &
!$omp   num_threads(4)
do idir = 1, N_DIRS
   ! step 1: general setup
   ...
   do ilim = 1, NLims
      ! step 2: other (more specific) setup
      ...
      ! HERE: main computation
      call rz%compute()
   enddo
   ! step 3: non relevant stuff..
   ...
enddo ! n dirs
!$omp end parallel do

核心计算（类型绑定）过程 #

如前一段代码注释所示，主计算在OMP threadprivate的派生类型变量rz的compute()过程中（注：示例中用proc_2()过程模拟）：

module subroutine compute(this)
   ...
   use omp_lib
   implicit none
   ...
   real(RDP), allocatable :: rres(:, :), intg(:)
   ...

   ! some computation based on *this instance's state
   ...

   allocate(intg(...))

   ! now invoking the main computing function pointer
   rres(:, 1) = getBFM_msh(...)
   
   ! for each of these calls, the local integral gets updated
   intg(:) = rres(:, 1) * somevar ! at first usage only assignment!!
   ...
   intg(:) = intg(:) + rres(:, 1) * somevar ! then, increment

   ! NOTE: this is done many times, even in some do loops !
   ...
   ...

   ! Once local integral is computed, global is updated and data dumped onto a global unit
   ! NOTE: critical to avoid data races and to guard file access to one thread only.
   !$omp critical
   m3mf_msh_ptr_ = m3mf_msh_ptr_ + intg   ! update main integral (shared OMP variable)

   call dumpData(this, rres)
   !$omp end critical
end subroutine

核心计算函数指针 #

实际核心计算在getBFM_msh()函数中，该函数接收输入并返回计算结果，后续会通过积分得到最终结果。函数内包含多层嵌套循环，最关键的是调用了LAPACK的dgesvd()例程：

module function procPtr_(...) result(rres)
   ...
   
   ! internal data allocation
   ...
   
   do itc = 1, NTCOMPS   ! first loop

      ...
      call dgesvd(&
              'O' &           ! min(M,N) columns of U are overwritten on array A (saves memory)
            , 'N' &           ! no rows of V are computed
            , NNODESL    &    ! n. of rows M
            , NNODESL    &    ! n. of cols N
            , S_uvw_w1   &    ! A matrix (overwritten with left-singular vectors)
            , NNODESL    &
            , D_S_uvw_w1 &    ! singular values
            , tmpv       &    ! U
            , 1          & 
            , tmpv       &    ! VT
            , 1          &
            , MSHR_SVD_WORK  &
            , MSHR_SVD_LWORK &
            , MSHR_SVD_INFO  &
         )

      ! Some other computations, which lead to  "rres"
      ...
   enddo   ! itc
end function

其中tmpv、S_uvw_w1和D_S_uvw_w1为局部变量，MSHR_SVD_WORK、MSHR_SVD_LWORK和MSHR_SVD_INFO为全局模块变量，已设为firstprivate（因每次调用dgesvd()都会修改这些变量）。

单线程运行时所有功能正常，但线程数设置为>1时结果错误。调试发现，调用dgesvd()后，作为输入输出的变量S_uvw_w1在单线程和多线程下存在差异，因此怀疑LAPACK例程自身已并行，导致并行嵌套出现问题。除此之外，我无法确定其他可能的错误原因。

module mod1

   implicit none
   integer :: sum_ = 0, ival_ = 0

   procedure(interf), pointer :: procPtr_ => null()
   abstract interface
      function interf(ival) result(ires)
         integer, intent(in) :: ival
         integer :: ires
      end function
   end interface

contains

   function ifunc1(ival) result(ires)
      integer, intent(in) :: ival
      integer :: ires

      if (mod(ival, 2) == 0) then
         ires = ival * 2
      else
         ires = ival
      endif
   end function

end module


program main

   use omp_lib
   use mod1
   implicit none

   procPtr_ => ifunc1
   call proc_1()
   print *, ' total sum = ', sum_

contains

   subroutine proc_1()
      integer :: i, j

      ! from some previous computations..
      sum_ = 10

      !$omp parallel do default(firstprivate) shared(sum_) &
      !$omp   num_threads(4)
      do i = 1, 4
         do j = 1, 10
            ival_ = i * j
            call proc_2()
         enddo
      enddo
      !$omp end parallel do
   end subroutine

   subroutine proc_2()
      integer :: j, temp_

      temp_ = procPtr_(ival_)

      !$omp critical
      sum_ = sum_ + temp_
      !$omp end critical
   end subroutine

end program

内容的提问来源于stack exchange，提问作者mEm