我的问题是：我有一个带有一定数量嵌套循环的fortran代码，首先我想知道是否可以优化(重新排列)它们以获得时间增益？其次，我想知道我是否可以使用OpenMP来优化它们？
我看过很多关于fortran中嵌套do循环以及如何优化它们的帖子，但我没有找到一个适合我的例子。我也在fortran中搜索过关于OpenMP的嵌套do循环，但我在OpenMP中是0级的，我很难知道如何在我的情况下使用它。
这里有两个非常相似的循环示例，第一个：
do p=1,N
do q=1,N
do ab=1,nVV
cd = 0
do c=nO+1,N
do d=c+1,N
cd = cd + 1
A(p,q,ab) = A(p,q,ab) + (B(p,q,c,d) - B(p,q,d,c))*C(cd,ab)
end do
end do
kl = 0
do k=1,nO
do l=k+1,nO
kl = kl + 1
A(p,q,ab) = A(p,q,ab) + (B(p,q,k,l) - B(p,q,l,k))*D(kl,ab)
end do
end do
end do
do ij=1,nOO
cd = 0
do c=nO+1,N
do d=c+1,N
cd = cd + 1
E(p,q,ij) = E(p,q,ij) + (B(p,q,c,d) - B(p,q,d,c))*F(cd,ij)
end do
end do
kl = 0
do k=1,nO
do l=k+1,nO
kl = kl + 1
E(p,q,ij) = E(p,q,ij) + (B(p,q,k,l) - B(p,q,l,k))*G(kl,ij)
end do
end do
end do
end do
end do
另一个是：
do p=1,N
do q=1,N
do ab=1,nVV
cd = 0
do c=nO+1,N
do d=nO+1,N
cd = cd + 1
A(p,q,ab) = A(p,q,ab) + B(p,q,c,d)*C(cd,ab)
end do
end do
kl = 0
do k=1,nO
do l=1,nO
kl = kl + 1
A(p,q,ab) = A(p,q,ab) + B(p,q,k,l)*D(kl,ab)
end do
end do
end do
do ij=1,nOO
cd = 0
do c=nO+1,N
do d=nO+1,N
cd = cd + 1
E(p,q,ij) = E(p,q,ij) + B(p,q,c,d)*F(cd,ij)
end do
end do
kl = 0
do k=1,nO
do l=1,nO
kl = kl + 1
E(p,q,ij) = E(p,q,ij) + B(p,q,k,l)*G(kl,ij)
end do
end do
end do
end do
end do
这两个例子之间的微小差异主要在于循环的索引。我不知道你是否需要更多关于循环中不同整数的信息，但你通常有：nO<nO<N<nV。因此，我不知道是否有可能优化这些循环和/或以便于使用OpenMP的方式放置它们(我还不知道我是否会使用OpenMP，这将取决于我在没有它的情况下通过优化循环可以获得多少收益)。
我已经尝试过以不同的方式重新排列循环，但没有任何成功(没有时间增益)，我也尝试过一点OpenMP，但我对此了解不多，所以再次没有成功。
从最初的注释中可以看出，至少在某些情况下，您使用的内存可能比可用RAM更多，这意味着您可能使用的是交换文件，所有这些都会对性能产生不良影响。要解决这个问题，如果可能的话，你必须安装更多的RAM，或者深度重组你的代码，使其不同时存储完整的B阵列(到目前为止最大的一个)(如果可能，再次)。
现在，让我们假设您有足够的RAM。正如我在评论中所写的，B数组上的访问模式远不是最佳的，因为内部循环对应于B的最后一个索引，这可能会导致许多缓存未命中(在给定B的大小的情况下更是如此)。如果可能的话，可以改变循环顺序。
看看你的第一个例子，我专注于数组A的计算(数组E的计算看起来完全独立于A，所以它可以单独处理)：
!! test it at first without OpenMP
!!$OMP PARALLEL DO PRIVATE(cd,c,d,kl,k,l)
do ab=1,nVV
cd = 0
do c=nO+1,N
do d=c+1,N
cd = cd + 1
A(:,:,ab) = A(:,:,ab) + (B(:,:,c,d) - B(:,:,d,c))*C(cd,ab)
end do
end do
kl = 0
do k=1,nO
do l=k+1,nO
kl = kl + 1
A(:,:,ab) = A(:,:,ab) + (B(:,:,k,l) - B(:,:,l,k))*D(kl,ab)
end do
end do
end do
!!$OMP END PARALLEL DO
我做了什么：
- 将p和q上的循环从外部位置移动到内部位置(这并不总是像这里那么容易)
- 用数组语法替换了它们(没有预期的性能提升，只是代码更易于阅读)
现在内部循环(由数组语法抽象)处理内存中的连续元素，这对性能有很大好处。该代码甚至可以在(现在)外循环上进行OpenMP多线程处理。
编辑/提示
Fortran将数组存储在"；列主要顺序"；，即当递增第一索引时，访问存储器中的连续元素。在C中，数组存储在"0"中；行主要顺序"；，即当递增最后一个索引时，访问存储器中的连续元素。因此，一个一般的规则是在第一个茚上有内环(在C中则相反)。
如果您可以使用张量表示法和爱因斯坦求和规则来描述您想要执行的操作，那将非常有用。我觉得使用NumPy中的np.einsum之类的东西可以更简洁地编写代码。
对于循环嵌套的第二块(在B的方形子部分而不是三角形中迭代)，您可以尝试引入一些子程序或基元，从中构建完整的解决方案。
从下到上，您从两个矩阵的简单和开始。
!
! a_ij := a_ij + beta * b_ij
!
pure subroutine apb(A,B,beta)
real(dp), intent(inout) :: A(:,:)
real(dp), intent(in) :: B(:,:)
real(dp), intent(in) :: beta
A = A + beta*B
end subroutine
(对于原始帖子中的第一个代码块，您可以用一个只更新矩阵的上/下三角形的基元来替换这个基元)
更高一步是张量收缩
!
! a_ij := a_ij + b_ijkl c_kl
!
pure subroutine reduce_b(A,B,C)
real(dp), intent(inout) :: A(:,:)
real(dp), intent(in) :: B(:,:,:,:)
real(dp), intent(in) :: C(:,:)
integer :: k, l
do l = 1, size(B,4)
do k = 1, size(B,3)
call apb( A, B(:,:,k,l), C(k,l) )
end do
end do
end subroutine
请注意，C的维度必须与B的最后两个维度匹配。(在上面的原始循环嵌套中，C的存储顺序被交换(即c_lk而不是c_kl。)
向上，我们得到了B的两个不同子块的收缩，此外，A、C和D还有一个额外的外部维度：
!
! A_n := A_n + B1_cd C_cdn + B2_kl D_kln
!
! The elements of A_n are a_ijn
! The elements of B1_cd are B1_ijcd
! The elements of B2_kl are B2_ijkl
!
subroutine abcd(A,B1,C,B2,D)
real(dp), intent(inout), contiguous :: A(:,:,:)
real(dp), intent(in) :: B1(:,:,:,:)
real(dp), intent(in) :: B2(:,:,:,:)
real(dp), intent(in), contiguous, target :: C(:,:), D(:,:)
real(dp), pointer :: p_C(:,:,:) => null()
real(dp), pointer :: p_D(:,:,:) => null()
integer :: k
integer :: nc, nd
nc = size(B1,3)*size(B1,4)
nd = size(B2,3)*size(B2,4)
if (nc /= size(C,1)) then
error stop "FATAL ERROR: Dimension mismatch between B1 and C"
end if
if (nd /= size(D,1)) then
error stop "FATAL ERROR: Dimension mismatch between B2 and D"
end if
! Pointer remapping of arrays C and D to rank-3
p_C(1:size(B1,3),1:size(B1,4),1:size(C,2)) => C
p_D(1:size(B2,3),1:size(B2,4),1:size(D,2)) => D
!$omp parallel do default(private) shared(A,B1,p_C,B2,p_D)
do k = 1, size(A,3)
call reduce_b( A(:,:,k), B1, p_C(:,:,k))
call reduce_b( A(:,:,k), B2, p_D(:,:,k))
end do
!$omp end parallel do
end subroutine
最后，我们到达了选择B的子块的主要级别
program doit
use transform, only: abcd, dp
implicit none
! n0 [2,10]
!
integer, parameter :: n0 = 6
integer, parameter :: n00 = n0*n0
integer, parameter :: N, nVV
real(dp), allocatable :: A(:,:,:), B(:,:,:,:), C(:,:), D(:,:)
! N [100,200]
!
read(*,*) N
nVV = (N - n0)**2
allocate(A(N,N,nVV))
allocate(B(N,N,N,N))
allocate(C(nVV,nVV))
allocate(D(n00,nVV))
print *, "Memory occupied (MB): ", &
real(sizeof(A) + sizeof(B) + sizeof(C) + sizeof(D),dp) / 1024._dp**2
A = 0
call random_number(B)
call random_number(C)
call random_number(D)
call abcd(A=A, &
B1=B(:,:,n0+1:N,n0+1:N), &
B2=B(:,:,1:n0,1:n0), &
C=C, &
D=D)
deallocate(A,B,C,D)
end program
与PierU的回答类似，并行化是最外层的循环。在我的电脑上，对于N=50，这个重新设计的例程在串行执行时大约快8倍。在4个线程上使用OpenMP时，系数为20。对于N=100，我厌倦了等待原始代码；在4个线程上重新设计的版本花费了大约3分钟。
我用于测试的完整代码，可通过环境变量(ORIG=<0|1> N=100 ./abcd)进行配置，可在此处获得：https://gist.github.com/ivan-pi/385b3ae241e517381eb5cf84f119423d
通过更多的微调，应该可以进一步降低数字。使用像cuTENSOR这样的专用库(也在Fortran内部函数的框架下使用，如将张量核心引入标准Fortran中所述)或像张量压缩引擎这样的工具，可以获得更好的性能。
最后一件我觉得奇怪的事是B的大部分都没有使用过。子部分CCD_ 12和CCD_。