ft:2D-omp outer loop parallelization

PoissonSolver/2D-omp/src/main.c

@@ -6,10 +6,10 @@
 #include <stdlib.h>
 
 #include "likwid-marker.h"
+#include "omp.h"
 #include "parameter.h"
 #include "solver.h"
 #include "timing.h"
-#include "omp.h"
 
 #define LIKWID_PROFILE(tag, call)                                                        \
     startTime = getTimeStamp();                                                          \
@@ -18,17 +18,13 @@
     LIKWID_MARKER_STOP(#tag);                                                            \
     endTime = getTimeStamp();
 
-enum VARIANT { SOR = 1, RB, RBA };
-
 int main(int argc, char** argv)
 {
     int volatile dummy = 0;
-    int variant = RB;
     double startTime, endTime;
     Parameter params;
     Solver solver;
     initParameter(&params);
-    LIKWID_MARKER_INIT;
 #pragma omp parallel
     {
         if (dummy == 1 || omp_get_thread_num() == 0)
@@ -40,37 +36,13 @@ int main(int argc, char** argv)
     }
 
     readParameter(&params, argv[1]);
-    // printParameter(&params);
-    if (argc == 3) {
-        variant = atoi(argv[2]);
-    }
-    if (argc == 4) {
-        sscanf("%lf", argv[3], &params.omg);
-    }
 
     initSolver(&solver, &params, 2);
-    writeResult(&solver, "p-0.dat");
+    startTime = getTimeStamp();
-
+    solveRB(&solver);
-    switch (variant) {
+    endTime = getTimeStamp();
-    case SOR:
-        printf("Plain SOR\n");
-        fflush(stdout);
-        LIKWID_PROFILE("SOR", solve);
-        break;
-    case RB:
-        printf("Red-black SOR\n");
-        fflush(stdout);
-        LIKWID_PROFILE("RB", solveRB);
-        break;
-    case RBA:
-        printf("Red-black SOR with acceleration\n");
-        fflush(stdout);
-        LIKWID_PROFILE("RBA", solveRBA);
-        break;
-    }
     printf(" %.2fs\n", endTime - startTime);
-    writeResult(&solver, "p-final.dat");
+    writeResult(&solver, "p.dat");
 
-    LIKWID_MARKER_CLOSE;
     return EXIT_SUCCESS;
 }

PoissonSolver/2D-omp/src/solver.c

@@ -9,10 +9,42 @@
 #include "allocate.h"
 #include "parameter.h"
 #include "solver.h"
+#include <omp.h>
+#include <stddef.h>
 
 #define PI        3.14159265358979323846
-#define P(i, j)   p[(j) * (imax + 2) + (i)]
+#define P(i, j)   p[(i) * (jmax + 2) + (j)]
-#define RHS(i, j) rhs[(j) * (imax + 2) + (i)]
+#define RHS(i, j) rhs[(i) * (jmax + 2) + (j)]
+
+void omp_create_dim(int num_threads, int* dim)
+{
+    int center_int = (int)sqrt(num_threads);
+    dim[0]         = num_threads;
+    dim[1]         = 1;
+    for (int num = center_int; num != 1; num--)
+        if (num_threads % num == 0) {
+            dim[0] = num;
+            dim[1] = num_threads / num;
+            break;
+        }
+}
+
+int distribute_dim(int dim_rank, int dim_comm_size, int dim_size)
+{
+    int base_size = dim_size / dim_comm_size;
+    return dim_rank < (dim_size % dim_comm_size) ? base_size + 1 : base_size;
+}
+
+int get_dim_start(int dim_rank, int dim_comm_size, int dim_size)
+{
+    int dim_start = 0;
+    for (int i = 0; i < dim_rank; i++)
+        dim_start += distribute_dim(i, dim_comm_size, dim_size);
+    return dim_start;
+}
+
+int get_x_choord(const int proc_num, const int const* dims) { return proc_num / dims[1]; }
+int get_y_choord(const int proc_num, const int const* dims) { return proc_num % dims[1]; }
 
 void initSolver(Solver* solver, Parameter* params, int problem)
 {
@@ -35,220 +67,146 @@ void initSolver(Solver* solver, Parameter* params, int problem)
     double dy       = solver->dy;
     double* p       = solver->p;
     double* rhs     = solver->rhs;
-    #pragma omp parallel for collapse(2)
+    int dim[2]      = { 0 };
-    for (int j = 0; j < jmax + 2; j++) {
+    int num_threads = 1;
-        for (int i = 0; i < imax + 2; i++) {
+#pragma omp parallel
+    {
+#pragma omp critical
+        num_threads = omp_get_num_threads();
+    }
+    omp_create_dim(num_threads, dim);
+    printf("%d: { %d, %d}\n", num_threads, dim[0], dim[1]);
+#pragma omp parallel
+    {
+        int jsw, isw;
+        double local_res = 0.0;
+        int li_start     = get_dim_start(get_x_choord(omp_get_thread_num(), dim),
+            dim[0],
+            solver->imax);
+        int lj_start     = get_dim_start(get_y_choord(omp_get_thread_num(), dim),
+            dim[1],
+            solver->jmax);
+        int limax = li_start + distribute_dim(get_x_choord(omp_get_thread_num(), dim),
+                                   dim[0],
+                                   solver->imax);
+        int ljmax = lj_start + distribute_dim(get_y_choord(omp_get_thread_num(), dim),
+                                   dim[1],
+                                   solver->jmax);
+
+        for (int i = li_start; i < limax + 2; i++) {
+            for (int j = lj_start; j < ljmax + 2; j++) {
                 P(i, j) = sin(2.0 * PI * i * dx * 2.0) + sin(2.0 * PI * j * dy * 2.0);
             }
         }
 
         if (problem == 2) {
-        #pragma omp parallel for collapse(2)
+            for (int i = li_start; i < limax + 2; i++) {
-        for (int j = 0; j < jmax + 2; j++) {
+                for (int j = lj_start; j < ljmax + 2; j++) {
-            for (int i = 0; i < imax + 2; i++) {
                     RHS(i, j) = sin(2.0 * PI * i * dx);
                 }
             }
         } else {
-        #pragma omp parallel for collapse(2)
+            for (int i = li_start; i < limax + 2; i++) {
-        for (int j = 0; j < jmax + 2; j++) {
+                for (int j = lj_start; j < ljmax + 2; j++) {
-            for (int i = 0; i < imax + 2; i++) {
                     RHS(i, j) = 0.0;
                 }
             }
         }
     }
-
-void solve(Solver* solver)
-{
-    int imax      = solver->imax;
-    int jmax      = solver->jmax;
-    double eps    = solver->eps;
-    int itermax   = solver->itermax;
-    double dx2    = solver->dx * solver->dx;
-    double dy2    = solver->dy * solver->dy;
-    double idx2   = 1.0 / dx2;
-    double idy2   = 1.0 / dy2;
-    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
-    double* p     = solver->p;
-    double* rhs   = solver->rhs;
-    double epssq  = eps * eps;
-    int it        = 0;
-    double res    = 1.0;
-    char filename[20];
-
-    while ((res >= epssq) && (it < itermax)) {
-        res = 0.0;
-
-        for (int j = 1; j < jmax + 1; j++) {
-            for (int i = 1; i < imax + 1; i++) {
-
-                double r = RHS(i, j) -
-                           ((P(i - 1, j) - 2.0 * P(i, j) + P(i + 1, j)) * idx2 +
-                               (P(i, j - 1) - 2.0 * P(i, j) + P(i, j + 1)) * idy2);
-
-                P(i, j) -= (factor * r);
-                res += (r * r);
-            }
-        }
-
-        for (int i = 1; i < imax + 1; i++) {
-            P(i, 0)        = P(i, 1);
-            P(i, jmax + 1) = P(i, jmax);
-        }
-
-        for (int j = 1; j < jmax + 1; j++) {
-            P(0, j)        = P(1, j);
-            P(imax + 1, j) = P(imax, j);
-        }
-
-        res = res / (double)(imax * jmax);
-#ifdef DEBUG
-        printf("%d Residuum: %e\n", it, res);
-#endif
-#ifdef ANIMATE
-        sprintf(filename, "p-%d.dat", it);
-        writeResult(solver, filename);
-#endif
-        it++;
-    }
-
-    printf("%d, %f\n", it, solver->omega);
 }
 
 void solveRB(Solver* solver)
 {
-    int imax      = solver->imax;
-    int jmax      = solver->jmax;
-    double eps    = solver->eps;
-    int itermax   = solver->itermax;
-    double dx2    = solver->dx * solver->dx;
-    double dy2    = solver->dy * solver->dy;
-    double idx2   = 1.0 / dx2;
-    double idy2   = 1.0 / dy2;
-    double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
-    double* p     = solver->p;
-    double* rhs   = solver->rhs;
-    double epssq  = eps * eps;
-    int it        = 0;
-    double res    = 1.0;
-    int pass, jsw, isw;
 
-    while ((res >= epssq) && (it < itermax)) {
+    const int imax     = solver->imax;
-        res = 0.0;
+    const int jmax     = solver->jmax;
-        jsw = 1;
+    const int itermax  = solver->itermax;
+    const double epssq = solver->eps * solver->eps;
 
-        for (pass = 0; pass < 2; pass++) {
+    const double dx2    = solver->dx * solver->dx;
-            isw = jsw;
+    const double dy2    = solver->dy * solver->dy;
-            #pragma omp parallel for firstprivate(isw)
+    const double idx2   = 1.0 / dx2;
-            for (int j = 1; j < jmax + 1; j++) {
+    const double idy2   = 1.0 / dy2;
-                for (int i = isw; i < imax + 1; i += 2) {
+    const double factor = solver->omega * 0.5 * (dx2 * dy2) / (dx2 + dy2);
+
+    double* __restrict p   = solver->p;
+    double* __restrict rhs = solver->rhs;
+
+    int dim[2] = { 0 };
+#pragma omp parallel
+#pragma omp single
+    {
+        omp_create_dim(omp_get_num_threads(), dim);
+    }
+
+    double res = 0.0;
+#pragma omp parallel shared(res)
+    {
+        const int tid      = omp_get_thread_num();
+        const int li_start = get_dim_start(get_x_choord(tid, dim), dim[0], imax);
+        const int lj_start = get_dim_start(get_y_choord(tid, dim), dim[1], jmax);
+
+        const int limax = li_start + distribute_dim(get_x_choord(tid, dim), dim[0], imax);
+        const int ljmax = lj_start + distribute_dim(get_y_choord(tid, dim), dim[1], jmax);
+
+        for (int it = 0; it < itermax; ++it) {
+
+#pragma omp single
+            {
+                res = 0;
+            }
+            double local_res = 0;
+            int jsw          = ((li_start) % 2 == 0) == ((lj_start) % 2 == 0) ? 1 : 2;
+
+            for (int pass = 0; pass < 2; ++pass) {
+                int isw = jsw;
+                for (int i = li_start + 1; i < limax + 1; ++i) {
+                    for (int j = lj_start + isw; j < ljmax + 1; j += 2) {
 
                         double r = RHS(i, j) -
                                    ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
-                                   (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
+                                       (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) *
+                                           idy2);
 
-                    P(i, j) -= (factor * r);
+                        P(i, j) -= factor * r;
-                    res += (r * r);
+                        local_res += r * r; /* reduction variable */
                     }
                     isw = 3 - isw;
                 }
+ #pragma omp barrier
                 jsw = 3 - jsw;
             }
-        #pragma omp parallel for 
+#pragma omp critical
-        for (int i = 1; i < imax + 1; i++) {
+            {
+                res += local_res;
+            }
+            if (lj_start == 0)
+                for (int i = li_start + 1; i < limax + 1; i++)
                     P(i, 0) = P(i, 1);
-            P(i, jmax + 1) = P(i, jmax);
+            if (ljmax == jmax)
-        }
+                for (int i = li_start + 1; i < limax + 1; i++)
-        #pragma omp parallel for
+                    P(i, ljmax + 1) = P(i, ljmax);
-        for (int j = 1; j < jmax + 1; j++) {
+            if (li_start == 0)
+                for (int j = lj_start + 1; j < ljmax + 1; j++)
                     P(0, j) = P(1, j);
-            P(imax + 1, j) = P(imax, j);
+            if (limax == imax)
-        }
+                for (int j = lj_start + 1; j < ljmax + 1; j++)
-
+                    P(limax + 1, j) = P(limax, j);
-        res = res / (double)(imax * jmax);
+#pragma omp single
+            {
+                res /= (double)(imax * jmax);
 #ifdef DEBUG
                 printf("%d Residuum: %e\n", it, res);
 #endif
-// #ifdef ANIMATE
-//         sprintf(filename, "p-%d.dat", it);
-//         writeResult(solver, filename);
-// #endif
-//         it++;
             }
-
+#pragma omp barrier
-    printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
+            if (res < epssq) {
-    printf("%d, %f\n", it, solver->omega);
+                break;
             }
-
+#pragma omp barrier
-void solveRBA(Solver* solver)
-{
-    int imax      = solver->imax;
-    int jmax      = solver->jmax;
-    double eps    = solver->eps;
-    int itermax   = solver->itermax;
-    double dx2    = solver->dx * solver->dx;
-    double dy2    = solver->dy * solver->dy;
-    double idx2   = 1.0 / dx2;
-    double idy2   = 1.0 / dy2;
-    double factor = 0.5 * (dx2 * dy2) / (dx2 + dy2);
-    double rho    = solver->rho;
-    double* p     = solver->p;
-    double* rhs   = solver->rhs;
-    double epssq  = eps * eps;
-    int it        = 0;
-    double res    = 1.0;
-    int pass, jsw, isw;
-    double omega = 1.0;
-
-    while ((res >= epssq) && (it < itermax)) {
-        res = 0.0;
-        jsw = 1;
-
-        for (pass = 0; pass < 2; pass++) {
-            isw = jsw;
-
-            for (int j = 1; j < jmax + 1; j++) {
-                for (int i = isw; i < imax + 1; i += 2) {
-
-                    double r = RHS(i, j) -
-                               ((P(i + 1, j) - 2.0 * P(i, j) + P(i - 1, j)) * idx2 +
-                                   (P(i, j + 1) - 2.0 * P(i, j) + P(i, j - 1)) * idy2);
-
-                    P(i, j) -= (omega * factor * r);
-                    res += (r * r);
         }
-                isw = 3 - isw;
     }
-            jsw   = 3 - jsw;
+    printf("Solver reached itermax (%d) with residual %e\n", itermax, sqrt(res));
-            omega = (it == 0 && pass == 0 ? 1.0 / (1.0 - 0.5 * rho * rho)
-                                          : 1.0 / (1.0 - 0.25 * rho * rho * omega));
-        }
-
-        for (int i = 1; i < imax + 1; i++) {
-            P(i, 0)        = P(i, 1);
-            P(i, jmax + 1) = P(i, jmax);
-        }
-
-        for (int j = 1; j < jmax + 1; j++) {
-            P(0, j)        = P(1, j);
-            P(imax + 1, j) = P(imax, j);
-        }
-
-        res = res / (double)(imax * jmax);
-#ifdef DEBUG
-        printf("%d Residuum: %e Omega: %e\n", it, res, omega);
-#endif
-#ifdef ANIMATE
-        sprintf(filename, "p-%d.dat", it);
-        writeResult(solver, filename);
-#endif
-        it++;
-    }
-
-    // printf("Final omega: %f\n", omega);
-    // printf("Solver took %d iterations to reach %f\n", it, sqrt(res));
-    printf("%d, %f\n", it, omega);
 }
 
 void writeResult(Solver* solver, char* filename)
@@ -265,8 +223,8 @@ void writeResult(Solver* solver, char* filename)
         exit(EXIT_FAILURE);
     }
 
-    for (int j = 0; j < jmax + 2; j++) {
+    for (int i = 1; i < imax + 1; i++) {
-        for (int i = 0; i < imax + 2; i++) {
+        for (int j = 1; j < jmax + 1; j++) {
             fprintf(fp, "%f ", P(i, j));
         }
         fprintf(fp, "\n");