Content

Aufteilung eines Indexbereiches

Es ist ein wiederkehrendes Problem, dass ein Indexbereich \([0, N)\) auf \(k\) Threads aufzuteilen ist. Insbesondere wenn die einzelnen Indizes mit Elementen eines im Speicher konsekutiv hintereinander angeordneten Arrays verbunden sind, dann es sinnvoll, die von den einzelnen Threads zu bearbeitenden Bereiche ebenfalls zusammenhängend zu belassen, damit Cache-Konflikte vermieden werden, bei denen zwei Threads auf die gleiche Cache-Line zugreifen, wobei mindestens einer der beiden Threads den Inhalt verändern möchte.

Dies könnte durchaus geschehen, wenn die Aufteilung anders erfolgen würde. Im folgenden Beispiel sind die Bereiche nicht zusammenhängend. Stattdessen erhält der erste Thread die Zeilen \(0, k, 2k, \dots\), der zweite Thread die Zeilen \(1, k+1, 2k+1, \dots\) usw.:

GeMatrix<double> A(N, N);
std::vector<std::thread> threads(nof_threads);
unsigned int job_size = A.numRows / nof_threads;
unsigned int remainder = A.numRows % nof_threads;
for (int index = 0; index < nof_threads; ++index) {
   GeMatrix<double>::Index num_rows = job_size;
   if (index < remainder) ++num_rows;
   auto A_ = GeMatrixView<double>(num_rows, A.numCols,
      &A(index, 0), A.incRow * nof_threads, A.incCol);
   threads[index] = std::thread([=]() mutable { randomInit(A_); });
}
for (int index = 0; index < nof_threads; ++index) {
   threads[index].join();
}

Da die Matrizen (per Voreinstellung) in ColMajor abgelegt sind, führt dies zu Konflikten und damit zu einer etwas längeren Zeit:

$shell> g++ -O3 -std=c++11 -o random_init7 -I/home/numerik/pub/hpc/session13 random_init7.cpp
$shell> ./random_init7
t1 = 0.191, t2 = 0.206
#include <thread>
#include <random>
#include <vector>
#include <hpc/matvec/gematrix.h>
#include <hpc/matvec/apply.h>
#include <hpc/matvec/print.h>
#include <hpc/aux/walltime.h>

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

   std::random_device random;
   std::mt19937 mt(random());
   std::uniform_real_distribution<ElementType> uniform(-100,100);

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

volatile double global;

int main() {
   using namespace hpc::matvec;
   using namespace hpc::aux;
   unsigned int nof_threads = std::thread::hardware_concurrency();

   WallTime<double> wall_time;

   constexpr unsigned int N = 5005;

   wall_time.tic();
   {
      GeMatrix<double> A(N, N);
      std::vector<std::thread> threads(nof_threads);
      unsigned int job_size = A.numRows / nof_threads;
      unsigned int remainder = A.numRows % nof_threads;
      GeMatrix<double>::Index start = 0;
      for (int index = 0; index < nof_threads; ++index) {
         GeMatrix<double>::Index num_rows = job_size;
         if (index < remainder) ++num_rows;
         auto A_ = A(start, 0, num_rows, A.numCols);
         threads[index] = std::thread([=]() mutable { randomInit(A_); });
         start += num_rows;
      }
      for (int index = 0; index < nof_threads; ++index) {
         threads[index].join();
      }
      global = A(7, 2);
   }
   double t1 = wall_time.toc();

   wall_time.tic();
   {
      GeMatrix<double> A(N, N);
      std::vector<std::thread> threads(nof_threads);
      unsigned int job_size = A.numRows / nof_threads;
      unsigned int remainder = A.numRows % nof_threads;
      for (int index = 0; index < nof_threads; ++index) {
         GeMatrix<double>::Index num_rows = job_size;
         if (index < remainder) ++num_rows;
         auto A_ = GeMatrixView<double>(num_rows, A.numCols,
            &A(index, 0), A.incRow * nof_threads, A.incCol);
         threads[index] = std::thread([=]() mutable { randomInit(A_); });
      }
      for (int index = 0; index < nof_threads; ++index) {
         threads[index].join();
      }
      global = A(7, 2);
   }
   double t2 = wall_time.toc();

   printf("t1 = %4.3lf, t2 = %4.3lf\n", t1, t2);
}

Aufgabe

Da es im Normalfall vorteilhaft erscheint, die Bereiche zusammenhängend zu belassen, lohnt es sich, hierfür eine kleine Hilfsklasse Slices anzulegen mit folgenden Anforderungen:

Dann könnte die Klasse folgendermaßen genutzt werden:

Slices slices(nof_threads, A.numRows);
for (int index = 0; index < nof_threads; ++index) {
   auto offset = slices.offset(index);
   auto size = slices.size(index);
   threads[index] = std::thread(
      /* lambda expression operating on [offset, offset+size) */
   );
}

Es könnte auch sinnvoll sein, Slices als Template-Klasse auszuführen mit einem wählbaren Typen für einen Wert aus dem Intervall.

Entwickeln Sie eine entsprechende Klasse und passen Sie die Initialisierung der Matrix so an, dass die Klasse hierfür verwendet wird.