=========================================
Boundaries and extent of an MPI data type                 [TOC]
=========================================

Assume that we are not that fortunate to have a matrix in
row-major order with `A.incCol() == 1`. In this case, we can
construct a data type in two steps:

 * Create a vector type for a row
 * Create a matrix type as vector of vectors

To do this, we need to understand how MPI interpretes
vectors in non-trivial cases. It is obvious if the
elements are stored consecutively in memory. In case of
a column-major matrix with `A.incRow() == 1` the rows
are not consecutively in memory. Instead they overlap
with offsets.

To understand how MPI sees this, we should return
to the definition of a data type as sequence of tuples:
$T = \left\{ (et_1, o_1), (et_2, o_2), \dots, (et_n, o_n) \right\}$.

 * We define $\mbox{extent}(T)$ for every data type $T$.
   In case of elementary types, this is the same value
   as returned by `sizeof`. Example:
   $\mbox{extent}(\mbox{MPI_DOUBLE}) = 8$.

 * Next we define a _lower bound_ $\mbox{lb}$:
   $\mbox{lb}(T) = \min\limits_{1 \le i \le n} \left\{ o_i \right\}$.

 * Likewise we have an _upper bound_ $\mbox{ub}$ which additionally
   considers the extent:
   $\mbox{ub}(T) = \max\limits_{1 \le i \le n} \left\{ o_i + \mbox{extent}(et_i) \right\} + \epsilon$.
   The variable $\epsilon$ is to be seen as a necessary round-up to
   the next alignment boundary. (In cases of vectors and matrices where
   we work with one elementary type only we can assume $\epsilon$ to be
   0.)

 * Now we can define $\mbox{extent}$ generally:
   $\mbox{extent}(T) = \mbox{ub}(T) - \mbox{lb}(T)$.

Assume we have a $2 \times 3$ matrix in column-major order.
The data type $R$ for the first row would then look as follows:
$R = \left\{ (DOUBLE, 0), (DOUBLE, 16), (DOUBLE, 32) \right\}$.
Then we have $\mbox{lb}(R) = 0$, $\mbox{ub}(R) = 40$, and
$\mbox{extent}(R) = \mbox{ub}(R) - \mbox{lb}(R) = 40$.
If we would combine two such rows to one vector, the absolute
offset of the first element of the second row vector would be 40
but actually it starts at 8.

The MPI library allows us to divert from the default
by using the function `MPI_Type_create_resized`:

---- CODE (type=cpp) ----------------------------------------------------------
int MPI_Type_create_resized(MPI_Datatype oldtype, MPI_Aint lb,
                            MPI_Aint extent, MPI_Datatype* newtype);
-------------------------------------------------------------------------------

Here `lb` specifies the new lower boundary $\mbox{lb}(T)$ and
`extent` specifies $\mbox{extent}(T)$ in bytes. From this
we get implicitly $\mbox{ub}(T)$. It is permitted to chose
an extent in a way where two "consecutively" following objects
actually overlap. In the example above of a column-major matrix
we could define the extent to be of 8 bytes as
the offset between `&A(0, 0)` and `&A(1, 0)` is just 8 bytes.

In the trivial case that all elements of a vector are consecutively
in memory, we can also use `MPI_Type_contiguous`:

---- CODE (type=cpp) ----------------------------------------------------------
int MPI_Type_contiguous(int count, MPI_Datatype oldtype,
                        MPI_Datatype* newtype);
-------------------------------------------------------------------------------

Exercise
========

 * Develop a function named `get_row_type` which constructs
   the MPI data type for the row of a matrix where the
   function `MPI_Type_create_resized` is to be used such
   that the data type can be used for consecutively following
   rows of a matrix. `lb` shall remain 0.

 * Develop a function named `get_type` for matrices
   where in case of `A.incCol() != 1` the function
   `get_row_type` is to be used to construct the data type
   as vector of row vectors.

   ---- CODE (type=cpp) -------------------------------------------------------
   template<typename T, template<typename> typename Matrix,
      Require<Ge<Matrix<T>>> = true>
   MPI_Datatype get_type(const Matrix<T>& A) {
      MPI_Datatype datatype;
      /* ... */
      MPI_Type_commit(&datatype);
      return datatype;
   }
   ----------------------------------------------------------------------------

 * Compare the extent of your matrix type with the size
   of the original matrix:

   ---- CODE (type=cpp) -------------------------------------------------------
   GeMatrix<double> A(m, n);
   MPI_Datatype datatype_A = get_type(A);

   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();

   assert(extent == size && true_extent == size);
   ----------------------------------------------------------------------------

   `MPI_Type_get_extent` returns the extent which has
   been possibly tampered with (as shown above) while
   `MPI_Type_get_true_extent` delivers the true extent.

 * Create a small test program that is to be executed with
   two processes that exchange
   matrices. Use it to test multiple configurations. Test,
   for example, the transfer of a row-major matrix to a
   column-major matrix and vice versa.

:navigate: up    -> doc:index
           back  -> doc:session06/page04
           next  -> doc:session06/page06