#include <iostream>
#include <memory>
#include <mpi.h>

#define PRINT_ARRAY(STR,PTR) \
  for (int r=0; r<nproc; r++){\
    if (rank==r){\
      for (int j=0; j<nproc; j++){\
        (STR) << (PTR)[j] << " ";\
      }\
      (STR) << std::endl << std::flush;\
    }\
    MPI_Barrier(MPI_COMM_WORLD);\
  }

int main(int argc, char* argv[])
{
int provided;
MPI_Init_thread(&argc,&argv, MPI_THREAD_MULTIPLE, &provided);
if (provided!=MPI_THREAD_MULTIPLE){
  std::cerr << "Your MPI installation does not provide threading level MPI_THREAD_MULTIPLE, required for this program." << std::endl;
  MPI_Abort(MPI_COMM_WORLD, -1);
}
int rank, nproc;
MPI_Comm_rank(MPI_COMM_WORLD,&rank);
MPI_Comm_size(MPI_COMM_WORLD,&nproc);

  std::unique_ptr<double[]> M1(new double[nproc]);
  std::unique_ptr<double[]> M2(new double[nproc]);

  for (int j=0; j<nproc; j++) M1[j] = rank*nproc+j+1;

  PRINT_ARRAY(std::cout,M1);

std::unique_ptr<MPI_Request[]> requests(new MPI_Request[nproc]);

#pragma omp parallel for schedule(static) shared(requests)
for (int j=0; j<nproc; j++){
  MPI_Igather(M1.get()+j,  1,  MPI_DOUBLE,
              M2.get(), 1, MPI_DOUBLE, j,
              MPI_COMM_WORLD, &requests[j]);
}

  MPI_Waitall(nproc,requests.get(),MPI_STATUSES_IGNORE);

  if (rank==0) std::cout << " => "<<std::endl;
  PRINT_ARRAY(std::cout,M2);
  MPI_Finalize();
  return 0;
}