wip: Exercise 01

Exercise_01/Makefile

@@ -0,0 +1,20 @@
+
+CPP=g++
+MPICPP=mpic++
+CPPFLAGS= -std=c++17 -O3 -Wall -Wpedantic -pedantic
+OUT=main
+
+seriel: main.cpp Matrix.h Solver.h
+	$(CPP) main.cpp $(CPPFLAGS) -o $(OUT)_seriel
+
+mpi: main.cpp Matrix.h Solver.h
+	$(MPICPP) main.cpp $(CPPFLAGS) -DUSE_MPI -o $(OUT)_mpi
+
+all: seriel mpi
+
+test: seriel mpi
+	./$(OUT)_seriel 120 10
+	mpirun -n 4 ./$(OUT)_mpi 120 10
+
+clean:
+	rm -f *.o main $(OUT)_seriel $(OUT)_mpi

Exercise_01/Matrix.h

@@ -3,49 +3,132 @@
 #include <stddef.h>
 #include <ostream>
 #include <vector>
+#define _USE_MATH_DEFINES   /* enables math constants from cmath */
+#include <cmath>
 
-template <typename T>
+inline size_t mod(int num, size_t module) {
-class MatrixView;
+    while (num < 0) { num += module; };
+
+    return static_cast<size_t>(num % module);
+}
+
+enum Norm { Frob, Max };
 
 template <typename T>
 class Matrix {
 public:
     Matrix(size_t nrow, size_t ncol) : _nrow{nrow}, _ncol{ncol} {
-        _elem.reserve(nrow * ncol);
+        _data.reserve(nrow * ncol);
+    };
+    Matrix(size_t nrow, size_t ncol, T elem) : Matrix(nrow, ncol) {
+        for (size_t i = 0; i < nrow * ncol; ++i) {
+            _data[i] = elem;
+        }
     };
 
     size_t nrow() const { return _nrow; };
     size_t ncol() const { return _ncol; };
+    size_t nx() const { return _ncol; };
+    size_t ny() const { return _nrow; };
     size_t size() const { return _nrow * _ncol; };
 
-          T& operator()(int i)       { return _elem[index(i)]; };
+          T* data()       { return _data.data(); };
-    const T& operator()(int i) const { return _elem[index(i)]; };
+    const T* data() const { return _data.data(); };
-          T& operator()(int i, int j)       { return _elem[index(i, j)]; };
+
-    const T& operator()(int i, int j) const { return _elem[index(i, j)]; };
+          T& operator()(int i)       { return _data[index(i)]; };
+    const T& operator()(int i) const { return _data[index(i)]; };
+
+          T& operator()(int i, int j)       { return _data[index(i, j)]; }
+    const T& operator()(int i, int j) const { return _data[index(i, j)]; }
+
+    template <enum Norm N = Norm::Frob>
+    double norm() const {
+        T accum = static_cast<T>(0);    /*< result accumulator */
+
+        // Enforce valid norm types (at compile time)
+        static_assert(N == Norm::Frob || N == Norm::Max);
+
+        if constexpr (N == Norm::Frob) {
+            // Sum of Squares
+            for (const T& elem : _data) {
+                accum += elem * elem;
+            }
+            // The Frobenius norm is the square root of the sum of squares
+            return std::sqrt(static_cast<double>(accum));
+        } else { // Maximum Norm
+            // Maximum absolute value
+            for (const T& elem : _data) {
+                if (accum < std::abs(elem)) {
+                    accum = std::abs(elem);
+                }
+            }
+            return static_cast<double>(accum);
+        }
+    }
+
+    /**
+     * Distance of this Matrix to given matrix A as || this - A ||_N
+     *
+     * Note: It there is a dimension missmatch, only the number of elements
+     * of the smaller matrix are considured.
+     */
+    template <enum Norm N = Norm::Frob>
+    double dist(const Matrix<T>& A) const {
+        T accum = static_cast<T>(0);    /*< result accomulator */
+        size_t nelem = size() < A.size() ? size() : A.size();
+
+        // Enforce valid norm types (at compile time)
+        static_assert(N == Norm::Frob || N == Norm::Max);
+
+        if constexpr (N == Norm::Frob) {
+            // Sum of Squared differences
+            for (size_t i = 0; i < nelem; ++i) {
+                T diff = this(i) - A(i);
+                accum += diff * diff;
+            }
+            return std::sqrt(static_cast<double>(accum));
+        } else { // Maximum Norm
+            for (size_t i = 0; i < nelem; ++i) {
+                T diff = std::abs(this(i) - A(i));
+                if (accum < diff) {
+                    accum = diff;
+                }
+            }
+            return static_cast<double>(accum);
+
+        }
+    }
+
+    /** overloaded standard swap for matrices (of same type)
+     *
+     * @example
+     * Matrix<int> A(1, 2);
+     * Matrix<int> B(3, 4);
+     * std::swap(A, B);
+     */
+    friend void swap(Matrix<T>& A, Matrix<T>& B) {
+        std::swap(A._nrow, B._nrow);
+        std::swap(A._ncol, B._ncol);
+        std::swap(A._data, B._data);
+    }
 
 private:
-    size_t _nrow;
+    size_t _nrow;           /*< Number of Rows */
-    size_t _ncol;
+    size_t _ncol;           /*< Number of Columns */
-    std::vector<T> _elem;
+    std::vector<T> _data;   /*< Matrix elements / Data */
 
     size_t index(int i) const {
-        int nelem = static_cast<int>(_nrow * _ncol);
+        const int nelem = static_cast<int>(_nrow * _ncol);
+
         while (i < 0)      { i += nelem; }
         while (i >= nelem) { i -= nelem; }
 
         return i;
     };
+
     size_t index(int i, int j) const {
-        int nrow = static_cast<int>(_nrow);
+        return static_cast<size_t>(mod(i, _nrow) + mod(j, _ncol) * _nrow);
-        int ncol = static_cast<int>(_ncol);
+    }
-
-        while (i < 0)     { i += nrow; }
-        while (i >= nrow) { i -= nrow; }
-        while (j < 0)     { j += ncol; }
-        while (j >= ncol) { j -= ncol; }
-
-        return static_cast<size_t>(i + j * _nrow);
-    };
 };
 
 template <typename T>
