diff --git a/lab10/jacobi/jacobi_mpi.cpp b/lab10/jacobi/jacobi_mpi.cpp index 4ad80e2..42e49f5 100644 --- a/lab10/jacobi/jacobi_mpi.cpp +++ b/lab10/jacobi/jacobi_mpi.cpp @@ -1,35 +1,19 @@ #include "jacobi.h" #include "jacobi_main.h" -std::vector getFilledBuffer(const Matrix &matrix, int row) -{ +std::vector getFilledBuffer(const Matrix &matrix, int row) { int numCols = matrix.cols(); std::vector rowBuffer(numCols); // Extract the row (matrix is column-major) - for (int j = 0; j < numCols; ++j) - { + for (int j = 0; j < numCols; ++j) { rowBuffer[j] = matrix(row, j); } return rowBuffer; } -/** - * Takes a matrix and an index indexRowLocal which defines which row of the local matrix has to be sent. - * - * Sends data of that row via MPI_Send to process with rank receiverRank. - * - * indexRowGlobal identifies row in the view of the whole matrix, used from receiving process - * to identify if lower or upper halo by comparing against rowHaloUpperIndex and rowHaloLowerIndex. - */ -void sendLocalRow(MPI_Request &request, std::vector &row, const Matrix &matrix, const int indexRowLocal, const int receiverRank, const int indexRowGlobal) -{ - row = getFilledBuffer(matrix, indexRowLocal); - MPI_Isend(row.data(), row.size(), MPI_DOUBLE, receiverRank, indexRowGlobal, MPI_COMM_WORLD, &request); -} - -Jacobi::Result JacobiMPI::run(const Matrix &init, double epsilon, int maxNumIter) -{ +Jacobi::Result JacobiMPI::run(const Matrix &init, double epsilon, + int maxNumIter) { int rank, numProc; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &numProc); @@ -56,28 +40,26 @@ Jacobi::Result JacobiMPI::run(const Matrix &init, double epsilon, int maxNumIter int t0 = 0; // array index of current timestep int t1 = 1; // array index of next timestep - while (dist > epsilon && nIter < maxNumIter) - { + while (dist > epsilon && nIter < maxNumIter) { dist = 0; - if (rank == 0) - { - std::vector row; - MPI_Request request; - sendLocalRow(request, row, phi[t0], numRows - 1, neighborLower, indexRowGlobalEnd); - MPI_Wait(&request, MPI_STATUS_IGNORE); + if (rank == 0) { + MPI_Request send_lower; + MPI_Isend(phi[t1].get_row(numRows - 1), numCols, MPI_DOUBLE, + neighborUpper, 0, MPI_COMM_WORLD, &send_lower); MPI_Request requestLower; - MPI_Irecv(haloLower.data(), haloLower.size(), MPI_DOUBLE, - neighborLower, rowHaloLowerIndex, MPI_COMM_WORLD, &requestLower); - MPI_Wait(&requestLower, MPI_STATUS_IGNORE); + MPI_Irecv(haloLower.data(), numCols, MPI_DOUBLE, neighborLower, 0, + MPI_COMM_WORLD, &requestLower); - for (int i = 1; i < numRows; ++i) - { - for (int j = 1; j < numCols - 1; ++j) - { + MPI_Wait(&send_lower, MPI_STATUS_IGNORE); + MPI_Wait(&requestUpper, MPI_STATUS_IGNORE); + + for (int i = 1; i < numRows; ++i) { + for (int j = 1; j < numCols - 1; ++j) { double valueBottom; - (i == numRows - 1) ? valueBottom = haloLower[j] : valueBottom = phi[t0](i + 1, j); + (i == numRows - 1) ? valueBottom = haloLower[j] + : valueBottom = phi[t0](i + 1, j); phi[t1](i, j) = .25 * (valueBottom + 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); @@ -86,23 +68,20 @@ Jacobi::Result JacobiMPI::run(const Matrix &init, double epsilon, int maxNumIter } } - else if (rank == (numProc - 1)) - { - std::vector row; + else if (rank == (numProc - 1)) { + MPI_Request send_upper; + MPI_Isend(phi[t1].get_row(rank * numRows), numCols, MPI_DOUBLE, + neighborLower, 0, MPI_COMM_WORLD, &send_upper); - MPI_Request request; - sendLocalRow(request, row, phi[t0], 0, neighborUpper, indexRowGlobalStart); - MPI_Wait(&request, MPI_STATUS_IGNORE); + MPI_Request request_upper; + MPI_Irecv(haloUpper.data(), numCols, MPI_DOUBLE, neighborUpper, 0, + MPI_COMM_WORLD, &request_upper); - MPI_Request requestUpper; - MPI_Irecv(haloUpper.data(), haloUpper.size(), MPI_DOUBLE, - neighborUpper, rowHaloUpperIndex, MPI_COMM_WORLD, &requestUpper); - MPI_Wait(&requestUpper, MPI_STATUS_IGNORE); + MPI_Wait(&send_upper, MPI_STATUS_IGNORE); + MPI_Wait(&request_upper, MPI_STATUS_IGNORE); - for (int i = 0; i < numRows - 1; ++i) - { - for (int j = 1; j < numCols - 1; ++j) - { + for (int i = 0; i < numRows - 1; ++i) { + for (int j = 1; j < numCols - 1; ++j) { double valueTop; (i == 0) ? valueTop = haloUpper[j] : valueTop = phi[t0](i - 1, j); phi[t1](i, j) = .25 * (phi[t0](i + 1, j) + valueTop + @@ -113,51 +92,41 @@ Jacobi::Result JacobiMPI::run(const Matrix &init, double epsilon, int maxNumIter } } - else - { - std::vector rowUpper; - std::vector rowLower; + else { + MPI_Request send_upper; + MPI_Request send_lower; + MPI_Isend(phi[t1].get_row(rank * numRows), numCols, MPI_DOUBLE, + neighborLower, 0, MPI_COMM_WORLD, &send_upper); + MPI_Isend(phi[t1].get_row(rank * numRows + numRows), numCols, MPI_DOUBLE, + neighborLower, 0, MPI_COMM_WORLD, &send_lower); - MPI_Request requestSendUpper; - MPI_Request requestSendLower; - sendLocalRow(requestSendUpper, rowUpper, phi[t0], 0, neighborUpper, indexRowGlobalStart); - sendLocalRow(requestSendLower, rowLower, phi[t0], numRows - 1, neighborLower, indexRowGlobalEnd); + MPI_Request request_upper; + MPI_Request request_lower; + MPI_Irecv(haloUpper.data(), numCols, MPI_DOUBLE, neighborUpper, 0, + MPI_COMM_WORLD, &request_upper); + MPI_Irecv(haloLower.data(), numCols, MPI_DOUBLE, neighborLower, 0, + MPI_COMM_WORLD, &request_lower); - MPI_Request sendRequests[2] = {requestSendUpper, requestSendLower}; - MPI_Waitall(2, sendRequests, MPI_STATUSES_IGNORE); + MPI_Wait(&send_upper, MPI_STATUS_IGNORE); + MPI_Wait(&send_lower, MPI_STATUS_IGNORE); + MPI_Wait(&request_upper, MPI_STATUS_IGNORE); + MPI_Wait(&request_lower, MPI_STATUS_IGNORE); - MPI_Request requestUpper; - MPI_Irecv(haloUpper.data(), haloUpper.size(), MPI_DOUBLE, - neighborUpper, rowHaloUpperIndex, MPI_COMM_WORLD, &requestUpper); - MPI_Request requestLower; - MPI_Irecv(haloLower.data(), haloLower.size(), MPI_DOUBLE, - neighborLower, rowHaloLowerIndex, MPI_COMM_WORLD, &requestLower); - - MPI_Request requests[2] = {requestUpper, requestLower}; - MPI_Waitall(2, requests, MPI_STATUSES_IGNORE); - - for (int i = 0; i < numRows; ++i) - { - for (int j = 1; j < numCols - 1; ++j) - { - if (i == 0) - { + for (int i = 0; i < numRows; ++i) { + for (int j = 1; j < numCols - 1; ++j) { + if (i == 0) { double valueTop = haloUpper[j]; phi[t1](i, j) = .25 * (phi[t0](i + 1, j) + valueTop + 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)); - } - else if (i == numRows - 1) - { + } else if (i == numRows - 1) { double valueBottom = haloLower[j]; phi[t1](i, j) = .25 * (valueBottom + 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)); - } - else - { + } else { 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); diff --git a/lab10/jacobi/matrix.h b/lab10/jacobi/matrix.h index fd3df1d..5cb1eb3 100644 --- a/lab10/jacobi/matrix.h +++ b/lab10/jacobi/matrix.h @@ -16,14 +16,13 @@ * limitations under the License. */ - #ifndef MATRIX_H #define MATRIX_H -#include -#include -#include #include +#include +#include +#include /** * Matrix class @@ -40,7 +39,7 @@ * */ class Matrix { - public: +public: /** * Create a matrix of size m x n and initialize all entries to zero. * @param rows number of rows @@ -79,19 +78,19 @@ class Matrix { */ static Matrix uninit(std::pair dim); - Matrix(const Matrix& other); - Matrix(Matrix&& other); + Matrix(const Matrix &other); + Matrix(Matrix &&other); ~Matrix(); - Matrix& operator=(const Matrix& other); - Matrix& operator=(Matrix&& other); + 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; + 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; + double *data(); + const double *data() const; // Getter functions for the dimensions std::pair dim() const; @@ -100,44 +99,48 @@ class Matrix { int numEntries() const; // Comparison operators - bool operator==(const Matrix& b); - bool operator!=(const Matrix& b); + bool operator==(const Matrix &b); + bool operator!=(const Matrix &b); // addition - Matrix& operator+=(const Matrix& b); + Matrix &operator+=(const Matrix &b); // subtraction - Matrix& operator-=(const Matrix& b); + Matrix &operator-=(const Matrix &b); // scalar multiplication - Matrix& operator*=(double x); + Matrix &operator*=(double x); // scalar division - Matrix& operator/=(double x); + Matrix &operator/=(double x); - private: + std::vector &get_row(int i) { + return std::vector slice(i * this->cols(), this->cols()); + } + +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 + 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: +public: static Vector zeros(int n) { Vector vec; vec.data_ = Matrix::zeros(n, 1); @@ -150,19 +153,19 @@ class Vector { 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(); } + 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) { + Vector operator+=(const Vector &b) { data_ += b.data_; return *this; } - Vector operator-=(const Vector& b) { + Vector operator-=(const Vector &b) { data_ -= b.data_; return *this; } @@ -177,7 +180,7 @@ class Vector { int size() const { return data_.rows(); } - private: +private: Matrix data_ = Matrix::zeros(0, 0); }; @@ -193,27 +196,21 @@ inline Matrix Matrix::zeros(int m, int n) { return mat; } -inline Matrix Matrix::zeros(int n) { - return zeros(n, n); -} +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 m, int n) { return Matrix(m, n); } -inline Matrix Matrix::uninit(int n) { - return uninit(n, 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) { +inline Matrix::Matrix(const Matrix &other) { numRows_ = other.numRows_; numCols_ = other.numCols_; if (numRows_ == 0 || numCols_ == 0) { @@ -228,7 +225,7 @@ inline Matrix::Matrix(const Matrix& other) { } } -inline Matrix::Matrix(Matrix&& other) { +inline Matrix::Matrix(Matrix &&other) { numRows_ = other.numRows_; numCols_ = other.numCols_; data_ = other.data_; @@ -244,7 +241,7 @@ inline Matrix::~Matrix() { } } -inline Matrix& Matrix::operator=(const Matrix& other) { +inline Matrix &Matrix::operator=(const Matrix &other) { if (this != &other) { if (data_ != nullptr) { delete[] data_; @@ -261,7 +258,7 @@ inline Matrix& Matrix::operator=(const Matrix& other) { return *this; } -inline Matrix& Matrix::operator=(Matrix&& other) { +inline Matrix &Matrix::operator=(Matrix &&other) { if (this != &other) { if (data_ != nullptr) { delete[] data_; @@ -276,39 +273,29 @@ inline Matrix& Matrix::operator=(Matrix&& other) { return *this; } -inline double& Matrix::operator()(int i, int j) { +inline double &Matrix::operator()(int i, int j) { return data_[i * numCols_ + j]; } -inline const double& Matrix::operator()(int i, int j) const { +inline const double &Matrix::operator()(int i, int j) const { return data_[i * numCols_ + j]; } -inline double* Matrix::data() { - return data_; -} +inline double *Matrix::data() { return data_; } -inline const double* Matrix::data() const { - 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::rows() const { return numRows_; } -inline int Matrix::cols() const { - return numCols_; -} +inline int Matrix::cols() const { return numCols_; } -inline int Matrix::numEntries() const { - return numRows_ * numCols_; -} +inline int Matrix::numEntries() const { return numRows_ * numCols_; } -inline bool Matrix::operator==(const Matrix& b) { +inline bool Matrix::operator==(const Matrix &b) { const double eps = 1e-12; if (numRows_ != b.numRows_ || numCols_ != b.numCols_) { return false; @@ -323,11 +310,9 @@ inline bool Matrix::operator==(const Matrix& b) { return true; } -inline bool Matrix::operator!=(const Matrix& b) { - return !operator==(b); -} +inline bool Matrix::operator!=(const Matrix &b) { return !operator==(b); } -inline Matrix& Matrix::operator+=(const Matrix& 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); @@ -336,7 +321,7 @@ inline Matrix& Matrix::operator+=(const Matrix& b) { return *this; } -inline Matrix& Matrix::operator-=(const Matrix& 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); @@ -345,7 +330,7 @@ inline Matrix& Matrix::operator-=(const Matrix& b) { return *this; } -inline Matrix& Matrix::operator*=(double x) { +inline Matrix &Matrix::operator*=(double x) { for (int i = 0; i < numRows_; ++i) { for (int j = 0; j < numCols_; ++j) { operator()(i, j) *= x; @@ -354,7 +339,7 @@ inline Matrix& Matrix::operator*=(double x) { return *this; } -inline Matrix& Matrix::operator/=(double x) { +inline Matrix &Matrix::operator/=(double x) { for (int i = 0; i < numRows_; ++i) { for (int j = 0; j < numCols_; ++j) { operator()(i, j) /= x; @@ -370,7 +355,7 @@ inline Matrix::Matrix(int m, int n) : numRows_(m), numCols_(n) { } } -inline std::ostream& operator<<(std::ostream& os, const Matrix& a) { +inline std::ostream &operator<<(std::ostream &os, const Matrix &a) { const int width = 10; const int precision = 4; @@ -389,8 +374,7 @@ inline std::ostream& operator<<(std::ostream& os, const Matrix& a) { return os; } -inline bool equalWithinRange(const Matrix& a, - const Matrix& b, +inline bool equalWithinRange(const Matrix &a, const Matrix &b, double eps = 1e-12) { if (a.rows() != b.rows() || a.cols() != b.cols()) return false; @@ -408,4 +392,4 @@ inline bool equalWithinRange(const Matrix& a, return true; } -#endif // MATRIX_H \ No newline at end of file +#endif // MATRIX_H