Single vs OpenMP vs MPI - Fortran

Question

我是MPI编程新手。我必须测试3个代码，例如顺序，OpenMP和MPI代码。这些3码（不是真正的代码，只是举例）分别给定为遵循

顺序码

program no_parallel 
implicit none 
integer, parameter       :: dp = selected_real_kind(15,307) 
integer         :: i, j 
real(kind = dp)       :: time1, time2 
real(kind = dp), dimension(1000000)  :: a 
    !Initialisation 
     do i = 1, 1000000 
      a(i) = sqrt(dble(i)/3.0d+0) 
     end do 
    call cpu_time(time1) 
     do j = 1, 1000 
      do i = 1, 1000000 
       a(i) = a(i) + sqrt(dble(i)) 
      end do 
     end do 
    call cpu_time(time2) 
    print *, a(1000000) 
    print *, 'Elapsed real time = ', time2 - time1, 'second(s)' 
end program no_parallel 

的OpenMP的代码

program openmp 
implicit none 
integer, parameter       :: dp = selected_real_kind(15,307) 
integer         :: i, j 
real(kind = dp)       :: time1, time2, omp_get_wtime 
real(kind = dp), dimension(1000000)  :: a 
    !Initialisation 
     do i = 1, 1000000 
      a(i) = sqrt(dble(i)/3.0d+0) 
     end do 
    time1 = omp_get_wtime() 
    !$omp parallel 
     do j = 1, 1000 
      !$omp do schedule(runtime) 
      do i = 1, 1000000 
       a(i) = a(i) + sqrt(dble(i)) 
      end do 
      !$omp end do 
     end do 
    !$omp end parallel 
    time2 = omp_get_wtime() 
    print *, a(1000000) 
    print *, 'Elapsed real time = ', time2 - time1, 'second(s)' 
end program openmp 

的MPI代码

program MPI 
implicit none 
include "mpif.h" 
integer, parameter       :: dp = selected_real_kind(15,307) 
integer         :: ierr, num_procs, my_id, destination, tag, source, stat, i, j 
real(kind = dp)       :: time1, time2 
real(kind = dp), dimension(1000000)  :: a 
    call MPI_INIT (ierr) 
    call MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr) 
    call MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr) 
    !Initialisation 
     do i = 1, 1000000 
      a(i) = sqrt(dble(i)/3.0d+0) 
     end do 
    destination = 0 
    tag = 999 
    source = 3 
    stat = MPI_STATUS_SIZE 
    time1 = MPI_Wtime() 
     do j = 1, 1000  
      do i = 1 + my_id, 1000000, num_procs 
       a(i) = a(i) + sqrt(dble(i)) 
      end do 
     end do 
    call MPI_BARRIER (MPI_COMM_WORLD, ierr) 
    if(my_id == source) then 
     call MPI_SEND (a(1000000), 1, MPI_DOUBLE_PRECISION, destination, tag, MPI_COMM_WORLD, ierr) 
    end if 
    if(my_id == destination) then 
     call MPI_RECV (a(1000000), 1, MPI_DOUBLE_PRECISION, source, tag, MPI_COMM_WORLD, stat, ierr) 
    end if 
    time2 = MPI_Wtime() 
    if(my_id == 0) then 
     print *, a(1000000) !, 'from ID =', my_id 
     print *, 'Elapsed real time = ', time2 - time1, 'second(s)' 
    end if 
    stop 
    call MPI_FINALIZE (ierr) 
end program MPI 

我使用Intel Fortran Compiler 17.0.3和-O0优化标志编译了这些代码。 OpenMP和MPI代码都是在4个核心Haswell Desktop上执行的。我分别获得了顺序OpenMP和MPI代码8.08s,2.1s和3.2s的CPU时间。实际上，我期待OpenMP和MPI代码之间的结果几乎相似;然而，事实并非如此。我的问题：

1. 关于MPI代码，如果我想打印出来的a(1000000)的结果，是有可能做到这一点在一个更聪明的方式，而不做这样call MPI_SEND和call MPI_RECV？

2. 您是否清楚MPI代码中哪些部分仍然可以优化？

3. 对于MPI代码中的source，是否可以自动定义它？在这种情况下，很容易对我来说，因为处理器的数量是4，所以a(1000000)必须被分配到线程3

预先感谢您。

Answer 1

最后，我得到了我的问题的解决方案。以前，我没有意识到parallelising的方式做循环串行代码：

do i = 1, 1000000 
    a(i) = a(i) + sqrt(dble(i)) 
end do 

是环状分布在MPI代码：

