0

I want to scatterv a 3D array in MPI Fortran using derived data types. I know how to do it by flattening an array as shown below. Is is possible to scatter without flattening the array (by using derived data types)

I have provided an example where I Scatter the (5,3,2) 3D array among let's say 2 procs. First by flattening the 3D array. Then scatter flattened amongst processors. Now I want to achieve the same without flattening the array. So that I avoid creating temporary arrays for storing flattened global and local arrays.

I have also included the use of MPI_Type_create_subarray to scatter sub-arrays. But it works only if sub-array sizes are equal. Can someone help me with how to do this when subarray sizes are different?

It works if I distribute (4,3,2) among two procs as (2,3,2) and (2,3,2). If I want to distribute (5,3,2) among two procs as (3,3,2) and (2,3,2) it doesn't work. Can someone help with that?

!Example code to Scatter 3D array using MPI subarray type 
program ex_scatterv
  use mpi  
  use iso_fortran_env, only : real64 
  implicit none
  
  real(real64), allocatable,dimension(:,:,:) :: array, array_local
  real(real64), allocatable,dimension(:) :: array_flat, array_local_flat
  integer :: rank, num_procs, mpierr, i, j, k
  integer :: nx, ny, nz, str_idx, end_idx, local_size, local_size_flat 
  integer, dimension(:), allocatable :: sendcounts, displacements
  integer :: sizes(3), sub_sizes(3), starts(3), recv_starts(3), recv_sizes(3), &
    send_type, resize_send_type, recv_type
  integer(kind=MPI_OFFSET_KIND) :: lb, extent
  call mpi_init(mpierr)
  call mpi_comm_size(mpi_comm_world, num_procs, mpierr)
  call mpi_comm_rank(mpi_comm_world, rank, mpierr)
  
  nx=5
  ny=3
  nz=2

  !allocate in the root rank
  if(rank==0) then
    allocate(array(nx,ny,nz))
    allocate(array_flat(nx*ny*nz))
  else !for other procs allocate with zero size
    allocate(array(0,0,0))
  endif
 
  !assign values to the array
  if(rank==0) then
    do k=1,nz
      do j=1,ny
        do i=1,nx
          array(i,j,k) = (i-1)+(j-1)*nx+(k-1)*nx*ny
        end do
      end do 
    end do 
    print*, "Before scattering..."
    print*, array   
    !flatten the 3D array
    forall(k=1:nz, j=1:ny, i=1:nx) array_flat(k+(j-1)*nz+(i-1)*ny*nz)=array(i,j,k)
  endif

  !distribute the 3d array among different procs 
  call distribute_points(nx, rank, num_procs, str_idx, end_idx)
  local_size = end_idx - str_idx + 1
  local_size_flat = local_size*ny*nz
  
  !allocate local(for each rank) arrays
  allocate(array_local_flat(local_size_flat))
  allocate(array_local(local_size, ny, nz))

  !allocate sendcoutns and displacements arrays for braodcasting
  allocate(sendcounts(num_procs), displacements(num_procs))

  !gather displacements and sendcounts for all ranks
  call MPI_Allgather(str_idx, 1, MPI_INTEGER, displacements, 1, MPI_INTEGER, &
     MPI_COMM_WORLD, mpierr)
  call MPI_Allgather(local_size, 1, MPI_INTEGER, sendcounts, 1, &
    MPI_INTEGER, MPI_COMM_WORLD, mpierr)
  
  !total sendcounts and displacements 
  sendcounts = sendcounts*ny*nz
  displacements = displacements - 1 !Array index starts with 0 in MPI (C)
  displacements = displacements*ny*nz

  !scatter the flattened array among procs 
  call MPI_Scatterv(array_flat, sendcounts, displacements, MPI_DOUBLE_PRECISION, &
    array_local_flat, local_size*ny*nz, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, &
    mpierr)
  
  !form 3D array from flattened local array 
  forall(k=1:nz, j=1:ny, i=1:local_size) array_local(i,j,k) = &
      array_local_flat(k+(j-1)*nz+(i-1)*ny*nz)
  
  print*, "Scattered array: ", rank
  print*, array_local 

  !Scatterning using subarray type 
  sizes = [nx, ny, nz]
  sub_sizes = [local_size, ny, nz]
  starts = [0, 0, 0] 
  recv_starts = [0, 0, 0] 

  !to get extent of MPI_DOUBLE_PRECISION
  call MPI_Type_get_extent(MPI_DOUBLE_PRECISION, lb, extent, mpierr)

  call MPI_Type_create_subarray(3, sizes, sub_sizes, starts, &
    MPI_ORDER_FORTRAN, MPI_DOUBLE_PRECISION, send_type, mpierr)

  call MPI_Type_create_resized(send_type, 0, local_size*extent, &
    resize_send_type, mpierr)
  call MPI_Type_commit(resize_send_type, mpierr)

  call MPI_Type_create_subarray(3, sub_sizes, sub_sizes, recv_starts, &
    MPI_ORDER_FORTRAN, MPI_DOUBLE_PRECISION, recv_type, mpierr)
  call MPI_Type_commit(recv_type, mpierr)
 
  sendcounts=1
 
  call MPI_Allgather(rank, 1, MPI_INTEGER, displacements, 1, MPI_INTEGER, &
     MPI_COMM_WORLD, mpierr)
  call MPI_Scatterv(array, sendcounts, displacements, resize_send_type, array_local, &
    sendcounts, recv_type, 0, MPI_COMM_WORLD, mpierr)

  print*, "Scattered array with subarray: ", rank
  print*, array_local 

  call MPI_Finalize(mpierr)


  contains 
  subroutine distribute_points(npts, rank, size, start_idx, end_idx)
    implicit none
    
    integer, intent(in) :: npts, size, rank
    integer, intent(out) :: start_idx, end_idx
    integer :: pts_per_proc

    pts_per_proc = npts/size
    
    if(rank < mod(npts, size)) then
      pts_per_proc=pts_per_proc + 1
    end if
   
    if(rank < mod(npts, size)) then
      start_idx = rank * pts_per_proc + 1
      end_idx = (rank + 1) * pts_per_proc
    else
      start_idx = mod(npts, size) + rank*pts_per_proc + 1
      end_idx   = mod(npts, size) + (rank + 1) * pts_per_proc
    end if

  end subroutine distribute_points

