使用 OpenMP 在 Matlab 的 Mex 函数 (Fortran) 中并行化循环答案

【问题标题】：Do loop to be parallelized in a Matlab's Mex function (Fortran) with OpenMP使用 OpenMP 在 Matlab 的 Mex 函数 (Fortran) 中并行化循环
【发布时间】：2021-09-11 02:01:36
【问题描述】：

我在 Matlab 的 Mex 函数（用 Fortran 编写）中有一个 do 循环，它为 FEM 网格的每个元素执行一些计算。我的网格由 250k 个元素组成，所以我认为值得并行化它。这是我第一次尝试使用 OpenMP 并行化此代码（我是编码初学者）。我使用reduction 命令来避免fintk(dofele) = fintk(dofele) + fintele 中的竞争条件。这是正确的吗？我可以在 Matlab 中毫无问题地编译它。但是，当我使用它（在 Matlab 中）时，它会为 12k 元素网格产生正确的结果，并且比序列化的更快，但是当我尝试将它用于 250k 元素网格时，Matlab 崩溃了。谢谢你帮助我

  subroutine loop_over_elements( &
  ! OUT  
  fintk,Sxyz,&
  ! IN                    
  Elem,Bemesh,Dofelemat,u,dt,NE,NDOF)

  use omp_lib
  
  implicit none
  
  mwSize NE, NDOF, ele
  integer,  parameter   :: dp = selected_real_kind(15,307)
  real(dp) :: fintk(NDOF), Sxyz(6,NE), Elemat(4,NE), Bemesh(6,12,NE), Dofelemat(12,NE)
  real(dp) :: u(NDOF)  
  
  real(dp) :: Bele(6,12), fintele(12), uele(12), si(6), dt
  integer*4  :: nodes(4), dofele(12) 
  
  fintk = 0.D0
  
  !$OMP PARALLEL DO REDUCTION(+:fintk(:)) PRIVATE(ele,nodes,Bele,dofele,uele,si,fintele) 
      DO ele = 1, NE 
        nodes   = Elemat(1:4,ele)
        Bele    = Bemesh(1:6,1:12,ele)
        dofele  = Dofelemat(1:12,ele)
        uele    = u(dofele)
        
        call comput_subroutine( &
!       IN
        Bele,uele,dt, &
!       OUT                                        
        si) 
        
        Sxyz(:,ele)   = si
        fintele       = MATMUL(TRANSPOSE(Bele),si)
        fintk(dofele) = fintk(dofele) + fintele
        
      END DO

  !$OMP END PARALLEL DO

  return
  end

【问题讨论】：

原始序列码是否适用于 250k 元素网格？
@CrisLuengo 是的，确实如此，新版本中唯一的变化是引入了!$OMP PARALLEL DO REDUCTION 命令。

标签： matlab parallel-processing openmp mex

【解决方案1】：

在调用loop_over_elements 子例程之前，我通过在通用mexFunction 子例程中添加这一行来解决我遇到的Matlab“崩溃”问题： call KMP_SET_STACKSIZE(100000000)。我想，因为当我使用带有非并行子程序的大型模型时，Matlab 没有崩溃，可能是内存问题。之后，我在使用 OpenMp 和大型数组时发现了众所周知的（不幸的是对我来说不是）分段错误问题（参见例如this）。我仍然对设置OMP_STACKSIZE（我不知道如何在Mex 函数中执行此操作）和KMP_SET_STACKSIZE 之间的区别感到困惑，但现在并行化代码适用于大型模型。

【讨论】：