Funktionen als Parameter

Bislang war die Funktion zur Initialisierung einer Matrix wenig flexibel:

template <typename T, typename Index>
void
initGeMatrix(Index m, Index n, T *A, Index incRowA, Index incColA)
{
    std::random_device                random;
    std::mt19937                      mt(random());
    std::uniform_real_distribution<T> uniform(-100, 100);

    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            A[i*incRowA+j*incColA] = uniform(mt);
        }
    }
}

In manchen Fällen wäre es angemessener, eine andere Verteilung zu wählen oder die Initialisierung deterministisch festzulegen. Gerade bei Tests und Benchmarks wird die Reproduzierbarkeit häufig vorausgesetzt. Somit liegt es nahe, die Initialisierung von der Bestimmung des zu initialisierenden Werts zu trennen.

Naheliegend wäre es hier, der initGeMatrix-Funktion eine Funktion mitzugeben, die als Parameter die Indizes i und j und die Dimensionen m und n erhält und \(A_{i,j}\) als Wert zurückliefert. Beginnen wir mit der ganz einfachen Initialisierung \(A_{i,j} \leftarrow j*n+i+1\):

template<typename T, typename Index>
T init_value(Index i, Index j, Index m, Index n) {
   return T(j)*T(n) + T(i) + T(1);
}

Der allgemeine Datentyp so einer Funktion ist dann T(Index, Index, Index, Index), d.h. sie erhält vier Parameter des Typs Index und liefert einen Wert des Typs T zurück. Eine Fassung von initGeMatrix, das so einen Funktionsparameter entgegennimmt, könnte so aussehen:

template <typename T, typename Index>
void
initGeMatrix(Index m, Index n, T *A, Index incRowA, Index incColA,
             T (*initValue)(Index, Index, Index, Index))
{
    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            A[i*incRowA+j*incColA] = initValue(i, j, m, n);
        }
    }
}

Zu beachten ist hier, dass initValue nicht die Funktion selbst, sondern ein Zeiger auf eine Funktion ist, so dass der Parametername mit dem Dereferenzierungsoperator * versehen ist.

Die entsprechende Funktion aus gematrix.h muss dann analog erweitert werden:

template <typename GeMatrix>
void
initGeMatrix(GeMatrix &A,
             typename GeMatrix::ElementType
                 (*initValue)(
                     typename GeMatrix::Index,
                     typename GeMatrix::Index,
                     typename GeMatrix::Index,
                     typename GeMatrix::Index))
{
    bench::initGeMatrix(A.m, A.n, A.data, A.incRow, A.incCol, initValue);
}

Der Aufruf könnte dann so aussehen:

GeMatrix<double> A(3, 7, StorageOrder::ColMajor);
initGeMatrix(A, init_value<double, std::size_t>);

Aufgabe

Wieso musste typename vor GeMatrix::Index angegeben werden?

Schreiben Sie analog zu init_value eine Variante init_complex, die für komplexe Zahlen funktioniert.

Die Parameterliste der beiden initGeMatrix-Funktionen wurde etwas unhandlich. Versuchen Sie dies mit Hilfe eines weiteren Template-Parameters zu vereinfachen.

Vorlage

session11/step01/bench.h

#ifndef HPC_BENCH_H
#define HPC_BENCH_H 1

#include <chrono>
#include <cmath>
#include <complex>
#include <random>

namespace bench {

template <typename T, typename Index>
T
asumDiffGeMatrix(Index m, Index n,
                 const T *A, Index incRowA, Index incColA,
                 T *B, Index incRowB, Index incColB)
{

    T asum = 0;

    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            asum += std::abs(B[i*incRowB+j*incColB] - A[i*incRowA+j*incColA]);
        }
    }
    return asum;
}

template <typename T, typename Index>
void
initGeMatrix(Index m, Index n, T *A, Index incRowA, Index incColA,
             T (*initValue)(Index, Index, Index, Index))
{
    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            A[i*incRowA+j*incColA] = initValue(i, j, m, n);
        }
    }
}

template <typename T>
struct WallTime
{
    void
    tic()
    {
        t0 = std::chrono::high_resolution_clock::now();
    }

    T
    toc()
    {
        using namespace std::chrono;

        elapsed = high_resolution_clock::now() - t0;
        return duration<T,seconds::period>(elapsed).count();
    }

    std::chrono::high_resolution_clock::time_point t0;
    std::chrono::high_resolution_clock::duration   elapsed;
};


} // namespace bench


#endif // HPC_BENCH_H