do i = 1 + my_id, 1000000, num_procs 
    a(i) = a(i) + sqrt(dble(i)) 
end do 

的问题。我认为这是因为更多缓存未命中发生。因此，与其环状分布，我申请嵌段分布到MPI代码，这是更有效的（对于这种情况!!!）。我现在写修订MPI代码：

program Revised_MPI 
use mpi 
implicit none 
integer, parameter       :: dp = selected_real_kind(15,307), array_size = 1000000 
integer          :: ierr, num_procs, my_id, ista, iend, i, j 
integer, dimension(:), allocatable   :: ista_idx, iend_idx 
real(kind = dp)        :: time1, time2 
real(kind = dp), dimension(:), allocatable :: a 

    call MPI_INIT (ierr) 
    call MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr) 
    call MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr) 

    !Distribute loop with block distribution 
    call para_range (1, array_size, num_procs, my_id, ista, iend) 
    allocate (a(ista : iend), ista_idx(num_procs), iend_idx(num_procs)) 

    !Initialisation and saving ista and iend 
     do i = ista, iend 
      a(i) = sqrt(dble(i)/3.0d+0) 
      ista_idx(my_id + 1) = ista 
      iend_idx(my_id + 1) = iend 
     end do 

    time1 = MPI_Wtime() 

    !Performing main calculation for all processors (including master and slaves) 
    do j = 1, 1000  
     do i = ista_idx(my_id + 1), iend_idx(my_id + 1) 
      a(i) = a(i) + sqrt(dble(i)) 
     end do 
    end do 
    call MPI_BARRIER (MPI_COMM_WORLD, ierr) 

    time2 = MPI_Wtime() 

    if(my_id == num_procs - 1) then 
     print *, a(array_size) 
     print *, 'Elapsed real time = ', time2 - time1, 'second(s)' 
    end if 

    call MPI_FINALIZE (ierr) 
    deallocate (a) 
end program Revised_MPI 

!----------------------------------------------------------------------------------------- 
subroutine para_range (n1, n2, num_procs, my_id, ista, iend) 

implicit none 

integer          :: n1, n2, num_procs, my_id, ista, iend, & 
                iwork1, iwork2 

    iwork1 = (n2 - n1 + 1)/num_procs 
    iwork2 = mod(n2 - n1 + 1, num_procs) 
    ista = my_id * iwork1 + n1 + min(my_id, iwork2) 
    iend = ista + iwork1 - 1 
    if(iwork2 > my_id) then 
     iend = iend + 1 
    end if 

end subroutine para_range 
!----------------------------------------------------------------------------------------- 

现在，MPI代码可以实现A（N）（几乎）类似的CPU时间与OpenMP的的。此外，它适用于优化标志-O3和-fast的使用。

谢谢大家的帮助。 :)

Answer 2

我发现在SUBROUTINE或FUNCTION中通常需要做更多的工作才能使并行性得到回报，所以在这个例子中，关注矢量化是最好的方法。
的绰号是“Vecorize内 - Parallelise外”（VIPO）
对于第二种情况我建议如下：

MODULE MyOMP_Funcs 
IMPLICIT NONE 
PRIVATE 

integer, parameter, PUBLIC   :: dp = selected_real_kind(15,307) 
real(kind = dp), dimension(1000000) :: a 

PUBLIC MyOMP_Init, MyOMP_Sum 

CONTAINS 

!================================= 
SUBROUTINE MyOMP_Init(N,A) 
IMPLICIT NONE 

integer      , INTENT(IN ) :: N 
real(kind = dp), dimension(n), INTENT(INOUT) :: A 

integer          :: I 

!Initialisation 
DO i = 1, n 
    A(i) = sqrt(dble(i)/3.0d+0) 
ENDDO 

RETURN 
END SUBROUTINE MyOMP_Init 


!================================= 
SUBROUTINE MyOMP_Sum(N,A,SumA) 
!$OMP DECLARE SIMD(MyOMP_Sum) UNIFORM(N,SumA) linear(ref(A)) 
USE OMPLIB 
IMPLICIT NONE 

integer      , INTENT(IN ) :: N 
!DIR$ ASSUME_ALIGNED A: 64     :: A 
real(kind = dp), dimension(n), INTENT(IN ) :: A 
real(kind = dp)    , INTENT( OUT) :: SumA 

integer          :: I 

SumA = 0.0 

!Maybe also try... !DIR$ VECTOR ALWAYS 

!$OMP SIMD REDUCTION(+:SumA) 
Sum_Loop: DO i = 1, N 
    SumA = SumA + A(i) + sqrt(dble(i)) 
