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