add lab 11

lab11/game_of_life/animate.py

@@ -0,0 +1,30 @@
+import numpy as np
+import glob
+import os
+import matplotlib.pyplot as plt
+import matplotlib.animation as animation
+
+def load_all_grids(folder):
+    i = 0
+    grids = []
+    while os.path.exists(f"{folder}/step_{i}.txt"):
+        grid = np.loadtxt(f"{folder}/step_{i}.txt", dtype=int)
+        grids.append(grid)
+        i += 1
+    print(f"Loaded {len(grids)} grids from {folder}")
+    return grids
+
+filename = "output"
+grids = load_all_grids(filename)
+
+fig, ax = plt.subplots()
+ax.set_aspect('equal')
+img = ax.imshow(grids[0], cmap='binary')
+
+def update(frame):
+    img.set_data(grids[frame])
+    ax.set_title(f"Generation {frame}")
+    return img
+
+ani = animation.FuncAnimation(fig, update, frames=len(grids), interval=150)
+plt.show()

lab11/game_of_life/common.h

@@ -0,0 +1,8 @@
+#ifndef COMMON_H
+#define COMMON_H
+
+#include <array>
+
+using MPIGridSize = std::array<int, 2>;
+
+#endif // COMMON_H

lab11/game_of_life/game_of_life.cpp

@@ -0,0 +1,56 @@
+
+#include "game_of_life.h"
+
+GameOfLife::GameOfLife(const Matrix& grid, MPIGridSize mpiProcs)
+    : grid_(grid), mpiProcs_(mpiProcs) {}
+
+void GameOfLife::step() {
+  Matrix next = Matrix::zeros(grid_.rows(), grid_.cols());
+  const int rows = grid_.rows();
+  const int cols = grid_.cols();
+  for (int i = 0; i < rows; ++i) {
+    for (int j = 0; j < cols; ++j) {
+      const int numLiveNeighbors = countLiveNeighbors(i, j);
+      next(i, j) = updateCell(grid_(i, j), numLiveNeighbors);
+    }
+  }
+
+  grid_ = next;
+}
+
+int GameOfLife::countLiveNeighbors(int row, int col) const {
+  int count = 0;
+  const int rows = grid_.rows();
+  const int cols = grid_.cols();
+  for (int i = -1; i <= 1; ++i) {
+    for (int j = -1; j <= 1; ++j) {
+      if (i == 0 && j == 0)
+        continue;                             // Skip the cell itself
+      int nextRow = (row + i + rows) % rows;  // Wrap around
+      int nextCol = (col + j + cols) % cols;  // Wrap around
+      if (grid_(nextRow, nextCol) == 1) {
+        count++;
+      }
+    }
+  }
+
+  return count;
+}
+
+int GameOfLife::updateCell(int currentState, int numLiveNeighbors) const {
+  if (numLiveNeighbors == 3) {
+    return 1;
+  } else if (numLiveNeighbors == 2) {
+    return currentState;
+  } else {
+    return 0;
+  }
+}
+
+Matrix GameOfLife::getGrid() const {
+  return grid_;
+}
+
+MPIGridSize GameOfLife::mpiProcs() const {
+  return mpiProcs_;
+}

lab11/game_of_life/game_of_life.h

@@ -0,0 +1,38 @@
+#ifndef GAME_OF_LIFE_H
+#define GAME_OF_LIFE_H
+
+#include <iostream>
+#include <vector>
+#include <fstream>
+#include <memory>
+#include <array>
+#include <mpi.h>
+
+#include "common.h"
+#include "matrix.h"
+
+class GameOfLife {
+ public:
+  GameOfLife(const Matrix& grid, MPIGridSize mpiProcs);
+
+  void step();
+
+  Matrix getGrid() const;
+
+  MPIGridSize mpiProcs() const;
+
+ private:
+  int countLiveNeighbors(int row, int col) const;
+
+  int updateCell(int currentState, int numLiveNeighbors) const;
+
+  Matrix grid_;
+  MPIGridSize mpiProcs_;
+  int myRank_ = 0;
+
+  std::array<std::array<int, 3>, 3> neighborRanks_;
+
+  friend class GameOfLifeTest;
+};
+
+#endif  // GAME_OF_LIFE_H

lab11/game_of_life/main.cpp

@@ -0,0 +1,80 @@
+#include <iostream>
+#include <fstream>
+#include <memory>
+#include <mpi.h>
+#include <string>
+#include <filesystem>
+
+#include "game_of_life.h"
+#include "patterns.h"
+#include "matrix_io.h"
+#include "utils.h"
+#include "common.h"
+
+/**
+ * Main function to run the simulation of the game of life
+ */
+void gameOfLife(MPIGridSize mpiProcs) {
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+  // Create a grid that results in some interesting patterns
+  const Matrix grid = Pattern(20, 40, mpiProcs)
+                          .glider(10, 17)
+                          .beeHive(7, 10)
+                          .octagon(6, 27)
+                          .octagon(12, 0)
+                          .getGrid();
+
+  GameOfLife game(grid, mpiProcs);
+
+  if (rank == 0)
+    std::cout << "Initial State:" << std::endl;
+
+  print(game);
+
+  for (int i = 0; i < 50; ++i) {
+    game.step();
+  }
+  if (rank == 0)
+    std::cout << "Final state" << std::endl;
+  print(game);
+
+  storeAnimation("output", grid, 150, mpiProcs);
+}
+
+/**
+ * Main entry point for the MPI program.
+ * Initializes MPI, checks command line arguments, and starts the game of life
+ * simulation.
+ */
+int main(int argc, char* argv[]) {
+  MPI_Init(&argc, &argv);
+  if (argc != 3) {
+    std::cout << "Specify number of processes in x and y as arguments\n";
+    std::cout << "jacobi <np0> <np1>\n";
+    return 1;
+  }
+
+  const int np0 = std::stoi(argv[1]);
+  const int np1 = std::stoi(argv[2]);
+  int numProc;
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  MPI_Comm_size(MPI_COMM_WORLD, &numProc);
+  if (np0 * np1 != numProc) {
+    std::cout << "Error: nproc != np0 x np1 (" << numProc << "!= " << np0 << "x"
+              << np1 << ")\n";
+    return 2;
+  }
+
+  try {
+    gameOfLife({np0, np1});
+  } catch (const std::exception& e) {
+    std::cerr << "Error: " << e.what() << std::endl;
+    MPI_Abort(MPI_COMM_WORLD, 1);
+    return 1;
+  }
+
+  MPI_Finalize();
+}

