Commit d4e98e78 authored by Alexandru Dura's avatar Alexandru Dura

Distributed the Diffusion problem using MPI

parent 07b8c90d
......@@ -6,73 +6,33 @@
#include <limits>
#include "Matrix.h"
class RowView {
double * const row;
RowView(double *row) : row(row) {}
friend class Field;
public:
double &operator[](unsigned i) {
return row[i];
}
double operator[](unsigned i) const {
return row[i];
}
double* raw() {
return row;
}
};
class Field {
class Field : public Matrix<double> {
const unsigned xdim, ydim;
const unsigned m, n;
double *field;
void init() {
field = new double[m * n];
}
public:
Field(int xdim, int ydim) : xdim(xdim), ydim(ydim), m(ydim + 2), n(xdim + 2) {
init();
}
~Field() {
delete[] field;
}
RowView operator[](unsigned i) {
return RowView(field + i * n);
}
RowView operator[](unsigned i) const {
return RowView(field + i * n);
Field(int xdim, int ydim) : Matrix<double>(ydim + 2, xdim + 2), xdim(xdim), ydim(ydim) {
}
RowView firstRow() {
return operator[](1);
Row<double> firstRow() {
return (*this)[1];
}
RowView lastRow() {
return operator[](m-2);
Row<double> lastRow() {
return (*this)[ydim];
}
RowView firstHalo() {
return operator[](0);
Row<double> firstHalo() {
return (*this)[0];
}
RowView lastHalo() {
return operator[](m-1);
Row<double> lastHalo() {
return (*this)[ydim + 1];
}
void print(std::ostream &os) {
for (int i = 0; i < m; ++i) {
for (int j = 0; j < n; ++j) {
os << *(field + i * n + j) << " ";
for (int i = 0; i < numRows(); ++i) {
for (int j = 0; j < numCols(); ++j) {
os << (*this)[i][j] << " ";
}
os << "\n";
}
......
......@@ -2,6 +2,7 @@
#include "Field.h"
#include "Util.h"
#include "Problem.h"
#include <cmath>
#define DEFAULT_COMM MPI_COMM_WORLD
......@@ -21,20 +22,59 @@ int main(int argc, char **argv) {
unsigned localXdim = xdim;
unsigned localYdim = ydim / nsize;
int upNeighbour = rank + 1;
int dnNeighbour = rank - 1;
double gateVal = 100;
unsigned nGates = 2;
Problem p(nsize, rank, localXdim, localYdim, nGates, nGates, gateVal);
Field &f = p.getField();
for (int i = 0; i < maxiter; ++i) {
MPI_Request request_up, request_dn;
if (dnNeighbour >= 0)
MPI_Isend(f.firstRow().raw(), f.numCols(), MPI_DOUBLE, dnNeighbour, 0, DEFAULT_COMM, &request_dn);
if (upNeighbour < nsize)
MPI_Isend(f.lastRow().raw(), f.numCols(), MPI_DOUBLE, upNeighbour, 0, DEFAULT_COMM, &request_up);
if (dnNeighbour >= 0)
MPI_Recv(f.firstHalo().raw(), f.numCols(), MPI_DOUBLE, dnNeighbour, 0, DEFAULT_COMM, MPI_STATUS_IGNORE);
if (upNeighbour < nsize)
MPI_Recv(f.lastHalo().raw(), f.numCols(), MPI_DOUBLE, upNeighbour, 0, DEFAULT_COMM, MPI_STATUS_IGNORE);
if (upNeighbour < nsize)
MPI_Wait(&request_up, MPI_STATUS_IGNORE);
if (dnNeighbour >= 0)
MPI_Wait(&request_dn, MPI_STATUS_IGNORE);
p.step(cx, cy);
}
double resNorm2 = p.residualSquared(cx, cy);
int bitdepth = 256;
auto img = p.getField().discretize(bitdepth);
write_pgm("out.pgm", img, bitdepth);
auto count = img.numCols() * img.numRows();
Matrix<int> resultImg(ydim, xdim);
MPI_Request gatherImg;
MPI_Igather(img.raw(), count, MPI_INT, resultImg.raw(), count, MPI_INT, 0, DEFAULT_COMM, &gatherImg);
double resNorm2 = p.residualSquared(cx, cy);
double gatherResNorm2;
MPI_Request gatherNorm;
MPI_Ireduce(&resNorm2, &gatherResNorm2, 1, MPI_DOUBLE, MPI_SUM, 0, DEFAULT_COMM, &gatherNorm);
MPI_Wait(&gatherImg, MPI_STATUS_IGNORE);
if (rank == 0) {
write_pgm("polution.pgm", resultImg, bitdepth);
}
MPI_Wait(&gatherNorm, MPI_STATUS_IGNORE);
if (rank == 0) {
printf("Norm of the residual: %f\n", std::sqrt(gatherResNorm2));
}
MPI_Finalize();
return 0;
......
#pragma once
template<typename> class Matrix;
template<typename> class MatrixView;
template<typename T>
class Row {
T * const v;
Row(T *v) : v(v) {}
friend class Matrix<T>;
friend class MatrixView<T>;
public:
T &operator[](unsigned i) {
return v[i];
......@@ -24,15 +25,17 @@ public:
};
template<typename T>
class Matrix {
T *data;
class MatrixView {
unsigned nRows, nCols;
public:
Matrix(unsigned nRows, unsigned nCols) : nRows(nRows), nCols(nCols) {
data = new T[nRows * nCols];
protected:
T *data;
MatrixView(unsigned nRows, unsigned nCols, T *data) : nRows(nRows), nCols(nCols), data(data) {
}
public:
Row<T> operator[](unsigned i) {
return Row<T>(data + i * nCols);
}
......@@ -41,10 +44,26 @@ public:
return Row<T>(data + i * nCols);
}
~Matrix() {
delete[] data;
MatrixView<T> extractRows(unsigned fromRow, unsigned numRows) {
return MatrixView(numRows, nCols, data + fromRow * nCols);
}
unsigned numCols() const { return nCols; };
unsigned numRows() const { return nRows; };
T *raw() { return data; }
};
template<typename T>
class Matrix : public MatrixView<T> {
public:
Matrix(unsigned nRows, unsigned nCols) : MatrixView<T>(nRows, nCols, new T[nRows * nCols]) {
}
virtual ~Matrix() {
delete[] MatrixView<T>::data;
MatrixView<T>::data = nullptr;
}
};
......@@ -31,9 +31,9 @@ void Problem::addGates(Field &f, unsigned xGates, unsigned yGates, double value)
// add the ygates
for (unsigned i = 0; i < yGates; ++i) {
unsigned gateCentre = (1 + 2*i) * ydim / (2 * yGates);
unsigned gateStart = gateCentre - 0.5 * GATE_WIDTH * ydim;
unsigned gateEnd = gateCentre + 0.5 * GATE_WIDTH * ydim;
unsigned gateCentre = (1 + 2*i) * nFields * ydim / (2 * yGates);
unsigned gateStart = gateCentre - 0.5 * GATE_WIDTH * nFields * ydim;
unsigned gateEnd = gateCentre + 0.5 * GATE_WIDTH * nFields * ydim;
if (gateStart >= starty && gateStart < endy) {
for (unsigned j = gateStart; j < std::min(gateEnd, endy); j++) {
......
......@@ -41,7 +41,7 @@ public:
double residualSquared(double cx, double cy) const;
const Field &getField() const {
Field &getField() {
return *currentf;
}
};
#include <cstdio>
#include <cstdlib>
#include <string>
#include "Matrix.h"
void read_param(int* pxdim, int* pydim, double* pcx, double* pcy, int* pmaxiter)
......@@ -42,9 +43,9 @@ void read_param(int* pxdim, int* pydim, double* pcx, double* pcy, int* pmaxiter)
}
void write_pgm(const char *filename, const Matrix<int> &image, int bitwidth) {
void write_pgm(const std::string &filename, const Matrix<int> &image, int bitwidth) {
int i, j;
FILE *fp = fopen( filename, "w");
FILE *fp = fopen(filename.c_str(), "w");
int width = image.numCols();
int height = image.numRows();
......
#pragma once
#include "Matrix.h"
#include <string>
void read_param(int* pxdim, int* pydim, double* pcx, double* pcy, int* pmaxiter);
void write_pgm(const char *filename, const Matrix<int> &image, int bitwidth);
void write_pgm(const std::string &filename, const Matrix<int> &image, int bitwidth);
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment