Vorbereitungen

Für diese Session gibt es eine neue Fassung der Vorlesungsbibliothek, die unter /home/numerik/pub/hpc/session13/hpc öffentlich auf unseren Rechnern zur Verfügung steht. Auf den Rechnern in unserem Pool in E.44 und auf der Thales können Sie die Bibliothek unmittelbar verwenden, wenn beim Übersetzen zusätzlich die Option -I/home/numerik/pub/hpc/session13 angegeben wird. Da die Vorlesungsbibliothek nur aus Header-Dateien besteht, muss keine Bibliothek beim Zusammenbau des Programms angegeben werden.

In dieser Fassung der Vorlesungsbibliothek stehen bereits einige Erweiterungen zur Verfügung. Ein Beispiel hierfür sind die apply-Template-Funktionen aus <hpc/matvec/apply.h>:

/home/numerik/pub/hpc/session13/hpc/matvec/apply.h

#ifndef HPC_MATVEC_APPLY_H
#define HPC_MATVEC_APPLY_H 1

#include <type_traits>
#include <hpc/matvec/isgematrix.h>

namespace hpc { namespace matvec {

template <typename MA, typename Func>
typename std::enable_if<IsGeMatrix<MA>::value,
         void>::type
apply(MA &A, Func func)
{
    typedef typename MA::Index    Index;

    if (A.incRow<A.incCol) {
        for (Index j=0; j<A.numCols; ++j) {
            for (Index i=0; i<A.numRows; ++i) {
                func(A(i,j), i, j);
            }
        }
    } else {
        for (Index i=0; i<A.numRows; ++i) {
            for (Index j=0; j<A.numCols; ++j) {
                func(A(i,j), i, j);
            }
        }
    }
}

template <typename MA, typename Func>
typename std::enable_if<IsGeMatrix<MA>::value,
         void>::type
apply(const MA &A, Func func)
{
    typedef typename MA::Index    Index;

    if (A.incRow<A.incCol) {
        for (Index j=0; j<A.numCols; ++j) {
            for (Index i=0; i<A.numRows; ++i) {
                func(A(i,j), i, j);
            }
        }
    } else {
        for (Index i=0; i<A.numRows; ++i) {
            for (Index j=0; j<A.numCols; ++j) {
                func(A(i,j), i, j);
            }
        }
    }
}

} } // namespace matvec, hpc

#endif // HPC_MATVEC_APPLY_H

Nun muss nicht jedes mal erneut entschieden werden, ob es für den Cache günstiger ist, zuerst die Zeilen oder die Spalten zu durchlaufen. Den beiden apply-Methoden wird nur eine Matrix übergeben und ein Funktionsobjekt und letzteres wird dann auf allen Elementen zusammen mit den jeweiligen Indizes aufgerufen. Die Template-Funktion apply kommt hier in zwei Varianten, die sich nur in der const-ness von des Matrix-Parameters unterscheiden. Wie üblich erfolgt die Auswahl automatisch durch den Übersetzer in Abhängigkeit der const-ness der übergebenen Matrix.

Aufgabe

In der beigefügten Vorlage nutzt randomInit noch nicht die apply-Template-Funktion. Stattdessen durchläuft es selbst die Matrix in der jeweils günstigeren Weise. Vereinfachen Sie die Funktion, indem Sie apply verwenden. So können Sie die Vorlage übersetzen:

$shell> g++ -O3 -std=c++11 -o random_init1 -I/home/numerik/pub/hpc/session13 random_init1.cpp
$shell> ./random_init1
A = 
       92.0      -65.2      -97.4      -87.5      -13.1       27.3      -62.3       93.9
       -5.5       47.1      -99.3       20.9       35.6      -22.3      -30.7        9.8
       66.4      -49.1       34.4       65.4        3.4      -55.5       68.8      -79.4
       98.1       14.9      -48.5      -59.6      -74.1      -33.2      -47.3       95.3
        5.4       37.0        2.5       32.4      -34.1        8.3       12.4       80.1
       59.1      -70.5       59.5       -4.3       63.5      -43.3        9.4      -58.3
      -57.0      -56.8      -40.5       61.3      -29.2       69.8       76.9       -1.2

Vorlage

session13/random_init1.cpp

#include <random>
#include <hpc/matvec/gematrix.h>
#include <hpc/matvec/print.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);

   if (A.incRow < A.incCol) {
      for (Index j = 0; j < A.numCols; ++j) {
         for (Index i = 0; i < A.numRows; ++i) {
            A(i, j) = uniform(mt);
         }
      }
   } else {
      for (Index i = 0; i < A.numRows; ++i) {
         for (Index j = 0; j < A.numCols; ++j) {
            A(i, j) = uniform(mt);
         }
      }
   }
}

int main() {
   using namespace hpc::matvec;
   GeMatrix<double> A(7, 8);
   randomInit(A);
   print(A, "A");
}