using MPI_Gatherv for Fortran

Question

This question follows an existing thread on MPI_type_create_subarray and MPI_Gather. My aim is to gather subarrays of a bigger array from all the slave processes (4 in number) into a larger array on a master process (rank=0) using MPI_Type_Create_Subarray and MPI_Gatherv in Fortran 90. This would help me understand MPI_Gatherv for my other projects. Following is my sample code:

    program main
    implicit none
    include "mpif.h"
    integer :: ierr, myRank, nProcs
    integer :: sendsubarray, recvsubarray, resizedrecvsubarray
    integer, dimension(2) :: starts,sizes,subsizes
    integer, dimension(:), allocatable :: counts, disps
    integer, parameter :: nx_glb=10, ny_glb=10, nx=5, ny=5
    integer, dimension(:,:), target, allocatable :: mat, matG
    integer, pointer :: sendPtr(:,:), recvPtr(:,:)
    integer :: i, j

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

    sizes(1)=nx+2; sizes(2)=ny+2
    subsizes(1)=nx; subsizes(2)=ny
    starts(1)=2; starts(2)=2
    call mpi_type_create_subarray(2, sizes, subsizes, starts, mpi_order_fortran, &
                                mpi_integer, sendsubarray, ierr)
    call mpi_type_commit(sendsubarray,ierr)

    allocate(mat(1:nx+2,1:ny+2))
    do j=1, ny+2
     do i=1, nx+2
      if(i.eq.1 .or. i.eq.nx+2 .or. j.eq.1 .or. j.eq.ny+2) then
       mat(i,j)=1000
      else
       mat(i,j) = myRank
      end if
     end do
    end do

    sendPtr=>mat
    if(myRank.eq.0) then
     allocate(matG(nx_glb,ny_glb))
     matG=1000
     sizes(1)=nx_glb; sizes(2)=ny_glb
     subsizes(1)=nx; subsizes(2)=ny
     starts(1)=1; starts(2)=1
     call mpi_type_create_subarray(2, sizes, subsizes, starts, mpi_order_fortran, &
                                   mpi_integer, recvsubarray, ierr)
     call mpi_type_commit(recvsubarray, ierr)
     call mpi_type_create_resized(recvsubarray, 1, sizeof(i), resizedrecvsubarray, ierr)
     call mpi_type_commit(resizedrecvsubarray,ierr)
     recvPtr=>matG
    end if

    counts(1:4) = (/1, 1, 1, 1/)
    disps(1:4) = (/0, 5, 50, 55/)
    call mpi_gatherv(sendPtr,1,sendsubarray,recvPtr,counts,disps,resizedrecvsubarray, &
                     0,mpi_comm_world,ierr)

    if(myRank.eq.0) then
     do i=1, nx_glb
      write(1000,*) (matG(i,j),j=1, ny_glb)
     end do
    end if

    call mpi_finalize(ierr)

    end program main

However, executing this code results in forrtl: severe(174): SIGSEGV, segmentation fault occurred.

It seems like I am trying to point to a variable/location of an array that has not be initialized or declared while gathering. I tried to debug in many ways, but in vain.

Many thanks in advance.

Answer 1

You'll kick yourself when you see the main problem here; you didn't allocate counts or disps.

As an aside, I strongly recommend using use mpi rather than include mpif.h; the use statement (before the implicit none) brings in the F90 interface which has much better typechecking. When you do that, you'll also see that for your type create resized, you'll need integers of kind mpi_address_kind.

Update:

Ok, so for the larger question of how to do the gatherv, you had things largely right, but you're right, the starts, disps etc have to be zero-indexed, not 1, because the actual MPI library is doing things from a C point of view, even with the FORTRAN bindings. So for the sendsubarray, starts have to be [1,1]; for the recv subarray it has to be [0,0], and the resize, start has to be 0 and extent has to be sizeof(type) (and both those have to be integers of kind mpi_address_kind).

I'm attaching a version of your code with those updates, and with the underlying arrays being of type character so it's easier to print out diagnostics and see what's going on:

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=10, ny_glb=10, nx=5, ny=5
    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)

    sizes(1)=nx+2; sizes(2)=ny+2
    subsizes(1)=nx; subsizes(2)=ny
    starts(1)=1; starts(2)=1
    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+2,1:ny+2))
    mat='.'
    forall (i=2:nx+1,j=2:ny+1) 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, 50, 55/)
    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+2
          print *, (mat(i,j),j=1,ny+2)
         end do
      endif
      call MPI_Barrier(MPI_COMM_WORLD,ierr)
    enddo
    if(myRank.eq.0) then
     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

With output:

 Local array for rank            0
 .......
 .00000.
 .00000.
 .00000.
 .00000.
 .00000.
 .......
 Local array for rank            1
 .......
 .11111.
 .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: 
 0000022222
 0000022222
 0000022222
 0000022222
 0000022222
 1111133333
 1111133333
 1111133333
 1111133333
 1111133333

...make sense? This is very similar to the C version of this question which is answered here ( MPI_Type_create_subarray and MPI_Gather ) but you had already figured things mostly out...

Oh, yeah, one more thing -- you don't actually need to set up pointers to the send/recv data in Fortran. In C, you need to explicitly pass pointers to arrays of data; in fortran, you can just pass arrays (and they are already passed "by reference", eg the equivalent of C's passing pointers to a variable). So you can just pass the array.