lab11/game_of_life/makefile

@@ -0,0 +1,33 @@
+PROGRAM_NAME = game_of_life
+SOURCE_FILES = game_of_life.cpp \
+							 utils.cpp \
+							 patterns.cpp \
+							 matrix_io.cpp 
+
+COMPILER_FLAGS = -std=c++17 -Wall -Wextra -Wall -Werror \
+								 -Wno-sign-compare \
+								 -Wno-unused-result \
+								 -Wno-cast-function-type \
+								 -Wno-unused-variable \
+								 -Wno-unused-parameter 
+
+default: debug
+
+release: $(SOURCE_FILES)
+	mpicxx $(SOURCE_FILES) main.cpp -O3 -o $(PROGRAM_NAME) ${COMPILER_FLAGS}
+
+debug: $(SOURCE_FILES)
+	mpicxx $(SOURCE_FILES) main.cpp -g -o $(PROGRAM_NAME) ${COMPILER_FLAGS}
+
+test:
+	mpicxx $(SOURCE_FILES) test.cpp -g -o game_of_life_test ${COMPILER_FLAGS}
+	mpirun -np 16 --oversubscribe ./game_of_life_test
+
+runserial:
+	mpirun -np 1 ./$(PROGRAM_NAME) 1 1
+
+run: release
+	mpirun -np 16 --oversubscribe ./$(PROGRAM_NAME) 4 4
+
+ani: animate.py
+	python3 animate.py

lab11/game_of_life/matrix.h

