add: utils,

wip: ex_01 (almost done)

Exercise_01/main.cpp

@@ -7,17 +7,13 @@
 #include <iomanip>
 #include <sstream>
 #include <functional>
+#include <chrono>
 #define _USE_MATH_DEFINES   /* enables math constants from cmath */
 #include <cmath>
 #ifdef USE_MPI
 #include <mpi.h>
 #endif
 
-
-#include <fstream>
-#include <iomanip>
-
-
 #include "Matrix.h"
 #include "Solver.h"
 #include "utils.h"
@@ -32,19 +28,24 @@ 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] << " <resolution> <iterations>" << 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") {
@@ -53,10 +54,13 @@ int main(int argn, char* argv[]) {
         std::cerr << "error: Parsing arg 1 <dim> failed" << std::endl;
         return -1;
     }
+#endif
+    size_t resolution;
     if (!(std::istringstream(argv[2]) >> resolution)) {
         std::cerr << "error: Parsing arg 2 <resolution> failed" << std::endl;
         return -1;
     }
+    size_t iterations;
     if (!(std::istringstream(argv[3]) >> iterations)) {
         std::cerr << "error: Parsing arg 3 <iterations> failed" << std::endl;
         return -1;
@@ -71,18 +75,22 @@ int main(int argn, char* argv[]) {
             << "out of domain [0, 65536]" << std::endl;
         return -1;
     }
-    }
 
     // Report configuration
+    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;
+    }
+
+    /************************* 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<Max>()
-        << std::endl;
-
     // Get solution
     Matrix<double>& solution = solver.solution();
 
@@ -314,51 +291,57 @@ int main(int argn, char* argv[]) {
         }
     }
 
-    // Calculate error
+#ifdef USE_MPI
+    // 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<Max>(), solution.dist<Max>(analytic)
+    };
+
+    // 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<Max>(), solution.dist<Max>(analytic)
+    };
+#endif
+
+    // calculate runtime time
+    auto stop = std::chrono::high_resolution_clock::now();
+    auto time = std::chrono::duration_cast<std::chrono::duration<double>>(stop - start)
+        .count();
+
+    // Report global stats (from "root" only)
+    if (mpi_world_rank == 0) {
         std::cout
-        <<   "|| error ||_2:       " << std::setw(15) << solution.dist(analytic)
+            <<   "Time [sec]:          " << std::setw(15) << time
-        << "\n|| error ||_inf:     " << std::setw(15) << solution.dist<Max>(analytic)
+            << "\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;
-
-    /*************************__________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());
-
-    out << "sol <- matrix(c(" << solution(0);
-    for (size_t i = 1; i < solution.size(); ++i) {
-        out << ',' << solution(i);
     }
-    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

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<int>(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;
+        }
+    }
+}