From e8e09b0637622ddac77ab9e3fe02c9547b5a5184 Mon Sep 17 00:00:00 2001 From: daniel Date: Sun, 13 Mar 2022 19:07:39 +0100 Subject: [PATCH] add: utils, wip: ex_01 (almost done) --- Exercise_01/main.cpp | 219 ++++++++++++++++++++----------------------- Exercise_01/utils.h | 105 +++++++++++++++++++++ 2 files changed, 206 insertions(+), 118 deletions(-) create mode 100644 Exercise_01/utils.h diff --git a/Exercise_01/main.cpp b/Exercise_01/main.cpp index b56025a..a7062ac 100644 --- a/Exercise_01/main.cpp +++ b/Exercise_01/main.cpp @@ -7,17 +7,13 @@ #include #include #include +#include #define _USE_MATH_DEFINES /* enables math constants from cmath */ #include #ifdef USE_MPI #include #endif - -#include -#include - - #include "Matrix.h" #include "Solver.h" #include "utils.h" @@ -32,57 +28,69 @@ int main(int argn, char* argv[]) { // Get MPI config int mpi_size; /*< MPI pool size (a.k.a. total number of processes) */ MPI_Comm_size(MPI_COMM_WORLD, &mpi_size); + + // Get MPI rank + int mpi_world_rank; /*< MPI world rank (a.k.a. grid process ID) */ + MPI_Comm_rank(MPI_COMM_WORLD, &mpi_world_rank); #else int mpi_size = 1; + int mpi_world_rank = 0; #endif /**************************** Parse Arguments *****************************/ - size_t dim, resolution, iterations; - { - if (0 < argn && argn < 4) { - std::cerr << "usage: " << argv[0] << " " << std::endl; - return -1; - } else if (argn > 4) { - std::cerr << "warning: " << "ignoring all but the first three params" << std::endl; - } - if (std::string(argv[1]) == "1D") { - dim = 1; - } else if (std::string(argv[1]) == "2D") { - dim = 2; - } else { - std::cerr << "error: Parsing arg 1 failed" << std::endl; - return -1; - } - if (!(std::istringstream(argv[2]) >> resolution)) { - std::cerr << "error: Parsing arg 2 failed" << std::endl; - return -1; - } - if (!(std::istringstream(argv[3]) >> iterations)) { - std::cerr << "error: Parsing arg 3 failed" << std::endl; - return -1; - } - if (resolution < 10 || resolution > 65536) { - std::cerr << "error: arg 1 " << resolution - << " out of domain [10, 65536]" << std::endl; - return -1; - } - if (iterations > 65536) { - std::cerr << "error: arg 2 " << iterations - << "out of domain [0, 65536]" << std::endl; - return -1; - } + if (0 < argn && argn < 4) { + std::cerr << "usage: " << argv[0] << " " << std::endl; + return -1; + } else if (argn > 4) { + std::cerr << "warning: " << "ignoring all but the first three params" << std::endl; + } +#ifdef USE_MPI + size_t dim; + if (std::string(argv[1]) == "1D") { + dim = 1; + } else if (std::string(argv[1]) == "2D") { + dim = 2; + } else { + std::cerr << "error: Parsing arg 1 failed" << std::endl; + return -1; + } +#endif + size_t resolution; + if (!(std::istringstream(argv[2]) >> resolution)) { + std::cerr << "error: Parsing arg 2 failed" << std::endl; + return -1; + } + size_t iterations; + if (!(std::istringstream(argv[3]) >> iterations)) { + std::cerr << "error: Parsing arg 3 failed" << std::endl; + return -1; + } + if (resolution < 10 || resolution > 65536) { + std::cerr << "error: arg 1 " << resolution + << " out of domain [10, 65536]" << std::endl; + return -1; + } + if (iterations > 65536) { + std::cerr << "error: arg 2 " << iterations + << "out of domain [0, 65536]" << std::endl; + return -1; } // Report configuration - std::cout << std::setfill(' ') + if (mpi_world_rank == 0) { + std::cout << std::setfill(' ') #ifdef USE_MPI - << "using MPI: true\n" + << "use MPI: YES\n" #else - << "using MPI: false\n" + << "use MPI: NO\n" #endif - << "nr. processes: " << std::setw(5) << mpi_size << '\n' - << "resolution: " << std::setw(5) << resolution << '\n' - << "iterations: " << std::setw(5) << iterations << std::endl; + << "nr. processes: " << std::setw(5) << mpi_size << '\n' + << "resolution: " << std::setw(5) << resolution << '\n' + << "iterations: " << std::setw(5) << iterations << std::endl; + } + + /************************* Start Time Measurement *************************/ + auto start = std::chrono::high_resolution_clock::now(); /**************** Setup (local) PDE + Boundary Conditions *****************/ const double k = M_PI; @@ -119,10 +127,6 @@ int main(int argn, char* argv[]) { &mpi_comm_grid ); - // Get MPI rank - int mpi_world_rank; /*< MPI world rank (a.k.a. grid process ID) */ - MPI_Comm_rank(MPI_COMM_WORLD, &mpi_world_rank); - // Get MPI rank with respect to the grid communicator int mpi_grid_rank; /*< MPI grid rank (a.k.a. grid process ID) */ MPI_Comm_rank(mpi_comm_grid, &mpi_grid_rank); @@ -169,28 +173,6 @@ int main(int argn, char* argv[]) { double ymax = (mpi_neighbours.north == MPI_PROC_NULL) ? 1.0 : h * (partition_sum(resolution, mpi_dims[1], mpi_coords[1]) + 1); - // /*************************__________DEBUG__________*************************/ - // std::cout << " (" << mpi_coords[0] << ", " << mpi_coords[1] << ") "; - - // std::cout << " -> " - // << "N: " << ((mpi_neighbours.north != MPI_PROC_NULL) ? '1' : '.') << ", " - // << "E: " << ((mpi_neighbours.east != MPI_PROC_NULL) ? '1' : '.') << ", " - // << "S: " << ((mpi_neighbours.south != MPI_PROC_NULL) ? '1' : '.') << ", " - // << "W: " << ((mpi_neighbours.west != MPI_PROC_NULL) ? '1' : '.'); - - // std::cout << " -> " - // << nx << " x " << ny; - - // std::cout << " -> " << std::setprecision(2) - // << "[" << std::setw(4) << std::left << xmin - // << ", " << std::setw(4) << std::left << xmax << "]" - // << " x " - // << "[" << std::setw(4) << std::left << ymin - // << ", " << std::setw(4) << std::left << ymax << "]"; - - // std::cout << std::endl; - // /*************************__________DEBUG__________*************************/ - // Create MPI vector Type (for boundary condition exchange) // Allows directly exchange of matrix rows (north/south bounds) since the // row elements are "sparce" in the sence that they are not directly aside @@ -296,11 +278,6 @@ int main(int argn, char* argv[]) { } /***************************** Solution Stats *****************************/ - std::cout - << "|| residuals ||_2: " << std::setw(15) << solver.resid_norm() - << "\n|| residuals ||_inf: " << std::setw(15) << solver.resid_norm() - << std::endl; - // Get solution Matrix& solution = solver.solution(); @@ -314,51 +291,57 @@ int main(int argn, char* argv[]) { } } - // Calculate error - std::cout - << "|| error ||_2: " << std::setw(15) << solution.dist(analytic) - << "\n|| error ||_inf: " << std::setw(15) << solution.dist(analytic) - << std::endl; - - /*************************__________DEBUG__________*************************/ - // Create R scripts for plotting the solution - std::ostringstream file_name; #ifdef USE_MPI - file_name << "Rplot_" << mpi_coords[0] << "_" << mpi_coords[1] << ".R"; -#else - file_name << "Rplot_seriel.R"; -#endif - std::ofstream out(file_name.str()); + // Frobenius Norm Accumulator + double frob_norm_sendbuf[2] = { + solver.resid_norm(), solution.dist(analytic) + }; + // Square the norms (Accumulation of partial results is the square root of + // the sum of the squared partial norms) + frob_norm_sendbuf[0] = frob_norm_sendbuf[0] * frob_norm_sendbuf[0]; + frob_norm_sendbuf[1] = frob_norm_sendbuf[1] * frob_norm_sendbuf[1]; + // Maximum Norm Accumulator + double max_norm_sendbuf[2] = { + solver.resid_norm(), solution.dist(analytic) + }; - out << "sol <- matrix(c(" << solution(0); - for (size_t i = 1; i < solution.size(); ++i) { - out << ',' << solution(i); + // Global result reduction buffers + double frob_norm_recvbuf[2], max_norm_recvbuf[2]; + + // Accumulate global stats by sum/max reduction of local stats. + MPI_Reduce(frob_norm_sendbuf, frob_norm_recvbuf, 2, MPI_DOUBLE, MPI_SUM, + 0, MPI_COMM_WORLD); + MPI_Reduce(max_norm_sendbuf, max_norm_recvbuf, 2, MPI_DOUBLE, MPI_MAX, + 0, MPI_COMM_WORLD); + + // Finish by taking the square root of the accumulated global + // squared frobenius norms + frob_norm_recvbuf[0] = std::sqrt(frob_norm_recvbuf[0]); + frob_norm_recvbuf[1] = std::sqrt(frob_norm_recvbuf[1]); +#else + double frob_norm_recvbuf[2] = { + solver.resid_norm(), solution.dist(analytic) + }; + double max_norm_recvbuf[2] = { + solver.resid_norm(), solution.dist(analytic) + }; +#endif + + // calculate runtime time + auto stop = std::chrono::high_resolution_clock::now(); + auto time = std::chrono::duration_cast>(stop - start) + .count(); + + // Report global stats (from "root" only) + if (mpi_world_rank == 0) { + std::cout + << "Time [sec]: " << std::setw(15) << time + << "\n|| residuals ||_2: " << std::setw(15) << frob_norm_recvbuf[0] + << "\n|| residuals ||_inf: " << std::setw(15) << max_norm_recvbuf[0] + << "\n|| error ||_2: " << std::setw(15) << frob_norm_recvbuf[1] + << "\n|| error ||_inf: " << std::setw(15) << max_norm_recvbuf[1] + << std::endl; } - out << "), " << nx << ", " << ny << ")\n" -#ifdef USE_MPI - << "png(filename = \"mpi_plot_" << mpi_coords[0] << - "_" << mpi_coords[1] << ".png\")\n" -#else - << "png(filename = \"seriel_plot.png\")\n" -#endif -#ifdef USE_MPI - << "image(sol, zlim = c(-270, 270), axes = FALSE, main = \"Cart Coords: " - << mpi_coords[0] << ", " << mpi_coords[1] << "\")\n" -#else - << "image(sol, zlim = c(-270, 270), axes = FALSE, main = \"Seriel\")\n" -#endif - << "axis(1, at = c(0, 1), labels = c(" - << "\"" << std::setprecision(2) << xmin << "\", " - << "\"" << std::setprecision(2) << xmax << "\"))\n" - << "axis(2, at = c(0, 1), labels = c(" - << "\"" << std::setprecision(2) << ymin << "\", " - << "\"" << std::setprecision(2) << ymax << "\"))\n" - // << "axis(2, at = seq(100, 800, by = 100))" - << "dev.off()" << std::endl; - out.close(); - /*************************__________DEBUG__________*************************/ - - // MPI shutdown/cleanup #ifdef USE_MPI diff --git a/Exercise_01/utils.h b/Exercise_01/utils.h new file mode 100644 index 0000000..28b4a23 --- /dev/null +++ b/Exercise_01/utils.h @@ -0,0 +1,105 @@ + +/** + * Partitions an integer `num` into `div` summands and returns the `i`th + * of the partition. + * + * @example + * num = 17 + * div = 5 + * + * partition(17, 5, 0) -> 4 + * partition(17, 5, 1) -> 4 + * partition(17, 5, 2) -> 3 + * partition(17, 5, 3) -> 3 + * partition(17, 5, 5) -> 3 + * + * 17 = num = div * num + num % div = 4 + 4 + 3 + 3 + 3 + * 1st 2nd 3rd 4th 5th + * i=0 i=1 i=2 i=3 i=4 + */ +int partition(int num, int div, int i) { + return num / div + static_cast(i < (num % div)); +} + +/** + * Computes the partial sum of the first `i` integer `num` partitiones into + * `div` parts. + * + * @example + * num = 17 + * div = 5 + * + * partition_sum(17, 5, 0) -> 4 + * partition_sum(17, 5, 1) -> 8 + * partition_sum(17, 5, 2) -> 11 + * partition_sum(17, 5, 3) -> 14 + * partition_sum(17, 5, 5) -> 17 + */ +int partition_sum(int num, int div, int i) { + int sum = 0; + for (int j = 0; j <= i; ++j) { + sum += partition(num, div, j); + } + return sum; +} + +/** + * Factorized a integer number `num` into two multiplicative factors. + * These factors are as close together as possible. This means for example for + * square numbers that the two factors are the square root of `num`. + * + * Assumes small numbers and therefore uses a simple linear search. + * + * The first few integers are factorized as follows: + * + * 0 = 1 * 0 + * 1 = 1 * 1 + * 2 = 1 * 2 (is prime) + * 3 = 1 * 3 (is prime) + * 4 = 2 * 2 (is square) + * 5 = 1 * 5 (is prime) + * 6 = 2 * 3 + * 7 = 1 * 7 (is prime) + * 8 = 2 * 4 + * 9 = 3 * 3 (is square) + * 10 = 2 * 5 + * 11 = 1 * 11 (is prime) + * 12 = 3 * 4 + * 13 = 1 * 13 (is prime) + * 14 = 2 * 7 + * 15 = 3 * 5 + * 16 = 4 * 4 (is square) + * + * @param num integer to be factorized + * @param factor [out] output parameter of length 2 where the two factors are + * written into. + * + * @example + * int factors[2]; + * two_factors(15, factors); // -> factors = {3, 5} + */ +void two_factors(int num, int* factors) { + // In case of zero, set both to zero + if (!num) { + factors[0] = factors[1] = 0; + } + + // Ensure `num` is positive + if (num < 0) { + num *= -1; + } + + // Set initial factorization (this always works) + factors[0] = 1; + factors[1] = num; + + // Check all numbers `i` untill the integer square-root + for (int i = 2; i * i <= num; ++i) { + // Check if `i` is a divisor + if (!(num % i)) { + // Update factors as `i` divides `num` + factors[0] = i; + factors[1] = num / i; + } + } +}