高效设置大型数组 Fortran 77答案

【问题标题】：Efficiently setting up large arrays Fortran 77高效设置大型数组 Fortran 77
【发布时间】：2014-02-22 22:09:25
【问题描述】：

我正在使用 SUNDIAL 求解一个大型耦合非线性 ODE 系统，但遇到了性能问题。一些分析表明，我在设置定义我的方程组的数组（在 Fortran 77 中）时遇到了严重的瓶颈。特别是，我需要定义一个数组，它是一个嵌套的三重 do 循环，在最里面的循环中执行求和。这些数组在每个时间步被填充（我使用隐式微分代数方程求解器，因为我的控制方程有一个非平凡的系数矩阵），我正在寻找有关加速这些计算的潜在方法的建议/参考。我将提供一个示例代码来显示我的问题，以及到目前为止我尝试过的一些细节。

首先，一些上下文：

我的有问题的数组采用以下形式（完整的测试代码将包含在下面）

  do n= 0,niv
     do l= 0,niv
        do k= 0,niv
         sumdum = 0
         do m = 1, niv     
         sumdum = sumdum+(p(-m,l)*qmat(m,k,n)+p(m,k)*qmat(-m,l,n))
         enddo
         sqmat(k,l,n) = sumdum
         enddo
     enddo
  enddo

其中 niv 是自变量的数量。 niv 最终将介于 100 和 300 之间。p(m,l) 和 qmat(m,k,n) 是先前定义的（复杂）数组（参见下面的代码）。

现在，让我们做一些基准测试（Mac 2.5 GHz i5、gfortran、gcc v 4.8.2 x86_64）。使用 OpenMP 时钟，我看到对于 niv = 100，我的机器上的总测试代码需要 1.49 秒（编辑：我已经更改了这次以反映包含 O2 标志）

  call omp_set_num_threads ( 1 )

我编译用

  gfortran -fopenmp -O2 -Wl,-stack_size,0x40000000 test.f

现在，如果我看到我的代码在填充 sqmat 数组之前需要多长时间，我发现 t = 0.21 秒，这意味着大约 87% 的时间用于定义这些数组。这并不奇怪，因为这些数组有 4*(101)^3*100 次计算，而下一个最大的数组有 (201)^3 个元素，也就是说它小了大约 100 倍。另请注意，对于 niv = 150，我发现代码需要约 7.01 秒，其中 90% 的时间由 sqmat 数组消耗。最后，我注意到在 SUNDIALS 和求解这个 ODE 系统的上下文中，这些矩阵在 niv =100 时消耗了大约 50% 的计算时间。

希望我已经充分定义了我的问题，现在我将简要讨论我为加快速度所做的尝试。

显然，根据我上面所写的内容，我已转向 OpenMP。这提供了一些加速 - 2 个线程给了我 1.63 倍的加速，而 4 个线程只提供了 4 线程情况下的 1.02 倍加速。我使用 OpenMP 的方式包含在下面的代码中。我可以访问具有许多内核的机器，因此我很想看看我使用 OpenMP 的方式是否缺少某些东西，这可能会带来更大的加速。从表面上看，存在不同 sqmat 组件调用 p(m,j) 和 qmat(i,j,k) 的相同组件的问题。我不确定这是否会降低这里的性能，我也没有发现能够找到一种方法来规避这个问题，所以任何关于这方面的建议都将不胜感激。

第二条攻击线是利用 p 和 qmat 的对称性来尝试减少这些循环中的计算次数。这使我能够将计算次数减少一半。此外，通过检查，这些数组似乎是“稀疏的”，但我找不到这种稀疏性的结构 - 也许这可以以某种方式被利用。

我最近一直在考虑转向现有的处理矩阵乘法的 fortran 套件（类似于 PBLAS），希望这些能够更有效地计算 sqmat 项（因为 sqmat 是 p 和 q 的张量积）。是否有人对 PBLAS 有任何经验 - 与 fortran 的规范数组结构相比，这会显着提高速度吗？我对此很好奇，尤其是在并行化的背景下。

所以，总而言之，我的问题是如何有效地填充这些类型数组的组件。任何建议或参考，将不胜感激。

尼克

