Suganthi Selvaraj
Suganthi Selvaraj

Reputation: 41

MPI_Gather 2-d arrays from all processors into bigger 2-d array on root with FORTRAN

Here's what I am trying to do with my code.

I have to calculate Fw and Fi on the 2-d grid that is of size nx by nz. I split the k loop between all processors so that each processor calculates nx by(nz/p) where p is the number of processors being used. After each processor is done I want to gather all the chunks, that is, each nx by nz/p Fw and Fi and put it into Fw and Fi in the root. I eventually want to use allgather i.e. Gather all the calculated Fw and Fi into all processors.

I have attached the code below.

I am not sure that I am specifying the sendcount and recvcount right or why my code is deadlocking. Any help is appreciated. Thanks!

PROGRAM gridtestpar
  IMPLICIT NONE

  INTEGER :: nx, nz, i, k, t
  INTEGER :: order, mx, mz
  INTEGER :: count
  INTEGER :: ierror, comm, p, rank, npr, s, f, np(2)

  REAL(KIND = 8) :: dx, dz, startx, startz, finishx, finishz
  REAL(KIND = 8) :: dt
  REAL(KIND = 8) :: cx, cz
  REAL(KIND = 8) :: cbx, cbz
  REAL(KIND = 8), ALLOCATABLE ::X(:), Z(:), Fw(:,:), Fi(:,:)
  REAL(KIND = 8), ALLOCATABLE :: Fn(:,:), Fnp1(:,:)



  include 'mpif.h'

  !----------------------------------------------------------
  !Parameters that can be changed
  !---------------------------------------------------------

  !Time step
  dt = 0.000000001d0
  !Number of points in x and z direction(i.e. streamwise and
  !spanwise) directions respectively
  nx = (400*5)
  nz = (400*5)

  !First and last grid point locations in x and z directions
  startx = 0.d0
  finishx = 60.d0*5.d0
  startz = 0.d0
  finishz = 60.d0*5.d0
  !Distance between grid points
  dx = (finishx-startx)/REAL(nx-1)
  dz = (finishz-startz)/REAL(nz-1)


  !Allocate
  ALLOCATE(X(nx),  Z(nz))
  ALLOCATE(Fw(nx,nz), Fi(nx,nz))
  ALLOCATE(Fn(nx,nz), Fnp1(nx,nz))


  ! Make Grid
  !--------------------------------------------------------------
  DO i = 1, nx
     X(i) = (i-1)*dx
  END DO

  DO k = 1, nz
     Z(k) = (k-1)*dx
  END DO

  CALL MPI_INIT(ierror)
  comm = MPI_COMM_WORLD
  !Get rank
  CALL MPI_COMM_RANK(comm, rank, ierror)
  !Get number of processors
  CALL MPI_COMM_SIZE(comm, p, ierror)


  !split job between all processors
  npr = INT((nz-1)/p)
  DO k = rank*npr+1, (rank+1)*npr
     DO i = 1, nx
        cx = 50.d0
        Fi(i,k) = 0.d0
        DO mx = 1,30
           cz = 0.d0;
           DO mz = 1,13*5
              Fi(i,k) = Fi(i,k) + EXP(-0.9d0*((X(i)-cx)**2+(Z(k)-cz)**2))
              cz = cz + 5.d0
           END DO
          cx = cx + 5.d0
        END DO
        cbz = 0.d0
        cbx = 30.d0
        DO mx = 1,4*5
           Fw(i,k) = Fw(i,k) + 0.05d0 + 7.d0*EXP(-0.1*((X(i)-cbx)**2 &
                + (Z(k)-cbz)**2)) + 0.1d0*Fi(i,k) 
           cbz = cbz + 20.d0
        END DO
     END DO
  END DO


  s = rank*npr+1
  f = (rank+1)*npr
  np(1) = nx
  np(2) = npr


  CALL MPI_GATHER(Fw(:,s:f), np , MPI_DOUBLE_PRECISION, &
       Fw,np , MPI_DOUBLE_PRECISION, 0,  comm, ierror)
  CALL MPI_GATHER(Fi(:,s:f), np , MPI_DOUBLE_PRECISION, &
       Fi,np , MPI_DOUBLE_PRECISION, 0, comm, ierror)

  Fn(:,:) = Fw(:,:) - Fi(:,:)
  Fnp1 = Fn

  WRITE(*,*) "I'm here"


  IF(rank == 0) THEN
     !Output initial condition
     !----------------------------------------------------------------
     OPEN(unit = 11, file = "Fiinitial.dat")
     WRITE(11,*) 'Variables = "X", "Z", "Fi"'
     WRITE(11,*) 'Zone I = ', nx, 'J = ', nz, 'F = POINT'
     DO k = 1, nz
        DO i = 1, nx
           WRITE(11,*) X(i), Z(k), Fi(i,k)
        END DO
     END DO
     WRITE(11,*) 'Zone I = ', nx, 'J = ', nz, 'F = POINT'
     DO k = 1, nz
        DO i = 1, nx
           WRITE(11,*) X(i), Z(k), Fw(i,k)
        END DO
     END DO
     CLOSE(11)
  END IF

  CALL MPI_FINALIZE(ierror)

END PROGRAM gridtestpar

Upvotes: 1

Views: 1981

Answers (1)

Bálint Aradi
Bálint Aradi

Reputation: 3812

You call the mpi_gather() subroutine incorrectly. You have to pass it the total nr. of elements which should be communicated as one integer for the send buffer and as an other integer for the receive buffer. You passed instead of each integer an array with two integers, which contained the number of elements along each dimension. Just multiply the numbers in your array, and pass the result as an integer instead:

program gridtestpar
  use mpi
  implicit none

  integer, parameter :: dp = kind(1.0d0)
  integer :: nx, nz
  integer :: ierror, comm, p, rank, npr, s, f, np(2)
  real(dp), allocatable :: Fw(:,:), Fi(:,:)

  nx = (400*5)
  nz = (400*5)

  allocate(Fw(nx,nz))
  allocate(Fi(nx,nz))
  Fw(:,:) = 0.0_dp
  Fi(:,:) = 0.0_dp

  call mpi_init(ierror)
  comm = MPI_COMM_WORLD
  call mpi_comm_rank(comm, rank, ierror)
  call mpi_comm_size(comm, p, ierror)

  s = rank * npr + 1
  f = (rank + 1) * npr

  call mpi_gather(Fw(:,s:f), nx * (f - s + 1), MPI_DOUBLE_PRECISION, &
       Fw, nx * npr, MPI_DOUBLE_PRECISION, 0, comm, ierror)
  call mpi_gather(Fi(:,s:f), nx * (f - s + 1), MPI_DOUBLE_PRECISION, &
       Fi, nx * npr, MPI_DOUBLE_PRECISION, 0, comm, ierror)
  write(*,*) "I'm here"
  call mpi_finalize(ierror)

end program gridtestpar

Maybe a few additional comments:

  • Please always post the shortest possible self containing code, which demonstrates the issue. Nobody likes to spend time on reading and trying to understand irrelevant code snippets. Leave away everything which is not essential for reproducing your problem. Maybe, this way you will even find the solution yourself.

  • Do not use kind = 8 when specifying precision. See the last part of this answer and some of the comments to it for alternatives.

  • You should use the mpi module instead of the include file.

Upvotes: 1

Related Questions