Studium_Angewandte_Informatik_Bachelor/hpc_laboratory_exercises
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity

add get row to matrix

Browse Source
This commit is contained in:
kai
2025-05-31 12:51:41 +02:00
parent d71c943b89
commit 5224b72924
2 changed files with 112 additions and 159 deletions
Download Patch File Download Diff File Expand all files Collapse all files

133
lab10/jacobi/jacobi_mpi.cpp
View File

@@ -1,35 +1,19 @@
#include "jacobi.h"
#include "jacobi_main.h"
std::vector<double> getFilledBuffer(const Matrix &matrix, int row)
{
std::vector<double> getFilledBuffer(const Matrix &matrix, int row) {
int numCols = matrix.cols();
std::vector<double> 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<double> &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<double> 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<double> 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<double> rowUpper;
std::vector<double> 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);

138
lab10/jacobi/matrix.h
View File

@@ -16,14 +16,13 @@
* limitations under the License.
*/
#ifndef MATRIX_H
#define MATRIX_H
#include <vector>
#include <iostream>
#include <iomanip>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <vector>
/**
* 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<int, int> 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<int, int> 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<double> &get_row(int i) {
return std::vector<double> 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<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 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<int, int> 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<int, int> Matrix::dim() const {
return std::pair<int, int>(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
#endif // MATRIX_H