session11/step01/gematrix.h

#ifndef HPC_GEMATRIX_H
#define HPC_GEMATRIX_H 1

#include <cstddef>
#include <cassert>
#include "bench.h"

namespace matvec {

enum class StorageOrder {
    ColMajor,
    RowMajor
};

template <typename T, typename I>
    struct GeMatrixView;

template <typename T, typename I=std::size_t>
struct GeMatrix
{
    typedef T                       ElementType;
    typedef I                       Index;
    typedef GeMatrix<T,Index>       NoView;
    typedef GeMatrixView<T,Index>   View;

    GeMatrix(Index m, Index n, StorageOrder order=StorageOrder::ColMajor)
        : m(m), n(n),
          incRow(order==StorageOrder::ColMajor ? 1: n),
          incCol(order==StorageOrder::RowMajor ? 1: m),
          data(new T[m*n])
    {
    }

    ~GeMatrix()
    {
        delete[] data;
    }

    const ElementType &
    operator()(Index i, Index j) const
    {
        assert(i<m && j<n);
        return data[i*incRow + j*incCol];
    }

    ElementType &
    operator()(Index i, Index j)
    {
        assert(i<m && j<n);
        return data[i*incRow + j*incCol];
    }

    View
    operator()(Index i, Index j, Index numRows, Index numCols)
    {
        assert(i+numRows<=m);
        assert(j+numCols<=n);
        return View(numRows, numCols, &(operator()(i,j)), incRow, incCol);
    }

    const Index     m, n, incRow, incCol;
    ElementType*    data;
};

template <typename T, typename I=std::size_t>
struct GeMatrixView
{
    typedef T                       ElementType;
    typedef I                       Index;
    typedef GeMatrix<T,Index>       NoView;
    typedef GeMatrixView<T,Index>   View;

    GeMatrixView(Index m, Index n, T *data, Index incRow, Index incCol)
        : m(m), n(n), incRow(incRow), incCol(incCol), data(data)
    {
    }

    GeMatrixView(const GeMatrixView &rhs)
        : m(rhs.m), n(rhs.n), incRow(rhs.incRow), incCol(rhs.incCol),
          data(rhs.data)
    {
    }

    const ElementType &
    operator()(Index i, Index j) const
    {
        assert(i<m && j<n);
        return data[i*incRow + j*incCol];
    }

    ElementType &
    operator()(Index i, Index j)
    {
        assert(i<m && j<n);
        return data[i*incRow + j*incCol];
    }

    View
    operator()(Index i, Index j, Index numRows, Index numCols)
    {
        assert(i+numRows<=m);
        assert(j+numCols<=n);
        return View(numRows, numCols, &(operator()(i,j)), incRow, incCol);
    }

    const Index     m, n, incRow, incCol;
    ElementType*    data;
};

//
//  Interface for bench
//
template <typename GeMatrix>
void
initGeMatrix(GeMatrix &A,
             typename GeMatrix::ElementType
                 (*initValue)(
                     typename GeMatrix::Index,
                     typename GeMatrix::Index,
                     typename GeMatrix::Index,
                     typename GeMatrix::Index))
{
    bench::initGeMatrix(A.m, A.n, A.data, A.incRow, A.incCol, initValue);
}

template <typename GeMatrix>
typename GeMatrix::ElementType
asumDiffGeMatrix(const GeMatrix &A, const GeMatrix &B)
{
    return bench::asumDiffGeMatrix(A.m, A.n,
                                   A.data, A.incRow, A.incCol,
                                   B.data, B.incRow, B.incCol);
}

} // namespace matvec

#endif // HPC_GEMATRIX_H

session11/step01/test_initmatrix.cpp

#include <cstdlib>
#include <cstdio>
#include "gematrix.h"

template<typename T, typename Index>
T init_value(Index i, Index j, Index m, Index n) {
   return T(j)*T(n) + T(i) + T(1);
}

void print_value(double value) {
   std::printf(" %4.1lf", value);
}

template<typename Matrix>
void print_matrix(const Matrix& A) {
   for (std::size_t i = 0; i < A.m; ++i) {
      std::printf("  ");
      for (std::size_t j = 0; j < A.n; ++j) {
         print_value(A(i, j));
      }
      std::printf("\n");
   }
}

int main() {
   using namespace matvec;
   GeMatrix<double> A(3, 7, StorageOrder::ColMajor);
   initGeMatrix(A, init_value<double, std::size_t>);
   printf("A:\n");
   print_matrix(A);
}