@@ -0,0 +1,417 @@
+/* matrix.h, an extrmely simple matrix class.
+ * Version 2.1
+ * Copyright (C) 2022-2025 Tobias Kreilos, Offenburg University of Applied
+ * Sciences
+ *
+ * Licensed under the Apache License, Version 2.0(the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+
+#ifndef MATRIX_H
+#define MATRIX_H
+
+#include <vector>
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+#include <cassert>
+
+/**
+ * Matrix class
+ *
+ * This class implements a matrix of size m x n. The matrix is stored in
+ * a one-dimensional array of size m*n. The matrix is stored in column-major
+ * order, i.e. the first column is stored first, then the second column, etc.
+ *
+ * Static functions are provided to create a matrix of zeros or an uninitialized
+ * matrix. The uninitialized matrix is not initialized, i.e. the entries are
+ * not set to zero. This is useful for performance reasons, e.g. to control
+ * the placement of matrix entries in locality-domain memory.
+ *
+ *
+ */
+class Matrix {
+ public:
+  /**
+   * Create a matrix of size m x n and initialize all entries to zero.
+   * @param rows number of rows
+   * @param cols number of columns
+   */
+  static Matrix zeros(int numRows, int numCols);
+
+  /**
+   * Create a square matrix of size n x n and initialize all entries to zero.
+   * @param n number of rows and columns
+   */
+  static Matrix zeros(int n);
+
+  /**
+   * Create a matrix of size m x n and initialize all entries to zero.
+   * @param dim number of rows and columns
+   */
+  static Matrix zeros(std::pair<int, int> dim);
+
+  /**
+   * Create a matrix of size m x n and do not initialize the entries.
+   * @param rows number of rows
+   * @param cols number of columns
+   */
+  static Matrix uninit(int m, int n);
+
+  /**
+   * Create a square matrix of size n x n and do not initialize the entries.
+   * @param n number of rows and columns
+   */
+  static Matrix uninit(int n);
+
+  /**
+   * Create a matrix of size m x n and do not initialize the entries.
+   * @param dim number of rows and columns
+   */
+  static Matrix uninit(std::pair<int, int> dim);
+
+  Matrix(const Matrix& other);
+  Matrix(Matrix&& other);
+  ~Matrix();
+  Matrix& operator=(const Matrix& other);
+  Matrix& operator=(Matrix&& other);
+
+  // Access element a_ij of the matrix
+  double& operator()(int i, int j);
+  const double& operator()(int i, int j) const;
+
+  // Obtain a pointer to the underlying data
+  double* data();
+  const double* data() const;
+
+  // Getter functions for the dimensions
+  std::pair<int, int> dim() const;
+  int rows() const;
+  int cols() const;
+  int numEntries() const;
+
+  // Comparison operators
+  bool operator==(const Matrix& b) const;
+  bool operator!=(const Matrix& b) const;
+
+  // addition
+  Matrix& operator+=(const Matrix& b);
+
+  // subtraction
+  Matrix& operator-=(const Matrix& b);
+
+  // scalar multiplication
+  Matrix& operator*=(double x);
+
+  // scalar division
+  Matrix& operator/=(double x);
+
+ private:
+  // Constructor is private to prevent creating an uninitialized matrix
+  // accidentally. Use Matrix::zeros() or Matrix::uninit() instead
+  Matrix(int m, int n);
+
+  int numRows_;   // number of rows
+  int numCols_;   // number of columns
+  double* data_;  // the matrix' entries
+};
+
+/**
+ * Vector class
+ * 
+ * This class implements a vector of size n. The vector is stored in a
+ * Matrix of size n x 1.
+ * 
+ * Constructors are provided to create a vector of zeros or an uninitialized
+ * vector. The uninitialized vector is not initialized, i.e. the entries are
+ * not set to zero. This is useful for performance reasons, e.g. to control
+ * the placement of vector entries in locality-domain memory.
+ */
+class Vector {
+ public:
+  static Vector zeros(int n) {
+    Vector vec;
+    vec.data_ = Matrix::zeros(n, 1);
+    return vec;
+  }
+
+  static Vector uninit(int n) {
+    Vector vec;
+    vec.data_ = Matrix::uninit(n, 1);
+    return vec;
+  }
+
+  bool operator==(const Vector& b) const { return data_ == b.data_; }
+  bool operator!=(const Vector& b) const { return !operator==(b); }
+  double& operator()(int i) { return data_(i, 0); }
+  const double& operator()(int i) const { return data_(i, 0); }
+  double* data() { return data_.data(); }
+  const double* data() const { return data_.data(); }
+
+  Vector operator+=(const Vector& b) {
+    data_ += b.data_;
+    return *this;
+  }
+
+  Vector operator-=(const Vector& b) {
+    data_ -= b.data_;
+    return *this;
+  }
+  Vector operator*=(double x) {
+    data_ *= x;
+    return *this;
+  }
+  Vector operator/=(double x) {
+    data_ /= x;
+    return *this;
+  }
+
+  int size() const { return data_.rows(); }
+
+ private:
+  Matrix data_ = Matrix::zeros(0, 0);
+};
+
+/********** Implementation below ********************/
+
+inline Matrix Matrix::zeros(int m, int n) {
+  Matrix mat(m, n);
+  for (int i = 0; i < m; ++i) {
+    for (int j = 0; j < n; ++j) {
+      mat(i, j) = 0;
+    }
+  }
+  return mat;
+}
+
+inline Matrix Matrix::zeros(int n) {
+  return zeros(n, n);
+}
+
+inline Matrix Matrix::zeros(std::pair<int, int> dim) {
+  return zeros(dim.first, dim.second);
+}
+
+inline Matrix Matrix::uninit(int m, int n) {
+  return Matrix(m, n);
+}
+
+inline Matrix Matrix::uninit(int n) {
+  return uninit(n, n);
+}
+
+inline Matrix Matrix::uninit(std::pair<int, int> dim) {
+  return uninit(dim.first, dim.second);
+}
+
+inline Matrix::Matrix(const Matrix& other) {
+  numRows_ = other.numRows_;
+  numCols_ = other.numCols_;
+  if (numRows_ == 0 || numCols_ == 0) {
+    data_ = nullptr;
+    return;
+  }
+  data_ = new double[numRows_ * numCols_];
+  for (int i = 0; i < numRows_; ++i) {
+    for (int j = 0; j < numCols_; ++j) {
+      operator()(i, j) = other(i, j);
+    }
+  }
+}
+
+inline Matrix::Matrix(Matrix&& other) {
+  numRows_ = other.numRows_;
+  numCols_ = other.numCols_;
+  data_ = other.data_;
+  other.data_ = nullptr;
+  other.numRows_ = 0;
+  other.numCols_ = 0;
+}
+
+inline Matrix::~Matrix() {
+  if (data_ != nullptr) {
+    delete[] data_;
+    data_ = nullptr;
+  }
+}
+
+inline Matrix& Matrix::operator=(const Matrix& other) {
+  if (this != &other) {
+    if (data_ != nullptr) {
+      delete[] data_;
+    }
+    numRows_ = other.numRows_;
+    numCols_ = other.numCols_;
+    data_ = new double[numRows_ * numCols_];
+    for (int i = 0; i < numRows_; ++i) {
+      for (int j = 0; j < numCols_; ++j) {
+        operator()(i, j) = other(i, j);
+      }
+    }
+  }
+  return *this;
+}
+
+inline Matrix& Matrix::operator=(Matrix&& other) {
+  if (this != &other) {
+    if (data_ != nullptr) {
+      delete[] data_;
+    }
+    numRows_ = other.numRows_;
+    numCols_ = other.numCols_;
+    data_ = other.data_;
+    other.data_ = nullptr;
+    other.numRows_ = 0;
+    other.numCols_ = 0;
+  }
+  return *this;
+}
+
+inline double& Matrix::operator()(int i, int j) {
+  assert(i >= 0 && i < numRows_);
+  assert(j >= 0 && j < numCols_);
+  return data_[i * numCols_ + j];
+}
+
+inline const double& Matrix::operator()(int i, int j) const {
+  assert(i >= 0 && i < numRows_);
+  assert(j >= 0 && j < numCols_);
+  return data_[i * numCols_ + j];
+}
+
+inline double* Matrix::data() {
+  return data_;
+}
+
+inline const double* Matrix::data() const {
+  return data_;
+}
+
+inline std::pair<int, int> Matrix::dim() const {
+  return std::pair<int, int>(numRows_, numCols_);
+}
+
+inline int Matrix::rows() const {
+  return numRows_;
+}
+
+inline int Matrix::cols() const {
+  return numCols_;
+}
+
+inline int Matrix::numEntries() const {
+  return numRows_ * numCols_;
+}
+
+inline bool Matrix::operator==(const Matrix& b) const {
+  const double eps = 1e-12;
+  if (numRows_ != b.numRows_ || numCols_ != b.numCols_) {
+    return false;
+  }
+  for (int i = 0; i < numRows_; ++i) {
+    for (int j = 0; j < numCols_; ++j) {
+      if (fabs(operator()(i, j) - b(i, j)) > eps) {
+        return false;
+      }
+    }
+  }
+  return true;
+}
+
+inline bool Matrix::operator!=(const Matrix& b) const {
+  return !operator==(b);
+}
+
+inline Matrix& Matrix::operator+=(const Matrix& b) {
+  for (int i = 0; i < numRows_; ++i) {
+    for (int j = 0; j < numCols_; ++j) {
+      operator()(i, j) += b(i, j);
+    }
+  }
+  return *this;
+}
+
+inline Matrix& Matrix::operator-=(const Matrix& b) {
+  for (int i = 0; i < numRows_; ++i) {
+    for (int j = 0; j < numCols_; ++j) {
+      operator()(i, j) -= b(i, j);
+    }
+  }
+  return *this;
+}
+
+inline Matrix& Matrix::operator*=(double x) {
+  for (int i = 0; i < numRows_; ++i) {
+    for (int j = 0; j < numCols_; ++j) {
+      operator()(i, j) *= x;
+    }
+  }
+  return *this;
+}
+
+inline Matrix& Matrix::operator/=(double x) {
+  for (int i = 0; i < numRows_; ++i) {
+    for (int j = 0; j < numCols_; ++j) {
+      operator()(i, j) /= x;
+    }
+  }
+  return *this;
+}
+
+inline Matrix::Matrix(int m, int n) : numRows_(m), numCols_(n) {
+  data_ = new double[numRows_ * numCols_];
+  if (data_ == nullptr) {
+    std::cout << "Error: not enough memory for matrix\n";
+    throw std::bad_alloc();
+  }
+}
+
+inline std::ostream& operator<<(std::ostream& os, const Matrix& a) {
+  const int width = 10;
+  const int precision = 4;
+
+  const auto originalPrecision = os.precision();
+  os << std::setprecision(precision);
+
+  for (int i = 0; i < a.rows(); ++i) {
+    for (int j = 0; j < a.cols(); ++j) {
+      os << std::setw(width) << a(i, j) << " ";
+    }
+    if (i != a.rows() - 1)
+      os << "\n";
+  }
+
+  os << std::setprecision(originalPrecision);
+  return os;
+}
+
+inline bool equalWithinRange(const Matrix& a,
+                             const Matrix& b,
+                             double eps = 1e-12) {
+  if (a.rows() != b.rows() || a.cols() != b.cols())
+    return false;
+
+  int m = a.rows();
+  int n = a.cols();
+  for (int i = 0; i < m; ++i) {
+    for (int j = 0; j < n; ++j) {
+      if (fabs(a(i, j) - b(i, j)) > eps) {
+        return false;
+      }
+    }
+  }
+
+  return true;
+}
+
+#endif  // MATRIX_H

