================
Lösungsvorschlag
================

Der Lösungsvorschlag wird hier in einzelnen Schritten erklärt, weiter
unten findet sich die gesamte Lösung.

Zuerst ordnen wir alle Prozesse in einem zwei-dimensionalen Gitter,
hier `grid` genannt:

---- CODE (type=cpp) ----------------------------------------------------------
/* create two-dimensional Cartesian grid for our prcesses */
int dims[2] = {0, 0}; int periods[2] = {false, false};
MPI_Dims_create(nof_processes, 2, dims);
MPI_Comm grid;
MPI_Cart_create(MPI_COMM_WORLD,
   2,        // number of dimensions
   dims,     // actual dimensions
   periods,  // both dimensions are non-periodical
   true,     // reorder is permitted
   &grid     // newly created communication domain
);
MPI_Comm_rank(MPI_COMM_WORLD, &rank); // update rank (could have changed)
-------------------------------------------------------------------------------

Dann wird die große Matrix, hier $100 \times 100$ angelegt. Um das Programm
einfach zu halten, wird diese Matrix von allen Prozessen deklariert, obwohl
nur der Prozess mit _rank_ 0 diese benötigt.

---- CODE (type=cpp) ----------------------------------------------------------
/* initialize the entire matrix, including its borders */
Matrix A(100, 100, StorageOrder::RowMajor);
using Index = GeMatrix<double>::Index;
if (rank == 0) {
   apply(A, [&](double& val, Index i, Index j) -> void {
      if (j == 0) {
	 val = std::sin(PI * ((double)i/(A.numRows-1)));
      } else if (j == A.numCols - 1) {
	 val = std::sin(PI * ((double)i/(A.numRows-1))) * E_POWER_MINUS_PI;
      } else {
	 val = 0;
      }
   });
}
-------------------------------------------------------------------------------

Als nächstes werden zwei Matrizen angelegt, die dem jeweiligen Prozess
zugeordneten Teil der Gesamtmatrix entsprechen. Wir benötigen zwei
Matrizen, so dass wir Jacobi-Schritte von `B1` nach `B2` und danach
von `B2` nach `B1` durchführen können. Die Variable `overlap` legt
hier fest, wie dick die Ränder sind bzw. wie weit sich die
Matrizen überlappen. Beim Jacobi-Verfahren genügt eine Überlappung von 1.

---- CODE (type=cpp) ----------------------------------------------------------
int overlap = 1;
/* get our position within the grid */
int coords[2];
MPI_Cart_coords(grid, rank, 2, coords);
Slices<int> rows(dims[0], A.numRows - 2*overlap);
Slices<int> columns(dims[1], A.numCols - 2*overlap);
Matrix B1(rows.size(coords[0]) + 2*overlap,
	  columns.size(coords[1]) + 2*overlap,
	  StorageOrder::RowMajor);
Matrix B2(rows.size(coords[0]) + 2*overlap,
	  columns.size(coords[1]) + 2*overlap,
	  StorageOrder::RowMajor);
-------------------------------------------------------------------------------

Nun können die einzelnen Teilblöcke mitsamt den jeweils überlappenden
Rändern an die einzelnen Prozesse verteilt werden. Das übernimmt die
bereits erarbeitete Funktion `scatter_by_block`. Dabei müssen wir
darauf achten, dass bei Prozessen in Randlage auch der globale Rand von
`B1` nach `B2` kopiert wird, da dieser wegen der fehlenden Nachbarn
später nicht mehr verändert wird. Der Einfachheit verwenden wir
hier die `copy`-Funktion.

---- CODE (type=cpp) ----------------------------------------------------------
/* distribute main body of A include left and right border */
scatter_by_block(A, B1, 0, grid, dims, coords, overlap);
copy(B1, B2); /* actually just the border needs to be copied */
-------------------------------------------------------------------------------

Dann bestimmen wir unsere Nachbarn, wie zuvor vorgestellt:

