Commit eb257440 authored by Alexandru Dura's avatar Alexandru Dura

Broacast paramaters and gather the results

parent d4e98e78
......@@ -4,6 +4,7 @@
#include <ostream>
#include <limits>
#include <tuple>
#include "Matrix.h"
class Field : public Matrix<double> {
......@@ -38,17 +39,20 @@ public:
}
}
Matrix<int> discretize(int bitdepth) const {
std::pair<double, double> minMax() const {
double maxval = std::numeric_limits<double>::min();
double minval = std::numeric_limits<double>::max();
Matrix<int> ret(ydim, xdim);
for (int y = 0; y < ydim; ++y) {
for (int x = 0; x < xdim; ++x) {
maxval = std::max(maxval, (*this)[y+1][x+1]);
minval = std::min(minval, (*this)[y+1][x+1]);
}
}
return std::make_pair(minval, maxval);
}
Matrix<int> discretize(int bitdepth, double minval, double maxval) const {
Matrix<int> ret(ydim, xdim);
double norm = (double)(bitdepth - 1) / (maxval - minval);
......@@ -60,4 +64,10 @@ public:
return ret;
}
Matrix<int> discretize(int bitdepth) const {
double minval, maxval;
std::tie(minval, maxval) = minMax();
return discretize(bitdepth, minval, maxval);
}
};
......@@ -3,9 +3,76 @@
#include "Util.h"
#include "Problem.h"
#include <cmath>
#include <cstddef>
#include <cassert>
#define DEFAULT_COMM MPI_COMM_WORLD
struct GlobalProblemConfig {
double cx;
double cy;
int xdim;
int ydim;
int maxiter;
static const int* block_lengths() {
static const int result[] = {
1,
1,
1,
1,
1
};
return result;
}
static const MPI_Aint* displacements() {
static const MPI_Aint result[] = {
offsetof(GlobalProblemConfig, cx),
offsetof(GlobalProblemConfig, cy),
offsetof(GlobalProblemConfig, xdim),
offsetof(GlobalProblemConfig, ydim),
offsetof(GlobalProblemConfig, maxiter)
};
return result;
}
static const MPI_Datatype* types() {
static const MPI_Datatype result[] = {
MPI_DOUBLE,
MPI_DOUBLE,
MPI_INT,
MPI_INT,
MPI_INT
};
return result;
}
};
void exchangeHalos(int rank, int nsize, Problem &p) {
MPI_Request request_up, request_dn;
Field &f = p.getField();
int upNeighbour = rank + 1;
int dnNeighbour = rank - 1;
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);
}
int main(int argc, char **argv) {
MPI_Init(nullptr, nullptr);
......@@ -14,13 +81,23 @@ int main(int argc, char **argv) {
MPI_Comm_size(DEFAULT_COMM, &nsize);
MPI_Comm_rank(DEFAULT_COMM, &rank);
int xdim, ydim, maxiter;
double cx, cy;
GlobalProblemConfig config;
MPI_Datatype GlobalProblemConfig_t;
MPI_Type_create_struct(5, GlobalProblemConfig::block_lengths(), GlobalProblemConfig::displacements(),
GlobalProblemConfig::types(), &GlobalProblemConfig_t);
MPI_Type_commit(&GlobalProblemConfig_t);
if (rank == 0) {
read_param(&config.xdim, &config.ydim, &config.cx, &config.cy, &config.maxiter);
}
MPI_Bcast(&config, 1, GlobalProblemConfig_t, 0, DEFAULT_COMM);
printf("%d: %d, %d\n", rank, config.xdim, config.ydim);
read_param(&xdim, &ydim, &cx, &cy, &maxiter);
unsigned localXdim = config.xdim;
unsigned localYdim = config.ydim / nsize;
unsigned localXdim = xdim;
unsigned localYdim = ydim / nsize;
assert(config.ydim % nsize == 0);
int upNeighbour = rank + 1;
int dnNeighbour = rank - 1;
......@@ -29,48 +106,51 @@ int main(int argc, char **argv) {
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);
for (int i = 0; i < config.maxiter; ++i) {
exchangeHalos(rank, nsize, p);
p.step(config.cx, config.cy);
}
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);
// Compute the global min and max
double localMin, localMax;
double globalMin, globalMax;
p.step(cx, cy);
}
std::tie(localMin, localMax) = p.getField().minMax();
MPI_Allreduce(&localMin, &globalMin, 1, MPI_DOUBLE, MPI_MIN, DEFAULT_COMM);
MPI_Allreduce(&localMax, &globalMax, 1, MPI_DOUBLE, MPI_MAX, DEFAULT_COMM);
// Discretize the local field
int bitdepth = 256;
auto img = p.getField().discretize(bitdepth);
auto img = p.getField().discretize(bitdepth, globalMin, globalMax);
auto count = img.numCols() * img.numRows();
Matrix<int> resultImg(ydim, xdim);
// Gather the discretize fields on rank 0
Matrix<int> resultImg(config.ydim, config.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);
// Compute the squred local residual; for this we need
// the up-to-date halos
exchangeHalos(rank, nsize, p);
double resNorm2 = p.residualSquared(config.cx, config.cy);
double gatherResNorm2;
// Sum up the residuals from all the nodes and put the result
// in rank 0
MPI_Request gatherNorm;
MPI_Ireduce(&resNorm2, &gatherResNorm2, 1, MPI_DOUBLE, MPI_SUM, 0, DEFAULT_COMM, &gatherNorm);
// We're done gathering the images, write it out
MPI_Wait(&gatherImg, MPI_STATUS_IGNORE);
if (rank == 0) {
write_pgm("polution.pgm", resultImg, bitdepth);
}
// We're done gathering the residuals and summing them up
MPI_Wait(&gatherNorm, MPI_STATUS_IGNORE);
if (rank == 0) {
printf("Norm of the residual: %f\n", std::sqrt(gatherResNorm2));
......
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