示例代码如下。请注意，对 openmp 的调用 c$omp 缩进以使代码全部包含在一个代码块中。

  program test

  implicit none

  include 'omp_lib.h'

  integer*4 iout(25), ipar, neq, niv, Mtot,ii
  parameter (Mtot = 10)
  double precision rout(10), rpar
  parameter (neq = 4*100)
  parameter (niv = neq/4)
  integer iatol, nout, jout, itask
  integer nst, kused, hused, No
  double precision t0, t1, rtol, tout, tret, tout1
  double precision y(neq), yp(neq), atol(neq), u(niv)
  data nst/3/, kused/9/, hused/2/
  integer reserr
  integer m, j, kk, l, n, l1, l2, l3, i, ll, indsum, indsum2
  integer indsum3, i4, k, k1, k2,k3, k4, k5, k6, l2p, l2pp
  integer k3a, k3b
  double precision  po, tcond, c, qout(1:niv,-niv:niv)
  double precision res(neq),res2(niv), res3(niv)
  double complex p(-niv:niv,-niv:niv), V(niv), resd2a,
 &  Vec(-niv:niv)
  double complex qmat(-niv:niv,-niv:niv,-niv:niv)
  double complex qmat2(-niv:niv,-niv:niv,-niv:niv)
  double complex sqmat(-niv:niv,-niv:niv,-niv:niv)
  double complex sqmat1(-niv:niv,-niv:niv,-niv:niv)
  double complex sqmat2(-niv:niv,-niv:niv,-niv:niv)
  double complex sqmat3(-niv:niv,-niv:niv,-niv:niv)
  double complex q(-niv:niv,-niv:niv), 
 & smat(-niv:niv,-niv:niv,-niv:niv)
  double complex sumdum, sum1, sum2(niv), sum3(niv), qsum,
 &  resd, resd2, sumdum2
  double complex ssum1, ssum2, forcing(niv)
  double precision seconds, seconds2, seconds3,seconds4,seconds5
  real :: start, finish
  call omp_set_num_threads ( 1 )
  seconds = omp_get_wtime ( )


  write ( *, '(a)' ) ' '
  write ( *, '(a,i8)' ) 
 &  '  Number of processors available = ', omp_get_num_procs ( )
  write ( * ,'(a,i8)' ) 
 &  '  Number of threads =              ', omp_get_max_threads ( )
  y = 0
  do i = 1, niv
  y(i) = 0.1d0*i
  enddo
  do i = niv+1,2*niv
  y(i) = -0.01d0*i
  enddo 
  c$omp parallel private(m,j,kk,l,k1,k2,k3a,k3b,k4,k5,k6) shared(p)
  c$omp do 
  do  j = -niv, niv
      do m = -niv, niv
      kk = m-j
      l = m-j
      k1 = m
      k2 = j
      k3a = m
      k3b = j
      k4 = m-j
      k5 = j
      k6 = m-j
      if (abs(kk) .gt. niv) then
      kk = 0
      else 
      kk = 1
      end if
      if (k1 .eq. 0) then
      k1 = 0
      else
      k1=1
      end if
      if (k2.eq.0) then
      k2 = 0
      else
      k2 = 1
      end if
      if (k3a .eq. 0) then
      k3a=1
      else 
      k3a=0
      end if
      if (k3b .eq. 0) then
      k3b=1
      else 
      k3b=0
      end if
      if (k4 .eq. 0) then
      k4 = 0
      else
      k4 = 1
      end if
      if (k5 .eq. 0) then
      k5=2
      else
      k5=1
      end if
      if (k6 .eq. 0) then
      k6=2
      else
      k6=1
      end if
  p(m,j)=(0.5*k5*k6*kk*abs(k4*sign(1,m-j)-k2*sign(1,j))*cmplx(k4
 & *y(abs(m-j))+1-k4,k4*y(niv+abs(m-j))*sign(1,m-j))-0.5*cmplx(k1*
 & y(abs(m))+(1-k1), k1*y(abs(m)+niv)*sign(1, m))*cmplx(k2*y(abs(j
 & ))+(1-k2),-k2*y(abs(j)+niv)*sign(1,j)))/((abs(m)+k3a)**(0.5)*(
 & k3b+abs(j)))
   enddo  
  enddo 
  c$omp end do
  c$omp end parallel 
  c$omp parallel reduction(+:qsum) private(m,j,n) shared(q)
  c$omp do 
  do j = -niv, niv
   do m = -niv, niv
    qsum = 0.0d0
    do n = 1, niv
    qsum = qsum + p(n,m)*p(-n,j)
    enddo
    q(m,j) = qsum
   enddo 
  enddo
  c$omp end do
  c$omp end parallel 
  q= 0.25d0*q


  c$omp parallel private(m,j,n,kk,k1,k2,k3a,k3b,k5,k6,l,indsum,
  c$omp& indsum2,i4,indsum3) shared(qmat)
  c$omp do  
  do m= -niv, niv
  do j= -niv, niv
  do n= -niv, niv
  kk = m-j
  k1 = m
  k2 = j
  k3a = m
  k3b = j
  k5 = m-j 
  k6 = j
  l = m-j
   if (abs(kk). gt. niv) then 
    kk = 0
   else
    kk = 1
   end if
   indsum = m-j-n
   if (indsum .eq. 0) then
   indsum = 1
   else
   indsum = 0
   end if
   indsum2 = m+j-n
   if (indsum2 .eq. 0) then
   indsum2 = 1
   else
   indsum2 = 0
   end if
   indsum3 = j+n
   if (indsum3 .eq. 0) then
   indsum3 = 1
   else
   indsum3 = 0
   end if
   i4 = m-n
   if (i4 .eq. 0) then
   i4 = 1
   else
   i4 = 0
   end if 
   if (k3a .eq. 0) then
   k3a=1
   else 
   k3a=0
   end if
   if (k3b .eq. 0) then
   k3b=1
   else 
   k3b=0
   end if
   if (k1 .eq. 0) then
   k1 = 0
   else
   k1 = 1
   end if
   if (k2 .eq. 0) then
   k2 = 0
   else
   k2 = 1
   end if
      if (k5 .eq. 0) then
      k5=2
      else
      k5=1
      end if
      if (k6 .eq. 0) then
      k6=2
      else
      k6=1
      end if
      if (l.eq.0) then
      l = 0
      else
      l = 1
      end if
   qmat(m,j,n) = (0.5*k5*k6*kk*abs(l*sign(1,m-j)-k2*sign(1,j))*
 &    cmplx(indsum,0)-0.5*indsum3*cmplx(y(abs(m))+1-k1,sign(1,
 &    m)*y(abs(m)+niv))-0.5*i4*cmplx(y(abs(j))+1-k2, -sign(1,
 &    j)*y(abs(j)+niv)))/((abs(m)+k3a)**(0.5)*(
 &    k3b+abs(j)))
  enddo 
  enddo 
  enddo 
  c$omp end do
  c$omp end parallel 


  c$omp parallel reduction(+:sumdum) private(k,l,n,m) shared(sqmat)
  c$omp do 
  do n= 0,niv
  do k= 0,niv
  do l= 0,niv
  sumdum = 0
  do m = 1, niv     
  sumdum = sumdum+(p(-m,l)*qmat(m,k,n)+p(m,k)*qmat(-m,l,n))
  enddo
  sqmat(k,l,n) = sumdum
  sqmat(-l,-k,-n) =conjg( sumdum)
  enddo
  enddo
  enddo
  c$omp end do 
  c$omp end parallel      
  c$omp parallel reduction(+:sumdum) private(k,l,n,m) shared(sqmat)
  c$omp do 
  do n= -niv, 0
  do l= 0,niv
  do k= 0,niv
  sumdum = 0
  do m = 1, niv    
  sumdum = sumdum+(p(-m,l)*qmat(m,k,n)+p(m,k)*qmat(-m,l,n))
  enddo
  sqmat(k,l,n) =sumdum
  sqmat(-l,-k,-n) = conjg(sumdum)
  enddo
  enddo
  enddo
  c$omp end do 
  c$omp end parallel 

  c$omp parallel reduction(+:sumdum) private(k,l,n,m) shared(sqmat)
  c$omp do         
  do n= -niv, 0
  do l= -niv, 0
  do k= 0, niv
  sumdum = 0
  do m = 1, niv     
  sumdum = sumdum+(p(-m,l)*qmat(m,k,n)+p(m,k)*qmat(-m,l,n))
  enddo
  sqmat(k,l,n) = sumdum
  sqmat(-l,-k,-n) = conjg(sumdum)
  enddo
  enddo
  enddo 
  c$omp end do
  c$omp end parallel 

  c$omp parallel reduction(+:sumdum) private(k,l,n,m) shared(sqmat)
  c$omp do SCHEDULE(DYNAMIC,1)
  do n= -niv, 0
  do k= -niv, 0
  do l=  0, niv      
  sumdum = 0
  do m = 1, niv    
  sumdum = sumdum+(p(-m,l)*qmat(m,k,n)+p(m,k)*qmat(-m,l,n))
  enddo
  sqmat(k,l,n) = sumdum
  sqmat(-l,-k,-n) = conjg(sumdum)
  enddo
  enddo
  enddo
  c$omp end do
  c$omp end parallel 


  sqmat = 0.25d0*sqmat   

  seconds = omp_get_wtime ( ) - seconds;
  print*, seconds
  return
  end

