CarpetLib: New parameter use_loopcontrol_in_operators

When this parameters is set, then (some) prolongation/restriction operators are parallelized via LoopControl instead of in dh.cc.

13 files changed, 295 insertions, 197 deletions

diff --git a/Carpet/CarpetLib/param.ccl b/Carpet/CarpetLib/param.ccl
index 3af3c401d..541ec2f49 100644
--- a/Carpet/CarpetLib/param.ccl
+++ b/Carpet/CarpetLib/param.ccl
@@ -39,6 +39,12 @@ BOOLEAN interpolate_from_buffer_zones "Use buffer points for interpolation" STEE
 
 
 
+BOOLEAN use_loopcontrol_in_operators "Use LoopControl to parallelize AMR operators" STEERABLE=always
+{
+} "no"
+
+
+
 restricted:
 
 BOOLEAN use_higher_order_restriction "Use third order cell centered restriction operators instead of first order"
diff --git a/Carpet/CarpetLib/src/data.cc b/Carpet/CarpetLib/src/data.cc
index 2aed5d914..fa1b3c92e 100644
--- a/Carpet/CarpetLib/src/data.cc
+++ b/Carpet/CarpetLib/src/data.cc
@@ -60,6 +60,13 @@ call_operator (void
                ibbox3 const & restrict dstregbbox,
                void * const extraargs)
 {
+  DECLARE_CCTK_PARAMETERS;
+  
+  if (use_loopcontrol_in_operators) {
+    (* the_operator)
+      (src, srcpadext, srcext, dst, dstpadext, dstext, srcbbox, dstbbox,
+       srcregbbox, dstregbbox, extraargs);
+  } else {
 #ifndef _OPENMP
   (* the_operator)
     (src, srcpadext, srcext, dst, dstpadext, dstext, srcbbox, dstbbox,
@@ -111,6 +118,7 @@ call_operator (void
   assert (alldstregbboxes == ibset (dstregbbox));
 #  endif
 #endif
+  }
 }
 
 template <typename T>
diff --git a/Carpet/CarpetLib/src/interpolate_3d_2tl.cc b/Carpet/CarpetLib/src/interpolate_3d_2tl.cc
index 33c00d8fd..bace93e5f 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_2tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_2tl.cc
@@ -1,13 +1,13 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <cstdlib>
 
-#include <loopcontrol.h>
-
 #include "operator_prototypes_3d.hh"
 #include "typeprops.hh"
 
diff --git a/Carpet/CarpetLib/src/interpolate_3d_3tl.cc b/Carpet/CarpetLib/src/interpolate_3d_3tl.cc
index f6dc7790d..6e39691f3 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_3tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_3tl.cc
@@ -1,13 +1,13 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <cstdlib>
 
-#include <loopcontrol.h>
-
 #include "operator_prototypes_3d.hh"
 #include "typeprops.hh"
 
diff --git a/Carpet/CarpetLib/src/interpolate_3d_4tl.cc b/Carpet/CarpetLib/src/interpolate_3d_4tl.cc
index 4f52f5723..f41d0b12b 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_4tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_4tl.cc
@@ -1,13 +1,13 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <cstdlib>
 
-#include <loopcontrol.h>
-
 #include "operator_prototypes_3d.hh"
 #include "typeprops.hh"
 
diff --git a/Carpet/CarpetLib/src/interpolate_3d_5tl.cc b/Carpet/CarpetLib/src/interpolate_3d_5tl.cc
index 0cfa144ca..c84f2c793 100644
--- a/Carpet/CarpetLib/src/interpolate_3d_5tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_3d_5tl.cc
@@ -1,13 +1,13 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <cstdlib>
 
-#include <loopcontrol.h>
-
 #include "operator_prototypes_3d.hh"
 #include "typeprops.hh"
 
diff --git a/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc b/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
index ac04121c0..abcee485a 100644
--- a/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
+++ b/Carpet/CarpetLib/src/interpolate_eno_3d_3tl.cc
@@ -1,13 +1,13 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
 #include <cstdlib>
 
-#include <loopcontrol.h>
-
 #include "operator_prototypes_3d.hh"
 #include "typeprops.hh"
 
diff --git a/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc b/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc
index 320635fd4..05aaad6ae 100644
--- a/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc
+++ b/Carpet/CarpetLib/src/prolongate_3d_cc_eno_rf2.cc
@@ -1,6 +1,8 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -966,6 +968,7 @@ namespace CarpetLib {
                             ibbox3 const & restrict regbbox,
                             void * extraargs)
   {
+    DECLARE_CCTK_PARAMETERS;
     assert (not extraargs);
     
     static_assert (ORDER>=0, "ORDER must be non-negative");
@@ -1091,6 +1094,8 @@ namespace CarpetLib {
     
     
     
+    if (not use_loopcontrol_in_operators) {
+
     // Loop over fine region
     // Label scheme: l 8 fk fj fi
     
@@ -1284,9 +1289,86 @@ namespace CarpetLib {
     
     // end k loop
   l9:;
+  
+    
+  
+  } else { // use_loopcontrol_in_operators
     
+    // Loop over fine region
+#pragma omp parallel
+    CCTK_LOOP3(prolongate_3d_cc_eno_rf2,
+               i,j,k, 0,0,0, regiext,regjext,regkext,
+               dstipadext,dstjpadext,dstkpadext)
+    {
+      const ptrdiff_t is = (srcioff + i) / 2;
+      const ptrdiff_t js = (srcjoff + j) / 2;
+      const ptrdiff_t ks = (srckoff + k) / 2;
+      const ptrdiff_t im = (srcioff + i) % 2;
+      const ptrdiff_t jm = (srcjoff + j) % 2;
+      const ptrdiff_t km = (srckoff + k) % 2;
+      const ptrdiff_t id = dstioff + i;
+      const ptrdiff_t jd = dstjoff + j;
+      const ptrdiff_t kd = dstkoff + k;
+      if (km == 0) {
+        if (jm == 0) {
+          if (im == 0) {
+            check_indices3<T,ORDER,0,0,0> (is,js,ks, srciext, srcjext, srckext);
+            dst[DSTIND3(id,jd,kd)] =
+              interp3<T,ORDER,0,0,0>
+              (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
+          } else {
+            check_indices3<T,ORDER,1,0,0> (is,js,ks, srciext, srcjext, srckext);
+            dst[DSTIND3(id,jd,kd)] =
+              interp3<T,ORDER,1,0,0>
+              (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
+          }
+        } else {
+          if (im == 0) {
+            check_indices3<T,ORDER,0,1,0> (is,js,ks, srciext, srcjext, srckext);
+            dst[DSTIND3(id,jd,kd)] =
+              interp3<T,ORDER,0,1,0>
+              (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
+          } else {
+            check_indices3<T,ORDER,1,1,0> (is,js,ks, srciext, srcjext, srckext);
+            dst[DSTIND3(id,jd,kd)] =
+              interp3<T,ORDER,1,1,0>
+              (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
+          }
+        }
+      } else {
+        if (jm == 0) {
+          if (im == 0) {
+            check_indices3<T,ORDER,0,0,1> (is,js,ks, srciext, srcjext, srckext);
+            dst[DSTIND3(id,jd,kd)] =
+              interp3<T,ORDER,0,0,1>
+              (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
+          } else {
+            check_indices3<T,ORDER,1,0,1> (is,js,ks, srciext, srcjext, srckext);
+            dst[DSTIND3(id,jd,kd)] =
+              interp3<T,ORDER,1,0,1>
+              (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
+          }
+        } else {
+          if (im == 0) {
+            check_indices3<T,ORDER,0,1,1> (is,js,ks, srciext, srcjext, srckext);
+            dst[DSTIND3(id,jd,kd)] =
+              interp3<T,ORDER,0,1,1>
+              (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
+          } else {
+            check_indices3<T,ORDER,1,1,1> (is,js,ks, srciext, srcjext, srckext);
+            dst[DSTIND3(id,jd,kd)] =
+              interp3<T,ORDER,1,1,1>
+              (& src[SRCIND3(is,js,ks)], srcdi, srcdj, srcdk);
+          }
+        }
+      }
+    } CCTK_ENDLOOP3(prolongate_3d_cc_eno_rf2);
+    
+    } // if use_loopcontrol_in_operators
+
   }
   
+  
 
   // Specialise for complex types
   
diff --git a/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc
index 12404ea25..3d1ce9acf 100644
--- a/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc
+++ b/Carpet/CarpetLib/src/restrict_3d_cc_o3_rf2.cc
@@ -1,6 +1,8 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -138,103 +140,102 @@ namespace CarpetLib {
     
     
     // Loop over coarse region
-#pragma omp parallel for collapse(3)
-    for (int k=0; k<regkext; ++k) {
-      for (int j=0; j<regjext; ++j) {
-        for (int i=0; i<regiext; ++i) {
-          
-#ifdef CARPET_DEBUG
-    if(not (2 * k + 2 + srckoff < srckext and
-            2 * j + 2 + srcjoff < srcjext and
-            2 * i + 2 + srcioff < srciext and
-            2 * k - 1 + srckoff >= 0 and
-            2 * j - 1 + srcjoff >= 0 and
-            2 * i - 1 + srcioff >= 0))
+#pragma omp parallel
+    CCTK_LOOP3(restrict_3d_cc_o3_rf2,
+               i,j,k, 0,0,0, regiext,regjext,regkext,
+               dstipadext,dstjpadext,dstkpadext)
     {
-      cout << "restrict_3d_cc_o3_rf2.cc\n";
-      cout << "regext " << regext << "\n";
-      cout << "srcext " << srcext << "\n";
-      cout << "srcbbox=" << srcbbox << "\n";
-      cout << "dstbbox=" << dstbbox << "\n";
-      cout << "regbbox=" << regbbox << "\n";
-      cout << "i,j,k=" << i << " " << j << " " << k << "\n";
-      assert(2 * k + 2 + srckoff < srckext);
-      assert(2 * j + 2 + srcjoff < srcjext);
-      assert(2 * i + 2 + srcioff < srciext);
-      assert(2 * k - 1 + srckoff >= 0);
-      assert(2 * j - 1 + srcjoff >= 0);
-      assert(2 * i - 1 + srcioff >= 0);
-    }
-#endif    
-          dst [DSTIND3(i, j, k)] =
-            - f1 * src [SRCIND3(2*i-1, 2*j-1, 2*k-1)]
-            + f2 * src [SRCIND3(2*i  , 2*j-1, 2*k-1)]
-            + f2 * src [SRCIND3(2*i+1, 2*j-1, 2*k-1)]
-            - f1 * src [SRCIND3(2*i+2, 2*j-1, 2*k-1)]
-            + f2 * src [SRCIND3(2*i-1, 2*j  , 2*k-1)]
-            - f3 * src [SRCIND3(2*i  , 2*j  , 2*k-1)]
-            - f3 * src [SRCIND3(2*i+1, 2*j  , 2*k-1)]
-            + f2 * src [SRCIND3(2*i+2, 2*j  , 2*k-1)]
-            + f2 * src [SRCIND3(2*i-1, 2*j+1, 2*k-1)]
-            - f3 * src [SRCIND3(2*i  , 2*j+1, 2*k-1)]
-            - f3 * src [SRCIND3(2*i+1, 2*j+1, 2*k-1)]
-            + f2 * src [SRCIND3(2*i+2, 2*j+1, 2*k-1)]
-            - f1 * src [SRCIND3(2*i-1, 2*j+2, 2*k-1)]
-            + f2 * src [SRCIND3(2*i  , 2*j+2, 2*k-1)]
-            + f2 * src [SRCIND3(2*i+1, 2*j+2, 2*k-1)]
-            - f1 * src [SRCIND3(2*i+2, 2*j+2, 2*k-1)]
-            + f2 * src [SRCIND3(2*i-1, 2*j-1, 2*k  )]
-            - f3 * src [SRCIND3(2*i  , 2*j-1, 2*k  )]
-            - f3 * src [SRCIND3(2*i+1, 2*j-1, 2*k  )]
-            + f2 * src [SRCIND3(2*i+2, 2*j-1, 2*k  )]
-            - f3 * src [SRCIND3(2*i-1, 2*j  , 2*k  )]
-            + f4 * src [SRCIND3(2*i  , 2*j  , 2*k  )]
-            + f4 * src [SRCIND3(2*i+1, 2*j  , 2*k  )]
-            - f3 * src [SRCIND3(2*i+2, 2*j  , 2*k  )]
-            - f3 * src [SRCIND3(2*i-1, 2*j+1, 2*k  )]
-            + f4 * src [SRCIND3(2*i  , 2*j+1, 2*k  )]
-            + f4 * src [SRCIND3(2*i+1, 2*j+1, 2*k  )]
-            - f3 * src [SRCIND3(2*i+2, 2*j+1, 2*k  )]
-            + f2 * src [SRCIND3(2*i-1, 2*j+2, 2*k  )]
-            - f3 * src [SRCIND3(2*i  , 2*j+2, 2*k  )]
-            - f3 * src [SRCIND3(2*i+1, 2*j+2, 2*k  )]
-            + f2 * src [SRCIND3(2*i+2, 2*j+2, 2*k  )]
-            + f2 * src [SRCIND3(2*i-1, 2*j-1, 2*k+1)]
-            - f3 * src [SRCIND3(2*i  , 2*j-1, 2*k+1)]
-            - f3 * src [SRCIND3(2*i+1, 2*j-1, 2*k+1)]
-            + f2 * src [SRCIND3(2*i+2, 2*j-1, 2*k+1)]
-            - f3 * src [SRCIND3(2*i-1, 2*j  , 2*k+1)]
-            + f4 * src [SRCIND3(2*i  , 2*j  , 2*k+1)]
-            + f4 * src [SRCIND3(2*i+1, 2*j  , 2*k+1)]
-            - f3 * src [SRCIND3(2*i+2, 2*j  , 2*k+1)]
-            - f3 * src [SRCIND3(2*i-1, 2*j+1, 2*k+1)]
-            + f4 * src [SRCIND3(2*i  , 2*j+1, 2*k+1)]
-            + f4 * src [SRCIND3(2*i+1, 2*j+1, 2*k+1)]
-            - f3 * src [SRCIND3(2*i+2, 2*j+1, 2*k+1)]
-            + f2 * src [SRCIND3(2*i-1, 2*j+2, 2*k+1)]
-            - f3 * src [SRCIND3(2*i  , 2*j+2, 2*k+1)]
-            - f3 * src [SRCIND3(2*i+1, 2*j+2, 2*k+1)]
-            + f2 * src [SRCIND3(2*i+2, 2*j+2, 2*k+1)]
-            - f1 * src [SRCIND3(2*i-1, 2*j-1, 2*k+2)]
-            + f2 * src [SRCIND3(2*i  , 2*j-1, 2*k+2)]
-            + f2 * src [SRCIND3(2*i+1, 2*j-1, 2*k+2)]
-            - f1 * src [SRCIND3(2*i+2, 2*j-1, 2*k+2)]
-            + f2 * src [SRCIND3(2*i-1, 2*j  , 2*k+2)]
-            - f3 * src [SRCIND3(2*i  , 2*j  , 2*k+2)]
-            - f3 * src [SRCIND3(2*i+1, 2*j  , 2*k+2)]
-            + f2 * src [SRCIND3(2*i+2, 2*j  , 2*k+2)]
-            + f2 * src [SRCIND3(2*i-1, 2*j+1, 2*k+2)]
-            - f3 * src [SRCIND3(2*i  , 2*j+1, 2*k+2)]
-            - f3 * src [SRCIND3(2*i+1, 2*j+1, 2*k+2)]
-            + f2 * src [SRCIND3(2*i+2, 2*j+1, 2*k+2)]
-            - f1 * src [SRCIND3(2*i-1, 2*j+2, 2*k+2)]
-            + f2 * src [SRCIND3(2*i  , 2*j+2, 2*k+2)]
-            + f2 * src [SRCIND3(2*i+1, 2*j+2, 2*k+2)]
-            - f1 * src [SRCIND3(2*i+2, 2*j+2, 2*k+2)];
-
-        }
+      
+#ifdef CARPET_DEBUG
+      if(not (2 * k + 2 + srckoff < srckext and
+              2 * j + 2 + srcjoff < srcjext and
+              2 * i + 2 + srcioff < srciext and
+              2 * k - 1 + srckoff >= 0 and
+              2 * j - 1 + srcjoff >= 0 and
+              2 * i - 1 + srcioff >= 0))
+      {
+        cout << "restrict_3d_cc_o3_rf2.cc\n";
+        cout << "regext " << regext << "\n";
+        cout << "srcext " << srcext << "\n";
+        cout << "srcbbox=" << srcbbox << "\n";
+        cout << "dstbbox=" << dstbbox << "\n";
+        cout << "regbbox=" << regbbox << "\n";
+        cout << "i,j,k=" << i << " " << j << " " << k << "\n";
+        assert(2 * k + 2 + srckoff < srckext);
+        assert(2 * j + 2 + srcjoff < srcjext);
+        assert(2 * i + 2 + srcioff < srciext);
+        assert(2 * k - 1 + srckoff >= 0);
+        assert(2 * j - 1 + srcjoff >= 0);
+        assert(2 * i - 1 + srcioff >= 0);
       }
-    }
+#endif    
+      dst [DSTIND3(i, j, k)] =
+        - f1 * src [SRCIND3(2*i-1, 2*j-1, 2*k-1)]
+        + f2 * src [SRCIND3(2*i  , 2*j-1, 2*k-1)]
+        + f2 * src [SRCIND3(2*i+1, 2*j-1, 2*k-1)]
+        - f1 * src [SRCIND3(2*i+2, 2*j-1, 2*k-1)]
+        + f2 * src [SRCIND3(2*i-1, 2*j  , 2*k-1)]
+        - f3 * src [SRCIND3(2*i  , 2*j  , 2*k-1)]
+        - f3 * src [SRCIND3(2*i+1, 2*j  , 2*k-1)]
+        + f2 * src [SRCIND3(2*i+2, 2*j  , 2*k-1)]
+        + f2 * src [SRCIND3(2*i-1, 2*j+1, 2*k-1)]
+        - f3 * src [SRCIND3(2*i  , 2*j+1, 2*k-1)]
+        - f3 * src [SRCIND3(2*i+1, 2*j+1, 2*k-1)]
+        + f2 * src [SRCIND3(2*i+2, 2*j+1, 2*k-1)]
+        - f1 * src [SRCIND3(2*i-1, 2*j+2, 2*k-1)]
+        + f2 * src [SRCIND3(2*i  , 2*j+2, 2*k-1)]
+        + f2 * src [SRCIND3(2*i+1, 2*j+2, 2*k-1)]
+        - f1 * src [SRCIND3(2*i+2, 2*j+2, 2*k-1)]
+        + f2 * src [SRCIND3(2*i-1, 2*j-1, 2*k  )]
+        - f3 * src [SRCIND3(2*i  , 2*j-1, 2*k  )]
+        - f3 * src [SRCIND3(2*i+1, 2*j-1, 2*k  )]
+        + f2 * src [SRCIND3(2*i+2, 2*j-1, 2*k  )]
+        - f3 * src [SRCIND3(2*i-1, 2*j  , 2*k  )]
+        + f4 * src [SRCIND3(2*i  , 2*j  , 2*k  )]
+        + f4 * src [SRCIND3(2*i+1, 2*j  , 2*k  )]
+        - f3 * src [SRCIND3(2*i+2, 2*j  , 2*k  )]
+        - f3 * src [SRCIND3(2*i-1, 2*j+1, 2*k  )]
+        + f4 * src [SRCIND3(2*i  , 2*j+1, 2*k  )]
+        + f4 * src [SRCIND3(2*i+1, 2*j+1, 2*k  )]
+        - f3 * src [SRCIND3(2*i+2, 2*j+1, 2*k  )]
+        + f2 * src [SRCIND3(2*i-1, 2*j+2, 2*k  )]
+        - f3 * src [SRCIND3(2*i  , 2*j+2, 2*k  )]
+        - f3 * src [SRCIND3(2*i+1, 2*j+2, 2*k  )]
+        + f2 * src [SRCIND3(2*i+2, 2*j+2, 2*k  )]
+        + f2 * src [SRCIND3(2*i-1, 2*j-1, 2*k+1)]
+        - f3 * src [SRCIND3(2*i  , 2*j-1, 2*k+1)]
+        - f3 * src [SRCIND3(2*i+1, 2*j-1, 2*k+1)]
+        + f2 * src [SRCIND3(2*i+2, 2*j-1, 2*k+1)]
+        - f3 * src [SRCIND3(2*i-1, 2*j  , 2*k+1)]
+        + f4 * src [SRCIND3(2*i  , 2*j  , 2*k+1)]
+        + f4 * src [SRCIND3(2*i+1, 2*j  , 2*k+1)]
+        - f3 * src [SRCIND3(2*i+2, 2*j  , 2*k+1)]
+        - f3 * src [SRCIND3(2*i-1, 2*j+1, 2*k+1)]
+        + f4 * src [SRCIND3(2*i  , 2*j+1, 2*k+1)]
+        + f4 * src [SRCIND3(2*i+1, 2*j+1, 2*k+1)]
+        - f3 * src [SRCIND3(2*i+2, 2*j+1, 2*k+1)]
+        + f2 * src [SRCIND3(2*i-1, 2*j+2, 2*k+1)]
+        - f3 * src [SRCIND3(2*i  , 2*j+2, 2*k+1)]
+        - f3 * src [SRCIND3(2*i+1, 2*j+2, 2*k+1)]
+        + f2 * src [SRCIND3(2*i+2, 2*j+2, 2*k+1)]
+        - f1 * src [SRCIND3(2*i-1, 2*j-1, 2*k+2)]
+        + f2 * src [SRCIND3(2*i  , 2*j-1, 2*k+2)]
+        + f2 * src [SRCIND3(2*i+1, 2*j-1, 2*k+2)]
+        - f1 * src [SRCIND3(2*i+2, 2*j-1, 2*k+2)]
+        + f2 * src [SRCIND3(2*i-1, 2*j  , 2*k+2)]
+        - f3 * src [SRCIND3(2*i  , 2*j  , 2*k+2)]
+        - f3 * src [SRCIND3(2*i+1, 2*j  , 2*k+2)]
+        + f2 * src [SRCIND3(2*i+2, 2*j  , 2*k+2)]
+        + f2 * src [SRCIND3(2*i-1, 2*j+1, 2*k+2)]
+        - f3 * src [SRCIND3(2*i  , 2*j+1, 2*k+2)]
+        - f3 * src [SRCIND3(2*i+1, 2*j+1, 2*k+2)]
+        + f2 * src [SRCIND3(2*i+2, 2*j+1, 2*k+2)]
+        - f1 * src [SRCIND3(2*i-1, 2*j+2, 2*k+2)]
+        + f2 * src [SRCIND3(2*i  , 2*j+2, 2*k+2)]
+        + f2 * src [SRCIND3(2*i+1, 2*j+2, 2*k+2)]
+        - f1 * src [SRCIND3(2*i+2, 2*j+2, 2*k+2)];
+      
+    } CCTK_ENDLOOP3(restrict_3d_cc_o3_rf2);
     
   }
   
diff --git a/Carpet/CarpetLib/src/restrict_3d_cc_o5_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_cc_o5_rf2.cc
index 2cd1b1f1f..5c2f88642 100644
--- a/Carpet/CarpetLib/src/restrict_3d_cc_o5_rf2.cc
+++ b/Carpet/CarpetLib/src/restrict_3d_cc_o5_rf2.cc
@@ -1,6 +1,8 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -124,8 +126,6 @@ namespace CarpetLib {
     
     DECLARE_CCTK_PARAMETERS;
     
-    typedef typename typeprops<T>::real RT;
-    
     
     
     if (any (srcbbox.stride() >= regbbox.stride() or
@@ -205,41 +205,40 @@ namespace CarpetLib {
     
     
     // Loop over coarse region
-#pragma omp parallel for collapse(3)
-    for (int k=0; k<regkext; ++k) {
-      for (int j=0; j<regjext; ++j) {
-        for (int i=0; i<regiext; ++i) {
+#pragma omp parallel
+    CCTK_LOOP3(restrict_3d_cc_o5_rf2,
+               i,j,k, 0,0,0, regiext,regjext,regkext,
+               dstipadext,dstjpadext,dstkpadext)
+    {
           
 #ifdef CARPET_DEBUG
-    if(not (2 * k + 3 + srckoff < srckext and
-            2 * j + 3 + srcjoff < srcjext and
-            2 * i + 3 + srcioff < srciext and
-            2 * k - 2 + srckoff >= 0 and
-            2 * j - 2 + srcjoff >= 0 and
-            2 * i - 2 + srcioff >= 0))
-    {
-      cout << "restrict_3d_cc_o3_rf2.cc\n";
-      cout << "regext " << regext << "\n";
-      cout << "srcext " << srcext << "\n";
-      cout << "srcbbox=" << srcbbox << "\n";
-      cout << "dstbbox=" << dstbbox << "\n";
-      cout << "regbbox=" << regbbox << "\n";
-      cout << "i,j,k=" << i << " " << j << " " << k << "\n";
-      assert(2 * k + 2 + srckoff < srckext);
-      assert(2 * j + 2 + srcjoff < srcjext);
-      assert(2 * i + 2 + srcioff < srciext);
-      assert(2 * k - 1 + srckoff >= 0);
-      assert(2 * j - 1 + srcjoff >= 0);
-      assert(2 * i - 1 + srcioff >= 0);
-    }
+      if(not (2 * k + 3 + srckoff < srckext and
+              2 * j + 3 + srcjoff < srcjext and
+              2 * i + 3 + srcioff < srciext and
+              2 * k - 2 + srckoff >= 0 and
+              2 * j - 2 + srcjoff >= 0 and
+              2 * i - 2 + srcioff >= 0))
+      {
+        cout << "restrict_3d_cc_o3_rf2.cc\n";
+        cout << "regext " << regext << "\n";
+        cout << "srcext " << srcext << "\n";
+        cout << "srcbbox=" << srcbbox << "\n";
+        cout << "dstbbox=" << dstbbox << "\n";
+        cout << "regbbox=" << regbbox << "\n";
+        cout << "i,j,k=" << i << " " << j << " " << k << "\n";
+        assert(2 * k + 2 + srckoff < srckext);
+        assert(2 * j + 2 + srcjoff < srcjext);
+        assert(2 * i + 2 + srcioff < srciext);
+        assert(2 * k - 1 + srckoff >= 0);
+        assert(2 * j - 1 + srcjoff >= 0);
+        assert(2 * i - 1 + srcioff >= 0);
+      }
 #endif    
-          dst [DSTIND3(i, j, k)] =
-            CarpetLib::restrict3
-            (& src[SRCIND3(2*i, 2*j, 2*k)], srcdi, srcdj, srcdk);
+      dst [DSTIND3(i, j, k)] =
+        CarpetLib::restrict3
+        (& src[SRCIND3(2*i, 2*j, 2*k)], srcdi, srcdj, srcdk);
 
-        }
-      }
-    }
+    } CCTK_ENDLOOP3(restrict_3d_cc_o5_rf2);
     
   }
   
diff --git a/Carpet/CarpetLib/src/restrict_3d_cc_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_cc_rf2.cc
index ee3fafc18..e186cbced 100644
--- a/Carpet/CarpetLib/src/restrict_3d_cc_rf2.cc
+++ b/Carpet/CarpetLib/src/restrict_3d_cc_rf2.cc
@@ -1,6 +1,8 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -127,26 +129,25 @@ namespace CarpetLib {
     
     
     // Loop over coarse region
-#pragma omp parallel for collapse(3)
-    for (int k=0; k<regkext; ++k) {
-      for (int j=0; j<regjext; ++j) {
-        for (int i=0; i<regiext; ++i) {
-          
-          // TODO: Introduce higher-order restriction operators (but
-          // don't use these for hydro!)
-          dst [DSTIND3(i, j, k)] =
-            + f1*f1*f1 * src [SRCIND3(2*i  , 2*j  , 2*k  )]
-            + f2*f1*f1 * src [SRCIND3(2*i+1, 2*j  , 2*k  )]
-            + f1*f2*f1 * src [SRCIND3(2*i  , 2*j+1, 2*k  )]
-            + f2*f2*f1 * src [SRCIND3(2*i+1, 2*j+1, 2*k  )]
-            + f1*f1*f2 * src [SRCIND3(2*i  , 2*j  , 2*k+1)]
-            + f2*f1*f2 * src [SRCIND3(2*i+1, 2*j  , 2*k+1)]
-            + f1*f2*f2 * src [SRCIND3(2*i  , 2*j+1, 2*k+1)]
-            + f2*f2*f2 * src [SRCIND3(2*i+1, 2*j+1, 2*k+1)];
-          
-        }
-      }
-    }
+#pragma omp parallel
+    CCTK_LOOP3(restrict_3d_cc_rf2,
+               i,j,k, 0,0,0, regiext,regjext,regkext,
+               dstipadext,dstjpadext,dstkpadext)
+    {
+      
+      // TODO: Introduce higher-order restriction operators (but
+      // don't use these for hydro!)
+      dst [DSTIND3(i, j, k)] =
+        + f1*f1*f1 * src [SRCIND3(2*i  , 2*j  , 2*k  )]
+        + f2*f1*f1 * src [SRCIND3(2*i+1, 2*j  , 2*k  )]
+        + f1*f2*f1 * src [SRCIND3(2*i  , 2*j+1, 2*k  )]
+        + f2*f2*f1 * src [SRCIND3(2*i+1, 2*j+1, 2*k  )]
+        + f1*f1*f2 * src [SRCIND3(2*i  , 2*j  , 2*k+1)]
+        + f2*f1*f2 * src [SRCIND3(2*i+1, 2*j  , 2*k+1)]
+        + f1*f2*f2 * src [SRCIND3(2*i  , 2*j+1, 2*k+1)]
+        + f2*f2*f2 * src [SRCIND3(2*i+1, 2*j+1, 2*k+1)];
+      
+    } CCTK_ENDLOOP3(restrict_3d_cc_rf2);
     
   }
   
diff --git a/Carpet/CarpetLib/src/restrict_3d_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_rf2.cc
index c1dade7fb..9f0936e1f 100644
--- a/Carpet/CarpetLib/src/restrict_3d_rf2.cc
+++ b/Carpet/CarpetLib/src/restrict_3d_rf2.cc
@@ -1,6 +1,8 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -113,15 +115,15 @@ namespace CarpetLib {
     
     
     // Loop over coarse region
-    for (int k=0; k<regkext; ++k) {
-      for (int j=0; j<regjext; ++j) {
-        for (int i=0; i<regiext; ++i) {
-          
-          dst [DSTIND3(i, j, k)] = src [SRCIND3(2*i, 2*j, 2*k)];
-          
-        }
-      }
-    }
+#pragma omp parallel
+    CCTK_LOOP3(restrict_3d_rf2,
+               i,j,k, 0,0,0, regiext,regjext,regkext,
+               dstipadext,dstjpadext,dstkpadext)
+    {
+      
+      dst [DSTIND3(i, j, k)] = src [SRCIND3(2*i, 2*j, 2*k)];
+      
+    } CCTK_ENDLOOP3(restrict_3d_rf2);
     
   }
   
diff --git a/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc b/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc
index e1539aa06..62450333c 100644
--- a/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc
+++ b/Carpet/CarpetLib/src/restrict_3d_vc_rf2.cc
@@ -1,6 +1,8 @@
 #include <cctk.h>
 #include <cctk_Parameters.h>
 
+#include <loopcontrol.h>
+
 #include <algorithm>
 #include <cassert>
 #include <cmath>
@@ -152,8 +154,6 @@ namespace CarpetLib {
     
     DECLARE_CCTK_PARAMETERS;
     
-    typedef typename typeprops<T>::real RT;
-    
     // false: vertex centered, true: cell centered
     ivect const icent (centi, centj, centk);
     assert (all (icent == 0 or icent == 1));
@@ -250,36 +250,35 @@ namespace CarpetLib {
     
     
     // Loop over coarse region
-#pragma omp parallel for collapse(3)
-    for (int k=0; k<regkext; ++k) {
-      for (int j=0; j<regjext; ++j) {
-        for (int i=0; i<regiext; ++i) {
-#ifdef CARPET_DEBUG
-    if(not (2 * k + centk < srckext and
-            2 * j + centj < srcjext and
-            2 * i + centi < srciext))
+#pragma omp parallel
+    CCTK_LOOP3(restrict_3d_vc_rf2,
+               i,j,k, 0,0,0, regiext,regjext,regkext,
+               dstipadext,dstjpadext,dstkpadext)
     {
-      cout << "restrict_3d_vc_rf2.cc\n";
-      cout << "regext " << regext << "\n";
-      cout << "srcext " << srcext << "\n";
-      cout << "srcbbox=" << srcbbox << "\n";
-      cout << "dstbbox=" << dstbbox << "\n";
-      cout << "regbbox=" << regbbox << "\n";
-      cout << "srcregbbox=" << srcregbbox << "\n";
-      cout << "icent=" << icent << "\n";
-    }
-    assert(2 * k + centk < srckext and
-           2 * j + centj < srcjext and
-           2 * i + centi < srciext);
-#endif    
-          
-          dst [DSTIND3(i, j, k)] =
-            restrict3<T,centi,centj,centk>::call
-            (& src[SRCIND3(2*i, 2*j, 2*k)], srcdi, srcdj, srcdk);
-          
-        }
+#ifdef CARPET_DEBUG
+      if(not (2 * k + centk < srckext and
+              2 * j + centj < srcjext and
+              2 * i + centi < srciext))
+      {
+        cout << "restrict_3d_vc_rf2.cc\n";
+        cout << "regext " << regext << "\n";
+        cout << "srcext " << srcext << "\n";
+        cout << "srcbbox=" << srcbbox << "\n";
+        cout << "dstbbox=" << dstbbox << "\n";
+        cout << "regbbox=" << regbbox << "\n";
+        cout << "srcregbbox=" << srcregbbox << "\n";
+        cout << "icent=" << icent << "\n";
       }
-    }
+      assert(2 * k + centk < srckext and
+             2 * j + centj < srcjext and
+             2 * i + centi < srciext);
+#endif    
+      
+      dst [DSTIND3(i, j, k)] =
+        restrict3<T,centi,centj,centk>::call
+        (& src[SRCIND3(2*i, 2*j, 2*k)], srcdi, srcdj, srcdk);
+      
+    } CCTK_ENDLOOP3(restrict_3d_vc_rf2);
     
   }