lab11/game_of_life/matrix_io.cpp

@@ -0,0 +1,141 @@
+#include <fstream>
+#include <mpi.h>
+#include <exception>
+#include "matrix.h"
+#include "matrix_io.h"
+
+MatrixIO::MatrixIO(MPIGridSize mpiProcs) : mpiProcs_(mpiProcs) {
+  std::array<int, 2> periods = {1, 1};
+  MPI_Cart_create(MPI_COMM_WORLD, 2, mpiProcs.data(), periods.data(), true,
+                  &comm_);
+  MPI_Comm_rank(comm_, &rank_);
+}
+
+void MatrixIO::saveDistributed(const Matrix& distributedMatrix,
+                               const std::string& filename) {
+  Matrix mat = gatherMatrixOnRoot(distributedMatrix);
+
+  if (rank() == 0) {
+    saveSerial(mat, filename);
+  }
+}
+
+void MatrixIO::saveSerial(const Matrix& m, const std::string& filename) {
+  std::ofstream fout(filename);
+  const int numRows = m.rows();
+  const int numCols = m.cols();
+  fout << "# " << numRows << "\t" << numCols << "\n";
+  for (int i = 0; i < numRows; ++i) {
+    for (int j = 0; j < numCols; ++j) {
+      fout << m(i, j) << "\t";
+    }
+    fout << "\n";
+  }
+}
+
+Matrix MatrixIO::load(const std::string& filename) {
+  Matrix matrixOnRoot = Matrix::zeros(0, 0);  // Initialize empty matrix
+  if (rank() == 0) {
+    matrixOnRoot = loadSerial(filename);
+  }
+
+  Matrix distributedMatrix = scatterMatrixFromRoot(matrixOnRoot);
+  return distributedMatrix;
+}
+
+Matrix MatrixIO::loadSerial(const std::string& filename) {
+  std::ifstream fin(filename);
+
+  // Check first character (has to be #)
+  std::string s;
+  fin >> s;
+  if (s != "#") {
+    throw std::runtime_error("Error, not reading expecte character #\n");
+  }
+
+  // Read in number of rows and cols and create matrix
+  int numRows, numCols;
+  fin >> numRows >> numCols;
+  Matrix mat = Matrix::zeros(numRows, numCols);
+
+  // Read in matrix contents
+  for (int i = 0; i < numRows; ++i) {
+    for (int j = 0; j < numCols; ++j) {
+      fin >> mat(i, j);
+    }
+  }
+
+  return mat;
+}
+
+Matrix MatrixIO::gatherMatrixOnRoot(const Matrix& distributedMatrix) {
+  const int numRowsLocal = distributedMatrix.rows();
+  const int numColsLocal = distributedMatrix.cols();
+  Matrix subMatrix(distributedMatrix);
+  Matrix matrixOnRoot = Matrix::zeros(0, 0);  // Initialize empty matrix
+  if (rank() == 0) {
+    const int numRowsTotal = numRowsLocal * mpiProcs_[0];
+    const int numColsTotal = numColsLocal * mpiProcs_[1];
+    matrixOnRoot = Matrix::zeros(numRowsTotal, numColsTotal);
+
+    for (int proc = 0; proc < nProc(); ++proc) {
+      if (proc != 0) {
+        MPI_Recv(&subMatrix(0, 0), subMatrix.numEntries(), MPI_DOUBLE, proc,
+                 proc, comm_, MPI_STATUS_IGNORE);
+      }
+
+      int coords[2];
+      MPI_Cart_coords(comm_, proc, 2, coords);
+      const int i0 = coords[0] * numRowsLocal;
+      const int j0 = coords[1] * numColsLocal;
+      for (int i = 0; i < numRowsLocal; ++i) {
+        for (int j = 0; j < numColsLocal; ++j) {
+          matrixOnRoot(i + i0, j + j0) = subMatrix(i, j);
+        }
+      }
+    }
+  } else {
+    MPI_Send(&subMatrix(0, 0), subMatrix.numEntries(), MPI_DOUBLE, 0, rank(),
+             comm_);
+  }
+
+  return matrixOnRoot;
+}
+
+Matrix MatrixIO::scatterMatrixFromRoot(const Matrix& matrixOnRoot) {
+  int numRowsTotal = matrixOnRoot.rows();
+  int numColsTotal = matrixOnRoot.cols();
+
+  MPI_Bcast(&numRowsTotal, 1, MPI_INT, 0, MPI_COMM_WORLD);
+  MPI_Bcast(&numColsTotal, 1, MPI_INT, 0, MPI_COMM_WORLD);
+
+  int numRowsLocal = numRowsTotal / mpiProcs_[0];
+  int numColsLocal = numColsTotal / mpiProcs_[1];
+
+  Matrix distributedMatrix = Matrix::zeros(numRowsLocal, numColsLocal);
+
+  if (rank() == 0) {
+    for (int proc = nProc() - 1; proc >= 0; --proc) {  // iterate backwards
+      int coords[2];
+      MPI_Cart_coords(comm_, proc, 2, coords);
+      const int i0 = coords[0] * numRowsLocal;
+      const int j0 = coords[1] * numColsLocal;
+      for (int i = 0; i < numRowsLocal; ++i) {
+        for (int j = 0; j < numColsLocal; ++j) {
+          distributedMatrix(i, j) = matrixOnRoot(i + i0, j + j0);
+        }
+      }
+
+      if (proc != 0) {
+        MPI_Send(&distributedMatrix(0, 0), distributedMatrix.numEntries(),
+                 MPI_DOUBLE, proc, proc, comm_);
+      }
+    }
+
+  } else {
+    MPI_Recv(&distributedMatrix(0, 0), distributedMatrix.numEntries(),
+             MPI_DOUBLE, 0, rank(), comm_, MPI_STATUS_IGNORE);
+  }
+
+  return distributedMatrix;
+}