ENDDO Sum_Loop 
!$omp end !<-- You probably do not need these 

RETURN 
END SUBROUTINE MyOMP_Sum 

!================================= 
SUBROUTINE My_NOVEC_Sum_Sum(N,A,SumA) 
IMPLICIT NONE 

integer      , INTENT(IN ) :: N 
!DIR$ ASSUME_ALIGNED A: 64     :: A 
real(kind = dp), dimension(n), INTENT(IN ) :: A 
real(kind = dp)    , INTENT( OUT) :: SumA 

integer          :: I 

SumA = 0.0 

!DIR$ NOVECTOR 
Sum_Loop: DO i = 1, N 
    SumA = SumA + A(i) + sqrt(dble(i)) 
ENDDO Sum_Loop 

RETURN 
END SUBROUTINE My_NOVEC_Sum 

!================================= 
END MODULE MyOMP_Funcs 
!================================= 


!================================= 
program openmp 
!USE OMP_LIB 
USE MyOMP_Funcs 
implicit none 

integer  , PARAMETER   :: OneM = 1000000 
integer  , PARAMETER   :: OneK = 1000 
integer        :: i, j 
real(kind = dp)      :: time1, time2, omp_get_wtime 
!DIR$ ATTRIBUTES ALIGNED:64   :: A, SumA 
real(kind = dp), dimension(OneM) :: A 
real(kind = dp)      :: SumA 
!Initialisation 

CALL MyOMP_Init(N,A) 
time1 = omp_get_wtime() 

! !$omp parallel 
! do j = 1, OneK 

CALL MyOMP_Sum(OneM, A, SumA) 

! end do 
! !$omp end parallel 
!!--> Put timing loops here 

time2 = omp_get_wtime() 
print *, a(1000000) 
print *, 'Elapsed real time = ', time2 - time1, 'second(s)' 

end program openmp 

一旦你有一个SIMD还原运行的版本，那么你可以尝试分层上并行性。
如果模块是库的一部分，那么编译器设置独立于程序。

Answer 3

其实，你的MPI程序对我来说没有多大意义。为什么所有的队伍都有相同的阵容？你为什么要复制整个阵列？为什么只在这个特定的源和目的地之间？

该程序没有计算任何有用的东西，所以很难判断什么是正确的程序（它没有正确计算任何有用的东西）。

在许多MPI程序中，您永远不会发送和接收整个阵列。甚至不是完整的本地数组，而只是它们之间的一些界限。

所以我想出了这个。请注意0​​，并且我从各处删除了magic number 1000000。

我也删除了stop。在end之前停止只是一个坏习惯，但它不是有害的。把它放在MPI_Finalize()之前是积极有害。

和最重要的是，我分配的工作不同。每个级别都有它自己的工作阵列的一部分。

program Test_MPI 

use mpi 

implicit none 

integer, parameter       :: dp = selected_real_kind(15,307) 
integer         :: ierr, num_procs, my_id, stat, i, j 
real(kind = dp)       :: time1, time2 
real(kind = dp), dimension(:), allocatable :: a 
integer, parameter       :: n = 1000000 
integer         :: my_n, ns 
    call MPI_INIT (ierr) 
    call MPI_COMM_RANK (MPI_COMM_WORLD, my_id, ierr) 
    call MPI_COMM_SIZE (MPI_COMM_WORLD, num_procs, ierr) 

    my_n = n/num_procs 
    ns = my_id * my_n 
    if (my_id == num_procs-1) my_n = n - ns 

    allocate(a(my_n)) 

    !Initialisation 
     do i = 1, my_n 
      a(i) = sqrt(real(i+ns, dp)/3.0d+0) 
     end do 

    stat = MPI_STATUS_SIZE 
    time1 = MPI_Wtime() 
     do j = 1, 1000  
      do i = 1 , my_n 
       a(i+my_id) = a(i) + sqrt(real(i+ns, dp)) 
      end do 
     end do 
    call MPI_BARRIER (MPI_COMM_WORLD, ierr) 

    time2 = MPI_Wtime() 
    if(my_id == 0) then 
     !!!! why??? print *, a(my_n) 
     print *, 'Elapsed real time = ', time2 - time1, 'second(s)' 
    end if 

    call MPI_FINALIZE (ierr) 
end program Test_MPI 

是的，那里没有沟通。我想不出为什么它应该在那里。如果它应该，你必须告诉我们为什么。它应该几乎完美地缩放。

也许你想收集最后一个排名？许多人都这样做，但通常根本不需要。目前尚不清楚为何需要您的情况。