Parallelisierung mit OpenMP

OpenMP ist ein seit 1997 bestehender Standard mit Pragma-basierten Spracherweiterungen zu Fortran, C und C++, die eine Parallelisierung mit Hilfe von Threads unterstützt. Dies kann so vorgestellt werden, dass ein globaler Thread-Pool zur Verfügung gestellt wird und einzelne Code-Bereiche parallelisiert werden können. Dabei ist immer ein implizites Fork-And-Join gegeben, d.h. der normale Programmtext wird immer nur von einem Thread ausgeführt. Eine Parallelisierung findet nur an isolierten Stellen statt.

Die Webseiten des zugehörigen Standardisierungsgremiums mitsamt dem aktuellen Standard finden sich unter http://www.openmp.org/.

Beim gcc bzw. g++ werden die OpenMP-Erweiterung mit der Option -fopenmp unterstützt. Diese Option ist sowohl beim Übersetzen als auch beim Zusammenbau anzugeben. Wenn Funktionen der OpenMP-Bibliothek verwendet werden, wird #include <omp.h> benötigt. Zur Laufzeit kann die Größe des Thread-Pools mit Hilfe der Umgebungsvariablen OMP_NUM_THREADS beeinflusst werden.

In der allereinfachsten Variante, die wir zunächst betrachten werden, kann mit Hilfe von OpenMP eine for-Schleife parallelisiert werden:

template<typename MA, typename RandomEnginePool>
typename std::enable_if<hpc::matvec::IsRealGeMatrix<MA>::value, void>::type
mt_randomInit(MA& A, RandomEnginePool& repool) {
   using ElementType = typename MA::ElementType;
   using Index = typename MA::Index;
   using namespace hpc::aux;

   #pragma omp parallel for
   for (Index row = 0; row < A.numRows; ++row) {
      auto A_ = A(row, 0, 1, A.numCols);
      randomInit(A_, repool);
   }
}

Die Pragma-Anweisung #pragma omp parallel for parallelisiert die folgende for-Schleife und sorgt für ein implizites Join am Ende der Schleife. OpenMP parallelisiert dabei die Schleife entsprechend der Zahl der zur Verfügung stehenden Threads, wobei der Bereich der for-Schleife automatisch entsprechend aufgeteilt wird.

Die Schleifenvariable (hier: row) ist dann in jedem Thread privat, alle anderen Variablen werden jedoch per Voreinstellung gemeinsam verwendet. Wie bisher ist hier auf Konflikte zu achten.

Eine solche automatische Parallelisierung ist nur möglich, wenn die for-Schleife einer der kanonischen Formen entspricht:

\[ \mbox{for (index = start; index } \left\{ \begin{array}{c} \mbox{<} \\ \mbox{<=} \\ \mbox{>=} \\ \mbox{>} \end{array} \right\} \mbox{end; } \left\{ \begin{array}{l} \mbox{index++} \\ \mbox{++index} \\ \mbox{index--} \\ \mbox{--index} \\ \mbox{index += inc} \\ \mbox{index -= inc} \\ \mbox{index = index + inc} \\ \mbox{index = inc + index} \\ \mbox{index = index - inc} \end{array} \right\} \mbox{)}\]

Die Schleifenvariable darf dabei natürlich auch innerhalb der for-Schleife deklariert werden.

Zusätzliche Optionen sind bei dieser Pragma-Anweisung möglich. Wenn eine Aggregierung (wie beispielsweise eine Summe) von Teilergebnissen gewünscht wird, lässt sich das beispielsweise mit

   double sum = 0;
   #pragma omp parallel for reduction(+:sum)
   for (...) {
      sum += /* ... */;
   }
   /* total sum now in sum */

erreichen. Hierbei erhält jeder Thread seine eigene Kopie der Variablen sum, die zu Beginn den vorherigen Wert der zuvor deklarierten Variablen hat (hier: 0). Bei der impliziten Join-Operation werden diese einzelnen sum-Werte entsprechend dem bei reduction angegebenen Operator aggregiert und das Result dann automatisch der zuvor deklarierten Variable sum zugewiesen.

Aufgabe

Erweitern Sie die Vorlage wiederum mit der Berechnung der auf asum basierenden Absolutbeträge der Elemente einer Matrix und parallelisieren Sie diese ebenfalls mit Hilfe von OpenMP. So kann das Programm übersetzt werden:

$shell> g++ -O3 -g -I/home/numerik/pub/hpc/session18 -std=c++11 -fopenmp -o random_init15 random_init15.cpp

Vorlage

session18/random_init15.cpp

#include <cstdio>
#include <hpc/matvec/gematrix.h>
#include <hpc/matvec/apply.h>
#include <hpc/matvec/print.h>
#include <hpc/matvec/asum.h>
#include <hpc/aux/slices.h>
#include <hpc/aux/repool.h>

template<typename MA, typename Pool>
typename std::enable_if<hpc::matvec::IsRealGeMatrix<MA>::value, void>::type
randomInit(MA& A, Pool& pool) {
   using ElementType = typename MA::ElementType;
   using Index = typename MA::Index;
   using EngineType = typename Pool::EngineType;

   std::uniform_real_distribution<double> uniform(-100, 100);
   hpc::aux::RandomEngineGuard<EngineType> guard(pool);
   auto& engine(guard.get());

   hpc::matvec::apply(A, [&](ElementType& val, Index i, Index j) -> void {
      val = uniform(engine);
   });
}

template<typename MA, typename RandomEnginePool>
typename std::enable_if<hpc::matvec::IsRealGeMatrix<MA>::value, void>::type
mt_randomInit(MA& A, RandomEnginePool& repool) {
   using ElementType = typename MA::ElementType;
   using Index = typename MA::Index;
   using namespace hpc::aux;

   #pragma omp parallel for
   for (Index row = 0; row < A.numRows; ++row) {
      auto A_ = A(row, 0, 1, A.numCols);
      randomInit(A_, repool);
   }
}

int main() {
   using namespace hpc::matvec;
   using namespace hpc::aux;
   RandomEnginePool<std::mt19937> repool(4);

   GeMatrix<double> A(51, 7);
   mt_randomInit(A, repool);
}