ft:2D-omp outer loop parallelization

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

@@ -45,6 +45,7 @@ int get_dim_start(int dim_rank, int dim_comm_size, int dim_size)
 
 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)
 {
     solver->imax    = params->imax;
@@ -139,25 +140,23 @@ void solveRB(Solver* solver)
     }
 
     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);
 
-    for (int it = 0; it < itermax; ++it) {
+        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);
 
-        res = 0.0; 
+        for (int it = 0; it < itermax; ++it) {
 
-#pragma omp parallel reduction(+ : 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);
-
-            
-            int jsw = ((li_start) % 2 == 0) == ((lj_start) % 2 == 0) ? 1 : 2;
+#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;
@@ -170,13 +169,17 @@ void solveRB(Solver* solver)
                                            idy2);
 
                         P(i, j) -= factor * r;
-                        res += r * r; /* reduction variable */
+                        local_res += r * r; /* reduction variable */
                     }
                     isw = 3 - isw;
                 }
-#pragma omp barrier
+ #pragma omp barrier
                 jsw = 3 - jsw;
             }
+#pragma omp critical
+            {
+                res += local_res;
+            }
             if (lj_start == 0)
                 for (int i = li_start + 1; i < limax + 1; i++)
                     P(i, 0) = P(i, 1);
@@ -189,21 +192,20 @@ void solveRB(Solver* solver)
             if (limax == imax)
                 for (int j = lj_start + 1; j < ljmax + 1; j++)
                     P(limax + 1, j) = P(limax, j);
-
-        } 
-
-        res /= (double)(imax * jmax); 
-
+#pragma omp single
+            {
+                res /= (double)(imax * jmax);
 #ifdef DEBUG
-        printf("%d Residuum: %e\n", it, res);
+                printf("%d Residuum: %e\n", it, res);
 #endif
-
-        if (res < epssq) {
-            printf("Solver took %d iterations to reach %e\n", it + 1, sqrt(res));
-            return; 
+            }
+#pragma omp barrier
+            if (res < epssq) {
+                break;
+            }
+#pragma omp barrier
         }
     }
-
     printf("Solver reached itermax (%d) with residual %e\n", itermax, sqrt(res));
 }