---- CODE (type=cpp) ----------------------------------------------------------
/* get the process numbers of our neighbors */
int left, right, upper, lower;
MPI_Cart_shift(grid, 0, 1, &upper, &lower);
MPI_Cart_shift(grid, 1, 1, &left, &right);
-------------------------------------------------------------------------------

Für den Austausch mit den Nachbarn benötigen wir zwei Typen,
einen für Zeilen und einen für Spalten:

---- CODE (type=cpp) ----------------------------------------------------------
/* compute type for inner rows and cols without the border */
auto B_inner_row = B1(0, 1, overlap, B1.numCols - 2*overlap);
MPI_Datatype rowtype = get_type(B_inner_row);
auto B_inner_col = B1(0, 1, B1.numRows - 2*overlap, overlap);
MPI_Datatype coltype = get_type(B_inner_col);
-------------------------------------------------------------------------------

Um die Ein- und Ausgabe zu vereinfachen, lohnt es sich,
eine Datenstruktur anzulegen. Eine I/O-Operation wird
definiert durch den Kommunikationspartner (`rank`), die
Anfangsadresse des Datenbereichs, der entweder versandt
oder gefüllt wird (`addr`) und den zugehörigen Datentyp
(`datatype`), der entweder `rowtype` oder `coltype` ist.

---- CODE (type=cpp) ----------------------------------------------------------
struct IOTransfer {
   int rank;
   void* addr;
   MPI_Datatype datatype;
};
-------------------------------------------------------------------------------

Die I/O-Operationen teilen wir in Versand- und Empfangs-Operationen
auf, die entsprechend mit `MPI_Isend` bzw. `MPI_Irecv` umgesetzt
werden. In `B1_in_vectors` finden sich die Spezifikationen für
all die Empfangsoperationen (in den äußeren Rand), die in `B1` landen;
in `B1_out_vectors` sind umgekehrt die Ausgabe-Operationen (in den
inneren Rand). Die gleichen Datenstrukturen gibt es auch noch einmal
für `B2`. Wir benötigen diese doppelt, da wir sowohl von `B1` nach
`B2` als auch von `B2` nach `B1` einen Jacobi-Schritt durchführen
möchten. Die `B2`-Variante ergibt sich einfach aus der `B1`-Variante,
wobei überall `B1` durch `B2` ersetzt wird. (In C++14 ließe sich
diese Copy&Paste-Aktion mit Hilfe von Template-Variablen vermeiden,
die es in C++11 noch nicht gibt.)

---- CODE (type=cpp) ----------------------------------------------------------
IOTransfer B1_in_vectors[] = {
   {left, &B1(1, 0), coltype},
   {upper, &B1(0, 1), rowtype},
   {right, &B1(1, B1.numCols-overlap), coltype},
   {lower, &B1(B1.numRows-overlap, 1), rowtype},
};
IOTransfer B1_out_vectors[] = {
   {left, &B1(1, 1), coltype},
   {upper, &B1(1, 1), rowtype},
   {right, &B1(1, B1.numCols-2*overlap), coltype},
   {lower, &B1(B1.numRows-2*overlap, 1), rowtype},
};
-------------------------------------------------------------------------------

So sieht dann die Hauptschleife aus. Da die globale Synchronisierung
bei `MPI_Reduce` und `MPI_Bcast` nicht ganz billig ist, führen wir
dies hier pragmatisch nur bei jedem 10. Schritt durch. Die
eigentliche Arbeit liegt in `jacobi_iteration`, das nicht nur
die Ein- und Ausgangsmatrix erhält, sondern auch die jeweils
passende Datenstruktur für die Kommunikation:

---- CODE (type=cpp) ----------------------------------------------------------
double eps = 1e-6; unsigned int iterations;
for (iterations = 0; ; ++iterations) {
   double maxdiff = jacobi_iteration(B1, B2, B2_in_vectors, B2_out_vectors);
   maxdiff = jacobi_iteration(B2, B1, B1_in_vectors, B1_out_vectors);
   if (iterations % 10 == 0) {
      double global_max;
      MPI_Reduce(&maxdiff, &global_max, 1, get_type(maxdiff),
	 MPI_MAX, 0, MPI_COMM_WORLD);
      MPI_Bcast(&global_max, 1, get_type(maxdiff), 0, MPI_COMM_WORLD);
      if (global_max < eps) break;
   }
}

