如何并行化嵌套循环

Question

一个与我的代码结构相同的小示例串行代码如下所示。

PROGRAM MAIN
IMPLICIT NONE
INTEGER          :: i, j
DOUBLE PRECISION :: en,ei,es
DOUBLE PRECISION :: ki(1000,2000), et(200),kn(2000)
OPEN(UNIT=3, FILE='output.dat', STATUS='UNKNOWN')
DO i = 1, 1000, 1
   DO j = 1, 2000, 1
      ki(i,j) = DBLE(i) + DBLE(j)
   END DO
END DO
DO i = 1, 200, 1
   en = 2.0d0/DBLE(200)*(i-1)-1.0d0
   et(i) = en
   es = 0.0d0
   DO j = 1, 1000, 1
      kn=ki(j,:)
      CALL CAL(en,kn,ei)
      es = es + ei
   END DO
   WRITE (UNIT=3, FMT=*) et(i), es
END DO
CLOSE(UNIT=3)
STOP
END PROGRAM MAIN

SUBROUTINE CAL (en,kn,ei)
IMPLICIT NONE
INTEGER          :: i
DOUBLE PRECISION :: en, ei, gf,p
DOUBLE PRECISION :: kn(2000)
p = 3.14d0
ei = 0.0d0
DO i = 1, 2000, 1
   gf = 1.0d0 / (en - kn(i) * p)
   ei = ei + gf
END DO
RETURN
END SUBROUTINE CAL

我在集群上运行我的代码，一个节点上有 32 个 CPU，一个节点上有 32 个 CPU 共享总共 250 GB 内存。我最多可以使用 32 个节点。

每完成一次内层循环，就会收集一个数据。完成所有外循环后，总共需要收集 200 条数据。如果只用一个 CPU 执行内层 Loop，则需要 3 天以上（72 小时以上）。

我想分别对内循环和外循环进行并行化？有人能建议如何并行化这段代码吗？

我可以分别对内循环和外循环使用 MPI 技术吗？如果是这样，如何区分执行不同循环（内循环和外循环）的不同CPU？

另一方面，我看到有人提到混合 MPI 和 OpenMP 方法的并行化。我可以对外部循环使用 MPI 技术，对内部循环使用 OpenMP 技术吗？如果是这样，如何在每次内循环完成后收集一个数据到CPU，在所有外循环完成后总共收集200个数据到CPU。如何区分分别执行inner Loop和outer Loop的不同CPU？

或者，有人会提供任何其他关于并行化代码和提高效率的建议吗？非常感谢您。

Answer 1

正如 cmets 中所述，一个好的答案需要更详细的问题。然而，乍一看，似乎并行化内部循环

DO j = 1, 1000, 1
  kn=ki(j,:)
  CALL CAL(en,kn,ei)
  es = es + ei
END DO

应该足以解决您的问题，或者至少它会是一个好的开始。首先我猜是循环有错误

DO i = 1, 1000, 1
  DO j = 1, 2000, 1
    ki(j,k) = DBLE(j) + DBLE(k)
  END DO
END Do

因为 k 设置为 0 并且没有地址对应于 0 的单元格（请参阅您的变量声明）。 ki 也被声明为 ki(1000,2000) 数组，而 ki(j,i) 是 (2000,1000) 数组。除了这些错误，我想 ki 应该计算为

ki(i,j) = DBLE(j) + DBLE(i)

如果属实，我建议你以下解决方案

PROGRAM MAIN
IMPLICIT NONE
INTEGER          :: i, j, k,icr,icr0,icr1
DOUBLE PRECISION :: en,ei,es,timerRate
DOUBLE PRECISION :: ki(1000,2000), et(200),kn(2000)
INTEGER,PARAMETER:: nthreads=1
call system_clock(count_rate=icr)
timerRate=real(icr)
call system_clock(icr0)
call omp_set_num_threads(nthreads)
OPEN(UNIT=3, FILE='output.dat', STATUS='UNKNOWN')
DO i = 1, 1000, 1
  DO j = 1, 2000, 1
    ki(i,j) = DBLE(j) + DBLE(i)
  END DO
END DO

DO i = 1, 200, 1
  en = 2.0d0/DBLE(200)*(i-1)-1.0d0
  et(i) = en
  es = 0.0d0
  !$OMP PARALLEL DO private(j,kn,ei) firstpribate(en) shared(ki) reduction(+:es)
  DO j = 1, 1000, 1
    kn=ki(j,:)
    CALL CAL(en,kn,ei)
    es = es + ei
  END DO
  !$OMP END PARALLEL DO 
  WRITE (UNIT=3, FMT=*) et(i), es
END DO
CLOSE(UNIT=3)
call system_clock(icr1)
write (*,*) (icr1-icr0)/timerRate ! return computing time 
STOP

END PROGRAM MAIN

SUBROUTINE CAL (en,kn,ei)
IMPLICIT NONE
INTEGER          :: i
DOUBLE PRECISION :: en, ei, gf,p
DOUBLE PRECISION :: kn(2000)
p = 3.14d0
ei = 0.0d0

DO i = 1, 2000, 1
  gf = 1.0d0 / (en - kn(i) * p)
  ei = ei + gf
END DO

RETURN
END SUBROUTINE CAL

我添加了一些变量来检查计算时间;-)。

对于 nthreads=1，此解的计算时间为 5.14 秒，对于 nthreads=2，计算时间为 2.75 秒。它不会将计算时间除以 2，但对于第一次尝试来说似乎很划算。不幸的是，在这台机器上我有一个核心 i3 proc。所以我不能比 nthreads=2 做得更好。但是，我想知道，代码在 nthreads=16 的情况下会如何表现？？？

请告诉我

希望对你有所帮助。

最后，我警告在实际代码中可能要仔细考虑的变量状态（私有、第一私有和共享）的选择。