Lösungsvorschlag

Zum ersten Punkt: Hier kommt die Variable last_value hinzu, die wir in der capture als Referenz aufnehmen müssen:

std::random_device random;
std::mt19937 mt(random());
std::uniform_real_distribution<double> uniform(1, 1<<16);
double last_value = 0;
initGeMatrix(A, [=,&last_value]
      (std::size_t i, std::size_t j) mutable -> double {
   last_value += uniform(mt);
   return last_value;
});

Die für die zweite Aufgabe entwickelte Template-Funktion applyGeMatrix ist recht ähnlich zu initGeMatrix. Nur die Zuweisung wird hier ersetzt durch den Aufruf von apply:

template <typename GeMatrix, typename ApplyOperator>
void
applyGeMatrix(GeMatrix &A, ApplyOperator apply)
{
    for (typename GeMatrix::Index i = 0; i < A.m; ++i) {
        for (typename GeMatrix::Index j = 0; j < A.n; ++j) {
            apply(A(i, j), i, j);
        }
    }
}

So könnte dann eine Initialisierung mit applyGeMatrix aussehen:

std::random_device random;
std::mt19937 mt(random());
std::uniform_real_distribution<double> uniform(1, 1<<16);
double last_value = 0;
applyGeMatrix(B, [=,&last_value]
      (double& Bij, std::size_t i, std::size_t j) mutable {
   last_value += uniform(mt);
   Bij = last_value;
});

Und die Berechnung der Summe aller Elemente:

double sum = 0;
applyGeMatrix(B, [&sum]
      (double Bij, std::size_t i, std::size_t j) mutable {
   sum += Bij;
});
printf("sum of all elements of B: %10.2lf\n", sum);

Und schließlich der Betrag des betragsgrößten Elements:

double amax = 0;
applyGeMatrix(B, [&amax]
      (double& Bij, std::size_t i, std::size_t j) mutable {
   if (abs(Bij) > amax) amax = abs(Bij);
});
printf("absolute max of B: %10.2lf\n", amax);

Für die letzte Aufgabe musste zunächst die Template-Klasse GeMatrixCombinedConstView erstellt werden:

template <typename Matrix1, typename Matrix2>
struct GeMatrixCombinedConstView
{
    typedef typename std::common_type<typename Matrix1::ElementType,
        typename Matrix2::ElementType>::type ElementType;
    typedef typename std::common_type<typename Matrix1::Index,
        typename Matrix2::Index>::type Index;

    template <typename ApplyOperator>
    GeMatrixCombinedConstView(const Matrix1& A, const Matrix2& B,
            ApplyOperator apply)
        : m(A.m), n(A.n), A(A), B(B), apply(apply)
    {
        assert(A.m == B.m && A.n == B.n);
    }

    const ElementType
    operator()(Index i, Index j) const
    {
        return apply(A(i, j), B(i, j));
    }

    const Index m;
    const Index n;
    const Matrix1& A;
    const Matrix2& B;
    std::function<ElementType(ElementType, ElementType)> apply;
};

Mit Hilfe von GeMatrixCombinedConstView lässt sich dann die Summe der Beträge der Differenzen zweier Matrizen leicht berechnen:

GeMatrixCombinedConstView<GeMatrix<double>, GeMatrix<double>>
   C(A, B, [](double a, double b) -> double {
      return abs(a - b);
   });
printf("C:\n");
print_matrix(C);
{
   double sum = 0;
   applyGeMatrix(C, [&sum]
         (double Cij, std::size_t i, std::size_t j) mutable {
      sum += Cij;
   });
   printf("sum of all elements of C: %10.2lf\n", sum);
}

Hier sind alle Dateien zum Lösungsvorschlag:

#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, typename InitValue>
void
initGeMatrix(Index m, Index n, T *A, Index incRowA, Index incColA,
             InitValue initValue)
{
    for (Index j=0; j<n; ++j) {
        for (Index i=0; i<m; ++i) {
            A[i*incRowA+j*incColA] = initValue(i, j);
        }
    }
}

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
#ifndef HPC_GEMATRIX_H
#define HPC_GEMATRIX_H 1

