Sample solution

#include <cassert>
#include <mpi.h>
#include <printf.hpp>
#include <hpc/matvec/gematrix.hpp>
#include <hpc/matvec/iterators.hpp>
#include <hpc/matvec/print.hpp>
#include <hpc/mpi/matrix.hpp>
#include <hpc/mpi/vector.hpp>

int main(int argc, char** argv) {
   MPI_Init(&argc, &argv);

   int nof_processes; MPI_Comm_size(MPI_COMM_WORLD, &nof_processes);
   int rank; MPI_Comm_rank(MPI_COMM_WORLD, &rank);
   assert(nof_processes == 2);

   using namespace hpc::matvec;
   using namespace hpc::mpi;
   using namespace hpc::aux;

   GeMatrix<double> A(3, 7, Order::ColMajor);
   GeMatrix<double> B(3, 7, Order::RowMajor);
   MPI_Datatype datatype_A = get_type(A);
   MPI_Datatype datatype_B = get_type(B);
   
   if (rank == 0) {
      MPI_Aint true_extent; MPI_Aint true_lb;
      MPI_Type_get_true_extent(datatype_A, &true_lb, &true_extent);
      MPI_Aint extent; MPI_Aint lb;
      MPI_Type_get_extent(datatype_A, &lb, &extent);
      auto size = sizeof(double) * A.numRows() * A.numCols();
      fmt::printf("true extent of A = %d\n", true_extent);
      fmt::printf("extent of A = %d\n", true_extent);
      fmt::printf("size = %zd\n", size);
      assert(true_extent == size);
   }

   if (rank == 0) {
      for (auto [i, j, Aij]: A) {
	 Aij = i * 100 + j;
      }
      MPI_Send(&A(0, 0), 1, datatype_A, 1, 0, MPI_COMM_WORLD);
      MPI_Status status;
      MPI_Recv(&B(0, 0), 1, datatype_B, 1, 0, MPI_COMM_WORLD, &status);
      for (auto [i, j, Aij]: A) {
	 if (Aij != B(i, j)) {
	    fmt::printf("verification failed for (%d,%d): %lg vs %lg\n",
	       i, j, Aij, B(i, j));
	 }
      }
   } else {
      MPI_Status status;
      MPI_Recv(&A(0, 0), 1, datatype_A, 0, 0, MPI_COMM_WORLD, &status);
      MPI_Send(&A(0, 0), 1, datatype_A, 0, 0, MPI_COMM_WORLD);
   }

   MPI_Finalize();
}

#ifndef HPC_MPI_MATRIX_HPP
#define HPC_MPI_MATRIX_HPP 1

#include <cassert>
#include <memory>
#include <type_traits>
#include <vector>
#include <mpi.h>
#include <hpc/aux/slices.hpp>
#include <hpc/matvec/gematrix.hpp>
#include <hpc/mpi/fundamental.hpp>

namespace hpc { namespace mpi {

/* construct MPI data type for a row of Matrix A
   where the extent is adapted such that consecutive
   rows can be combined even if they overlap */
template<typename T, template<typename> class Matrix,
   Require<Ge<Matrix<T>>> = true>
MPI_Datatype get_row_type(const Matrix<T>& A) {
   MPI_Datatype rowtype;
   MPI_Type_vector(
      /* count = */ A.numCols(),
      /* blocklength = */ 1,
      /* stride = */ A.incCol(),
      /* element type = */ get_type(A(0, 0)),
      /* newly created type = */ &rowtype);

   /* in case of row-major we are finished */
   if (A.incRow() == A.numCols()) {
      MPI_Type_commit(&rowtype);
      return rowtype;
   }

   /* the extent of the MPI data type does not match
      the offset of subsequent rows -- this is a problem
      whenever we want to handle more than one row;
      to fix this we need to use the resize function
      which allows us to adapt the extent to A.incRow() */
   MPI_Datatype resized_rowtype;
   MPI_Type_create_resized(rowtype, 0, /* lb remains unchanged */
      A.incRow() * sizeof(T), &resized_rowtype);
   MPI_Type_commit(&resized_rowtype);
   MPI_Type_free(&rowtype);
   return resized_rowtype;
}

/* create MPI data type for matrix A */
template<typename T, template<typename> class Matrix,
   Require<Ge<Matrix<T>>> = true>
MPI_Datatype get_type(const Matrix<T>& A) {
   MPI_Datatype datatype;
   if (A.incCol() == 1) {
      MPI_Type_vector(
	 /* count = */ A.numRows(),
	 /* blocklength = */ A.numCols(),
	 /* stride = */ A.incRow(),
	 /* element type = */ get_type(A(0, 0)),
	 /* newly created type = */ &datatype);
   } else {
      /* vector of row vectors */
      MPI_Datatype rowtype = get_row_type(A);
      MPI_Type_contiguous(A.numRows(), rowtype, &datatype);
      MPI_Type_free(&rowtype);
   }
   MPI_Type_commit(&datatype);
   return datatype;
}

} } // namespaces mpi, hpc

#endif // HPC_MPI_MATRIX_H