lab11/game_of_life/matrix_io.h

@@ -0,0 +1,77 @@
+#ifndef MATRIX_IO_H
+#define MATRIX_IO_H
+
+#include <fstream>
+#include <mpi.h>
+#include <exception>
+#include <memory>
+#include "matrix.h"
+#include "common.h"
+
+/**
+ * File i/o for distributed matrices.
+ * Matrices is distributed along 2 dimensions (i.e. rows and columns are
+ * distributed).
+ *
+ * File format is compatible with python numpy loadtxt
+ * First line: # numRows numCols
+ * One line of the matrix per line in the file, separated by whitespace
+ */
+class MatrixIO {
+ public:
+  MatrixIO(MPIGridSize mpiProcs);
+
+  /**
+   * Store a matrix on disk, that is not distribued. Should be called by a
+   * single rank only.
+   */
+  void saveSerial(const Matrix& m, const std::string& filename);
+
+  /**
+   * Load a matrix from disk, that is not distribued. Should be called by a
+   * single rank only.
+   */
+  Matrix loadSerial(const std::string& filename);
+
+  /**
+   * Store a distributed matrix on disk. Should be called collectively by
+   * all ranks.
+   */
+  void saveDistributed(const Matrix& m, const std::string& filename);
+
+  /**
+   * Load a matrix from disk and distribute it along the first axis.
+   * Should be called collectively by all ranks.
+   */
+  Matrix load(const std::string& filename);
+
+  /**
+   * Gather a distributed matrix on the root rank.
+   * The matrix is gathered in a single matrix on the root rank.
+   * Other ranks return an empty matrix.
+   * Should be called collectively by all ranks.
+   */
+  Matrix gatherMatrixOnRoot(const Matrix& distributedMatrix);
+
+  /**
+   * Scatter a matrix from the root rank to all other ranks.
+   * The matrix is scattered along the first axis.
+   * Should be called collectively by all ranks.
+   */
+  Matrix scatterMatrixFromRoot(const Matrix& matrixOnRoot);
+
+ private:
+  int rank() const { return rank_; };
+  int nProc() const {
+    int size;
+    MPI_Comm_size(comm_, &size);
+    return size;
+  }
+
+  int rank_ = 0;
+  MPIGridSize mpiProcs_ = {0, 0};  // Number of processes in each dimension
+
+  MPI_Comm comm_ = MPI_COMM_NULL;  // The communicator for the matrix operations
+};
+
+#endif  // MATRIX_IO_H

lab11/game_of_life/patterns.cpp