end program ex_scatterv

I want to do it without flattening and creating temporary arrays.

nhm
  • 104
  • 7
  • I'm not totally sure I understand what you are trying to do - showing your "flattened" version would help a lot in making it clear. – Ian Bush Jun 03 '23 at 20:57
  • The easiest way for us to understand is for you to edit this question a **self-tested** [mcve] (leave blanks for the derived datatype creation and `MPI_Scatterv()` invokation). – Gilles Gouaillardet Jun 04 '23 at 03:03
  • @GillesGouaillardet I have added a working example where I scatter by flattening the array. – nhm Jun 04 '23 at 06:27
  • 1
    Do not use `real(kind=8)`, it is even longer than `kind(rk)` where `rk` is some named constant that you set appropriately. The value 8 is not portable, will not be available in some compilers and `MPI_REAL8` is not guaranteed to be compatible with it. – Vladimir F Героям слава Jun 04 '23 at 06:40
  • 2
    Did you try `MPI_type_create_subarray` (https://stackoverflow.com/questions/13345147/how-to-use-mpi-type-create-subarray) at all? Are there some particular points that you do not understand? It is really best to try it yourself or at least read the manual and ask for clarifications. Otherwise we would just repeat the work of authors of those manuals or of many tutorials available and still we would not know if the information will help you. – Vladimir F Героям слава Jun 04 '23 at 06:43
  • 1
    yeah, better use `DOUBLE_PRECISION` in Fortran! I did the homework and had the same idea of using `MPI_Type_create_subarray()`. You can create derived datatype describing a subarray for `(1,1:3,1:2)` and then resize it so the next type starts at `(2,1,1)` – Gilles Gouaillardet Jun 04 '23 at 07:01
  • I have tried using MPI_Type_create_subarray() but it works only for equal-sized subarrays, i.e., I can distribute (4,3,2) among 2 procs as (2,3,2) and (2,3,2) but if I have (5,3,2), I can't distribute as (3,3,2) and (2,3,2). – nhm Jun 04 '23 at 19:34
  • You can achieve the expected result with sendcounts=[3,2] – Gilles Gouaillardet Jun 04 '23 at 22:40
  • I tried. It doesn't work. I have also added how I have used MPI_tyep_create_subarray. Can someone help me if I am making some mistakes? – nhm Jun 05 '23 at 05:23
  • 1
    you really need to simplify your program and make it self-tested first. For example, initialize the `(5,3,2)` send array with `0` and `1` so rank `0` should receive a `(3,3,2)` array full of `0`, and rank `1` should receive a `(2,3,2)` array full of `1`. I do not understand why `nx=9` also things will be much easier if the first sub-size is `1`. note `recv_type` will be differents on ranks `0` and `1` – Gilles Gouaillardet Jun 05 '23 at 05:48
  • 1
    And as explained by Vladimir, use `DOUBLE PRECISION` and `MPI_DOUBLE_PRECISION` – Gilles Gouaillardet Jun 05 '23 at 05:48
  • 1
    I do not think it is necessary to use `double precision`, `real(real64)` should be portably compliant with `MPI_REAL8` as well. There must some misunderstanding with how you use the derived types. I hope someone will help you in the meantime, I have to go to the airport. – Vladimir F Героям слава Jun 05 '23 at 06:07
  • 2
    `MPI_REAL8` is an optional datatype, that is why I generally strongly discourage its use. My apologies for putting my words into Vladimir's mouth. – Gilles Gouaillardet Jun 05 '23 at 06:40
  • 1
    OK, I see. double precision is fine semantically, but I would prefer it in the kind syntax using a named constant with the right kind value, e.g. `kind(1.d0)`. – Vladimir F Героям слава Jun 05 '23 at 08:18
  • 1
    Regarding the main problem, I must admit I did not think it fully and for the Scatter operation. it will be tricky. Thank you for showing your subarray attempt, it is really helpful. I will think a bit more. The subarray is good for point to point, but for this collective less so. – Vladimir F Героям слава Jun 05 '23 at 08:26
  • @VladimirFГероямслава What's the better solution for this problem? I mean is there a way to scatter without flattening arrays? – nhm Jun 05 '23 at 08:32
  • yes there is, read my comments and if you are still stuck, share your best effort – Gilles Gouaillardet Jun 05 '23 at 09:56
  • Hi, @GillesGouaillardet the code that I shared is the best I could do. It will be of great help if you can help me – nhm Jun 06 '23 at 07:06
  • 1
    In order to build `sendtype`, you can `call mpi_type_create_subarray(3, [5,3,2], [1,3,2], [0,0,0], mpi_order_fortran, MPI_DOUBLE_PRECISION, tmptype, ierr)` and then resize `tmptype` with `extent=8`. And do similar thing for `recvtype` – Gilles Gouaillardet Jun 07 '23 at 03:47
  • Sorry, I didn't understand. Right now in the code that I have I created sub_array_sizes as [3,3,2] and [2,3,2] with sendcounts=[1,1] and displacements=[0,1] in terms of extent. But it only works for equal chunks. I am not sure, I understood your comments correctly. – nhm Jun 07 '23 at 19:20
  • 1
    It is easier if you use `sendcounts = [3,2]` – Gilles Gouaillardet Jun 08 '23 at 01:09
  • 1
    @GillesGouaillardet It worked with your suggestions and help. Thank you so much! – nhm Jun 08 '23 at 06:41

1 Answers1

0

Here is the full code with Gatherv

!Example code to Scatterv and Gatherv 3D array using MPI subarray 
!derived data type.

!For example, a 3D array of size (5,3,2) can be distributed
!to rank 0 and rank 1 as (3,3,2) and (2,3,2) or to 3 ranks 
!(2,3,2), (2,3,2), (1,3,2), etc. 

program ex_scatterv
  use mpi  
  use iso_fortran_env, only : real64 
  implicit none
  
  !allocate arrays
  real(real64), allocatable,dimension(:,:,:) :: array, array_local
  real(real64), allocatable,dimension(:) :: array_flat, array_local_flat
  integer :: rank, num_procs, i, j, k
  integer :: nx, ny, nz, str_idx, end_idx, local_size, local_size_flat 
  integer, dimension(:), allocatable :: sendcounts, displacements
  integer :: sizes(3), sub_sizes(3), starts(3), recv_starts(3), recv_sizes(3), &
    send_type, resize_send_type, recv_type, resize_recv_type
  integer(kind=8) :: lb, extent, lb_resize
  real(real64) :: start_time
  integer :: mpierr 
  call mpi_init(mpierr)
  call mpi_comm_size(mpi_comm_world, num_procs, mpierr)
  call mpi_comm_rank(mpi_comm_world, rank, mpierr)
  
  !size of array
  nx=5
  ny=3
  nz=2
  

  if(rank==0) then 
    if(num_procs>nx) then 
      print*, "Number of procs should be less than or equal to first dimension of the array"
      call MPI_Abort(mpi_comm_world, 1, mpierr)
    endif
  endif

  start_time=MPI_Wtime()
  !allocate in the root rank
  if(rank==0) then
    allocate(array(nx,ny,nz))
    allocate(array_flat(nx*ny*nz))
  else !for other procs allocate with zero size
    allocate(array(0,0,0))
  endif
 
  !assign values to the array
  if(rank==0) then
    do k=1,nz
      do j=1,ny
        do i=1,nx
          array(i,j,k) = (i-1)+(j-1)*nx+(k-1)*nx*ny
        end do
      end do    
    end do 
    print*, "Before scattering..."
    print*, array   
    !flatten the 3D array
    forall(k=1:nz, j=1:ny, i=1:nx) array_flat(k+(j-1)*nz+(i-1)*ny*nz)=array(i,j,k)
  endif

  !distribute the 3d array among different procs 
  call distribute_points(nx, rank, num_procs, str_idx, end_idx)
  local_size = end_idx - str_idx + 1
  local_size_flat = local_size*ny*nz
  
  !allocate local(for each rank) arrays
  allocate(array_local_flat(local_size_flat))
  allocate(array_local(0:local_size+1, ny, nz))

  !allocate sendcoutns and displacements arrays for braodcasting
  allocate(sendcounts(num_procs), displacements(num_procs))

  !gather displacements and sendcounts for all ranks
  call MPI_Allgather(str_idx, 1, MPI_INTEGER, displacements, 1, MPI_INTEGER, &
     MPI_COMM_WORLD, mpierr)
  call MPI_Allgather(local_size, 1, MPI_INTEGER, sendcounts, 1, &
    MPI_INTEGER, MPI_COMM_WORLD, mpierr)
  
  !total sendcounts and displacements 
  sendcounts = sendcounts*ny*nz
  displacements = displacements - 1 !Array index starts with 0 in MPI (C)
  displacements = displacements*ny*nz

  !scatter the flattened array among procs 
  call MPI_Scatterv(array_flat, sendcounts, displacements, MPI_DOUBLE_PRECISION, &
    array_local_flat, local_size*ny*nz, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, &
    mpierr)
  
  !form 3D array from flattened local array 
  forall(k=1:nz, j=1:ny, i=1:local_size) array_local(i,j,k) = &
      array_local_flat(k+(j-1)*nz+(i-1)*ny*nz)
  
  print*, "Scattered array: ", rank
  print*, array_local 
  !if(rank==0) then 
  !  print*, "Time taken by flatten and scatter: ", MPI_Wtime()-start_time 
  !endif

  call MPI_Barrier(mpi_comm_world, mpierr)
  !deallocate(array_flat, array_local_flat)

  start_time=MPI_Wtime()
  !Scatterning using subarray type 
  sizes = [nx, ny, nz]
  recv_sizes=[local_size+2, ny, nz]
  sub_sizes = [1, ny, nz]
  starts = [0, 0, 0] 
  recv_starts = [1, 0, 0] 

  !to get extent of MPI_DOUBLE_PRECISION
  call MPI_Type_get_extent(MPI_DOUBLE_PRECISION, lb, extent, mpierr)
  
  !create a mpi subarray data type for sending data
  call MPI_Type_create_subarray(3, sizes, sub_sizes, starts, &
    MPI_ORDER_FORTRAN, MPI_DOUBLE_PRECISION, send_type, mpierr)
  lb_resize=0
  !resize the send subarray for starting at correct location for next send
  call MPI_Type_create_resized(send_type, lb_resize, extent, &
    resize_send_type, mpierr)
  call MPI_Type_commit(resize_send_type, mpierr)

  !create a mpi subarray data type for receiving data
  call MPI_Type_create_subarray(3, recv_sizes, sub_sizes, recv_starts, &
    MPI_ORDER_FORTRAN, MPI_DOUBLE_PRECISION, recv_type, mpierr)
  
  !resize the receive subarray for starting at correct location for next receive
  call MPI_Type_create_resized(recv_type, lb_resize, extent, &
    resize_recv_type, mpierr)
  call MPI_Type_commit(resize_recv_type, mpierr)
 
  !sendcounts and displacement for sending and receiving subarrays 
  sendcounts=sendcounts/(ny*nz)
  displacements = displacements/(ny*nz)
 
  !if(rank==0) then 
  !  print*, "Time taken for creating MPI type subarrays: ", MPI_Wtime()-start_time 
  !endif

  call MPI_Barrier(mpi_comm_world, mpierr)
  start_time=MPI_Wtime()
  !scatter the subarrays 
   call MPI_Scatterv(array, sendcounts, displacements, resize_send_type, &
    array_local, sendcounts, resize_recv_type, 0, MPI_COMM_WORLD, mpierr)

  !if(rank==0) then 
  !  print*, "Time taken for scattering using MPI type subarrays: ", MPI_Wtime()-start_time 
  !endif
  call MPI_Barrier(mpi_comm_world, mpierr)
  !print the scattered array 
  print*, "Scattered array with subarray: ", rank
  print*, array_local 
  
  !do some computations on the scattered local arrays
  array_local = array_local+1
 
  call MPI_Barrier(mpi_comm_world, mpierr)
  start_time=MPI_Wtime()
  !Gather the local arrays to global (array) using the same subarrays 
  call MPI_Gatherv(array_local, local_size, resize_recv_type, array, &
    sendcounts, displacements, resize_send_type, 0, MPI_COMM_WORLD, mpierr)
  
  !if(rank==0) then 
  !  print*, "Time taken by MPI_Type_create_subarray Gathering: ", MPI_Wtime()-start_time 
  !endif
  
  if(rank==0) then 
    print*, "Gathered array: ------------------"
    print*, array
  endif
  call MPI_Finalize(mpierr)
end program ex_scatterv
nhm
  • 104
  • 7