@@ -68,6 +151,11 @@ public:
 
     const size_t size() const { return _nelem; };
 
+    const size_t stride() const { return _stride; };
+
+          T* data()       { return _matrix.data() + _index; };
+    const T* data() const { return _matrix.data() + _index; };
+
           T& operator()(int i)       { return _matrix(_index + i * _stride); };
     const T& operator()(int i) const { return _matrix(_index + i * _stride); };
 
@@ -90,15 +178,17 @@ std::ostream& operator<<(std::ostream& out, const MatrixView<T>& view) {
 template <typename T>
 class Row : public MatrixView<T> {
 public:
-    Row(Matrix<T>& matrix, size_t index) :
+    Row(Matrix<T>& matrix, int index) :
-        MatrixView<T>(matrix, index * matrix.ncol(), 1, matrix.nrow()) { };
+        MatrixView<T>(matrix, mod(index, matrix.nrow()),
+            matrix.nrow(), matrix.ncol()) { };
 };
 
 template <typename T>
 class Col : public MatrixView<T> {
 public:
-    Col(Matrix<T>& matrix, size_t index) :
+    Col(Matrix<T>& matrix, int index) :
-        MatrixView<T>(matrix, index, matrix.nrow(), matrix.ncol()) { };
+        MatrixView<T>(matrix, mod(index, matrix.ncol()) * matrix.nrow(),
+            1, matrix.nrow()) { };
 };
 
 template <typename T>

Exercise_01/Solver.h

@@ -0,0 +1,114 @@
+#pragma once
+
+#include <functional>
+
+#include "Matrix.h"
+
+/** Stencil
+ */
+struct Stencil {
+    const double C; /*< Center */
+    const double B; /*< Bottom */
+    const double L; /*< Left */
+    const double U; /*< Up */
+    const double R; /*< Right */
+};
+
+class Solver {
+public:
+    Solver(const size_t nx, const size_t ny, const double h, const double k,
+        std::function<double(double, double)> fun,
+        std::function<double(double)> gN,
+        std::function<double(double)> gE,
+        std::function<double(double)> gS,
+        std::function<double(double)> gW
+    ) :
+        _iter{0},
+        _nx{nx},
+        _ny{ny},
+        _stencil{ 4.0 / (h * h) + 4.0 * k * k,  /* Center */
+                 -1.0 / (h * h),                /* Bottom */
+                 -1.0 / (h * h),                /* Left */
+                 -1.0 / (h * h),                /* Up */
+                 -1.0 / (h * h)},               /* Right */
+        _rhs(nx, ny),
+        _sol(nx, ny, 0.0),
+        _tmp(nx, ny)
+    {
+        // Initialize Right Hand Size _rhs(x, y) = f(X(x), Y(y))
+        // Note that the correspondence usual matrix indexing sceeme as
+        // row/colums indices lead to an missplaces representation if the matrix
+        // is printed in the usual format (x <-> rows, y <-> negative columns).
+        // If this in entierly ignored and only considured in the case of
+        // printing the matrix, everything is fine.
+        for (size_t x = 0; x < nx; ++x) {
+            for (size_t y = 0; y < ny; ++y) {
+                _rhs(x, y) = fun(X(x), Y(y));
+            }
+        }
+
+        // Set Dirichlet boundary conditions
+        // North
+        for (size_t x = 0; x < nx; ++x) { _sol(x, -1) = _tmp(x, -1) = gN(X(x)); }
+        // East
+        for (size_t y = 0; y < ny; ++y) { _sol(-1, y) = _tmp(-1, y) = gE(Y(y)); }
+        // South
+        for (size_t x = 0; x < nx; ++x) { _sol(x, 0) = _tmp(x, 0) = gS(X(x)); }
+        // West
+        for (size_t y = 0; y < ny; ++y) { _sol(0, y) = _tmp(0, y) = gW(Y(y)); }
+    };
+
+    /* Maps x/y index to x/y position in the [0, 1] x [0, 1] domain */
+    double X(size_t x) const {
+        return static_cast<double>(x) / static_cast<double>(_nx - 1);
+    }
+    double Y(size_t y) const {
+        return static_cast<double>(y) / static_cast<double>(_ny - 1);
+    }
+
+    enum class Dir {
+        North, East, South, West,
+        N = North, E = East, S = South, W = West
+    };
+
+    MatrixView<double> boundary(enum Dir dir) {
+        switch (dir) {
+            case Dir::North:
+                return Col<double>(_sol, -1);
+            case Dir::East:
+                return Row<double>(_sol, -1);
+            case Dir::South:
+                return Col<double>(_sol, 0);
+            default:    // Dir::West
+                return Row<double>(_sol, 0);
+        }
+    };
+
+    /**
+     * Performs a single Jacobian Iteration
+     */
+    void iterate() {
+        double s = 1.0 / _stencil.C;
+
+        for (size_t y = 1; y < _ny - 1; ++y) {
+            for (size_t x = 1; x < _nx - 1; ++x) {
+                _tmp(x, y) = s * (_rhs(x, y) - (_stencil.B * _sol(x, y - 1) +
+                                                _stencil.L * _sol(x - 1, y) +
+                                                _stencil.U * _sol(x, y + 1) +
+                                                _stencil.R * _sol(x + 1, y)));
+            }
+        }
+
+        std::swap(_tmp, _sol);
+        _iter++;
+    };
+
+private:
+    size_t _iter;           /*< Iteration count */
+    const size_t _nx;       /*< Number of X-axis grid points */
+    const size_t _ny;       /*< Number of Y-axis grid points */
+    const Stencil _stencil; /*< Simple, + shaped stencil */
+    Matrix<double> _rhs;    /*< Grid evaluated RHS of the PDE = f(x, y) */
+    Matrix<double> _sol;    /*< Solution after _iter iterations */
+    Matrix<double> _tmp;    /*< Temp. datablock, used in iterate() */
+};

Exercise_01/examples/MPI_Send_Recv.cpp

@@ -0,0 +1,78 @@
+/**
+ * Compile
+ *  mpic++ MPI_Send_Recv.cpp -Wall -pedantic -Wpedantic -o MPI_Send_Recv
+ * Usage
+ *  mpirun -n 4 ./MPI_Send_Recv
+ * 
+ * Note: An easy way to finde the location of the `mpi.h` file is
+ *  mpic++ --showme:compile
+ * which might be usefull to configure the intellicense.
+ */
+
+#include <iostream>
+#include <sstream>
+#include <mpi.h>
+
+/** min of two ints (fine for an example) */
+int min(int a, int b) { return a < b ? a : b; }
+
+int main(int argn, char* argv[]) {
+
+    // Initialize MPI (always required)
+    MPI_Init(nullptr, nullptr);
+
+    // Allocate MPI Settings
+    int mpi_size;       /*< Number of processes */
+    int mpi_rank;       /*< This process rank (a.k.a. the MPI process ID) */
+    // Set/Get MPI Settings
+    MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
+    MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
+
+    // Write MPI info to stdout
+    std::cout << "Hello World, i am rank " << mpi_rank << " of " << mpi_size
+        << " processes." << std::endl;
+
+    // setting this to small truncates the message (Try it, works fine.)
+    constexpr int max_size = 128;   /*< Maximum message length */
+    // Distinguish between rank 0 and the rest
+    if (mpi_rank == 0) {
+        // Send a different message to all other ranks
+        for (int rank = 1; rank < mpi_size; ++rank) {
+            // Build message
+            std::ostringstream message;
+            message << "Hello rank nr. " << rank;
+            std::string sendbuf = message.str();
+
+            // Send message
+            MPI_Send(sendbuf.c_str(), min(max_size, sendbuf.size()),
+                     MPI_CHAR, rank, sendbuf.size(), MPI_COMM_WORLD);
+            // Note: the tag beeing set to the string length (hacky, but :-})
+        }
+
+        // Report your done
+        std::cout << "Done: rank 0 send all messages." << std::endl;
+    } else {
+        char recvbuf[max_size]; /*< MPI_Recv buffer, e.g. memory to write to */
+        MPI_Status status;      /*< MPI Status object, e.g. communication info */
+
+        // Each rank listens for a single message from rank 0
+        MPI_Recv(recvbuf, max_size,
+                 MPI_CHAR, 0, MPI_ANY_TAG, MPI_COMM_WORLD, &status);
+
+        if (status.MPI_ERROR) {
+            std::cerr << "Error: Recv reported an Error?!" << std::endl;
+        } else {
+            // Ensure C-style string is terminated and of propper size as
+            // the tag encodes the string length.
+            // Should be null terminated, but better save than sorry.
+            recvbuf[min(max_size - 1, status.MPI_TAG)] = '\0';
+            // Print receved message to console
+            std::cout << "Recv: " << recvbuf << std::endl;
+        }
+    }
+
+    // Shutdown MPI (always required)
+    MPI_Finalize();
+
+    return 0;
+}

Exercise_01/main.cpp

@@ -2,29 +2,98 @@
  * g++ main.cpp -std=c++17 -Wall -Wpedantic -pedantic -o main; ./main
  */
 
+#include <stddef.h>
 #include <iostream>
+#include <functional>
+#define _USE_MATH_DEFINES   /* enables math constants from cmath */
+#include <cmath>
+#ifdef USE_MPI
+#include <mpi.h>
+#endif
 
 #include "Matrix.h"
+#include "Solver.h"
 
 int main(int argn, char* argv[]) {
 
-    Matrix<int> mat(3, 4);
+    /******************************* MPI Setup ********************************/
-    for (size_t i = 0; i < mat.size(); ++i) {
+#ifdef USE_MPI
-        mat(i) = i;
+    // Initialize MPI
+    MPI_Init(nullptr, nullptr);
+
+    // Get MPI config
+    int mpi_size; /*< MPI pool size (a.k.a. total number of processes) */
+    int mpi_rank; /*< MPI rank (a.k.a. process ID in the context of MPI) */
+    MPI_Comm_size(MPI_COMM_WORLD, &mpi_size);
+    MPI_Comm_rank(MPI_COMM_WORLD, &mpi_rank);
+#endif
+
+    /**************************** Parse Arguments *****************************/
+    if (argn < 3) {
+        std::cerr << "usage: " << argv[0] << " <resolution> <iterations>" << std::endl;
+        return -1;
+    } else if (argn > 3) {
+        std::cerr << "warning: " << "ignoring all but the first two params" << std::endl;
+    }
+    // TODO: make this proper!!!
+    size_t resolution = atol(argv[1]);
+    size_t iterations = atol(argv[2]);
+    if (resolution < 1 || resolution > 65536
+    ||  iterations < 1 || iterations > 65536) {
+        std::cerr << "error: parsing arguments failed" << std::endl;
     }
 
-    // std::cout << "Matrix:\n" << mat(2, 1) << std::endl;
-    std::cout << "Matrix:\n" << mat.nrow() << std::endl;
 
-    std::cout << "Matrix:\n" << mat << std::endl;
+    /************************* Initialize PDE Solver **************************/
+    size_t nx = resolution;
+    size_t ny = resolution;
+    const double k = M_PI;
+    const double h = 1.0 / static_cast<double>(resolution - 1);
 
-    Col<int> colView(mat, 2);
+    // Declare right hand side function f(x, y) = k^2 sin(2 pi x) sinh(2 pi y)
-    Row<int> rowView(mat, 1);
+    std::function<double(double, double)> fun = [k](double x, double y) {
-    Diag<int> diagView(mat);
+        return k * k * sin(M_2_PI * x) * sinh(M_2_PI * y);
+    };
+    // Boundary conditions
+    /** North boundary condition g(x) = k^2 sin(2 pi x) sinh(2 pi) */
+    std::function<double(double)> gN = [k](double x) {
+        return k * k * sin(M_2_PI * x) * sinh(M_2_PI);
+    };
+    /** East, South and West boundary conditions are simply = 0 */
+    std::function<double(double)> g0 = [k](double) { return 0.0; };
 
-    std::cout << "col:\n" << colView << std::endl;
+    // Instanciate solver
-    std::cout << "row:\n" << rowView << std::endl;
+    Solver solver(nx, ny, h, k, fun, gN, g0, g0, g0);
-    std::cout << "diag:\n" << diagView << std::endl;
+
+    // Set Diriclet boundary conditions (East, South and West to 0)
+    for (auto dir : { Solver::Dir::E, Solver::Dir::S, Solver::Dir::W }) {
+        MatrixView<double> boundary = solver.boundary(dir);
+        for (size_t i = 0; i < boundary.size(); ++i) {
+            boundary(i) = 0.0;
+        }
+    }
+    // The North boundary condition is g(x) = sin(2 pi x) sinh(2 pi)
+    {
+        MatrixView<double> boundary = solver.boundary(Solver::Dir::North);
+        for (size_t i = 0; i < boundary.size(); ++i) {
+            double x = static_cast<double>(i) / static_cast<double>(nx - 1);
+            boundary(i) = sin(M_2_PI * x) * sinh(M_2_PI);
+        }
+    }
+
+    /******************************* Solve PDE ********************************/
+    // Run solver iterations
+    for (size_t iter = 0; iter < iterations; ++iter) {
+        solver.iterate();
+    }
+
+    /****************************** Tests/Report ******************************/
+
+
+    // MPI shutdown/cleanup
+#ifdef USE_MPI
+    MPI_Finalize();
+#endif
 
     return 0;
 }