Fortran并行计算中OpenMP，MPI及Coarrays之单节点多线程应用比较

风平老涡

如今的数值计算可供选择的方法有很多，一般的计算Python足以应付。而在高性能计算中，为了缩短计算时间及提高效率，通常采用并行计算，这是Fortran的强项。在Fortran并行计算时可有多种并行技术供选择，比如共享内存的OpenMP，分布式的MPI， 以及Fortran语言标准的Coarrays。

OpenMP通常用于单节点多线程进程（多线程并行），其特点是多线程进程可随时产生或合并（Fork & Merge），并且所产生的多线程进程可共享同一内存的变量（Shared Memory），避免了数据传输。OpenMP使用相对容易，只要在原先的单进程程序上，对所需并行部分语句加上OpenMP特有的注释行及少许改动即可。

MPI通常用于多节点分布型多进程进程（多进程并行），主要用于超算（Supercomputer）。其特点是多进程进程不可随时产生或合并，每个进程拥有本地内存（Local Memory）及变量。不同进程间不共享变量，进程间变量交换需通过数据传输。MPI应用是基于对通讯库的访问来实现，使用相对困难，其构架与原先的单进程程序上可能完全不一样。

Coarrays是Fortran2008语言标准，是属于一种叫PGAS（Partitioned Global Address Space）并行编程模型，可用于多节点多核多线程的单一程序多数据（Single Program Multiple Data， SPMD）类并行编程。由于Coarrays是语言标准的一部分，避免了用户直接使用类似MPI对通讯库的访问，简化了进程间变量交换的编程。 并且语言标准提供了内在的同步功能（Synchronization），使竞态（Race Condition， RC）和死锁（Deadlock）得以避免。


计算\pi有很多种方法，这里通过对Gregory-Leibniz级数\pi = 4 \sum_{n=1}^{\infty}\frac{(-1)^{(n-1)}}{2n-1}计算来比较几种不同的并行技术。


使用的软硬件如下：Intel i7-4710HQ CPU @ 2.50GHz（4核）， GNU Gfortran V10.1.0， OpenMP V4.5.0， OpenMPI V4.0.4， Opencoarray V2.9.0


一  OpenMP法

[Fortran] 纯文本查看 复制代码

program omp_parallel
  use omp_lib
  implicit none
  integer, parameter :: rk = 8
  integer :: n_threads, i, n_limit, c1, c2, c_rate
  real(kind=rk) :: pi

  print *, "n_limit="
  read(*,*)  n_limit
  print *, 'num_threads='
  read(*,*) n_threads
  pi = 0.0

  call system_clock(c1, c_rate)  
!$OMP PARALLEL DO DEFAULT(NONE) PRIVATE(i) NUM_THREADS(n_threads) REDUCTION(+:pi)  
  do i = 1, n_limit
    pi = pi + (-1)**(i+1) / real( 2*i-1, kind=rk )
  end do
!$OMP END PARALLEL DO
  pi = pi * 4.0_rk
  call system_clock(c2, c_rate)

  write(*,*) pi, real(c2-c1)/real(c_rate), n_threads
 
end program omp_parallel

二  MPI法

[Fortran] 纯文本查看 复制代码

