收集大小不等的数组时在 Fortran 中使用 MPI_GATHERV

Question

此题与this one, with the complication that the size of the matrices being collected is not equal in row length, but they are equal in column length. (I have also looked at this question and this one类似，但没弄明白）。

背景

我正在执行一个计算，在计算结束之前我不知道结果矩阵的行数。连续地，我分配了一个非常大的矩阵，它被填满，在计算结束时（当我知道行的限制时）我砍掉这个大数组的末尾，剩下我想要的结果。使用 MPI，我应用相同的逻辑：

在每个进程中，我都有一个大数组被填满
在计算结束时，我需要在每个数组各自的限制下截断它们（每个进程的限制都不同）
然后我需要将每个进程的结果数组收集到一个数组中，然后将其提供给根进程，以便继续执行程序的其余部分。

目前已尝试

为了理解 MPI_GATHERV 的工作原理以及如何在我的案例中使用它，我编辑了 this question 的答案中给出的代码，以便从每个进程中接受可变大小的数组.

program main
use mpi
implicit none
integer :: ierr, myRank, nProcs
integer :: sendsubarray, recvsubarray, resizedrecvsubarray
integer, dimension(2) :: starts,sizes,subsizes
integer, dimension(:), allocatable :: counts, disps
integer, parameter :: nx_glb=20, ny_glb=5, ny=5
integer :: nx
integer, dimension(:), allocatable :: nx_all  
character, dimension(:,:), target, allocatable :: mat, matG
character :: c
integer :: i, j
integer(kind=mpi_address_kind) :: start, extent

call mpi_init(ierr)
call mpi_comm_rank(mpi_comm_world, myRank, ierr)
call mpi_comm_size(mpi_comm_world, nProcs, ierr)
allocate(nx_all(nProcs))
nx_all = (/5, 4, 5, 5/)
nx = nx_all(myRank+1)  
sizes(1)=nx; sizes(2)=ny
subsizes(1)=nx; subsizes(2)=ny
starts(1)=0; starts(2)=0
call mpi_type_create_subarray(2, sizes, subsizes, starts,    mpi_order_fortran, &
                            mpi_character, sendsubarray, ierr)
call mpi_type_commit(sendsubarray,ierr)
allocate(mat(1:nx,1:ny))
mat='.'
forall (i=1:nx,j=1:ny) mat(i,j)=ACHAR(ICHAR('0')+myRank)

if(myRank.eq.0) then
 allocate(matG(nx_glb,ny_glb))
 matG='.'
 sizes(1)=nx_glb; sizes(2)=ny_glb
 subsizes(1)=nx; subsizes(2)=ny
 starts(1)=0; starts(2)=0
 call mpi_type_create_subarray(2, sizes, subsizes, starts, mpi_order_fortran, &
                               mpi_character, recvsubarray, ierr)
 call mpi_type_commit(recvsubarray, ierr)
 extent = sizeof(c)
 start = 0
 call mpi_type_create_resized(recvsubarray, start, extent, resizedrecvsubarray, ierr)
 call mpi_type_commit(resizedrecvsubarray,ierr)
end if

allocate(counts(4),disps(4))
counts(1:4) = (/1, 1, 1, 1/)
disps(1:4) = (/0, 5, 10, 15/)
call mpi_barrier(mpi_comm_world,ierr)
print *, mat, "process", myRank    
call mpi_gatherv(mat,1,sendsubarray,matG,counts,disps,resizedrecvsubarray, &
                 0,mpi_comm_world,ierr)
do p=0,nProcs
  if (myRank == p) then
     print *, 'Local array for rank ', myRank
     do i=1, nx
      print *, (mat(i,j),j=1,ny)
     end do
  endif
enddo
call MPI_Barrier(MPI_COMM_WORLD,ierr)
if(myRank.eq.0) then
 print * , matG, "process", myRank    
 print *, 'Global array: '
 do i=1, nx_glb
  print *, (matG(i,j),j=1, ny_glb)
 end do
end if

call mpi_finalize(ierr)

end program main

理想的结果（注意排名 1 如何少了一行）：

Local array for rank 0
00000
00000
00000
00000
00000
Local array for rank 1
11111
11111
11111
11111
Local array for rank 2
22222
22222
22222
22222
22222
Local array for rank 3
33333
33333
33333
33333
33333

Global array: 
00000
00000
00000
00000
00000
11111
11111
11111
11111
22222
22222
22222
22222
22222
33333
33333
33333
33333
33333

实际结果（注意在局部我有预期的行为，但在全局矩阵中，额外的行 1s 在那里，它在末尾有点）：

Local array for rank 0
00000
00000
00000
00000
00000
Local array for rank 1
11111
11111
11111
11111
Local array for rank 2
22222
22222
22222
22222
22222
Local array for rank 3
33333
33333
33333
33333
33333

Global array: 
00000
00000
00000
00000
00000
1111.
1111.
1111.
1111.
1111.
22222
22222
22222
22222
22222
33333
33333
33333
33333
33333

我了解到在内存中，矩阵被保存为数组，所以我得到的全局数组是这样的：

0000011111222223333300000111112222233333000001111122222333330000011111222223333300000.....2222233333

问题

如何删除点（代表排名 1 的空行）？如何让它显示为具有正确行数的矩阵？

编辑全局数组中出现额外行的原因是因为在根进程中创建的 recvsubarray 的维度为 5x5，尽管来自进程 1 的 sendsubarray 的维度为 4x5。现在的问题是，如何定义具有可变维度的 recvsubarray，具体取决于它从哪个等级接收信息？

Answer 1

您将全局矩阵定义为在第一维 (nx) 而不是第二维 (ny) 更大，这让您的生活变得非常艰难。 Fortran 存储数组的方式，具有较大的 ny 更自然，因为它对应于将所有子矩阵按顺序存储在内存中。

如果您乐于交换 nx 和 ny，那么您就不需要使用任何复杂的派生类型。事实上，我怀疑您是否可以使用 Scatterv 执行此模式，因为该函数需要单一接收类型，但您对每个传入矩阵都有不同的模式（由于您选择了 nx 和 ny 的顺序）。

这段代码，交换了 nx 和 ny，似乎工作正常。虚线在最后 - 我猜你总会有一些点，因为你分配的 space 比子矩阵占用的多。

program main
use mpi
implicit none
integer :: ierr, myRank, nProcs
integer :: sendsubarray, recvsubarray, resizedrecvsubarray
integer, dimension(2) :: starts,sizes,subsizes
integer, dimension(:), allocatable :: counts, disps
integer, parameter :: ny_glb=20, nx_glb=5, nx=5
integer :: ny
integer, dimension(:), allocatable :: ny_all
character, dimension(:,:), target, allocatable :: mat, matG
character :: c
integer :: i, j, p
integer(kind=mpi_address_kind) :: start, extent

call mpi_init(ierr)
call mpi_comm_rank(mpi_comm_world, myRank, ierr)
call mpi_comm_size(mpi_comm_world, nProcs, ierr)
allocate(ny_all(nProcs))
ny_all = (/5, 4, 5, 5/)
ny = ny_all(myRank+1)
allocate(mat(1:nx,1:ny))
mat='.'
forall (i=1:nx,j=1:ny) mat(i,j)=ACHAR(ICHAR('0')+myRank)

if(myRank.eq.0) then
 allocate(matG(nx_glb,ny_glb))
 matG='.'
end if

allocate(counts(4),disps(4))
counts(:) = nx*ny_all(:)

disps(1)=0
do i = 2, 4
  disps(i) = disps(i-1)+counts(i-1)
end do

call mpi_barrier(mpi_comm_world,ierr)
print *, mat, "process", myRank
call mpi_gatherv(mat,nx*ny,MPI_CHARACTER,matG,counts,disps,MPI_CHARACTER, &
                 0,mpi_comm_world,ierr)
do p=0,nProcs
  if (myRank == p) then
     print *, 'Local array for rank ', myRank
     do i=1, nx
      print *, (mat(i,j),j=1,ny)
     end do
  endif
enddo
call MPI_Barrier(MPI_COMM_WORLD,ierr)
if(myRank.eq.0) then
 print * , matG, "process", myRank
 print *, 'Global array: '
 do i=1, nx_glb
  print *, (matG(i,j),j=1, ny_glb)
 end do
end if

call mpi_finalize(ierr)

end program main

这是一些输出：

 Global array: 
 0000011112222233333.
 0000011112222233333.
 0000011112222233333.
 0000011112222233333.
 0000011112222233333.

000000000000000000000000011111111111111111111222222222222222222222222233333333333333333333333333.....

希望这有用。

Answer 2

我就是这样解决了上面的问题（答案是同事给我的，所以都归功于他）：

由于最终矩阵的维度之一是固定的，并且由于矩阵无论如何都存储在数组中，因此最好将 mpi_gather 与 mpi_type_vector 一起使用而不是 mpi_type_create_subarray .因此，程序结构如下：a) 确定每个秩中子矩阵的感兴趣大小，b) 将其转换为向量，c) 从每个秩收集向量，d) 将向量整形为最终矩阵。这样，就不需要收集不需要的信息（如上面问题中的点所示），因此不需要在使用 mpi_gather.

