CarpetReduce: Add support for OpenMP

darcs-hash:20070821190958-dae7b-1ed180793dd5d6482f48aa0bbf13bd5523a4b40b.gz

5 files changed, 37 insertions, 13 deletions

diff --git a/Carpet/CarpetReduce/interface.ccl b/Carpet/CarpetReduce/interface.ccl
index 5c55b526b..71b96c3a6 100644
--- a/Carpet/CarpetReduce/interface.ccl
+++ b/Carpet/CarpetReduce/interface.ccl
@@ -2,6 +2,7 @@
 
 IMPLEMENTS: reduce
 
+uses include header: defs.hh
 uses include header: dist.hh
 uses include header: vect.hh
 
diff --git a/Carpet/CarpetReduce/src/mask_carpet.cc b/Carpet/CarpetReduce/src/mask_carpet.cc
index 11c2f6642..5593dfbfa 100644
--- a/Carpet/CarpetReduce/src/mask_carpet.cc
+++ b/Carpet/CarpetReduce/src/mask_carpet.cc
@@ -133,6 +133,7 @@ namespace CarpetMask {
                 
                 // Set weight on the boundary to 1/2
                 assert (dim == 3);
+#pragma omp parallel for
                 for (int k=imin[2]; k<imax[2]; ++k) {
                   for (int j=imin[1]; j<imax[1]; ++j) {
                     for (int i=imin[0]; i<imax[0]; ++i) {
@@ -195,6 +196,7 @@ namespace CarpetMask {
               
               // Set weight in the restricted region to 0
               assert (dim == 3);
+#pragma omp parallel for
               for (int k=imin[2]; k<imax[2]; ++k) {
                 for (int j=imin[1]; j<imax[1]; ++j) {
                   for (int i=imin[0]; i<imax[0]; ++i) {
@@ -212,6 +214,7 @@ namespace CarpetMask {
           vector<int> mask (cctk_lsh[0] * cctk_lsh[1] * cctk_lsh[2]);
           
           assert (dim == 3);
+#pragma omp parallel for
           for (int k=0; k<cctk_lsh[2]; ++k) {
             for (int j=0; j<cctk_lsh[1]; ++j) {
               for (int i=0; i<cctk_lsh[0]; ++i) {
@@ -248,6 +251,7 @@ namespace CarpetMask {
                 
                 // Set weight on the boundary to 1/2
                 assert (dim == 3);
+#pragma omp parallel for
                 for (int k=imin[2]; k<imax[2]; ++k) {
                   for (int j=imin[1]; j<imax[1]; ++j) {
                     for (int i=imin[0]; i<imax[0]; ++i) {
@@ -266,6 +270,7 @@ namespace CarpetMask {
           } // for d
           
           assert (dim == 3);
+#pragma omp parallel for
           for (int k=0; k<cctk_lsh[2]; ++k) {
             for (int j=0; j<cctk_lsh[1]; ++j) {
               for (int i=0; i<cctk_lsh[0]; ++i) {
diff --git a/Carpet/CarpetReduce/src/mask_coords.c b/Carpet/CarpetReduce/src/mask_coords.c
index e148d56bb..ccb931d8e 100644
--- a/Carpet/CarpetReduce/src/mask_coords.c
+++ b/Carpet/CarpetReduce/src/mask_coords.c
@@ -81,6 +81,7 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
             CCTK_VInfo (CCTK_THORNSTRING,
                         "Setting boundary points in direction %d face %d to weight 0", d, f);
           }
+#pragma omp parallel for
           for (k=bmin[2]; k<bmax[2]; ++k) {
             for (j=bmin[1]; j<bmax[1]; ++j) {
               for (i=bmin[0]; i<bmax[0]; ++i) {
@@ -127,6 +128,7 @@ CoordBase_SetupMask (CCTK_ARGUMENTS)
                 CCTK_VInfo (CCTK_THORNSTRING,
                             "Setting non-staggered boundary points in direction %d face %d to weight 1/2", d, f);
               }
+#pragma omp parallel for
               for (k=bmin[2]; k<bmax[2]; ++k) {
                 for (j=bmin[1]; j<bmax[1]; ++j) {
                   for (i=bmin[0]; i<bmax[0]; ++i) {
diff --git a/Carpet/CarpetReduce/src/mask_init.c b/Carpet/CarpetReduce/src/mask_init.c
index ed1297871..7365b0637 100644
--- a/Carpet/CarpetReduce/src/mask_init.c
+++ b/Carpet/CarpetReduce/src/mask_init.c
@@ -35,6 +35,7 @@ MaskBase_InitMask (CCTK_ARGUMENTS)
     CCTK_INFO ("Initialising to weight 1");
   }
   
+#pragma omp parallel for
   for (k=0; k<cctk_lsh[2]; ++k) {
     for (j=0; j<cctk_lsh[1]; ++j) {
       for (i=0; i<cctk_lsh[0]; ++i) {
@@ -74,6 +75,7 @@ MaskBase_InitMask (CCTK_ARGUMENTS)
         }
         
         /* Loop over the boundary slab */
+#pragma omp parallel for
         for (k=imin[2]; k<imax[2]; ++k) {
           for (j=imin[1]; j<imax[1]; ++j) {
             for (i=imin[0]; i<imax[0]; ++i) {
@@ -142,6 +144,7 @@ MaskBase_InitMask (CCTK_ARGUMENTS)
           }
           
           /* Loop over the boundary slab */
+#pragma omp parallel for
           for (k=imin[2]; k<imax[2]; ++k) {
             for (j=imin[1]; j<imax[1]; ++j) {
               for (i=imin[0]; i<imax[0]; ++i) {
diff --git a/Carpet/CarpetReduce/src/reduce.cc b/Carpet/CarpetReduce/src/reduce.cc
index cfdf61b85..f57810d54 100644
--- a/Carpet/CarpetReduce/src/reduce.cc
+++ b/Carpet/CarpetReduce/src/reduce.cc
@@ -15,6 +15,7 @@
 #include <util_ErrorCodes.h>
 #include <util_Table.h>
 
+#include <defs.hh>
 #include <dist.hh>
 #include <vect.hh>
 
@@ -492,21 +493,33 @@ namespace CarpetReduce {
       imax[d] = lsh[d] - (bbox[2*d+1] ? 0 : nghostzones[d]);
     }
     static_assert (dim==3, "Only 3 dimensions are currently supported");
-    for (int k=imin[2]; k<imax[2]; ++k) {
-      for (int j=imin[1]; j<imax[1]; ++j) {
-        for (int i=imin[0]; i<imax[0]; ++i) {
-	  const int index = i + lsh[0] * (j + lsh[1] * k);
-          CCTK_REAL const w = (weight ? weight[index] : 1.0) * levfac;
-          T myinval = T(0);
-          for (size_t tl=0; tl<inarrays.size(); ++tl) {
-            myinval += (static_cast<const T*>(inarrays.at(tl))[index]
-                        * tfacs.at(tl));
+#pragma omp parallel
+    {
+      T myoutval_local;
+      OP::initialise (myoutval_local);
+      CCTK_REAL mycnt_local = 0;
+#pragma omp for
+      for (int k=imin[2]; k<imax[2]; ++k) {
+        for (int j=imin[1]; j<imax[1]; ++j) {
+          for (int i=imin[0]; i<imax[0]; ++i) {
+            const int index = i + lsh[0] * (j + lsh[1] * k);
+            CCTK_REAL const w = weight ? weight[index] * levfac : levfac;
+            T myinval = T(0);
+            for (size_t tl=0; tl<inarrays.size(); ++tl) {
+              myinval +=
+                static_cast<const T*>(inarrays.AT(tl))[index] * tfacs.AT(tl);
+            }
+            OP::reduce (myoutval_local, myinval, w);
+            mycnt_local += w;
           }
-	  OP::reduce (myoutval, myinval, w);
-	  mycnt += w;
-	}
+        }
       }
-    }
+#pragma omp critical
+      {
+        OP::reduce (myoutval, myoutval_local, mycnt_local);
+        mycnt += mycnt_local;
+      }
+    } // end omp parallel
     *(T*)outval = myoutval;
     *(T*)cnt    = mycnt;
   }