program MPI_parallel
  use mpi
  implicit none
  integer, parameter :: rk = 8
  integer :: i, n_limit, ierr, numprocs, myid, c1, c2, c_rate
  real(kind=rk) :: pi, picalc

  call MPI_INIT(ierr)
  call MPI_COMM_RANK(MPI_COMM_WORLD, myid, ierr)
  call MPI_COMM_SIZE(MPI_COMM_WORLD, numprocs, ierr)

  if(myid == 0) then
     print *, "n_limit="
     read(*,*) n_limit
  end if
  call MPI_BCAST(n_limit, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
  pi = 0.0

  call system_clock(c1, c_rate)

  do i = myid + 1, n_limit, numprocs
     pi = pi + (-1)**(i+1) / real( 2*i-1, kind=rk )
  end do
  call MPI_REDUCE( pi, picalc, 1, MPI_DOUBLE_PRECISION, &
                 MPI_SUM, 0, MPI_COMM_WORLD, ierr )

  picalc = picalc * 4.0_rk
  call system_clock(c2, c_rate)
  if(myid == 0) then
     write(*,*) picalc, real(c2-c1)/real(c_rate), numprocs
  end if
  call MPI_FINALIZE(ierr)
end program MPI_parallel

三 Coarray法

[Fortran] 纯文本查看 复制代码

program coarray_parallel
  implicit none
  integer, parameter :: rk = 8
  integer :: n_images, i, c1, c2, c_rate, n_limit[*]
  real(kind=rk) :: pi[*]

  n_images = num_images()
  if(this_image() == 1) then
     print *, "n_limit="
     read(*,*) n_limit
  end if
  call co_broadcast(n_limit, 1)

  pi = 0.0
  
  call system_clock(c1, c_rate)
  do i = this_image(), n_limit, n_images
     pi = pi + (-1)**(i+1) / real( 2*i-1, kind=rk )
  end do
  
  call co_sum(pi)
  pi = pi * 4.0_rk
  call system_clock(c2, c_rate)

  if(this_image() == 1) then
     write(*,*) pi, real(c2-c1)/real(c_rate), n_images
  end if
  
end program coarray_parallel

四 比较
当N上限取值为2000000000，所需计算时间（秒）如下:
线程／进程    OpenMP       MPI       Coarrays
       1                   8.27           8.23          8.17
       2                   4.18           4.19          4.17
       4                   2.15           2.15          2.13
从上列数据可看出，三种方法在性能上差不多。这里使用的Opencoarray库是基于MPI3.0标准。因为Fortran语言标准没有指定采用那一种通讯技术，如MPI，OpenMP，SHMEM，GASnet，ARMCI，DMAPP等，所以对Coarrays的性能影响是多方面的。不管如何，Coarrays的通用性和容易使用是显而易见的。


五 总结

                  OpenMP      MPI        Coarrays
语言标准       否             否              是
共享内存       是             是              是
分布式           否             是              是
易使用           是             否              是
可塑性         一般           好              好
性能             一般           高              高

weixing1531

OpenMP法第一个注释语句结尾加入shared(n_limit)

风平老涡

唐汉 发表于 2021-7-19 19:26
这个帖子有些时间了，不知道答主有没有后续检查。我运行了一下coarray的例子也并没有发现问题，提供一点可 ...

多谢你的建议。加了sync all确实是比较保险的办法。因为co_broadcast语句出现在Fortran2018标准，还比较新，并不是所有的编译器能正确解释标准。新标准对co_broadcast之类的语句同步要求比较模糊，一边说执行该语句时没有自动同步，另一边确又说执行该语句时系统会提供同步。所以加sync all是比较保险的办法。

唐汉

这个帖子有些时间了，不知道答主有没有后续检查。我运行了一下coarray的例子也并没有发现问题，提供一点可能的思路吧。在coarray例子下，在co_broadcast之前（11行12行之间）你可能要加一个sync call。如果有一个image先于image 1赋值之前co_broadcast了image 1的值，就有可能得不到正确的n_limit的值。我看了一个opencoarrays里的例子（https://github.com/sourceryinstitute/OpenCoarrays/blob/ef8069d94ecd3a0d93e0d4413a63c2366b894130/src/tests/unit/collectives/co_broadcast.F90），其中

[Fortran] 纯文本查看 复制代码

! Verify broadcasting of character data from image 1
  c_char_co_broadcast: block
    character(kind=c_char,len=14), save :: string_received[*]
    character(kind=c_char,len=*), parameter :: string_sent=c_char_"Hello, world!"! Character test message
    if (me==1) string_received=string_sent
    sync all
    call co_broadcast(string_received,source_image=1)
    if (string_received/=string_sent) then
      write(error_unit,*) "Incorrect co_broadcast(",string_received,") on image",me
    else
      c_char_test_passes=.true.
    end if
  end block c_char_co_broadcast

就是在co_broadcast之前加了sync all。
一点小建议哈~

风平老涡

weixing1531 发表于 2020-11-4 18:38
Coarray法得不到PI的正确结果

不知道你用的是什么编译器，我用的是gfortran 10.1和opencoarray2.9, 计算结果在三种方法上都一样。

weixing1531

Coarray法得不到PI的正确结果

chiangtp

[Fortran] 纯文本查看 复制代码

!---------------------------------------
! Intel Core2 Quad @2.83GHz, Windows 10
! Compiler default optimizations
!
! Intel         9 (sec)
!
! GNU          14
! Absooft      21
! G95          21
! Lahey        24
!
! NAG         121
! PGI         127
!
! Silverfrost 331
! Compaq      546
!---------------------------------------
program no_parallel
  implicit none

  integer, parameter :: rk = SELECTED_REAL_KIND(P=15)
  real(kind=rk) :: pi

  integer :: i, n_limit, c1, c2, c_rate, c_max

  !-----------------

  n_limit = 2000000000

  call system_clock(c1, c_rate, c_max)

  pi = 0.0_rk
  do i = 1, n_limit
    pi = pi + (-1)**(i+1) / real( 2*i-1, kind=rk )
  end do
  pi = pi * 4.0_rk

  call system_clock(c2, c_rate)
  IF( c2 < c1 ) c2 = c2 + c_max

  write(*,*) pi, real(c2-c1)/real(c_rate)

end program no_parallel                                 

Intel的效率確實沒得比, GNU/Absoft/G95/Lahey持平, NAG/PGI令人跌破眼鏡, 還是我的測試有問題?