后摆脱它们

因此，为了将不同长度但宽度恒定的子数组收集到一个全局矩阵中，以下代码可以解决问题：

program main
  use mpi
  implicit none
  integer :: ierr, myRank, iProcs, nProcs, master
  integer :: ix, iy, ip
  integer :: nx, nxSum, offset, newtype
  integer, parameter :: ny=5
  integer, allocatable:: vec(:), vecG(:), nxAll(:), displs(:), rcounts(:), matG(:,:)

  call mpi_init(ierr)
  call mpi_comm_rank(mpi_comm_world, myRank, ierr)
  call mpi_comm_size(mpi_comm_world, nProcs, ierr)

  master = 0

  nx = myRank+1
  allocate(vec(nx*ny))
  do ix = 1,nx
     do iy = 1,ny
        ip = (ix-1)*ny + iy
        vec(ip) = myRank
     enddo
  enddo

  call mpi_barrier(mpi_comm_world,ierr)

  allocate(nxAll(nProcs))
  call mpi_gather(nx, 1, mpi_integer, nxAll, 1, mpi_integer, &
       master, mpi_comm_world, ierr)

  if (myRank == master) then
     ! print *, 'nxAll = ', nxAll, 'sum(nxAll) = ',sum(nxAll)
     nxSum = sum(nxAll)        
     allocate(vecG(nxSum*ny))

     allocate(displs(nProcs),rcounts(nProcs))
     offset = 0
     do iProcs = 1,nProcs
        displs(iProcs) = offset
        rcounts(iProcs) = nxAll(iProcs)*ny
        offset = offset + rcounts(iProcs)
        ! print *,'iProcs',iProcs,'displs = ',displs(iProcs),'rcounts',rcounts(iProcs)
     enddo

  endif

  call mpi_type_vector(nx*ny, 1, 1, mpi_integer,newtype,ierr)
  call mpi_type_commit(newtype,ierr)

  call mpi_gatherv(vec,1,newtype,vecG,rcounts,displs,mpi_integer, &
       master,mpi_comm_world,ierr)    


  if (myRank == master) then
     print *, 'Global vector, vecG = ',vecG

     ! Reshape into matrix
     print *, 'Global matrix'
     allocate(matG(nxSum,ny))
     do ix = 1,nxSum
        do iy = 1,ny
           ip = (ix-1)*ny + iy
           matG(ix,iy) = vecG(ip)
        enddo
        print *, (matG(ix,iy),iy=1,ny)
     enddo

  endif

  call mpi_finalize(ierr)

end program main

收集大小不等的数组时在 Fortran 中使用 MPI_GATHERV

Using MPI_GATHERV in Fortran when gathering unequal sized arrays

parallel-processing

fortran

mpi