使用 MPI 求解泊松方程时代码卡住答案

【问题标题】：Code stuck when solving Poisson equation with MPI使用 MPI 求解泊松方程时代码卡住
【发布时间】：2016-12-10 05:12:02
【问题描述】：

我使用 MPI 实现了一个简单的 1D Poisson 方程并行求解器，因此请熟悉 MPI 库。我将代码设计为使用未指定数量的处理器（包括只有 1 个）运行。

代码在 1 或 2 个处理器上运行并产生良好的结果。但是，它会卡在带有 4 个处理器的 mpi_send 和 mpi_recv 调用上。因此，我预计我的鬼点交换的实现是错误的。

由于代码太大，这里无法包含，我只包含了 Jacobi 方案和数据交换：

do iter=1,max_iter                                                                                    
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~             
    ! Initial guess, on interior points only                                                          
    Ujacob(min_x+1:max_x-1) = 0._dp                                                                   
    Ujacob_all(0:grid_nx-1) = 0._dp                                                                   

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~             
    ! Store solution vector from last iteration                                                       
    Uold    (:) = Ujacob    (:)                                                                       
    Uold_all(:) = Ujacob_all(:)                                                                       

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~             
    ! Jacobi scheme                                                                                   
    do ii=min_x+1,max_x-1                                                                             
        !Ujacob(ii) = 0.5_dp * (Uold  (ii-1) + Uold  (ii+1) - grid_delta_x**2 * Urhs(ii))             
    end do                                                                                            

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~             
    ! Gather Ujacob vector                                                                            
    call mpi_allgather(Ujacob(0:proc_nx-1), proc_nx, mpi_float,                         &             
                    &  Ujacob_all,          proc_nx, mpi_float, mpi_comm_world, ierror)               

    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~             
    ! Compute error and check if less than tolerance value                                            
    error = sqrt((sum(Ujacob_all - Uold_all)**2) / dble(grid_nx))                                     

    if(error < error_tol) return                                                                      
    !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~             
    ! Exchange data points                                                                            

    ! Interior processors                                                                         
    if(Xsrc /= -1 .AND. Xdes /= -1) then
        call mpi_send(Ujacob(        0), 1, mpi_float, Xsrc, 200, mpi_comm_world, ierror)         
        call mpi_send(Ujacob(proc_nx-1), 1, mpi_float, Xdes, 100, mpi_comm_world, ierror)         

        call mpi_recv(Ujacob(     -1), 1, mpi_float, Xsrc, 100, mpi_comm_world, stat, ierror)     
        call mpi_recv(Ujacob(proc_nx), 1, mpi_float, Xdes, 200, mpi_comm_world, stat, ierror)
    ! First processor                                                                             
    elseif(Xsrc == -1) then                                                                    
        call mpi_send(Ujacob(proc_nx-1), 1, mpi_float, Xdes, 100, mpi_comm_world, ierror)         
        call mpi_recv(Ujacob(proc_nx  ), 1, mpi_float, Xdes, 200, mpi_comm_world, stat, ierror)
    ! Last processor                                                                              
    elseif(Xdes == -1) then                                                                   
        call mpi_send(Ujacob( 0), 1, mpi_float, Xsrc, 200, mpi_comm_world, ierror)                
        call mpi_recv(Ujacob(-1), 1, mpi_float, Xsrc, 100, mpi_comm_world, stat, ierror)          
    end if                                                                                                           
end do

Xsrc 和 Xdes 设置如下：

   !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
    ! Setting the source and destination neighbors of each processor                              
    if(myid == 0) then                                                                            
        Xsrc = -1                                                                                 
        Xdes = myid + 1                                                                           
    elseif(myid == nprocs-1) then                                                                 
        Xsrc = myid -1                                                                            
        Xdes = -1                                                                                 
    else                                                                                          
        Xsrc = myid - 1                                                                           
        Xsrc = myid + 1                                                                           
    end if

另外，我检查了处理器等级 0 和 nprocs-1 确实对应于左右边界处理器。

我已检查标签是否设置正确。另外，请随时评论您认为可以改进的任何内容。

【问题讨论】：

我认为看看你如何设置Xsrc 和Xdes 可能会有用。
@d_1999 是的，你是对的，请查看编辑
您的光环交换在概念上存在缺陷，因为您依赖 MPI_SEND 被缓冲，但情况可能并非总是如此。使用MPI_SENDRECV可以同时发送和接收，不会阻塞。另外，不要使用如此复杂的逻辑进行发送和接收。只需使用MPI_PROC_NULL 而不是-1 用于边界等级的不存在的邻居，并始终在两个方向上执行sendrecv。发送到MPI_PROC_NULL 或从它接收是无操作的。
我认为您在编辑的最后一个分支中有错字。
计算速度很快，但通信速度却不是这样。您不会看到小型数组的任何加速。

标签： parallel-processing fortran mpi

【解决方案1】：

@Hristo 是正确的，您的代码在原则上存在概念上的缺陷。但是，几乎每个 MPI 实现都会为包含单个实数的消息缓冲 MPI_Send（当然这不能保证），所以这不是您的代码的问题。

我认为您的标签不匹配 - 边缘情况应该将标签颠倒：

  elseif(Xsrc == -1) then                                                                    
        call mpi_send(Ujacob(proc_nx-1), 1, mpi_float, Xdes, 200, mpi_comm_world, ierror)         
        call mpi_recv(Ujacob(proc_nx  ), 1, mpi_float, Xdes, 100, mpi_comm_world, stat, ierror)
    ! Last processor                                                                              
    elseif(Xdes == -1) then                                                                   
        call mpi_send(Ujacob( 0), 1, mpi_float, Xsrc, 100, mpi_comm_world, ierror)                
        call mpi_recv(Ujacob(-1), 1, mpi_float, Xsrc, 200, mpi_comm_world, stat, ierror)          
    end if

代码上的其他几个cmets：

使用 allgather 计算误差项非常低效：您应该只对局部元素求和，然后使用 MPI_Allreduce 计算全局误差；
对于 Fortran 代码，您应该使用 MPI_REAL 而不是 MPI_FLOAT；
我看不出我们的代码如何在单个进程上运行 - 这里进程将执行第一个 elseif 子句，然后尝试发送到不存在的等级。

一旦你检查了你的标签现在是正确的，你就应该修复@Hristo指出的问题。

【讨论】：

其实标签是对的。主要问题是由于定义Xsrc 和Xdes 时的拼写错误。我修复了它以及我遇到的其他非高效问题 - 请参阅上面的 cmets。不过我确实有一个问题。 “几乎每个 MPI 实现都会为包含单个实数的消息缓冲 MPI_SEND 是什么意思？
问题是发送是同步的（即在匹配的接收发布之前不会完成）还是异步的（无论是否有匹配的接收都会完成）。 Ssend 保证同步，Bsend 保证异步。但是，Send 可以。这就是@HristoIliev 提出的观点——如果 Send 是同步的，那么你的代码将死锁，因为所有进程都在发送而没有接收。异步发送需要稍后进行复制和交付。对于小消息，MPI 通常会复制一份；对于大消息没有空间，所以它使用同步发送。