diff --git a/Serie06/poisson/solution-asynchronous/AUTHORS b/Serie06/poisson/solution-asynchronous/AUTHORS new file mode 100644 index 0000000..9253635 --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/AUTHORS @@ -0,0 +1,4 @@ +Authors: + - Nicolas Richart + - Vincent Keller + . Vittoria Rezzonico diff --git a/Serie06/poisson/solution-asynchronous/COPYING b/Serie06/poisson/solution-asynchronous/COPYING new file mode 100644 index 0000000..e551330 --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/COPYING @@ -0,0 +1,24 @@ +Copyright (c) 2017, Ecole Polytechnique Federal de Lausanne - Scientific IT and +Application Support (SCITAS) All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Serie06/poisson/solution-asynchronous/Makefile b/Serie06/poisson/solution-asynchronous/Makefile new file mode 100644 index 0000000..3692c48 --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/Makefile @@ -0,0 +1,14 @@ +CXX=mpicxx +LD=${CXX} +CXXFLAGS+=-Wall -Wextra -Werror -pedantic -std=c++11 -O3 -march=native +LDFLAGS+=-lm $(CXXFLAGS) + +OBJS=poisson.o simulation.o double_buffer.o grid.o dumpers.o + +all: poisson + +poisson: $(OBJS) + $(LD) -o $@ $(OBJS) $(LDFLAGS) + +clean: + rm -f hello poisson *.o *~ diff --git a/Serie06/poisson/solution-asynchronous/double_buffer.cc b/Serie06/poisson/solution-asynchronous/double_buffer.cc new file mode 100644 index 0000000..e6ecf0f --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/double_buffer.cc @@ -0,0 +1,19 @@ +/* -------------------------------------------------------------------------- */ +#include "double_buffer.hh" +#include "grid.hh" +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +DoubleBuffer::DoubleBuffer(int m, int n) + : m_current(new Grid(m, n)), m_old(new Grid(m, n)) {} + +/* -------------------------------------------------------------------------- */ +Grid & DoubleBuffer::current() { return *m_current; } + +/* -------------------------------------------------------------------------- */ +Grid & DoubleBuffer::old() { return *m_old; } + +/* -------------------------------------------------------------------------- */ +void DoubleBuffer::swap() { + m_current.swap(m_old); +} diff --git a/Serie06/poisson/solution-asynchronous/double_buffer.hh b/Serie06/poisson/solution-asynchronous/double_buffer.hh new file mode 100644 index 0000000..67a1405 --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/double_buffer.hh @@ -0,0 +1,24 @@ +#ifndef DOUBLE_BUFFER +#define DOUBLE_BUFFER + +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ +#include "grid.hh" +/* -------------------------------------------------------------------------- */ + +class DoubleBuffer { +public: + DoubleBuffer(int m, int n); + + Grid & current(); + Grid & old(); + + void swap(); +private: + + std::unique_ptr m_current; + std::unique_ptr m_old; +}; + +#endif /* DOUBLE_BUFFER */ diff --git a/Serie06/poisson/solution-asynchronous/dumpers.cc b/Serie06/poisson/solution-asynchronous/dumpers.cc new file mode 100644 index 0000000..429400a --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/dumpers.cc @@ -0,0 +1,145 @@ +/* -------------------------------------------------------------------------- */ +#include "dumpers.hh" +#include "grid.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +void Dumper::set_min(float min) { + m_min = min; +} + +void Dumper::set_max(float max) { + m_max = max; +} + +float Dumper::min() const { + return m_min; +} + +float Dumper::max() const { + return m_max; +} + +/* -------------------------------------------------------------------------- */ +void DumperASCII::dump(int step) { + std::ofstream fout; + std::stringstream sfilename; + + sfilename << "output/out_" << std::setfill('0') << std::setw(5) << step << ".pgm"; + fout.open(sfilename.str()); + + int m = m_grid.m(); + int n = m_grid.n(); + + fout << "P2" << std::endl << "# CREATOR: Poisson program" << std::endl; + fout << m << " " << n << std::endl; + fout << 255 << std::endl; + + for (int i = 0; i < m; i++) { + for (int j = 0; j < n; j++) { + int v = 255. * (m_grid(i, j) - m_min) / (m_max - m_min); + v = std::min(v, 255); + fout << v << std::endl; + } + } +} + +/* -------------------------------------------------------------------------- */ +void DumperBinary::dump(int step) { + std::ofstream fout; + std::stringstream sfilename; + + sfilename << "output/out_" << std::setfill('0') << std::setw(5) << step << ".bmp"; + fout.open(sfilename.str(), std::ios_base::binary); + + int prank, psize; + MPI_Comm_rank(MPI_COMM_WORLD, &prank); + MPI_Comm_size(MPI_COMM_WORLD, &psize); + + int nghosts = ((psize == 1) ? 0 : + (prank == 0 or prank == psize - 1) ? 1 : 2); + int h = m_grid.m() - nghosts; + int w = m_grid.n(); + + // remove ghosts at beginning if rank not 0 + int i_start = (prank == 0 ? 0 : 1); + + std::vector buffer; + if (prank == 0) { + buffer.resize(w * w); + MPI_Gather(&m_grid(i_start, 0), h * w, MPI_FLOAT, + buffer.data(), h*w, MPI_FLOAT, 0, MPI_COMM_WORLD); + } else { + MPI_Gather(&m_grid(i_start, 0), h * w, MPI_FLOAT, + NULL, 0, MPI_FLOAT, 0, MPI_COMM_WORLD); + + // it prank is not stop here + return; + } + + + int total_h = h * psize; + int row_size = 3 * w; + // if the file width (3*w) is not a multiple of 4 adds enough bytes to make it + // a multiple of 4 + int padding = (4 - (row_size) % 4) % 4; + row_size += padding; + + int filesize = 54 + (row_size)*total_h; + + std::vector img(row_size*total_h); + std::fill(img.begin(), img.end(), 0); + + + for (int i = 0; i < total_h; i++) { + for (int j = 0; j < w; j++) { + float v = (buffer[i * w + j] - m_min) / (m_max - m_min); + + float r = v * 255; // Red channel + float g = v * 255; // Green channel + float b = v * 255; // Red channel + + r = std::min(r, 255.f); + g = std::min(g, 255.f); + b = std::min(b, 255.f); + + img[row_size * i + 3 * j + 2] = r; + img[row_size * i + 3 * j + 1] = g; + img[row_size * i + 3 * j + 0] = b; + } + } + + std::array bmpfileheader = {'B', 'M', 0, 0, 0, 0, 0, + 0, 0, 0, 54, 0, 0, 0}; + std::array bmpinfoheader = {40, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 1, 0, 24, 0}; + + bmpfileheader[2] = filesize; + bmpfileheader[3] = filesize >> 8; + bmpfileheader[4] = filesize >> 16; + bmpfileheader[5] = filesize >> 24; + + bmpinfoheader[4] = w; + bmpinfoheader[5] = w >> 8; + bmpinfoheader[6] = w >> 16; + bmpinfoheader[7] = w >> 24; + bmpinfoheader[8] = total_h; + bmpinfoheader[9] = total_h >> 8; + bmpinfoheader[10] = total_h >> 16; + bmpinfoheader[11] = total_h >> 24; + bmpinfoheader[20] = (filesize - 54); + bmpinfoheader[21] = (filesize - 54) >> 8; + bmpinfoheader[22] = (filesize - 54) >> 16; + bmpinfoheader[23] = (filesize - 54) >> 24; + + fout.write(bmpfileheader.data(), 14); + fout.write(bmpinfoheader.data(), 40); + + fout.write(img.data(), total_h * row_size); +} diff --git a/Serie06/poisson/solution-asynchronous/dumpers.hh b/Serie06/poisson/solution-asynchronous/dumpers.hh new file mode 100644 index 0000000..03a78c5 --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/dumpers.hh @@ -0,0 +1,41 @@ +#ifndef DUMPERS_HH +#define DUMPERS_HH + +/* -------------------------------------------------------------------------- */ +class Grid; + +/* -------------------------------------------------------------------------- */ +class Dumper { +public: + explicit Dumper(const Grid & grid) : m_grid(grid), m_min(-1.), m_max(1.) {} + + virtual void dump(int step) = 0; + + void set_min(float min); + void set_max(float max); + + float min() const; + float max() const; + +protected: + const Grid & m_grid; + float m_min, m_max; +}; + +/* -------------------------------------------------------------------------- */ +class DumperASCII : public Dumper { +public: + explicit DumperASCII(const Grid & grid) : Dumper(grid) {} + + virtual void dump(int step); +}; + +/* -------------------------------------------------------------------------- */ +class DumperBinary : public Dumper { +public: + explicit DumperBinary(const Grid & grid) : Dumper(grid) {} + + virtual void dump(int step); +}; + +#endif /* DUMPERS_HH */ diff --git a/Serie06/poisson/solution-asynchronous/efficiency.sh b/Serie06/poisson/solution-asynchronous/efficiency.sh new file mode 100644 index 0000000..1184cbd --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/efficiency.sh @@ -0,0 +1,14 @@ +#!/usr/bin/bash -l +#SBATCH -N 2 +#SBATCH -n 56 +#SBATCH -C E5v4 + + + +np=(1 2 4 7 14 28 42 56) +size=(1024 1448 2048 2709 3836 5404 6636 7672) + +for i in 0 1 2 3 4 5 6 7 +do + srun -n ${np[$i]} ./poisson $(( 4 * ${size[$i]} )) >> eff +done diff --git a/Serie06/poisson/solution-asynchronous/grid.cc b/Serie06/poisson/solution-asynchronous/grid.cc new file mode 100644 index 0000000..4e9fb84 --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/grid.cc @@ -0,0 +1,15 @@ +/* -------------------------------------------------------------------------- */ +#include "grid.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +Grid::Grid(int m, int n) : m_m(m), m_n(n), m_storage(m * n) { clear(); } + +/* -------------------------------------------------------------------------- */ +void Grid::clear() { std::fill(m_storage.begin(), m_storage.end(), 0.); } + +/* -------------------------------------------------------------------------- */ +int Grid::m() const { return m_m; } +int Grid::n() const { return m_n; } diff --git a/Serie06/poisson/solution-asynchronous/grid.hh b/Serie06/poisson/solution-asynchronous/grid.hh new file mode 100644 index 0000000..ad872cb --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/grid.hh @@ -0,0 +1,29 @@ +#ifndef GRID_HH +#define GRID_HH + +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +class Grid { +public: + Grid(int m, int n); + + /// access the value [i][j] of the grid + inline float & operator()(int i, int j) { return m_storage[i * m_n + j]; } + inline const float & operator()(int i, int j) const { + return m_storage[i * m_n + j]; + } + + /// set the grid to 0 + void clear(); + + int m() const; + int n() const; + +private: + int m_m, m_n; + std::vector m_storage; +}; + +#endif /* GRID_HH */ diff --git a/Serie06/poisson/solution-asynchronous/output/.to-keep b/Serie06/poisson/solution-asynchronous/output/.to-keep new file mode 100644 index 0000000..e69de29 diff --git a/Serie06/poisson/solution-asynchronous/poisson.cc b/Serie06/poisson/solution-asynchronous/poisson.cc new file mode 100644 index 0000000..692b50e --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/poisson.cc @@ -0,0 +1,70 @@ +/* -------------------------------------------------------------------------- */ +#include "simulation.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +#define EPSILON 0.005 + +typedef std::chrono::high_resolution_clock clk; +typedef std::chrono::duration second; + +static void usage(const std::string & prog_name) { + std::cerr << prog_name << " " << std::endl; + exit(0); +} + +int main(int argc, char * argv[]) { + MPI_Init(&argc, &argv); + int prank, psize; + + MPI_Comm_rank(MPI_COMM_WORLD, &prank); + MPI_Comm_size(MPI_COMM_WORLD, &psize); + + if (argc != 2) usage(argv[0]); + + std::stringstream args(argv[1]); + int N; + args >> N; + + if(args.fail()) usage(argv[0]); + + + // Divide the total size by the number of processes + int m_local = N / psize; + + // Add 1 or 2 ghosts depending on the prank + m_local += + (psize == 1) ? 0 : + (prank == 0 or prank == psize -1) ? 1 : 2; + + Simulation simu(m_local, N); + + simu.set_initial_conditions(); + + simu.set_epsilon(EPSILON); + + float l2; + int k; + + auto start = clk::now(); + std::tie(l2, k) = simu.compute(); + auto end = clk::now(); + + second time = end - start; + + if(prank == 0) + std::cout << psize << " " << N << " " + << k << " " << std::scientific << l2 << " " + << time.count() / k << std::endl; + + MPI_Finalize(); + + return 0; +} diff --git a/Serie06/poisson/solution-asynchronous/scalability.sh b/Serie06/poisson/solution-asynchronous/scalability.sh new file mode 100644 index 0000000..376c6db --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/scalability.sh @@ -0,0 +1,12 @@ +#!/usr/bin/bash -l +#SBATCH -N 2 +#SBATCH -n 56 +#SBATCH -C E5v4 + +for size in 1024 2048 4096 8192 16384 +do + for i in 1 2 4 7 14 28 42 56 + do + srun -n $i ./poisson $size >> async_$size + done +done diff --git a/Serie06/poisson/solution-asynchronous/simulation.cc b/Serie06/poisson/solution-asynchronous/simulation.cc new file mode 100644 index 0000000..e500fa7 --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/simulation.cc @@ -0,0 +1,128 @@ +/* -------------------------------------------------------------------------- */ +#include "simulation.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +Simulation::Simulation(int m, int n) + // Compute h_m and h_n with n since m is different on all processors + : m_local_m(m), m_local_n(n), m_epsilon(1e-7), m_h_m(1. / n), + m_h_n(1. / n), m_grids(m, n), m_f(m, n), + m_dumper(new DumperBinary(m_grids.old())) { + + // retrieving the number of proc and the rank in the proc pool + MPI_Comm_rank(MPI_COMM_WORLD, &m_prank); + MPI_Comm_size(MPI_COMM_WORLD, &m_psize); +} + +/* -------------------------------------------------------------------------- */ +void Simulation::set_initial_conditions() { + // remove ghosts at beginning if rank not 0 + int i_start = (m_prank == 0 ? 0 : 1); + +// remove ghost at the end if prank is not the last + int i_end = (m_prank == m_psize - 1 ? m_local_m : m_local_m - 1); + + int nghosts = + (m_psize == 1) ? 0 : + (m_prank == 0 or m_prank == m_psize - 1 ? 1 : 2); + // computing the offsets of the local grid in the global one + int offset_m = + (m_local_m - nghosts) * m_prank; + int offset_n = 0; + + for (int i = i_start; i < i_end; i++) { + for (int j = 0; j < m_local_n; j++) { + m_f(i, j) = -2. * 100. * M_PI * M_PI * + std::sin(10. * M_PI * (offset_m + i - i_start) * m_h_m) * + std::sin(10. * M_PI * (offset_n + j) * m_h_n); + } + } +} + +/* ------------------------------------------------------------------------65;6203;1c-- */ +std::tuple Simulation::compute() { + int s = 0; + float l2 = 0; + //do { + for(s = 0; s < 100; ++s) { + l2 = compute_step(); + + m_grids.swap(); + } + //++s; + //} while (l2 > m_epsilon); + + //m_dumper->dump(s); + + return std::make_tuple(l2, s); +} + +/* -------------------------------------------------------------------------- */ +void Simulation::set_epsilon(float epsilon) { m_epsilon = epsilon; } + +/* -------------------------------------------------------------------------- */ +float Simulation::epsilon() const { return m_epsilon; } + +/* -------------------------------------------------------------------------- */ +float Simulation::compute_step() { + float l2 = 0.f; + + Grid & u = m_grids.current(); + Grid & uo = m_grids.old(); + + std::array requests; + + // determining the rank of the neighbors + int north_prank = (m_prank == 0 ? MPI_PROC_NULL : m_prank - 1); + int south_prank = (m_prank == (m_psize - 1) ? MPI_PROC_NULL : m_prank + 1); + + // MPI_PROC_NULL allows you to do a send en receive that does nothing + + MPI_Irecv(&uo(0, 0), m_local_n, MPI_FLOAT, north_prank, 0, + MPI_COMM_WORLD, &requests[0]); + MPI_Irecv(&uo(m_local_m - 1, 0), m_local_n, MPI_FLOAT, south_prank, 0, + MPI_COMM_WORLD, &requests[1]); + + MPI_Isend(&uo(1, 0), m_local_n, MPI_FLOAT, north_prank, 0, + MPI_COMM_WORLD, &requests[2]); + MPI_Isend(&uo(m_local_m - 2, 0), m_local_n, MPI_FLOAT, south_prank, 0, + MPI_COMM_WORLD, &requests[3]); + + auto compute_row = [&](int i) { + for (int j = 1; j < m_local_n - 1; j++) { + // computation of the new step + u(i, j) = 0.25 * (uo(i - 1, j) + uo(i + 1, j) + uo(i, j - 1) + + uo(i, j + 1) - m_f(i, j) * m_h_m * m_h_n); + + // L2 norm + l2 += (uo(i, j) - u(i, j)) * (uo(i, j) - u(i, j)); + } + }; + + // compute inner most part + for (int i = 2; i < m_local_m - 2; i++) { + compute_row(i); + } + + // Wait for receives to finish + for (int r = 0; r <2; ++r) { + int i; + MPI_Waitany(2, requests.data(), &i, MPI_STATUSES_IGNORE); + + // compute first or last line depending on with receive is finished + compute_row(i == 0 ? 1 : m_local_m - 2); + } + + // Wait for send to finish + MPI_Waitall(2, requests.data() + 2, MPI_STATUSES_IGNORE); + + // Summing the value of all the processors together + MPI_Allreduce(MPI_IN_PLACE, &l2, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD); + + return l2; +} diff --git a/Serie06/poisson/solution-asynchronous/simulation.hh b/Serie06/poisson/solution-asynchronous/simulation.hh new file mode 100644 index 0000000..be1879e --- /dev/null +++ b/Serie06/poisson/solution-asynchronous/simulation.hh @@ -0,0 +1,54 @@ +#ifndef SIMULATION_HH +#define SIMULATION_HH + +/* -------------------------------------------------------------------------- */ +#include "double_buffer.hh" +#include "dumpers.hh" +#include "grid.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +class Simulation { +public: + Simulation(int m, int n); + + /// set the initial conditions, Dirichlet and source term + virtual void set_initial_conditions(); + + /// perform the simulation + std::tuple compute(); + + /// access the precision + void set_epsilon(float epsilon); + float epsilon() const; + +protected: + /// compute one step and return an error + virtual float compute_step(); + +private: + /// Global problem size + int m_local_m, m_local_n; + + /// Precision to achieve + float m_epsilon; + + /// grid spacing + float m_h_m; + float m_h_n; + + /// Grids storage + DoubleBuffer m_grids; + + /// source term + Grid m_f; + + int m_prank{0}, m_psize{1}; + + /// Dumper to use for outputs + std::unique_ptr m_dumper; +}; + +#endif /* SIMULATION_HH */ diff --git a/Serie06/poisson/solution-sendrecv/AUTHORS b/Serie06/poisson/solution-sendrecv/AUTHORS new file mode 100644 index 0000000..9253635 --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/AUTHORS @@ -0,0 +1,4 @@ +Authors: + - Nicolas Richart + - Vincent Keller + . Vittoria Rezzonico diff --git a/Serie06/poisson/solution-sendrecv/COPYING b/Serie06/poisson/solution-sendrecv/COPYING new file mode 100644 index 0000000..e551330 --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/COPYING @@ -0,0 +1,24 @@ +Copyright (c) 2017, Ecole Polytechnique Federal de Lausanne - Scientific IT and +Application Support (SCITAS) All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Serie06/poisson/solution-sendrecv/Makefile b/Serie06/poisson/solution-sendrecv/Makefile new file mode 100644 index 0000000..3692c48 --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/Makefile @@ -0,0 +1,14 @@ +CXX=mpicxx +LD=${CXX} +CXXFLAGS+=-Wall -Wextra -Werror -pedantic -std=c++11 -O3 -march=native +LDFLAGS+=-lm $(CXXFLAGS) + +OBJS=poisson.o simulation.o double_buffer.o grid.o dumpers.o + +all: poisson + +poisson: $(OBJS) + $(LD) -o $@ $(OBJS) $(LDFLAGS) + +clean: + rm -f hello poisson *.o *~ diff --git a/Serie06/poisson/solution-sendrecv/double_buffer.cc b/Serie06/poisson/solution-sendrecv/double_buffer.cc new file mode 100644 index 0000000..e6ecf0f --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/double_buffer.cc @@ -0,0 +1,19 @@ +/* -------------------------------------------------------------------------- */ +#include "double_buffer.hh" +#include "grid.hh" +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +DoubleBuffer::DoubleBuffer(int m, int n) + : m_current(new Grid(m, n)), m_old(new Grid(m, n)) {} + +/* -------------------------------------------------------------------------- */ +Grid & DoubleBuffer::current() { return *m_current; } + +/* -------------------------------------------------------------------------- */ +Grid & DoubleBuffer::old() { return *m_old; } + +/* -------------------------------------------------------------------------- */ +void DoubleBuffer::swap() { + m_current.swap(m_old); +} diff --git a/Serie06/poisson/solution-sendrecv/double_buffer.hh b/Serie06/poisson/solution-sendrecv/double_buffer.hh new file mode 100644 index 0000000..67a1405 --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/double_buffer.hh @@ -0,0 +1,24 @@ +#ifndef DOUBLE_BUFFER +#define DOUBLE_BUFFER + +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ +#include "grid.hh" +/* -------------------------------------------------------------------------- */ + +class DoubleBuffer { +public: + DoubleBuffer(int m, int n); + + Grid & current(); + Grid & old(); + + void swap(); +private: + + std::unique_ptr m_current; + std::unique_ptr m_old; +}; + +#endif /* DOUBLE_BUFFER */ diff --git a/Serie06/poisson/solution-sendrecv/dumpers.cc b/Serie06/poisson/solution-sendrecv/dumpers.cc new file mode 100644 index 0000000..429400a --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/dumpers.cc @@ -0,0 +1,145 @@ +/* -------------------------------------------------------------------------- */ +#include "dumpers.hh" +#include "grid.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +void Dumper::set_min(float min) { + m_min = min; +} + +void Dumper::set_max(float max) { + m_max = max; +} + +float Dumper::min() const { + return m_min; +} + +float Dumper::max() const { + return m_max; +} + +/* -------------------------------------------------------------------------- */ +void DumperASCII::dump(int step) { + std::ofstream fout; + std::stringstream sfilename; + + sfilename << "output/out_" << std::setfill('0') << std::setw(5) << step << ".pgm"; + fout.open(sfilename.str()); + + int m = m_grid.m(); + int n = m_grid.n(); + + fout << "P2" << std::endl << "# CREATOR: Poisson program" << std::endl; + fout << m << " " << n << std::endl; + fout << 255 << std::endl; + + for (int i = 0; i < m; i++) { + for (int j = 0; j < n; j++) { + int v = 255. * (m_grid(i, j) - m_min) / (m_max - m_min); + v = std::min(v, 255); + fout << v << std::endl; + } + } +} + +/* -------------------------------------------------------------------------- */ +void DumperBinary::dump(int step) { + std::ofstream fout; + std::stringstream sfilename; + + sfilename << "output/out_" << std::setfill('0') << std::setw(5) << step << ".bmp"; + fout.open(sfilename.str(), std::ios_base::binary); + + int prank, psize; + MPI_Comm_rank(MPI_COMM_WORLD, &prank); + MPI_Comm_size(MPI_COMM_WORLD, &psize); + + int nghosts = ((psize == 1) ? 0 : + (prank == 0 or prank == psize - 1) ? 1 : 2); + int h = m_grid.m() - nghosts; + int w = m_grid.n(); + + // remove ghosts at beginning if rank not 0 + int i_start = (prank == 0 ? 0 : 1); + + std::vector buffer; + if (prank == 0) { + buffer.resize(w * w); + MPI_Gather(&m_grid(i_start, 0), h * w, MPI_FLOAT, + buffer.data(), h*w, MPI_FLOAT, 0, MPI_COMM_WORLD); + } else { + MPI_Gather(&m_grid(i_start, 0), h * w, MPI_FLOAT, + NULL, 0, MPI_FLOAT, 0, MPI_COMM_WORLD); + + // it prank is not stop here + return; + } + + + int total_h = h * psize; + int row_size = 3 * w; + // if the file width (3*w) is not a multiple of 4 adds enough bytes to make it + // a multiple of 4 + int padding = (4 - (row_size) % 4) % 4; + row_size += padding; + + int filesize = 54 + (row_size)*total_h; + + std::vector img(row_size*total_h); + std::fill(img.begin(), img.end(), 0); + + + for (int i = 0; i < total_h; i++) { + for (int j = 0; j < w; j++) { + float v = (buffer[i * w + j] - m_min) / (m_max - m_min); + + float r = v * 255; // Red channel + float g = v * 255; // Green channel + float b = v * 255; // Red channel + + r = std::min(r, 255.f); + g = std::min(g, 255.f); + b = std::min(b, 255.f); + + img[row_size * i + 3 * j + 2] = r; + img[row_size * i + 3 * j + 1] = g; + img[row_size * i + 3 * j + 0] = b; + } + } + + std::array bmpfileheader = {'B', 'M', 0, 0, 0, 0, 0, + 0, 0, 0, 54, 0, 0, 0}; + std::array bmpinfoheader = {40, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 1, 0, 24, 0}; + + bmpfileheader[2] = filesize; + bmpfileheader[3] = filesize >> 8; + bmpfileheader[4] = filesize >> 16; + bmpfileheader[5] = filesize >> 24; + + bmpinfoheader[4] = w; + bmpinfoheader[5] = w >> 8; + bmpinfoheader[6] = w >> 16; + bmpinfoheader[7] = w >> 24; + bmpinfoheader[8] = total_h; + bmpinfoheader[9] = total_h >> 8; + bmpinfoheader[10] = total_h >> 16; + bmpinfoheader[11] = total_h >> 24; + bmpinfoheader[20] = (filesize - 54); + bmpinfoheader[21] = (filesize - 54) >> 8; + bmpinfoheader[22] = (filesize - 54) >> 16; + bmpinfoheader[23] = (filesize - 54) >> 24; + + fout.write(bmpfileheader.data(), 14); + fout.write(bmpinfoheader.data(), 40); + + fout.write(img.data(), total_h * row_size); +} diff --git a/Serie06/poisson/solution-sendrecv/dumpers.hh b/Serie06/poisson/solution-sendrecv/dumpers.hh new file mode 100644 index 0000000..03a78c5 --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/dumpers.hh @@ -0,0 +1,41 @@ +#ifndef DUMPERS_HH +#define DUMPERS_HH + +/* -------------------------------------------------------------------------- */ +class Grid; + +/* -------------------------------------------------------------------------- */ +class Dumper { +public: + explicit Dumper(const Grid & grid) : m_grid(grid), m_min(-1.), m_max(1.) {} + + virtual void dump(int step) = 0; + + void set_min(float min); + void set_max(float max); + + float min() const; + float max() const; + +protected: + const Grid & m_grid; + float m_min, m_max; +}; + +/* -------------------------------------------------------------------------- */ +class DumperASCII : public Dumper { +public: + explicit DumperASCII(const Grid & grid) : Dumper(grid) {} + + virtual void dump(int step); +}; + +/* -------------------------------------------------------------------------- */ +class DumperBinary : public Dumper { +public: + explicit DumperBinary(const Grid & grid) : Dumper(grid) {} + + virtual void dump(int step); +}; + +#endif /* DUMPERS_HH */ diff --git a/Serie06/poisson/solution-sendrecv/efficiency.sh b/Serie06/poisson/solution-sendrecv/efficiency.sh new file mode 100644 index 0000000..f3a52e1 --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/efficiency.sh @@ -0,0 +1,14 @@ +#!/usr/bin/bash -l +#SBATCH -N 2 +#SBATCH -n 56 +#SBATCH -C E5v4 + +module load gcc mvapich2 + +np=(1 2 4 7 14 28 42 56) +size=(1024 1448 2048 2709 3836 5404 6636 7672) + +for i in 0 1 2 3 4 5 6 7 +do + srun -n ${np[$i]} ./poisson $(( 4 * ${size[$i]} )) >> eff +done diff --git a/Serie06/poisson/solution-sendrecv/grid.cc b/Serie06/poisson/solution-sendrecv/grid.cc new file mode 100644 index 0000000..4e9fb84 --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/grid.cc @@ -0,0 +1,15 @@ +/* -------------------------------------------------------------------------- */ +#include "grid.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +Grid::Grid(int m, int n) : m_m(m), m_n(n), m_storage(m * n) { clear(); } + +/* -------------------------------------------------------------------------- */ +void Grid::clear() { std::fill(m_storage.begin(), m_storage.end(), 0.); } + +/* -------------------------------------------------------------------------- */ +int Grid::m() const { return m_m; } +int Grid::n() const { return m_n; } diff --git a/Serie06/poisson/solution-sendrecv/grid.hh b/Serie06/poisson/solution-sendrecv/grid.hh new file mode 100644 index 0000000..ad872cb --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/grid.hh @@ -0,0 +1,29 @@ +#ifndef GRID_HH +#define GRID_HH + +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +class Grid { +public: + Grid(int m, int n); + + /// access the value [i][j] of the grid + inline float & operator()(int i, int j) { return m_storage[i * m_n + j]; } + inline const float & operator()(int i, int j) const { + return m_storage[i * m_n + j]; + } + + /// set the grid to 0 + void clear(); + + int m() const; + int n() const; + +private: + int m_m, m_n; + std::vector m_storage; +}; + +#endif /* GRID_HH */ diff --git a/Serie06/poisson/solution-sendrecv/output/.to-keep b/Serie06/poisson/solution-sendrecv/output/.to-keep new file mode 100644 index 0000000..e69de29 diff --git a/Serie06/poisson/solution-sendrecv/poisson b/Serie06/poisson/solution-sendrecv/poisson new file mode 100755 index 0000000..c13d331 Binary files /dev/null and b/Serie06/poisson/solution-sendrecv/poisson differ diff --git a/Serie06/poisson/solution-sendrecv/poisson.cc b/Serie06/poisson/solution-sendrecv/poisson.cc new file mode 100644 index 0000000..692b50e --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/poisson.cc @@ -0,0 +1,70 @@ +/* -------------------------------------------------------------------------- */ +#include "simulation.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +#define EPSILON 0.005 + +typedef std::chrono::high_resolution_clock clk; +typedef std::chrono::duration second; + +static void usage(const std::string & prog_name) { + std::cerr << prog_name << " " << std::endl; + exit(0); +} + +int main(int argc, char * argv[]) { + MPI_Init(&argc, &argv); + int prank, psize; + + MPI_Comm_rank(MPI_COMM_WORLD, &prank); + MPI_Comm_size(MPI_COMM_WORLD, &psize); + + if (argc != 2) usage(argv[0]); + + std::stringstream args(argv[1]); + int N; + args >> N; + + if(args.fail()) usage(argv[0]); + + + // Divide the total size by the number of processes + int m_local = N / psize; + + // Add 1 or 2 ghosts depending on the prank + m_local += + (psize == 1) ? 0 : + (prank == 0 or prank == psize -1) ? 1 : 2; + + Simulation simu(m_local, N); + + simu.set_initial_conditions(); + + simu.set_epsilon(EPSILON); + + float l2; + int k; + + auto start = clk::now(); + std::tie(l2, k) = simu.compute(); + auto end = clk::now(); + + second time = end - start; + + if(prank == 0) + std::cout << psize << " " << N << " " + << k << " " << std::scientific << l2 << " " + << time.count() / k << std::endl; + + MPI_Finalize(); + + return 0; +} diff --git a/Serie06/poisson/solution-sendrecv/scalability.sh b/Serie06/poisson/solution-sendrecv/scalability.sh new file mode 100644 index 0000000..f9ae444 --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/scalability.sh @@ -0,0 +1,12 @@ +#!/usr/bin/bash -l +#SBATCH -N 2 +#SBATCH -n 56 +#SBATCH -C E5v4 + +for size in 1024 2048 4096 8192 16384 +do + for i in 1 2 4 7 14 28 42 56 + do + srun -n $i ./poisson $size >> sendrecv_$size + done +done diff --git a/Serie06/poisson/solution-sendrecv/simulation.cc b/Serie06/poisson/solution-sendrecv/simulation.cc new file mode 100644 index 0000000..86b5265 --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/simulation.cc @@ -0,0 +1,111 @@ +/* -------------------------------------------------------------------------- */ +#include "simulation.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +Simulation::Simulation(int m, int n) + // Compute h_m and h_n with n since m is different on all processors + : m_local_m(m), m_local_n(n), m_epsilon(1e-7), m_h_m(1. / n), + m_h_n(1. / n), m_grids(m, n), m_f(m, n), + m_dumper(new DumperBinary(m_grids.old())) { + + // retrieving the number of proc and the rank in the proc pool + MPI_Comm_rank(MPI_COMM_WORLD, &m_prank); + MPI_Comm_size(MPI_COMM_WORLD, &m_psize); +} + +/* -------------------------------------------------------------------------- */ +void Simulation::set_initial_conditions() { + // remove ghosts at beginning if rank not 0 + int i_start = (m_prank == 0 ? 0 : 1); + +// remove ghost at the end if prank is not the last + int i_end = (m_prank == m_psize - 1 ? m_local_m : m_local_m - 1); + + int nghosts = + (m_psize == 1) ? 0 : + (m_prank == 0 or m_prank == m_psize - 1 ? 1 : 2); + // computing the offsets of the local grid in the global one + int offset_m = + (m_local_m - nghosts) * m_prank; + int offset_n = 0; + + for (int i = i_start; i < i_end; i++) { + for (int j = 0; j < m_local_n; j++) { + m_f(i, j) = -2. * 100. * M_PI * M_PI * + std::sin(10. * M_PI * (offset_m + i - i_start) * m_h_m) * + std::sin(10. * M_PI * (offset_n + j) * m_h_n); + } + } +} + +/* -------------------------------------------------------------------------- */ +std::tuple Simulation::compute() { + int s = 0; + float l2 = 0; + //do { + for(s = 0; s < 100; ++s) { + l2 = compute_step(); + + m_grids.swap(); + + // ++s; + //} while (l2 > m_epsilon); + } + + //m_dumper->dump(s); + + return std::make_tuple(l2, s); +} + +/* -------------------------------------------------------------------------- */ +void Simulation::set_epsilon(float epsilon) { m_epsilon = epsilon; } + +/* -------------------------------------------------------------------------- */ +float Simulation::epsilon() const { return m_epsilon; } + +/* -------------------------------------------------------------------------- */ +float Simulation::compute_step() { + float l2 = 0.; + + Grid & u = m_grids.current(); + Grid & uo = m_grids.old(); + + // determining the rank of the neighbors + int north_prank = (m_prank == 0 ? MPI_PROC_NULL : m_prank - 1); + int south_prank = (m_prank == (m_psize - 1) ? MPI_PROC_NULL : m_prank + 1); + + // MPI_PROC_NULL allows you to do a send en receive that does nothing + + // Taking care of communications going up (so send up, receiving from bottom) + MPI_Sendrecv(&uo(1, 0), m_local_n, MPI_FLOAT, north_prank, 0, + &uo(m_local_m - 1, 0), m_local_n, MPI_FLOAT, south_prank, 0, + MPI_COMM_WORLD, MPI_STATUS_IGNORE); + + // Taking care of communications going down (so sending down, + // receiving from top) + MPI_Sendrecv(&uo(m_local_m - 2, 0), m_local_n, MPI_FLOAT, south_prank, 0, + &uo(0, 0), m_local_n, MPI_FLOAT, north_prank, 0, + MPI_COMM_WORLD, MPI_STATUS_IGNORE); + + + for (int i = 1; i < m_local_m - 1; i++) { + for (int j = 1; j < m_local_n - 1; j++) { + // computation of the new step + u(i, j) = 0.25 * (uo(i - 1, j) + uo(i + 1, j) + uo(i, j - 1) + + uo(i, j + 1) - m_f(i, j) * m_h_m * m_h_n); + + // L2 norm + l2 += (uo(i, j) - u(i, j)) * (uo(i, j) - u(i, j)); + } + } + + // Summing the value of all the processors together + MPI_Allreduce(MPI_IN_PLACE, &l2, 1, MPI_FLOAT, MPI_SUM, MPI_COMM_WORLD); + + return l2; +} diff --git a/Serie06/poisson/solution-sendrecv/simulation.hh b/Serie06/poisson/solution-sendrecv/simulation.hh new file mode 100644 index 0000000..be1879e --- /dev/null +++ b/Serie06/poisson/solution-sendrecv/simulation.hh @@ -0,0 +1,54 @@ +#ifndef SIMULATION_HH +#define SIMULATION_HH + +/* -------------------------------------------------------------------------- */ +#include "double_buffer.hh" +#include "dumpers.hh" +#include "grid.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +class Simulation { +public: + Simulation(int m, int n); + + /// set the initial conditions, Dirichlet and source term + virtual void set_initial_conditions(); + + /// perform the simulation + std::tuple compute(); + + /// access the precision + void set_epsilon(float epsilon); + float epsilon() const; + +protected: + /// compute one step and return an error + virtual float compute_step(); + +private: + /// Global problem size + int m_local_m, m_local_n; + + /// Precision to achieve + float m_epsilon; + + /// grid spacing + float m_h_m; + float m_h_n; + + /// Grids storage + DoubleBuffer m_grids; + + /// source term + Grid m_f; + + int m_prank{0}, m_psize{1}; + + /// Dumper to use for outputs + std::unique_ptr m_dumper; +}; + +#endif /* SIMULATION_HH */ diff --git a/Serie07/pi/solution/Makefile b/Serie07/pi/solution/Makefile new file mode 100644 index 0000000..30f109f --- /dev/null +++ b/Serie07/pi/solution/Makefile @@ -0,0 +1,14 @@ +OPTIM+=-O3 + +CXX=mpicxx +CC=mpicxx +LD=${CXX} +CXXFLAGS+=-Wall -Wextra -std=c++11 $(OPTIM) +LDFLAGS+=$(OPTIM) -lm + +EXECS=pi_derived_datatype pi_p2p_persistent_ring + +all: clean $(EXECS) + +clean: + rm -f $(EXECS) *.o *~ diff --git a/Serie07/pi/solution/pi_derived_datatype.cc b/Serie07/pi/solution/pi_derived_datatype.cc new file mode 100644 index 0000000..d744f6a --- /dev/null +++ b/Serie07/pi/solution/pi_derived_datatype.cc @@ -0,0 +1,111 @@ +/* + This exercise is taken from the class Parallel Programming Workshop (MPI, + OpenMP and Advanced Topics) at HLRS given by Rolf Rabenseifner + */ + +#include +#include +#include +#include +#include + +using clk = std::chrono::high_resolution_clock; +using second = std::chrono::duration; +using time_point = std::chrono::time_point; + +inline int digit(double x, int n) { + return std::trunc(x * std::pow(10., n)) - std::trunc(x * std::pow(10., n - 1)) *10.; +} + +inline double f(double a) { return (4. / (1. + a * a)); } + +const int n = 10000000; + +int main(int /* argc */ , char ** /* argv */) { + int i; + double dx, x, l_sum, pi; + int psize, prank; + + MPI_Init(NULL, NULL); + + MPI_Comm_size(MPI_COMM_WORLD, &psize); + MPI_Comm_rank(MPI_COMM_WORLD, &prank); + + struct Sum { + double sum; + int rank; + }; + + Sum sum{0., 0}; + + int blk_length[2] = {1, 1}; + + MPI_Aint zero_address, first_address, second_address; + MPI_Get_address(&sum, &zero_address); + MPI_Get_address(&sum.sum, &first_address); + MPI_Get_address(&sum.rank, &second_address); + + MPI_Aint displs[2]; + displs[0] = MPI_Aint_diff(first_address, zero_address);; + displs[1] = MPI_Aint_diff(second_address, first_address); + + MPI_Datatype types[2] = {MPI_DOUBLE, MPI_INT}; + MPI_Datatype sum_t; + MPI_Type_create_struct(2, blk_length, displs, types, &sum_t); + MPI_Type_commit(&sum_t); + + auto mpi_t1 = MPI_Wtime(); + auto t1 = clk::now(); + + int nlocal = n / psize; + int istart = 1 + nlocal * prank; + int iend = nlocal * (prank + 1); + + /* calculate pi = integral [0..1] 4 / (1 + x**2) dx */ + dx = 1. / n; + l_sum = 0.0; + for (i = istart; i <= iend; i++) { + x = (1. * i - 0.5) * dx; + l_sum = l_sum + f(x); + } + + sum.sum = l_sum; + sum.rank = prank; + + if(prank == 0) { + std::vector sums(psize); + MPI_Gather(&sum, 1, sum_t, sums.data(), 1, sum_t, 0, MPI_COMM_WORLD); + + l_sum = 0.; + for(auto s : sums) { + l_sum += s.sum; + } + + sum.sum = l_sum; + } else { + MPI_Gather(&sum, 1, sum_t, NULL, 0, sum_t, 0, MPI_COMM_WORLD); + } + + + MPI_Bcast(&sum, 1, sum_t, 0, MPI_COMM_WORLD); + + pi = dx * sum.sum; + + auto mpi_elapsed = MPI_Wtime() - mpi_t1; + second elapsed = clk::now() - t1; + + if(prank == 0) { + std::printf("computed pi = %.16g\n", pi); + std::printf("wall clock time (mpi_wtime) = %.4gs with %d process\n", mpi_elapsed, psize); + std::printf("wall clock time (chrono) = %.4gs\n", elapsed.count()); + + for(int d = 1; d <= 15; ++d) { + std::printf("%d", digit(pi, d)); + } + } + + MPI_Type_free(&sum_t); + MPI_Finalize(); + + return 0; +} diff --git a/Serie07/pi/solution/pi_p2p_persistent_ring.cc b/Serie07/pi/solution/pi_p2p_persistent_ring.cc new file mode 100644 index 0000000..a2c6742 --- /dev/null +++ b/Serie07/pi/solution/pi_p2p_persistent_ring.cc @@ -0,0 +1,93 @@ +/* + This exercise is taken from the class Parallel Programming Workshop (MPI, + OpenMP and Advanced Topics) at HLRS given by Rolf Rabenseifner + */ + +#include +#include +#include +#include + +using clk = std::chrono::high_resolution_clock; +using second = std::chrono::duration; +using time_point = std::chrono::time_point; + +inline int digit(double x, int n) { + return std::trunc(x * std::pow(10., n)) - std::trunc(x * std::pow(10., n - 1)) *10.; +} + +inline double f(double a) { return (4. / (1. + a * a)); } + +const int n = 10000000; + +int main(int /* argc */ , char ** /* argv */) { + int i; + double dx, x, sum, pi; + int psize, prank; + + MPI_Init(NULL, NULL); + + MPI_Comm_size(MPI_COMM_WORLD, &psize); + MPI_Comm_rank(MPI_COMM_WORLD, &prank); + + auto mpi_t1 = MPI_Wtime(); + auto t1 = clk::now(); + + int nlocal = n / psize; + int istart = 1 + nlocal * prank; + int iend = nlocal * (prank + 1); + + /* calculate pi = integral [0..1] 4 / (1 + x**2) dx */ + dx = 1. / n; + sum = 0.0; + for (i = istart; i <= iend; i++) { + x = (1. * i - 0.5) * dx; + sum = sum + f(x); + } + + int next = (prank + 1) % psize; + int prev = (prank - 1 + psize) % psize; + + double send, recv; + MPI_Request request[2]; + + send = sum; + MPI_Send_init(&send, 1, MPI_DOUBLE, next, 13, MPI_COMM_WORLD, request); + MPI_Recv_init(&recv, 1, MPI_DOUBLE, prev, 13, MPI_COMM_WORLD, request + 1); + + for(int s = 1; s < psize; ++s) { + MPI_Startall(2, request); + + // ensure that receive is finished before using recv + MPI_Wait(request + 1, MPI_STATUS_IGNORE); + + sum += recv; + + // ensure that send is finished + MPI_Wait(request, MPI_STATUS_IGNORE); + + send = recv; + } + + pi = dx * sum; + + MPI_Request_free(request); + MPI_Request_free(request + 1); + + auto mpi_elapsed = MPI_Wtime() - mpi_t1; + second elapsed = clk::now() - t1; + + if(prank == 0) { + std::printf("computed pi = %.16g\n", pi); + std::printf("wall clock time (mpi_wtime) = %.4gs with %d process\n", mpi_elapsed, psize); + std::printf("wall clock time (chrono) = %.4gs\n", elapsed.count()); + + for(int d = 1; d <= 15; ++d) { + std::printf("%d", digit(pi, d)); + } + } + + MPI_Finalize(); + + return 0; +} diff --git a/Serie07/write_bmp/solution/AUTHORS b/Serie07/write_bmp/solution/AUTHORS new file mode 100644 index 0000000..9253635 --- /dev/null +++ b/Serie07/write_bmp/solution/AUTHORS @@ -0,0 +1,4 @@ +Authors: + - Nicolas Richart + - Vincent Keller + . Vittoria Rezzonico diff --git a/Serie07/write_bmp/solution/COPYING b/Serie07/write_bmp/solution/COPYING new file mode 100644 index 0000000..e551330 --- /dev/null +++ b/Serie07/write_bmp/solution/COPYING @@ -0,0 +1,24 @@ +Copyright (c) 2017, Ecole Polytechnique Federal de Lausanne - Scientific IT and +Application Support (SCITAS) All rights reserved. + +Redistribution and use in source and binary forms, with or without +modification, are permitted provided that the following conditions are met: + * Redistributions of source code must retain the above copyright + notice, this list of conditions and the following disclaimer. + * Redistributions in binary form must reproduce the above copyright + notice, this list of conditions and the following disclaimer in the + documentation and/or other materials provided with the distribution. + * Neither the name of the nor the + names of its contributors may be used to endorse or promote products + derived from this software without specific prior written permission. + +THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND +ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED +WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE +DISCLAIMED. IN NO EVENT SHALL BE LIABLE FOR ANY +DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES +(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; +LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND +ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT +(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS +SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. diff --git a/Serie07/write_bmp/solution/Makefile b/Serie07/write_bmp/solution/Makefile new file mode 100644 index 0000000..6f64d80 --- /dev/null +++ b/Serie07/write_bmp/solution/Makefile @@ -0,0 +1,14 @@ +CXX=mpic++ +LD=${CXX} +CXXFLAGS+=-Wall -Wextra -pedantic -std=c++11 -O3 +LDFLAGS+=-lm $(CXXFLAGS) + +OBJS=grid.o dumpers.o main.o + +all: write_bmp + +write_bmp: $(OBJS) + $(LD) -o $@ $(OBJS) $(LDFLAGS) + +clean: + rm -f write_bmp *.o *~ diff --git a/Serie07/write_bmp/solution/dumpers.cc b/Serie07/write_bmp/solution/dumpers.cc new file mode 100644 index 0000000..d9a034f --- /dev/null +++ b/Serie07/write_bmp/solution/dumpers.cc @@ -0,0 +1,142 @@ +/* -------------------------------------------------------------------------- */ +#include "dumpers.hh" +#include "grid.hh" +/* -------------------------------------------------------------------------- */ +#include +#include +#include +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +void Dumper::set_min(float min) { + m_min = min; +} + +void Dumper::set_max(float max) { + m_max = max; +} + +float Dumper::min() const { + return m_min; +} + +float Dumper::max() const { + return m_max; +} + +/* -------------------------------------------------------------------------- */ +void DumperASCII::dump(int step) { + std::ofstream fout; + std::stringstream sfilename; + + sfilename << "output/out_" << std::setfill('0') << std::setw(5) << step << ".pgm"; + fout.open(sfilename.str()); + + int m = m_grid.m(); + int n = m_grid.n(); + + fout << "P2" << std::endl << "# CREATOR: Poisson program" << std::endl; + fout << m << " " << n << std::endl; + fout << 255 << std::endl; + + for (int i = 0; i < m; i++) { + for (int j = 0; j < n; j++) { + int v = 255. * (m_grid(i, j) - m_min) / (m_max - m_min); + v = std::min(v, 255); + fout << v << std::endl; + } + } +} + +/* -------------------------------------------------------------------------- */ +void DumperBinary::dump(int step) { + std::ofstream fout; + std::stringstream sfilename; + + sfilename << "out_" << std::setfill('0') << std::setw(5) << step << ".bmp"; + fout.open(sfilename.str(), std::ios_base::binary); + + int h = m_grid.m(); + int w = m_grid.n(); + + int row_size = 3 * w; + // if the file width (3*w) is not a multiple of 4 adds enough bytes to make it + // a multiple of 4 + int padding = (4 - (row_size) % 4) % 4; + row_size += padding; + + int psize; + MPI_Comm_size(MPI_COMM_WORLD, &psize); + + int filesize = 54 + (row_size)* h * psize; + + std::vector img(row_size*h); + std::fill(img.begin(), img.end(), 0); + + + for (int i = 0; i < h; i++) { + for (int j = 0; j < w; j++) { + float v = ((m_grid(i, j) - m_min) / (m_max - m_min)); + + float r = v * 255; // Red channel + float g = v * 255; // Green channel + float b = v * 255; // Red channel + + r = std::min(r, 255.f); + g = std::min(g, 255.f); + b = std::min(b, 255.f); + + img[row_size * i + 3 * j + 2] = r; + img[row_size * i + 3 * j + 1] = g; + img[row_size * i + 3 * j + 0] = b; + } + } + + std::array bmpfileheader = {'B', 'M', 0, 0, 0, 0, 0, + 0, 0, 0, 54, 0, 0, 0}; + std::array bmpinfoheader = {40, 0, 0, 0, 0, 0, 0, 0, + 0, 0, 0, 0, 1, 0, 24, 0}; + + bmpfileheader[2] = filesize; + bmpfileheader[3] = filesize >> 8; + bmpfileheader[4] = filesize >> 16; + bmpfileheader[5] = filesize >> 24; + + bmpinfoheader[4] = w; + bmpinfoheader[5] = w >> 8; + bmpinfoheader[6] = w >> 16; + bmpinfoheader[7] = w >> 24; + bmpinfoheader[8] = (h * psize); + bmpinfoheader[9] = (h * psize) >> 8; + bmpinfoheader[10] = (h * psize) >> 16; + bmpinfoheader[11] = (h * psize) >> 24; + bmpinfoheader[20] = (filesize - 54); + bmpinfoheader[21] = (filesize - 54) >> 8; + bmpinfoheader[22] = (filesize - 54) >> 16; + bmpinfoheader[23] = (filesize - 54) >> 24; + + MPI_File fh; + MPI_Status status; + + // opening a file in write and create mode + MPI_File_open(MPI_COMM_WORLD, sfilename.str().c_str(), + MPI_MODE_WRONLY | MPI_MODE_CREATE, + MPI_INFO_NULL, &fh); + // defining the size of the file + MPI_File_set_size(fh, filesize); + + int prank; + MPI_Comm_rank(MPI_COMM_WORLD, &prank); + + // rank 0 writes the header + if (prank == 0) { + MPI_File_write_at(fh, 0, bmpfileheader.data(), 14, MPI_CHAR, &status); + MPI_File_write_at(fh, 14, bmpinfoheader.data(), 40, MPI_CHAR, &status); + } + + int offset = 54 + row_size * h * prank; + + MPI_File_write_at(fh, offset, img.data(), img.size(), MPI_CHAR, &status); + MPI_File_close(&fh); +} diff --git a/Serie07/write_bmp/solution/dumpers.hh b/Serie07/write_bmp/solution/dumpers.hh new file mode 100644 index 0000000..27a5f12 --- /dev/null +++ b/Serie07/write_bmp/solution/dumpers.hh @@ -0,0 +1,48 @@ +#include + +#ifndef DUMPERS_HH +#define DUMPERS_HH + +/* -------------------------------------------------------------------------- */ +class Grid; + +/* -------------------------------------------------------------------------- */ +class Dumper { +public: + explicit Dumper(const Grid & grid, MPI_Comm communicator) + : m_grid(grid), m_min(-1.), m_max(1.), m_communicator(communicator) {} + + virtual void dump(int step) = 0; + + void set_min(float min); + void set_max(float max); + + float min() const; + float max() const; + +protected: + const Grid & m_grid; + float m_min, m_max; + + MPI_Comm m_communicator; +}; + +/* -------------------------------------------------------------------------- */ +class DumperASCII : public Dumper { +public: + explicit DumperASCII(const Grid & grid, MPI_Comm communicator) + : Dumper(grid, communicator) {} + + virtual void dump(int step); +}; + +/* -------------------------------------------------------------------------- */ +class DumperBinary : public Dumper { +public: + explicit DumperBinary(const Grid & grid, MPI_Comm communicator) + : Dumper(grid, communicator) {} + + virtual void dump(int step); +}; + +#endif /* DUMPERS_HH */ diff --git a/Serie07/write_bmp/solution/grid.cc b/Serie07/write_bmp/solution/grid.cc new file mode 100644 index 0000000..bd635e9 --- /dev/null +++ b/Serie07/write_bmp/solution/grid.cc @@ -0,0 +1,24 @@ +/* -------------------------------------------------------------------------- */ +#include "grid.hh" +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +Grid::Grid(int m, int n) : m_m(m), m_n(n), m_storage(m * n) { clear(); } + +/* -------------------------------------------------------------------------- */ +void Grid::clear() { std::fill(m_storage.begin(), m_storage.end(), 0.); } + +/* -------------------------------------------------------------------------- */ +void Grid::resize(int m, int n) { + m_m = m; + m_n = n; + m_storage.resize(m*n); +} + +/* -------------------------------------------------------------------------- */ +int Grid::m() const { return m_m; } +int Grid::n() const { return m_n; } + +/* -------------------------------------------------------------------------- */ diff --git a/Serie07/write_bmp/solution/grid.hh b/Serie07/write_bmp/solution/grid.hh new file mode 100644 index 0000000..5f0c367 --- /dev/null +++ b/Serie07/write_bmp/solution/grid.hh @@ -0,0 +1,31 @@ +#ifndef GRID_HH +#define GRID_HH + +/* -------------------------------------------------------------------------- */ +#include +/* -------------------------------------------------------------------------- */ + +class Grid { +public: + Grid(int m = 0, int n = 0); + + /// access the value [i][j] of the grid + inline float & operator()(int i, int j) { return m_storage[i * m_n + j]; } + inline const float & operator()(int i, int j) const { + return m_storage[i * m_n + j]; + } + + void resize(int m, int n); + + /// set the grid to 0 + void clear(); + + int m() const; + int n() const; + +private: + int m_m, m_n; + std::vector m_storage; +}; + +#endif /* GRID_HH */ diff --git a/Serie07/write_bmp/solution/main.cc b/Serie07/write_bmp/solution/main.cc new file mode 100644 index 0000000..8acd32b --- /dev/null +++ b/Serie07/write_bmp/solution/main.cc @@ -0,0 +1,53 @@ +/* -------------------------------------------------------------------------- */ +#include +#include +#include +/* -------------------------------------------------------------------------- */ +#include "grid.hh" +#include "dumpers.hh" +/* -------------------------------------------------------------------------- */ + +/* -------------------------------------------------------------------------- */ +static void usage(const std::string & prog_name) { + std::cerr << prog_name << " " << std::endl; + exit(0); +} + +/* -------------------------------------------------------------------------- */ +int main(int argc, char *argv[]) { + MPI_Init(&argc, &argv); + int prank, psize; + + MPI_Comm_rank(MPI_COMM_WORLD, &prank); + MPI_Comm_size(MPI_COMM_WORLD, &psize); + + if (argc != 2) usage(argv[0]); + + std::stringstream args(argv[1]); + int n; + args >> n; + + if(args.fail()) usage(argv[0]); + + int m_local = n / psize; + + Grid f(m_local, n); + + float h = 1./n; + int offset_m = m_local * prank; + int offset_n = 0; + + for (int i = 0; i < m_local; i++) { + for (int j = 0; j < n; j++) { + f(i, j) = std::sin(10. * M_PI * (offset_m + i) * h) * + std::cos(10. * M_PI * (offset_n + j) * h); + } + } + + DumperBinary dumper(f, MPI_COMM_WORLD); + dumper.dump(0); + + MPI_Finalize(); + +return 0; +}