@@ -0,0 +1,126 @@
+#include "patterns.h"
+
+Pattern::Pattern(int rows, int cols, MPIGridSize mpiProcs)
+    : mpiProcs_(mpiProcs),
+      grid_(Matrix::zeros(rows / np0(), cols / np1())) {
+  if (rows <= 0 || cols <= 0) {
+    throw std::invalid_argument("Rows and columns must be positive");
+  }
+  if (rows % np0() != 0) {
+    throw std::invalid_argument(
+        "Rows must be divisible by the number of processes in the first "
+        "dimension");
+  }
+  if (cols % np1() != 0) {
+    throw std::invalid_argument(
+        "Columns must be divisible by the number of processes in the second "
+        "dimension");
+  }
+
+  std::array<int, 2> periods = {1, 1};
+  MPI_Cart_create(MPI_COMM_WORLD, 2, mpiProcs.data(), periods.data(), true,
+                  &comm_);
+}
+
+int Pattern::np0() const {
+  return mpiProcs_[0];
+}
+
+int Pattern::np1() const {
+  return mpiProcs_[1];
+}
+
+int Pattern::numRowsLocal() const {
+  return grid_.rows();
+}
+
+int Pattern::numRowsTotal() const {
+  return grid_.rows() * np0();
+}
+
+int Pattern::numColsLocal() const {
+  return grid_.cols();
+}
+
+int Pattern::numColsTotal() const {
+  return grid_.cols() * np1();
+}
+
+Pattern Pattern::beeHive(int row, int col) {
+  if (row < 0 || col < 0 || row + 3 >= numRowsTotal() ||
+      col + 3 >= numColsTotal()) {
+    throw std::out_of_range("Bee hive pattern exceeds grid bounds");
+  }
+  setCell(row, col + 1);
+  setCell(row, col + 2);
+  setCell(row + 1, col);
+  setCell(row + 1, col + 3);
+  setCell(row + 2, col);
+  setCell(row + 2, col + 3);
+  setCell(row + 3, col + 1);
+  setCell(row + 3, col + 2);
+  return *this;
+}
+
+Pattern Pattern::glider(int row, int col) {
+  if (row < 0 || col < 0 || row + 3 >= numRowsTotal() ||
+      col + 3 >= numColsTotal()) {
+    throw std::out_of_range("Glider pattern exceeds grid bounds");
+  }
+  setCell(row, col + 1);
+  setCell(row + 1, col + 2);
+  setCell(row + 2, col);
+  setCell(row + 2, col + 1);
+  setCell(row + 2, col + 2);
+  return *this;
+}
+
+Pattern Pattern::octagon(int row, int col) {
+  if (row < 0 || col < 0 || row + 7 >= numRowsTotal() ||
+      col + 7 >= numColsTotal()) {
+    throw std::out_of_range("Octagon pattern exceeds grid bounds");
+  }
+  setCell(row + 0, col + 3);
+  setCell(row + 0, col + 4);
+  setCell(row + 1, col + 2);
+  setCell(row + 1, col + 5);
+  setCell(row + 2, col + 1);
+  setCell(row + 2, col + 6);
+  setCell(row + 3, col + 0);
+  setCell(row + 3, col + 7);
+  setCell(row + 4, col + 0);
+  setCell(row + 4, col + 7);
+  setCell(row + 5, col + 1);
+  setCell(row + 5, col + 6);
+  setCell(row + 6, col + 2);
+  setCell(row + 6, col + 5);
+  setCell(row + 7, col + 3);
+  setCell(row + 7, col + 4);
+
+  return *this;
+}
+
+Matrix Pattern::getGrid() const {
+  return grid_;
+}
+
+Pattern Pattern::setCell(int globalRow, int globalCol) {
+  if (globalRow < 0 || globalCol < 0 || globalRow >= numRowsTotal() ||
+      globalCol >= numColsTotal()) {
+    throw std::out_of_range("Cell indices are out of bounds");
+  }
+
+  int rank;
+  MPI_Comm_rank(comm_, &rank);
+
+  std::array<int, 2> coords;
+  MPI_Cart_coords(comm_, rank, 2, coords.data());
+  int localRow = globalRow - coords[0] * numRowsLocal();
+  int localCol = globalCol - coords[1] * numColsLocal();
+  if (localRow < 0 || localCol < 0 || localRow >= grid_.rows() ||
+      localCol >= grid_.cols()) {
+    return *this;  // Ignore out-of-bounds indices
+  }
+  grid_(localRow, localCol) = 1;
+  return *this;
+}

lab11/game_of_life/patterns.h

@@ -0,0 +1,51 @@
+#ifndef PATTERNS_H
+#define PATTERNS_H
+
+#include <memory>
+#include <array>
+#include <mpi.h>
+#include "matrix.h"
+#include "common.h"
+
+/**
+ * Class to create some patterns for the Game of Life.
+ *
+ * This class allows you to create predefined patterns like bee hive, glider,
+ * octagon, etc., on a grid that is distributed across multiple processes.
+ *
+ * Pattern creation methods return a reference to the Pattern object itself,
+ * allowing for method chaining.
+ *
+ * Example usage for creating a single glider on a 20x40 grid:
+ * Pattern pattern(20, 40, comm).glider(10, 17).getGrid();
+ */
+class Pattern {
+ public:
+  Pattern(int rows, int cols, MPIGridSize mpiProcs);
+
+  Matrix getGrid() const;
+
+  Pattern setCell(int globalRow, int globalCol);
+
+  Pattern beeHive(int row, int col);
+
+  Pattern glider(int row, int col);
+
+  Pattern octagon(int row, int col);
+
+ private:
+  int numRowsLocal() const;
+  int numRowsTotal() const;
+  int numColsLocal() const;
+  int numColsTotal() const;
+
+  int np0() const;
+  int np1() const;
+
+  MPIGridSize mpiProcs_;
+  Matrix grid_;
+
+  MPI_Comm comm_; 
+};
+
+#endif  // PATTERNS_H

lab11/game_of_life/test.h