【问题讨论】：

较低代码的哪一部分是重要的（最后一个并行块看起来就像它们做的完全一样）？顺便说一句，我的 i5 在 10 秒（1 个线程）到 3 秒（4 个线程）内运行此代码。
他们正在填写数组 sqmat(i,j,k) 但在不同的索引范围内。我一开始让 i, j, k 在 -niv,niv 上运行 - 但可以利用 sqmat 的对称性将迭代总数减少一半。我是我的机器，我发现 1 个线程 = 8.6s，2 个线程是 5.9s，4 个线程是 4.09s。您使用哪些标志来编译？您是否在做诸如将线程绑定到内核之类的事情？我对加速的比例比实际数字更感兴趣，如果包含 -O2 标志，实际数字实际上要小得多，但在我的机器上加速保持不变
我已编辑时间以反映包含 O2 标志。时间要小得多，但多处理器提供的速度提升大致相当
我只使用gfortran -O2 -fopenmp test.f。可以通过发出ulimit -s unlimited 来规避堆栈大小的问题。也许，这在一定程度上会有所帮助。
关于我的问题：我没有认识到循环中l 和k 相对于之前的块的更改顺序。为了加快速度，我尝试使用 strides 和 SUM 将这个总和重新表述为向量运算，以在第一个循环中得到 sqmat(k,l,n) = SUM(p(-1:-niv:-1,l)*qmat(1:niv,k,n)+p(1:niv,k)*qmat(-1:-niv:-1,l,n))，但我看不到任何明显的区别。不过你可以试一试。

