Clever/fast数组乘法求和的方式

Question

我要解一个二重积分

在我的程序中，可以将其转换为以下最小工作示例中的 i-,j- 循环：

program test

  implicit none

  integer :: i, j, n
  double precision, allocatable :: y(:), res(:), C(:,:,:)

  n=200

  allocate(y(n), res(n), C(n,n,n))

  call random_number(y)
  call random_number(C)

  res = 0.d0
  do i=1, n
    do j=1, n
      res(:) = res(:) + y(i) * y(j) * C(:, j, i)
    end do
  end do

  deallocate(y, res, C)

end program test

每次执行我都必须多次求解这个积分，分析告诉我这是我计算的瓶颈，消耗了超过 95% 的执行时间。

我想知道是否有可能以更聪明的方式解决这个问题，即，快速的方法，也许可以摆脱一个或两个循环。

我的问题不是使用编译器标志或并行化来优化代码，而是双循环是否是解决给定问题的最佳实践。通常循环很慢，我尽量避免它们。我当时在想，可以通过重塑或展开数组来避免循环。但是就是看不出来。

Answer 1

如果你用矩阵表示法写双循环，y(i)*y(j) 变成双元数 YY^t，Y 是一个 n x 1 矩阵。有了这个，您可以将循环重写为（伪代码）

do n=1,size(C,1)
  res(n) = sum( YY^t * C_n )
enddo

其中 C_n = C(n,:,:) 和 * 是逐元素乘法。除了您已经进行的逐元素计算之外，这还为您提供了另外两种计算结果的方法：

res(n) = sum( (YY^t) * C_n )
res(n) = sum( Y * (Y^t C_n) )

在这两种情况下，拥有连续的数据并对数组重新排序是有益的C：

  do i=1,n
    C2(:,:,i) = C(i,:,:)
  enddo !i

两种方法的浮点运算次数相同，略低于原始方法。那么让我们来衡量一下他们所有人的时间...

以下是使用 LAPACK 进行矩阵运算的实现（并在适用的情况下使用点积）：

1. sum( (YY^t) * C_n )

  call system_clock(iTime1)
  call dgemm('N','N',n,n,1,1.d0,y,n,y,1,0.d0,mat,n)

  nn=n*n
  do i=1,n
    res(i) = ddot( nn, mat, 1, C2(:,:,i), 1 ) 
  enddo !i

2。 sum( Y * (Y^t C_n) )

  do i=1,n
    call dgemm('N','N',1,n,n,1.d0,y,1,C2(:,:,i),n,0.d0,vec,1)
    res(i) = ddot( n, y, 1, vec, 1 ) 
  enddo !i

结果如下：

 Orig:   0.111000001    
 sum((YY^t)C):   0.116999999    
 sum(Y(Y^tC)):   0.187000006

你原来的实现是最快的！为什么？很可能是由于 CPU 上缓存的理想使用。 Fortran 编译器通常在优化循环方面非常聪明，在逐元素计算中，您只需添加和缩放向量，无需任何矩阵运算。这可以非常有效地利用。

那么，还有改进的余地吗？当然 :) 您在循环内执行的操作通常称为 axpy：y = a*x + y。这是一个常用的 BLAS 子程序——通常是高度优化的。利用这个会导致

  res = 0.d0
  do i=1, n
    do j=1, n
      call daxpy(n, y(i)*y(j), C(:,j,i), 1, res, 1)
    end do
  end do

并接受

 Orig (DAXPY):   0.101000004

大约快 10%。

这是完整的代码，所有的测量都是用 OpenBLAS 和 n=500 进行的（为了更好地查看影响）

program test

  implicit none

  integer :: i, j, n, nn
  double precision, allocatable, target :: y(:), res(:), resC(:), C(:,:,:), C2(:,:,:), mat(:,:), vec(:)
  integer           :: count_rate, iTime1, iTime2

  double precision :: ddot
  n=500

  allocate(y(n), res(n), resC(n), C(n,n,n), C2(n,n,n), mat(n,n), vec(n))

  call random_number(y)
  call random_number(C)

  ! Get the count rate
  call system_clock(count_rate=count_rate)

  ! Original Aproach
  call system_clock(iTime1)
  res = 0.d0
  do i=1, n
    do j=1, n
      res(:) = res(:) + y(i) * y(j) * C(:, j, i)
    end do
  end do
  call system_clock(iTime2)
  print *,'Orig: ',real(iTime2-iTime1)/real(count_rate)

  ! Original Aproach, DAXPY
  call system_clock(iTime1)
  resC = 0.d0
  do i=1, n
    do j=1, n
      call daxpy(n, y(i)*y(j), C(:,j,i), 1, resC, 1)
    end do
  end do
  call system_clock(iTime2)
  print *,'Orig (DAXPY): ',real(iTime2-iTime1)/real(count_rate)
!  print *,maxval( abs(resC-res) )

  ! Re-order
  do i=1,n
    C2(:,:,i) = C(i,:,:)
  enddo !i

  ! sum((YY^t)C)
  call system_clock(iTime1)
  call dgemm('N','N',n,n,1,1.d0,y,n,y,1,0.d0,mat,n)

  nn=n*n
  do i=1,n
    resC(i) = ddot( nn, mat, 1, C2(:,:,i), 1 ) 
  enddo !i
  call system_clock(iTime2)
  print *,'sum((YY^t)C): ',real(iTime2-iTime1)/real(count_rate)
!  print *,maxval( abs(resC-res) )

  ! sum(Y(Y^tC))
  call system_clock(iTime1)
  do i=1,n
    call dgemm('N','N',1,n,n,1.d0,y,1,C2(:,:,i),n,0.d0,vec,1)
    resC(i) = ddot( n, y, 1, vec, 1 ) 
  enddo !i
  call system_clock(iTime2)
  print *,'sum(Y(Y^tC)): ',real(iTime2-iTime1)/real(count_rate)
!  print *,maxval( abs(resC-res) )
end program test

Clever/fast数组乘法求和的方式

Clever/fast way of array multiplication and summation

arrays

integration

fortran

multiplication

do-loops