@@ -0,0 +1,319 @@
+/*
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ *     http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+/** test.h, an extremly simple test framework.
+ * Version 1.7
+ * Copyright (C) 2022-2024 Tobias Kreilos, Offenburg University of Applied
+ * Sciences
+ */
+
+/**
+ * The framework defines a function check(a,b) that can be called with
+ * parameters of different types. The function asserts
+ * that the two paramters are equal (within a certain, predefined range for
+ * floating point numbers) and prints the result of the comparison on the
+ * command line. Additionally a summary of all tests is printed at the end of
+ * the program.
+ * There is a TEST macro, which you can place outside main to group
+ * tests together. Code in the macro is automatically executed at the beginning
+ * of the program.
+ * The file also defines a class InstanceCount, that can be used to
+ * count how many instances of an object are still alive at the end of a
+ * program. To use it, derive your class from InstanceCount<ClassName> and the
+ * message is automatically printed at the end of the program.
+ *
+ * The functions are thread- and reentrant-safe. Support for OpenMP is included.
+ * Execution with MPI is supported, but no collection of the results occurs. All
+ * tests are executed locally, results are printed for every node separately.
+ *
+ * Caution: the TEST macro uses static storage of objects, so be aware of the
+ * static initialization order fiasco when using multiple source files.
+ *
+ * Example usage:
+ *
+ * #include "test.h"
+ * TEST(MyTest) {
+ *   check(1, 1);
+ * }
+ *
+ * int main() {
+ *   const std::string s = "Hi";
+ *   check(s, "Hi");
+ * }
+ */
+
+#ifndef VERY_SIMPLE_TEST_H
+#define VERY_SIMPLE_TEST_H
+
+#include <atomic>
+#include <cmath>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+
+#ifdef _OPENMP
+#include <omp.h>
+#endif
+
+/** Simple macro to execute the code that follows the macro (without call from
+ * main)
+ *
+ * Define a class, that is directly instantiated
+ * and contains the test code in the constructor.
+ *
+ * Usage:
+ * TEST(MyTest) {
+ *    // test code
+ * }
+ */
+#define TEST(name)              \
+  struct _TestClass##name {     \
+    _TestClass##name();         \
+  } _TestClass##name##Instance; \
+  _TestClass##name::_TestClass##name()
+
+// Use a namespace to hide implementation details
+namespace Test::Detail {
+
+/**
+ * Make it possible to print the underlying value of class enums with ostream
+ *
+ * The expression typename std::enable_if<std::is_enum<T>::value,
+ * std::ostream>::type decays to ostream if the type T is an enum. Otherwise,
+ * the function is not generated.
+ */
+template <typename T>
+std::ostream& operator<<(
+    typename std::enable_if<std::is_enum<T>::value, std::ostream>::type& stream,
+    const T& e) {
+  return stream << static_cast<typename std::underlying_type<T>::type>(e);
+}
+
+/**
+ * Convert anything to a string.
+ */
+template <typename T>
+std::string toString(const T& t) {
+  std::ostringstream ss;
+  ss << std::setprecision(10);
+  ss << t;
+  return "\"" + ss.str() + "\"";
+}
+
+/**
+ * Convert bools to string "true" or "false" instead of 0 and 1
+ */
+template <>
+inline std::string toString<bool>(const bool& b) {
+  return b ? "\"true\"" : "\"false\"";
+}
+
+/**
+ * Comparison function for different types
+ */
+template <typename T>
+bool isEqual(const T& t1, const T& t2) {
+  return t1 == t2;
+}
+
+/**
+ * Double values are equal if they differ no more than 1e-8
+ */
+template <>
+inline bool isEqual<double>(const double& expectedValue,
+                            const double& actualValue) {
+  const double epsilon = 1e-4;
+  const double distance = fabs(actualValue - expectedValue);
+  return (distance < epsilon);
+}
+
+/**
+ * Float values are equal if they differ no more than 1e-4
+ */
+template <>
+inline bool isEqual<float>(const float& expectedValue,
+                           const float& actualValue) {
+  const double epsilon = 1e-4;
+  const double distance = fabs(actualValue - expectedValue);
+  return (distance < epsilon);
+}
+
+/**
+ * This class realizes some basics of the test framework.
+ * Test summary is printed in the destructor.
+ * Apart from that, the class implements counting of total and failed tests,
+ * comparison of floating point numbers within sensible boundaries and prints
+ * the result of each test on the command line.
+ */
+class Test {
+ public:
+  /**
+   * Test class is a Singleton
+   */
+  static Test& instance() {
+    static Test test;
+    return test;
+  }
+
+  /**
+   * the main entry point for tests. Test two values for equality and output the
+   * result.
+   */
+  template <typename T>
+  bool check(const T& expectedValue, const T& actualValue) {
+    bool testResult = isEqual(expectedValue, actualValue);
+    if (testResult == true) {
+      registerPassingTest();
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+      std::cout << "Test successful! Expected value == actual value (="
+                << toString(expectedValue) << ")" << std::endl;
+    } else {
+      registerFailingTest();
+#ifdef _OPENMP
+#pragma omp critical
+#endif
+      std::cout << "Error in test: expected value " << toString(expectedValue)
+                << ", but actual value was " << toString(actualValue)
+                << std::endl;
+    }
+
+    return testResult;
+  }
+
+ private:
+  /**
+   * Print a summary of all tests at the end of program execution.
+   *
+   * Since the Test class is a static Singleton, destruction happens when the
+   * program terminates, so this is a good place to print the summary.
+   */
+  ~Test() {
+    std::cout << "\n--------------------------------------" << std::endl;
+    std::cout << "Test summary:" << std::endl;
+    std::cout << "Executed tests: " << numTests_ << std::endl;
+    std::cout << "Failed tests: " << numFailedTests_ << std::endl;
+  }
+
+  void registerPassingTest() { numTests_++; }
+
+  void registerFailingTest() {
+    numTests_++;
+    numFailedTests_++;
+  }
+
+  /**
+   * For statistics
+   */
+  std::atomic<int> numTests_ = 0;
+
+  /**
+   * For statistics
+   */
+  std::atomic<int> numFailedTests_ = 0;
+};
+
+template <typename T>
+class InstanceCounterHelper {
+ public:
+  ~InstanceCounterHelper() {
+    std::cout << "The remaining number of objects of type " << typeid(T).name()
+              << " at the end of the program is " << count;
+    if (count > 0)
+      std::cout << " (NOT zero!)";
+    std::cout << "\nThe total number of objects created was " << total
+              << std::endl;
+  }
+
+  void increment() {
+    count++;
+    total++;
+  }
+
+  void decrement() { count--; }
+
+ private:
+  std::atomic<int> count = 0;
+  std::atomic<int> total = 0;
+};
+
+}  // namespace Test::Detail
+
+/**
+ * Count the instances of a class T.
+ * Result gets printed automatically at the end of the program.
+ * To use it, inherit T from InstanceCounter<T>, e.g.
+ * class MyClass : InstanceCounter<MyClass>
+ */
+template <typename T>
+class InstanceCounter {
+ public:
+  InstanceCounter() { counter().increment(); }
+
+  InstanceCounter(const InstanceCounter&) { counter().increment(); }
+
+  InstanceCounter(const InstanceCounter&&) { counter().increment(); }
+
+  virtual ~InstanceCounter() { counter().decrement(); }
+
+  Test::Detail::InstanceCounterHelper<T>& counter() {
+    static Test::Detail::InstanceCounterHelper<T> c;
+    return c;
+  }
+};
+
+/**
+ * Check if the expected value is equal to the actual value.
+ * Result is printed on the command line and at the end of the program, a
+ * summary of all tests is printed.
+ */
+template <typename T1, typename T2>
+void check(const T1& actualValue, const T2& expectedValue) {
+  const T1& expectedValueCasted{
+      expectedValue};  // allows conversion in general, but avoids narrowing
+                       // conversion
+  Test::Detail::Test::instance().check(expectedValueCasted, actualValue);
+}
+
+// allow conversion from int to double explicitely
+template <>
+inline void check(const double& actualValue, const int& expectedValue) {
+  Test::Detail::Test::instance().check(static_cast<double>(expectedValue),
+                                       actualValue);
+}
+
+/**
+ * Check if the entered value is true.
+ * Result is printed on the command line and at the end of the program, a
+ * summary of all tests is printed.
+ */
+inline void check(bool a) {
+  Test::Detail::Test::instance().check(true, a);
+}
+
+#endif  // VERY_SIMPLE_TEST_H
+
+/**
+ * V1.0: Creation of framework
+ * V1.1: make check(bool) inline, automatically convert expected value type to
+ * actual value type
+ * V1.2: added possibilty to count constructions and destructions of some type
+ * V1.3: tweaks on check for int and double types
+ * V1.4: Adding thread safety in OpenMP programs
+ * V1.5: reduce accuraccy in comparing double and float to 1e-8
+ * V1.6: Increase precision for printing floating point values
+ * V1.7: Put #ifdef _OPENMP around pragmas to avoid warnings when compiling
+ *       without -fopenmp
+ */