标签： arrays fortran openmp

【解决方案1】：

阅读您的帖子后，我的第一个想法是您应该进行更认真的分析。你给的时序信息很简陋，并没有对浮点运算性能或内存带宽使用情况进行实际分析。

看看你声称是代码中最慢的部分，我可以看到本质上是很多小矩阵乘法。您将它们作为每个元素的点积来处理，这肯定会比具有高效重用缓存的优化 BLAS 库更糟糕。提高性能的最佳方法也是尝试仅使用大型矩阵乘法（组合某些维度）来编写代码，这也将增加使用线程 BLAS 库的有用性。这样你应该得到更好的性能。

【讨论】：

抱歉回复缓慢，但我一直在慢慢解决这个问题。我最终使用 OpenBLAS 来加速这个（张量）矩阵乘法，我最终发现计算速度接近 10 阶。感谢您的建议。

【解决方案2】：

我无法将其放在评论中，但这更像是一个问题而不是答案。您的问题似乎与代码的最后一部分有关：

  c$omp parallel reduction(+:sumdum) private(k,l,n,m) shared(sqmat)
  c$omp do SCHEDULE(DYNAMIC,1)
  do n= -niv, 0
  do k= -niv, 0
  do l=  0, niv      
    sumdum = 0
    do m = 1, niv    
      sumdum = sumdum+(p(-m,l)*qmat(m,k,n)+p(m,k)*qmat(-m,l,n))
    enddo
    sqmat(k,l,n) = sumdum
    sqmat(-l,-k,-n) = conjg(sumdum)
  enddo
  enddo
  enddo
  c$omp end do
  c$omp end parallel

您在此处采用的并行方法有些令人困惑。我不明白，为什么你根本需要减少项，sumdum 可能是一个完全私有的变量，你说，k，l，n 都应该是独立的。在您的代码中，求和发生在最里面的循环中，该循环甚至没有并行化，除非您忘记在示例代码中包含折叠子句。现在只有你最外层的循环被并行化，如果我理解正确的话，你有 100-300 次循环迭代来分布在线程上。对于少于 10 个线程，这实际上应该没问题，但是我猜想减少和动态调度在这里完全杀死你。

在我看来，以下仍然是有效的并行实现：

  !$omp parallel private(k,l,n,m,sumdum) shared(sqmat)
  !$omp do collapse(3)
  do n= -niv, 0
  do k= -niv, 0
  do l=  0, niv      
    sumdum = 0
    do m = 1, niv    
      sumdum = sumdum+(p(-m,l)*qmat(m,k,n)+p(m,k)*qmat(-m,l,n))
    enddo
    sqmat(k,l,n) = sumdum
    sqmat(-l,-k,-n) = conjg(sumdum)
  enddo
  enddo
  enddo
  !$omp end do
  !$omp end parallel

您可以通过手动展开求和循环并使用条带挖掘方法对求和操作进行矢量化来进一步加快这一速度。

【讨论】：