Commit 07b8c90d authored by Alexandru Dura's avatar Alexandru Dura

Translated the diffusion problem to C++

parent 45b7cb8a
......@@ -75,3 +75,26 @@ Question: analysis tools for avoiding race conditions in OpenMP, or avoiding dea
- ARM DDT parallel debugger
- ~module load arm-forge~
- ~ddt~
* <2020-03-05 Thu>
** MPI Collective Communication
- send information from a single node to an entire communicator
** Working with communicators
- duplicating communicators
- cartesian greeds
- MPI shines with regular grid codes
* <2020-03-10 Tue>
** MPI messages
- use derived data types to build big messages
- prefer large messages to many short ones
** MPI IO
- MPI_file_open
- MPI_file_set_view
* <2020-03-12 Thu>
- last lecture
- don't use mpirun --bind-to-core in the project
- https://lunarc-documentation.readthedocs.io/en/latest/batch_system/#hybrid-jobs-using-threads-within-an-mpi-framework
- 4-5 weeks
......@@ -3,6 +3,8 @@
#pragma once
#include <ostream>
#include <limits>
#include "Matrix.h"
class RowView {
double * const row;
......@@ -75,4 +77,27 @@ public:
os << "\n";
}
}
Matrix<int> discretize(int bitdepth) 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]);
}
}
double norm = (double)(bitdepth - 1) / (maxval - minval);
for (int y = 0; y < ydim; ++y) {
for (int x = 0; x < xdim; ++x) {
ret[y][x] = norm * ((*this)[y+1][x+1] - minval);
}
}
return ret;
}
};
......@@ -6,14 +6,13 @@
#define DEFAULT_COMM MPI_COMM_WORLD
int main(int argc, char **argv) {
MPI::Init();
MPI_Init(nullptr, nullptr);
int nsize, rank;
MPI_Comm_size(DEFAULT_COMM, &nsize);
MPI_Comm_rank(DEFAULT_COMM, &rank);
int xdim, ydim, maxiter;
double cx, cy;
......@@ -27,9 +26,16 @@ int main(int argc, char **argv) {
Problem p(nsize, rank, localXdim, localYdim, nGates, nGates, gateVal);
for (int i = 0; i < maxiter; ++i) {
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);
MPI::Finalize();
MPI_Finalize();
return 0;
}
ENTRY_POINTS=Main.cpp Test.cpp
CPP_SOURCE_FILES = $(filter-out $(ENTRY_POINTS), $(wildcard *.cpp))
C_SOURCE_FILES = $(wildcard *.c)
......@@ -11,6 +10,8 @@ CXX=g++
MPICXX=mpic++
CXX_FLAGS=-O3 -I/usr/lib/x86_64-linux-gnu/openmpi/include/
main-debug : Main.cpp $(CPP_SOURCE_FILES) $(HEADER_FILES)
$(CXX) $(MPIFLAGS) -O0 -g -Wall Main.cpp $(CPP_SOURCE_FILES) $(MPI_LD_FLAGS) -o $@
main : Main.cpp $(CPP_SOURCE_FILES) $(HEADER_FILES)
$(CXX) $(MPIFLAGS) -O3 Main.cpp $(CPP_SOURCE_FILES) $(MPI_LD_FLAGS) -o main
......
#pragma once
template<typename> class Matrix;
template<typename T>
class Row {
T * const v;
Row(T *v) : v(v) {}
friend class Matrix<T>;
public:
T &operator[](unsigned i) {
return v[i];
}
T &operator[](unsigned i) const {
return v[i];
}
T *raw() {
return v;
}
};
template<typename T>
class Matrix {
T *data;
unsigned nRows, nCols;
public:
Matrix(unsigned nRows, unsigned nCols) : nRows(nRows), nCols(nCols) {
data = new T[nRows * nCols];
}
Row<T> operator[](unsigned i) {
return Row<T>(data + i * nCols);
}
Row<T> operator[](unsigned i) const {
return Row<T>(data + i * nCols);
}
~Matrix() {
delete[] data;
}
unsigned numCols() const { return nCols; };
unsigned numRows() const { return nRows; };
};
......@@ -37,7 +37,7 @@ void Problem::addGates(Field &f, unsigned xGates, unsigned yGates, double value)
if (gateStart >= starty && gateStart < endy) {
for (unsigned j = gateStart; j < std::min(gateEnd, endy); j++) {
f[j - starty + 1][xdim] = value;
f[j - starty + 1][xdim + 1] = value;
}
}
}
......@@ -58,7 +58,7 @@ void Problem::step(double cx, double cy) {
std::swap(this->nextf, this->currentf);
}
double Problem::residual(double cx, double cy) const {
double Problem::residualSquared(double cx, double cy) const {
double res_sqr = 0.0;
const auto &solution = *this->currentf;
......
......@@ -39,6 +39,9 @@ public:
void step(double cx, double cy);
double residual(double cx, double cy) const;
double residualSquared(double cx, double cy) const;
const Field &getField() const {
return *currentf;
}
};
#include <cstdio>
#include <cstdlib>
#include "Matrix.h"
void read_param(int* pxdim, int* pydim, double* pcx, double* pcy, int* pmaxiter)
{
......@@ -40,3 +41,19 @@ void read_param(int* pxdim, int* pydim, double* pcx, double* pcy, int* pmaxiter)
fclose(fp);
}
void write_pgm(const char *filename, const Matrix<int> &image, int bitwidth) {
int i, j;
FILE *fp = fopen( filename, "w");
int width = image.numCols();
int height = image.numRows();
fprintf( fp, "P2\n");
fprintf( fp, "%i %i\n", width, height );
fprintf( fp, "%i \n", bitwidth-1 );
for (i=0; i<height; i++)
for (j=0; j<width; j++)
fprintf( fp, "%i ", image[i][j]);
fclose(fp);
}
#pragma once
#include "Matrix.h"
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);
64 80
0.4 -0.4
1000
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