diff --git a/lab10/jacobi/Makefile b/lab10/jacobi/Makefile new file mode 100644 index 0000000..e09bc13 --- /dev/null +++ b/lab10/jacobi/Makefile @@ -0,0 +1,23 @@ +PROGRAM_NAME = jacobi +SOURCE_FILES = jacobi_main.cpp jacobi_serial.cpp jacobi_mpi.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) -O3 -o $(PROGRAM_NAME) ${COMPILER_FLAGS} + +debug: $(SOURCE_FILES) + mpicxx $(SOURCE_FILES) -g -o $(PROGRAM_NAME) ${COMPILER_FLAGS} + +run: release + mpirun -np 4 ./$(PROGRAM_NAME) + +plot: plot.py + python3 plot.py \ No newline at end of file diff --git a/lab10/jacobi/jacobi.h b/lab10/jacobi/jacobi.h new file mode 100644 index 0000000..69dde9e --- /dev/null +++ b/lab10/jacobi/jacobi.h @@ -0,0 +1,30 @@ +#ifndef JACOBI_H +#define JACOBI_H + +#include "matrix.h" + +/** + * Base class for Jacobi algorithms. + * Defines a struct result which is returned by the algorithm. + */ +class Jacobi { + public: + /** + * Result of the computation + */ + struct Result { + Matrix phi = Matrix::zeros(0,0); // The converged solution matrix + double finalDist = 0; // The final distance + int numIter = 0; // The number of iterations made + }; + + /** + * Run the Jacobi algorithm on the initial matrix init. + * The algorithm terminates if either the maximum number of iterations + * exceeds maxNumIter or if the difference between two steps is smaller + * than eps. + */ + virtual Result run(const Matrix& init, double eps, int maxNumIter) = 0; +}; + +#endif // JACOBI_H \ No newline at end of file diff --git a/lab10/jacobi/jacobi_main.cpp b/lab10/jacobi/jacobi_main.cpp new file mode 100644 index 0000000..3eb3f2d --- /dev/null +++ b/lab10/jacobi/jacobi_main.cpp @@ -0,0 +1,182 @@ +#include "jacobi_main.h" + +Matrix getInitialCondition(int n, int rank, int numProc) { + if (n % numProc != 0) { + throw std::runtime_error( + "Matrix dimension not divisible by number of processes."); + } + + const int numRowsLocal = n / numProc; + const int numRowsTotal = n; + const int numCols = n; + + // const int numCols = n; + const int i0 = rank * numRowsLocal; + Matrix m = Matrix::zeros(numRowsLocal, numCols); + for (int i = 0; i < numRowsLocal; ++i) { + const double x = 1. * (i + i0) / numRowsTotal; + for (int j = 0; j < numCols; ++j) { + const double y = 1. * j / numCols; + + // The lower left corner is hot + m(i, j) += exp(-0.5 * (x * x + y * y) * 10) * exp(-100 * x * y); + + // The right side is cold + m(i, j) += -1.0 * exp(-20 * (1 - y)); + } + } + + return m; +} + +Matrix getDistributedInitialCondition(int n) { + int rank, numProc; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &numProc); + + return getInitialCondition(n, rank, numProc); +} + +Matrix getCompleteInitialCondition(int n) { + return getInitialCondition(n, 0, 1); +} + +Matrix run(Jacobi& jacobi, const Matrix& init, double eps, int maxNumIter) { + auto [m, dist, nIter] = jacobi.run(init, eps, maxNumIter); + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + if (rank == 0) { + std::cout << "Finished running Jacobi, nIter=" << nIter << ", dist=" << dist + << std::endl; + } + return m; +} + +double benchmark(Jacobi& jacobi, + const Matrix& init, + double eps, + int maxNumIter) { + auto start = std::chrono::system_clock::now(); + jacobi.run(init, eps, maxNumIter); + auto end = std::chrono::system_clock::now(); + + std::chrono::duration diff = end - start; + return diff.count(); +} + +double serialBenchmark(const int n, double eps, int maxNumIter) { + // Benchmark the serial version + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + double serialTime = 0; + if (rank == 0) { + std::cout << "Starting serial benchmark" << std::endl; + JacobiSerial jacobiSerial; + serialTime = benchmark(jacobiSerial, getCompleteInitialCondition(n), eps, + maxNumIter); + + std::cout << "Serial benchmark finished, time=" << serialTime << "ms" + << std::endl; + } + MPI_Bcast(&serialTime, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD); + return serialTime; +} + +double parallelBenchmark(const int n, double eps, int maxNumIter) { + JacobiMPI jacobi; + Matrix init = getDistributedInitialCondition(n); + double time = benchmark(jacobi, init, eps, maxNumIter); + + int rank; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + + if (rank == 0) { + std::cout << "Parallel benchmark finished, time=" << time << "ms" + << std::endl; + } + return time; +} + +void runSerial(int n, double eps, int maxNumIter) { + int rank, numProc; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &numProc); + + // To run the serial version, we just execute the code on rank 0, + // all other ranks are idle. + if (rank == 0) { + std::cout << "Starting serial Jacobi run" << std::endl; + Matrix init = getCompleteInitialCondition(n); + const std::string initFilename = "serial_init.asc"; + MatrixIO().saveSerial(init, initFilename); + JacobiSerial jacobi; + Matrix phi = run(jacobi, init, eps, maxNumIter); + const std::string resultFilename = "result_serial.asc"; + MatrixIO().saveSerial(phi, resultFilename); + std::cout << "Finished serial Jacobi run" << std::endl << std::endl; + } +} + +void runParallel(int n, double eps, int maxNumIter) { + int rank, numProc; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &numProc); + + if (rank == 0) { + std::cout << "Starting parallel Jacobi run" << std::endl; + } + + Matrix init = getDistributedInitialCondition(n); + + const std::string initFilename = "parallel_init.asc"; + MatrixIO().saveDistributed(init, initFilename); + + JacobiMPI jacobi; + Matrix phi = run(jacobi, init, eps, maxNumIter); + const std::string resultFilename = "result_parallel.asc"; + MatrixIO().saveDistributed(phi, resultFilename); + + if (rank == 0) { + std::cout << "Finished parallel Jacobi run" << std::endl; + } +} + +int main(int argc, char* argv[]) { + MPI_Init(&argc, &argv); + + int rank, numProc; + MPI_Comm_rank(MPI_COMM_WORLD, &rank); + MPI_Comm_size(MPI_COMM_WORLD, &numProc); + + // Set the Jacobi parameters + const int n = 9 * 64; + const double eps = 1e-7; + const int maxNumIter = 10000; + + // Run the serial and parallel versions + // Result is saved to file. + // Use this to graphically verify the correctness of the parallel + // implementation. + runSerial(n, eps, maxNumIter); + + runParallel(n, eps, maxNumIter); + + + // Run the benchmark + double serialTime = 0; + double parallelTime = 0; + + serialTime = serialBenchmark(n, eps, maxNumIter); + parallelTime = parallelBenchmark(n, eps, maxNumIter); + + if (rank == 0) { + std::cout << "Serial time: " << serialTime << "ms" << std::endl; + std::cout << "Parallel time: " << parallelTime << "ms" << std::endl; + std::cout << "Speedup: " << serialTime / parallelTime << std::endl; + std::ofstream fout("benchmark.txt", std::ios::app); + fout << numProc << "\t" << parallelTime << "\n"; + } + + MPI_Finalize(); +} \ No newline at end of file diff --git a/lab10/jacobi/jacobi_main.h b/lab10/jacobi/jacobi_main.h new file mode 100644 index 0000000..c65f548 --- /dev/null +++ b/lab10/jacobi/jacobi_main.h @@ -0,0 +1,74 @@ +#ifndef JACOBI_MAIN_H +#define JACOBI_MAIN_H + +#include +#include + +#include +#include + +#include "matrix.h" +#include "jacobi.h" +#include "jacobi_serial.h" +#include "jacobi_mpi.h" +#include "matrix_io.h" + +/** + * Generate initial condition for Jacobi iteration. + * Data is already appropriately distributed. + */ +Matrix getInitialCondition(int n, int rank, int numProc); + +/** + * Generate initial condition for Jacobi iteration. + * Data is distriburted along the first axis. + */ +Matrix getDistributedInitialCondition(int n); + +/** + * Generate initial condition for Jacobi iteration. + * Data is not distributed, i.e. all ranks have the same data. + */ +Matrix getCompleteInitialCondition(int n); + +/** + * Create an initial condition and run Jacobi algorithm until convergence. + * @return The converged steady state. + */ +Matrix run(Jacobi& jacobi, + const Matrix& init, + double eps, + int maxNumIter); + +/** + * Perform a single benchmark run: create a initial condition and measure + * the time until Jacobi converges. + * + * Allows to compare the performance of the serial and parallel + * Jacobi implementations by taking an instance of the Jacobi base class. + * + * @return time for jacobi convergence in ms + */ +double benchmark(Jacobi& jacobi, + const Matrix& init, + double eps, + int maxNumIter); +/** + * Run the serial Jacobi algorithm. + * + * Create the initial condition, run the Jacobi algorithm until convergence, + * and save the result to a file. + */ +void runSerial(int n, double eps, int maxNumIter); + +/** + * Run the parallel Jacobi algorithm. + * + * Create the initial condition, run the Jacobi algorithm until convergence, + * and save the result to a file. + */ +void runParallel(int n, double eps, int maxNumIter); + +int main(int argc, char* argv[]); + +#endif // JACOBI_MAIN_H \ No newline at end of file diff --git a/lab10/jacobi/jacobi_mpi.cpp b/lab10/jacobi/jacobi_mpi.cpp new file mode 100644 index 0000000..e69de29 diff --git a/lab10/jacobi/jacobi_mpi.h b/lab10/jacobi/jacobi_mpi.h new file mode 100644 index 0000000..6939da8 --- /dev/null +++ b/lab10/jacobi/jacobi_mpi.h @@ -0,0 +1,13 @@ +#ifndef JACOBI_MPI_H +#define JACOBI_MPI_H + +#include "jacobi.h" + +class JacobiMPI : public Jacobi { + public: + Result run(const Matrix& init, double epsilon, int maxNumIter) { + return Result{}; + } +}; + +#endif // JABOBI_MPI_H diff --git a/lab10/jacobi/jacobi_serial.cpp b/lab10/jacobi/jacobi_serial.cpp new file mode 100644 index 0000000..fa7f470 --- /dev/null +++ b/lab10/jacobi/jacobi_serial.cpp @@ -0,0 +1,41 @@ +#include +#include +#include +#include +#include +#include "matrix.h" +#include "jacobi_serial.h" + +Jacobi::Result JacobiSerial::run(const Matrix& init, + double eps, + int maxNumIter) { + std::vector phi(2, init); + const int numRows = phi[0].rows(); + const int numCols = phi[0].cols(); + + double dist = std::numeric_limits::max(); + int nIter = 0; + + int t0 = 0; // array index of current timestep + int t1 = 1; // array index of next timestep + while (dist > eps && nIter < maxNumIter) { + dist = 0; + for (int i = 1; i < numRows - 1; ++i) { + for (int j = 1; j < numCols - 1; ++j) { + phi[t1](i, j) = .25 * (phi[t0](i + 1, j) + phi[t0](i - 1, j) + + phi[t0](i, j + 1) + phi[t0](i, j - 1)); + const double diff = phi[t1](i, j) - phi[t0](i, j); + dist = std::max(dist, std::abs(diff)); + } + } + + // if (nIter % 1000 == 0) { + // std::cout << "Iteration " << nIter << ", dist=" << dist << "\n"; + // } + + nIter++; + std::swap(t0, t1); + } + + return Result{phi[t0], dist, nIter}; +} diff --git a/lab10/jacobi/jacobi_serial.h b/lab10/jacobi/jacobi_serial.h new file mode 100644 index 0000000..f7cae3e --- /dev/null +++ b/lab10/jacobi/jacobi_serial.h @@ -0,0 +1,12 @@ +#ifndef JACOBI_SERIAL_H +#define JACOBI_SERIAL_H + + +#include "jacobi.h" + +class JacobiSerial : public Jacobi { + public: + Result run(const Matrix& phi, double epsilon, int maxNumIter); +}; + +#endif // JABOBI_SERIAL_H diff --git a/lab10/jacobi/matrix.h b/lab10/jacobi/matrix.h new file mode 100644 index 0000000..fd3df1d --- /dev/null +++ b/lab10/jacobi/matrix.h @@ -0,0 +1,411 @@ +/* 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 +#include +#include +#include + +/** + * 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 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 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 dim() const; + int rows() const; + int cols() const; + int numEntries() const; + + // Comparison operators + bool operator==(const Matrix& b); + bool operator!=(const Matrix& b); + + // 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) { return data_ == b.data_; } + bool operator!=(const Vector& b) { 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 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 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) { + return data_[i * numCols_ + j]; +} + +inline const double& Matrix::operator()(int i, int j) const { + return data_[i * numCols_ + j]; +} + +inline double* Matrix::data() { + return data_; +} + +inline const double* Matrix::data() const { + return data_; +} + +inline std::pair Matrix::dim() const { + return std::pair(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 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) { + 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) { + 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 \ No newline at end of file diff --git a/lab10/jacobi/matrix_io.cpp b/lab10/jacobi/matrix_io.cpp new file mode 100644 index 0000000..b5d9708 --- /dev/null +++ b/lab10/jacobi/matrix_io.cpp @@ -0,0 +1,116 @@ +#include +#include +#include +#include "matrix.h" +#include "matrix_io.h" + +MatrixIO::MatrixIO() { + MPI_Comm_rank(MPI_COMM_WORLD, &rank_); + MPI_Comm_size(MPI_COMM_WORLD, &numProc_); +} + +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::loadDistributed(const std::string& filename) { + Matrix matrixOnRoot = Matrix::zeros(0, 0); + 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 numCols = distributedMatrix.cols(); + const int numRowsTotal = numRowsLocal * numProc_; + + Matrix matrixOnRoot = Matrix::zeros(0, 0); + if (rank_ == 0) { + matrixOnRoot = Matrix::zeros(numRowsTotal, numCols); + } + + const int sendSize = distributedMatrix.numEntries(); + + MPI_Gather(&distributedMatrix(0, 0), sendSize, MPI_DOUBLE, + &matrixOnRoot(0, 0), sendSize, MPI_DOUBLE, 0, MPI_COMM_WORLD); + + return matrixOnRoot; +} + +Matrix MatrixIO::scatterMatrixFromRoot(const Matrix& matrixOnRoot) { + int numRowsTotal = matrixOnRoot.rows(); + int numCols = matrixOnRoot.cols(); + + MPI_Bcast(&numRowsTotal, 1, MPI_INT, 0, MPI_COMM_WORLD); + MPI_Bcast(&numCols, 1, MPI_INT, 0, MPI_COMM_WORLD); + + if (numRowsTotal % numProc_ != 0) + throw std::runtime_error( + "Number of rows not divisible by number of processes."); + + const int numRowsLocal = numRowsTotal / numProc_; + + Matrix distributedMatrix = Matrix::zeros(numRowsLocal, numCols); + + const int sendSize = distributedMatrix.numEntries(); + + MPI_Scatter(matrixOnRoot.data(), sendSize, MPI_DOUBLE, + distributedMatrix.data(), sendSize, MPI_DOUBLE, 0, + MPI_COMM_WORLD); + + const int i0 = rank_ * numRowsLocal; + + for (int i = 0; i < numRowsLocal; ++i) { + std::cout << "Scattered Line " << i0 + i << ": " + << distributedMatrix(i, numCols - 1) << "\n"; + } + + return distributedMatrix; +} diff --git a/lab10/jacobi/matrix_io.h b/lab10/jacobi/matrix_io.h new file mode 100644 index 0000000..1560e93 --- /dev/null +++ b/lab10/jacobi/matrix_io.h @@ -0,0 +1,51 @@ +#include "matrix.h" +#include +#include +#include + +/** + * File i/o for distributed matrices. + * Matrices may be distributed along the first dimension (i.e. rows are + * distributed). Only works if number of rows is divisible by number of + * processes. + * + * 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(); + + /** + * 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 loadDistributed(const std::string& filename); + + private: + Matrix gatherMatrixOnRoot(const Matrix& distributedMatrix); + + Matrix scatterMatrixFromRoot(const Matrix& matrixOnRoot); + + int rank_ = MPI_PROC_NULL; + int numProc_ = 0; +}; diff --git a/lab10/jacobi/plot.py b/lab10/jacobi/plot.py new file mode 100644 index 0000000..0a007b5 --- /dev/null +++ b/lab10/jacobi/plot.py @@ -0,0 +1,24 @@ +import pylab as plt + + +ser_init = plt.loadtxt("serial_init.asc") +par_init = plt.loadtxt("parallel_init.asc") +ser_data = plt.loadtxt("result_serial.asc") +par_data = plt.loadtxt("result_parallel.asc") +# diff = par_data - ser_data + +plt.figure() +plt.subplot(221) +plt.title("Serial initial condition") +plt.contourf(ser_init, cmap='jet', levels=100) +plt.subplot(222) +plt.title("Parallel initial condition") +plt.contourf(par_init, cmap='jet', levels=100) +plt.subplot(223) +plt.title("Serial result") +plt.contourf(ser_data, cmap='jet', levels=100) +plt.subplot(224) +plt.title("Parallel result") +plt.contourf(par_data, cmap="jet", levels=100) + +plt.show()