if (rank == 0) printf("%d iterations\n", iterations);
-------------------------------------------------------------------------------

So sieht die Signatur für eine Jacobi-Iteration aus. Die `in_vectors`
und `out_vectors` werden als Referenz übergeben, weil dies die
elegante Schleifenform ermöglicht. (Sonst würde es als Zeiger übergeben
werden, wobei die Information über die Länge verloren geht.)

---- CODE (type=cpp) ----------------------------------------------------------
template<typename Matrix>
typename Matrix::ElementType jacobi_iteration(const Matrix& A, Matrix& B,
      const IOTransfer (&in_vectors)[4],
      const IOTransfer (&out_vectors)[4]) {
   using ElementType = typename Matrix::ElementType;
   using Index = typename Matrix::Index;
   assert(A.numRows > 2 && A.numCols > 2);
   ElementType maxdiff = 0;

   /* ... */

   return maxdiff;
}
-------------------------------------------------------------------------------

Um den Programmtext für die Berechnung nicht wiederholt verwenden
zu müssen, wird eine kleine lokale Funktion hierfür definiert:

---- CODE (type=cpp) ----------------------------------------------------------
auto jacobi = [&](Index i, Index j) {
   B(i, j) = 0.25 *
      (A(i - 1, j) + A(i + 1, j) + A(i, j - 1) + A(i, j + 1));
   double diff = std::fabs(A(i, j) - B(i, j));
   if (diff > maxdiff) maxdiff = diff;
};
-------------------------------------------------------------------------------

Dann erfolgt der erste Schritt, bei dem wir unseren Rand berechnen:

---- CODE (type=cpp) ----------------------------------------------------------
/* compute the border first which is sent in advance
   to our neighbors */
for (Index i = 1; i + 1 < B.numRows; ++i) {
   jacobi(i, 1);
   jacobi(i, B.numCols-2);
}
for (Index j = 2; j + 2 < B.numCols; ++j) {
   jacobi(1, j);
   jacobi(B.numRows-2, j);
}
-------------------------------------------------------------------------------

Im zweiten Schritt initiieren wir die Kommunikation mit unseren
Nachbarn unter Verwendung der übergebenen Datenstrukturen:

---- CODE (type=cpp) ----------------------------------------------------------
/* send border to our neighbors and
   initiate the receive operations */ 
MPI_Request requests[8]; int ri = 0;
for (auto& in_vector: in_vectors) {
   MPI_Irecv(in_vector.addr, 1, in_vector.datatype,
      in_vector.rank, 0, MPI_COMM_WORLD, &requests[ri++]);
}
for (auto& out_vector: out_vectors) {
   MPI_Isend(out_vector.addr, 1, out_vector.datatype,
      out_vector.rank, 0, MPI_COMM_WORLD, &requests[ri++]);
}
-------------------------------------------------------------------------------

Im dritten Schritt berechnen wir den inneren Teil,
während parallel dazu im Hintergrund die Kommunikation
abgewickelt wird.

---- CODE (type=cpp) ----------------------------------------------------------
/* compute the inner block in parallel
   to the initiated communication */
for (Index i = 2; i + 2 < B.numRows; ++i) {
   for (Index j = 2; j + 2 < B.numCols; ++j) {
      jacobi(i, j);
   }
}
-------------------------------------------------------------------------------

Schließlich warten wir auf die Beendigung der
Kommunikation:

---- CODE (type=cpp) ----------------------------------------------------------
for (auto& request: requests) {
   MPI_Status status;
   MPI_Wait(&request, &status);
}
-------------------------------------------------------------------------------

Gesamte Lösung
==============

:import:session22/jacobi-2d.cpp

:navigate: up    -> doc:index
           back  -> doc:session22/page11