Fortran OpenMP向量归约失效问题：求并行化替代方案

简化后的代码如下：

!$omp parallel do private(e, b0_vek) reduction(+:b_vek) schedule(static, chunk_elem)
   do e = 1, n_elem
      call AddSource(elem(e), b0_vek)
      slice=elem(e)%node(1:6) !vector with indices in the global node list
      b_vek(slice)=b_vek(slice)+b0_vek(1:6)
   end do  ! e = 1, n_elem

计算需遍历所有单元（n_elem），每个单元由6个节点定义。通过AddSource子程序计算每个单元的b0_vek(1:6)向量，需将其添加到全局向量b_vek的对应片段（slice）中，但相邻单元可能包含相同节点，导致片段重叠。当前reduction操作无法正常工作，请问使用OpenMP还有其他可行的实现方式吗？

可行的实现方式 #

1. 原子操作 #

对每个节点的更新单独加原子锁，避免多个线程同时修改同一节点的值：

!$omp parallel do private(e, b0_vek, slice, i) schedule(static, chunk_elem)
   do e = 1, n_elem
      call AddSource(elem(e), b0_vek)
      slice = elem(e)%node(1:6)
      do i = 1, 6
         !$omp atomic update
         b_vek(slice(i)) = b_vek(slice(i)) + b0_vek(i)
      end do
   end do

这种方式实现简单，但如果单元和节点数量极大，原子操作的锁开销会影响性能。

2. 私有临时数组+全局归约 #

每个线程先把更新累积到自己的私有临时数组，最后再合并到全局数组，从根源上避免竞争：

!$omp parallel private(e, b0_vek, slice, i) shared(b_vek, elem, n_elem)
      ! 每个线程初始化私有临时数组
      integer, allocatable :: b_private(:)
      allocate(b_private(size(b_vek)))
      b_private = 0

      !$omp do schedule(static, chunk_elem)
      do e = 1, n_elem
         call AddSource(elem(e), b0_vek)
         slice = elem(e)%node(1:6)
         do i = 1, 6
            b_private(slice(i)) = b_private(slice(i)) + b0_vek(i)
         end do
      end do
      !$omp end do

      ! 临界区合并结果
      !$omp critical
         b_vek = b_vek + b_private
      !$omp end critical

      deallocate(b_private)
   !$omp end parallel

如果编译器支持OpenMP 3.0及以上的数组归约，可以把临界区换成reduction子句，进一步优化：

!$omp parallel private(e, b0_vek, slice, i) reduction(+:b_vek)
      !$omp do schedule(static, chunk_elem)
      do e = 1, n_elem
         call AddSource(elem(e), b0_vek)
         slice = elem(e)%node(1:6)
         do i = 1, 6
            b_vek(slice(i)) = b_vek(slice(i)) + b0_vek(i)
         end do
      end do
      !$omp end do
   !$omp end parallel

注意：数组归约的支持度依赖编译器，比如GCC、Intel Fortran都支持，但部分小众编译器可能不兼容。

3. 按节点分区（依赖网格结构） #

如果节点的索引有明确的分组规律（比如按空间区域划分），可以把b_vek拆分到不同线程负责的区域，每个线程只更新自己区域内的节点，彻底消除竞争。但这种方式需要结合具体的网格结构，通用性不强。

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