#include <cstddef>
#include <cassert>
#include <functional>
#include <type_traits>
#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;
};

template <typename Matrix1, typename Matrix2>
struct GeMatrixCombinedConstView
{
    typedef typename std::common_type<typename Matrix1::ElementType,
        typename Matrix2::ElementType>::type ElementType;
    typedef typename std::common_type<typename Matrix1::Index,
        typename Matrix2::Index>::type Index;

    template <typename ApplyOperator>
    GeMatrixCombinedConstView(const Matrix1& A, const Matrix2& B,
            ApplyOperator apply)
        : m(A.m), n(A.n), A(A), B(B), apply(apply)
    {
        assert(A.m == B.m && A.n == B.n);
    }

    const ElementType
    operator()(Index i, Index j) const
    {
        return apply(A(i, j), B(i, j));
    }

    const Index m;
    const Index n;
    const Matrix1& A;
    const Matrix2& B;
    std::function<ElementType(ElementType, ElementType)> apply;
};

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

template <typename GeMatrix, typename ApplyOperator>
void
applyGeMatrix(GeMatrix &A, ApplyOperator apply)
{
    for (typename GeMatrix::Index i = 0; i < A.m; ++i) {
        for (typename GeMatrix::Index j = 0; j < A.n; ++j) {
            apply(A(i, j), i, j);
        }
    }
}

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
#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);
}

template<typename T, typename Index>
std::complex<T> init_complex_value(Index i, Index j, Index m, Index n) {
   return std::complex<T>(T(i), T(j));
}

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

void print_value(std::complex<double> value) {
   std::printf(" (%4.1lf, %4.1lf)", value.real(), value.imag());
}

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");
   }
}

template<typename T, typename Index>
struct IncreasingRandomValues {
   std::mt19937 mt;
   std::uniform_real_distribution<T> uniform;
   T last_value;

   IncreasingRandomValues() :
      mt(std::random_device()()), uniform(0, 1<<16), last_value(0) {
   }
   IncreasingRandomValues(unsigned int seed) :
      mt(seed), uniform(1, 1<<16) {
   }
   T operator()(Index i, Index j) {
      last_value += uniform(mt);
      return last_value;
   }
};

int main() {
   using namespace matvec;
   GeMatrix<double> A(3, 7, StorageOrder::ColMajor);
   {
      std::random_device random;
      std::mt19937 mt(random());
      std::uniform_real_distribution<double> uniform(1, 1<<16);
      double last_value = 0;
      initGeMatrix(A, [=,&last_value]
            (std::size_t i, std::size_t j) mutable -> double {
         last_value += uniform(mt);
         return last_value;
      });
   }
   printf("A:\n");
   print_matrix(A);

   GeMatrix<double> B(3, 7, StorageOrder::ColMajor);
   {
      std::random_device random;
      std::mt19937 mt(random());
      std::uniform_real_distribution<double> uniform(1, 1<<16);
      double last_value = 0;
      applyGeMatrix(B, [=,&last_value]
            (double& Bij, std::size_t i, std::size_t j) mutable {
         last_value += uniform(mt);
         Bij = last_value;
      });
   }
   printf("B:\n");
   print_matrix(B);
   {
      double sum = 0;
      applyGeMatrix(B, [&sum]
            (double Bij, std::size_t i, std::size_t j) mutable {
         sum += Bij;
      });
      printf("sum of all elements of B: %10.2lf\n", sum);
   }
   {
      double amax = 0;
      applyGeMatrix(B, [&amax]
            (double& Bij, std::size_t i, std::size_t j) mutable {
         if (abs(Bij) > amax) amax = abs(Bij);
      });
      printf("absolute max of B: %10.2lf\n", amax);
   }
   GeMatrixCombinedConstView<GeMatrix<double>, GeMatrix<double>>
      C(A, B, [](double a, double b) -> double {
         return abs(a - b);
      });
   printf("C:\n");
   print_matrix(C);
   {
      double sum = 0;
      applyGeMatrix(C, [&sum]
            (double Cij, std::size_t i, std::size_t j) mutable {
         sum += Cij;
      });
      printf("sum of all elements of C: %10.2lf\n", sum);
   }
}