lab11/game_of_life/utils.cpp

@@ -0,0 +1,84 @@
+#include "utils.h"
+
+void clearOrCreateFolder(const std::string& foldername) {
+  namespace fs = std::filesystem;
+  fs::path folder(foldername);
+  if (fs::exists(folder)) {
+    if (!fs::is_directory(folder)) {
+      throw std::runtime_error("Path exists but is not a directory: " +
+                               foldername);
+    }
+    if (!fs::is_empty(folder)) {
+      char choice;
+      std::cout << "Folder '" << foldername
+                << "' already exists and is not empty.\n";
+      std::cout << "Do you want to clear it? (y/N): ";
+      std::cin >> choice;
+      if (choice != 'y' && choice != 'Y') {
+        throw std::runtime_error(
+            "Folder is not empty and user chose not to clear it.");
+      }
+
+      std::cout << "Clearing folder '" << foldername << "'...\n";
+      // Clear the folder
+      for (const auto& entry : fs::directory_iterator(folder)) {
+        if (fs::is_directory(entry)) {
+          fs::remove_all(entry);
+        } else if (fs::is_regular_file(entry)) {
+          fs::remove(entry);
+        }
+      }
+    }
+  } else {
+    // Folder does not exist, create it
+    fs::create_directories(folder);
+  }
+}
+
+void storeAnimation(const std::string& foldername,
+                    const Matrix& initstate,
+                    int numSteps,
+                    MPIGridSize mpiProcs) {
+  GameOfLife game(initstate, mpiProcs);
+
+  MatrixIO io(mpiProcs);
+
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+  // Only rank 0 handles file system operations
+  if (rank == 0) {
+    clearOrCreateFolder(foldername);
+  }
+
+  if (rank == 0) {
+    std::cout << "Storing animation in folder: " << foldername << std::endl;
+  }
+  for (int step = 0; step < numSteps; ++step) {
+    if (rank == 0) {
+      std::string filename =
+          foldername + "/step_" + std::to_string(step) + ".txt";
+      io.saveDistributed(game.getGrid(), filename);
+    }
+    game.step();
+  }
+  if (rank == 0) {
+    std::cout << "Animation finished" << std::endl;
+  }
+}
+
+void print(const GameOfLife& game) {
+  MatrixIO io(game.mpiProcs());
+  Matrix grid = io.gatherMatrixOnRoot(game.getGrid());
+
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  if (rank == 0) {
+    for (int i = 0; i < grid.rows(); ++i) {
+      for (int j = 0; j < grid.cols(); ++j) {
+        std::cout << ((grid(i, j) == 1) ? "X " : ". ");
+      }
+      std::cout << std::endl;
+    }
+  }
+}

lab11/game_of_life/utils.h

@@ -0,0 +1,36 @@
+#ifndef UTILS_H
+#define UTILS_H
+
+#include <iostream>
+#include <memory>
+#include <fstream>
+#include <filesystem>
+
+#include "game_of_life.h"
+#include "matrix_io.h"
+#include "matrix.h"
+#include "common.h"
+
+/**
+ * Clears the contents of a folder or create it if it does not exist.
+ * If the folder exists, the user is asked if he wants to proceed and if yes,
+ * all folder contents are removed.
+ * If the folder does not exist, it is created.
+ */
+void clearOrCreateFolder(const std::string& foldername);
+
+/**
+ * Stores the animation of the Game of Life in a bunch of text files.
+ * Each step of the game is saved in a separate file named "step_<n>.txt",
+ * where <n> is the step number.
+ * The files are stored in the specified folder.
+ */
+void storeAnimation(const std::string& foldername,
+                    const Matrix& initstate,
+                    int numSteps,
+                    MPIGridSize mpiProcs);
+
+/** Print the current grid of the game on the console */
+void print(const GameOfLife& game);
+
+#endif  // UTILS_H