[test-suite] Corrections for MiniGMG

[homerdin]

Summary

The openMP loops were being macro'd out along with the pragmas. This patch allows the application to run fully.

The problem size was decreased as the test's runtime increased substantially. Additional verification output has been enable using FP_TOLERANCE.

[hfinkel]

mg.c

103

Why are you changing the number of periods in the printf strings (here and below)? Is this change necessary?

[homerdin]

I had removed the periods after getting a FP Comparison failed, not a numeric difference between '.' and '.' " error. They were not the issue and I have added them back in.

[thegameg]

This seems to break here: http://green.lab.llvm.org/green/job/lnt-test-suite-x86_64-Os-flto/3288/consoleFull

/Users/buildslave/jenkins/workspace/lnt-test-suite-x86_64-Os-flto/lnt-sandbox/build/tools/fpcmp: Compared: 4.440892e-16 and -1.110223e-15
abs. diff = 1.554312e-15 rel.diff = 1.400000e+00
Out of tolerance: rel/abs: 1.000000e-05/0.000000e+00

Can you please take a look @homerdin or revert it? Thanks!

[thegameg]

I reverted both commits in r335009 and r335010 to fix the bots.

[thegameg]

Ah, sorry, I just saw that your second commit fixed the issue. Relanded in:

https://reviews.llvm.org/rL335012

and

https://reviews.llvm.org/rL335011.

CMakeLists.txt

 set(PROG miniGMG)
+list(APPEND CXXFLAGS -D__PRINT_NORM)
 list(APPEND LDFLAGS -lm)
-set(RUN_OPTIONS 6  2 2 2  1 1 1)
+set(FP_TOLERANCE 0.00001)
+set(RUN_OPTIONS 5  2 2 2  1 1 1)
 llvm_multisource()

apply_op.inc

#ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k,s;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double h2inv = 1.0/(domain->h[level]*domain->h[level]);
     double * __restrict__ x      = domain->subdomains[box].levels[level].grids[     x_id] + ghosts*(1+pencil+plane); // i.e. [0] = first non ghost zone point
     double * __restrict__ Ax     = domain->subdomains[box].levels[level].grids[    Ax_id] + ghosts*(1+pencil+plane); 
     double * __restrict__ alpha  = domain->subdomains[box].levels[level].grids[ __alpha ] + ghosts*(1+pencil+plane);
     double * __restrict__ beta_i = domain->subdomains[box].levels[level].grids[ __beta_i] + ghosts*(1+pencil+plane);
     double * __restrict__ beta_j = domain->subdomains[box].levels[level].grids[ __beta_j] + ghosts*(1+pencil+plane);
     double * __restrict__ beta_k = domain->subdomains[box].levels[level].grids[ __beta_k] + ghosts*(1+pencil+plane);
     double * __restrict__ lambda = domain->subdomains[box].levels[level].grids[ __lambda] + ghosts*(1+pencil+plane);
 
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<dim_k;k++){
     for(j=0;j<dim_j;j++){
     for(i=0;i<dim_i;i++){
       int ijk = i + j*pencil + k*plane;
       double helmholtz =  a*alpha[ijk]*x[ijk] - b*h2inv*(
                              beta_i[ijk+1     ]*( x[ijk+1     ]-x[ijk       ] )
                             -beta_i[ijk       ]*( x[ijk       ]-x[ijk-1     ] )
                             +beta_j[ijk+pencil]*( x[ijk+pencil]-x[ijk       ] )
                             -beta_j[ijk       ]*( x[ijk       ]-x[ijk-pencil] )
                             +beta_k[ijk+plane ]*( x[ijk+plane ]-x[ijk       ] )
                             -beta_k[ijk       ]*( x[ijk       ]-x[ijk-plane ] )
                           );
       Ax[ijk] = helmholtz;
     }}}
   }
-#endif
+
   domain->cycles.apply_op[level] += (uint64_t)(CycleTime()-_timeStart);
 }
 //------------------------------------------------------------------------------------------------------------------------------

exchange_boundary.inc

void exchange_boundary(domain_type *domain, int level, int grid_id, int exchange_faces, int exchange_edges, int exchange_corners){
 
   // pack MPI send buffers...
   _timeStart = CycleTime();
   #pragma omp parallel for schedule(static,1)
   for(buffer=domain->bufferCopy_Pack_Start;buffer<domain->bufferCopy_Pack_End;buffer++){
     if( (domain->bufferCopies[level][buffer].isFace   && exchange_faces  ) || 
         (domain->bufferCopies[level][buffer].isEdge   && exchange_edges  ) ||
         (domain->bufferCopies[level][buffer].isCorner && exchange_corners) ){
           DoBufferCopy(domain,level,grid_id,buffer);
   }}
   _timeEnd = CycleTime();
+  domain->cycles.pack[level] += (_timeEnd-_timeStart);
+ 
+  // loop through MPI send buffers and post Isend's...
+  _timeStart = CycleTime();
+  #ifdef __MPI_THREAD_MULTIPLE
+  #pragma omp parallel for schedule(dynamic,1)
+  #endif
+  for(n=nMessages/2;n<nMessages;n++){
+    MPI_Isend(buffers_packed[n],sizes_packed[n],MPI_DOUBLE,ranks_packed[n],tags_packed[n],MPI_COMM_WORLD,&requests_packed[n]);
+  }
+  _timeEnd = CycleTime();
+  domain->cycles.send[level] += (_timeEnd-_timeStart);
+  #endif
+
+
+  // exchange locally... try and hide within Isend latency... 
+  _timeStart = CycleTime();
+#ifdef OMP
+  #pragma omp parallel for schedule(static,1)
+#endif
+  for(buffer=domain->bufferCopy_Local_Start;buffer<domain->bufferCopy_Local_End;buffer++){
+    if( (domain->bufferCopies[level][buffer].isFace   && exchange_faces  ) || 
+        (domain->bufferCopies[level][buffer].isEdge   && exchange_edges  ) ||
+        (domain->bufferCopies[level][buffer].isCorner && exchange_corners) ){
+          DoBufferCopy(domain,level,grid_id,buffer);
+  }}
+  _timeEnd = CycleTime();
   domain->cycles.grid2grid[level] += (_timeEnd-_timeStart);
 
 
   // wait for MPI to finish...
   #ifdef __MPI 
   _timeStart = CycleTime();
   MPI_Waitall(nMessages,requests_packed,status_packed);
   _timeEnd = CycleTime();
   domain->cycles.wait[level] += (_timeEnd-_timeStart);
-  #endif
+
   // unpack MPI receive buffers 
   _timeStart = CycleTime();
   #pragma omp parallel for schedule(static,1)
   for(buffer=domain->bufferCopy_Unpack_Start;buffer<domain->bufferCopy_Unpack_End;buffer++){
     if( (domain->bufferCopies[level][buffer].isFace   && exchange_faces  ) || 
         (domain->bufferCopies[level][buffer].isEdge   && exchange_edges  ) ||
         (domain->bufferCopies[level][buffer].isCorner && exchange_corners) ){
           DoBufferCopy(domain,level,grid_id,buffer);
   }}
   _timeEnd = CycleTime();
   domain->cycles.unpack[level] += (_timeEnd-_timeStart);
   #endif
  
  
   domain->cycles.communication[level] += (uint64_t)(CycleTime()-_timeCommunicationStart);
 }

mg.h

int create_domain(domain_type * domain,
                    int subdomains_per_rank_in_i, int subdomains_per_rank_in_j, int subdomains_per_rank_in_k,
                    int ranks_in_i, int ranks_in_j, int ranks_in_k,
                    int rank, 
                    int *boundary_conditions,
                    int numGrids, int ghosts, int numLevels);
 void      MGBuild(domain_type * domain, double a, double b, double h0);
 void      MGSolve(domain_type * domain, int u_id, int F_id, double a, double b, double desired_mg_norm);
 void print_timing(domain_type *domain);
+void DoBufferCopy(domain_type *domain, int level, int grid_id, int buffer);
 //------------------------------------------------------------------------------------------------------------------------------
 #endif

mg.c

int create_domain(domain_type * domain,
               int *boundary_conditions,
               int numGrids, int ghosts, int numLevels
              ){
   int level;
   int  i, j, k;
   int di,dj,dk;
   // FIX ghosts = array (varies with level)
   domain->rank = rank;
-  if(domain->rank==0){printf("creating domain...       ");fflush(stdout);}
+  if(domain->rank==0){printf("creating domain.       ");fflush(stdout);}
   //FIX, limit ghosts based on coarsest problem size
   if(ghosts>(subdomain_dim_i>>(numLevels-1))){if(domain->rank==0)printf("#ghosts(%d)>bottom grid size(%d)\n",ghosts,subdomain_dim_i>>(numLevels-1));exit(0);}
 
   if( (subdomain_dim_i!=subdomain_dim_j)||(subdomain_dim_j!=subdomain_dim_k)||(subdomain_dim_i!=subdomain_dim_k) ){
   if(domain->rank==0)printf("subdomain_dim's must be equal\n");exit(0);
   }
   uint64_t memory_allocated =0;
   // processes are laid out in x, then y, then z
int create_domain(domain_type * domain,
   printf("  allocated %llu MB\n",memory_allocated>>20);
   fflush(stdout);
   }
   return(memory_allocated);
 }
 
 
 void destroy_domain(domain_type * domain){
-  if(domain->rank==0){printf("deallocating domain...   ");fflush(stdout);}
+  if(domain->rank==0){printf("deallocating domain.   ");fflush(stdout);}
   int box;for(box=0;box<domain->subdomains_per_rank;box++){
     destroy_subdomain(&domain->subdomains[box]);
   }
   free(domain->subdomains);
   // FIX, free buffer_copies, etc...
   if(domain->rank==0){printf("done\n");fflush(stdout);}
 }
 
void MGResetTimers(domain_type * domain){
   domain->MGSolves_performed          = 0;
   domain->Krylov_iterations           = 0;
   domain->CAKrylov_formations_of_G    = 0;
 }
 //------------------------------------------------------------------------------------------------------------------------------
 void MGBuild(domain_type * domain, double a, double b, double h0){
   int box,level,numLevels = domain->numLevels;
 
-  if(domain->rank==0){printf("MGBuild...                      ");fflush(stdout);}
+  if(domain->rank==0){printf("MGBuild.                      ");fflush(stdout);}
 
   // initialize timers...
   MGResetTimers(domain);
   domain->cycles.MGBuild = 0;
 
   uint64_t _timeStartMGBuild = CycleTime();
 
 
void MGSolve(domain_type * domain, int u_id, int F_id, double a, double b, double desired_mg_norm){
   int numLevels  = domain->numLevels;
   #ifdef __TEST_MG_CONVERGENCE
   int maxVCycles       = 20;
   #else
   int maxVCycles       = 10;
   #endif
 
   //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
-  if(domain->rank==0){printf("MGSolve...                      ");fflush(stdout);}
+  if(domain->rank==0){printf("MGSolve.                      ");fflush(stdout);}
   uint64_t _timeStartMGSolve = CycleTime();
 
   //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   // make initial guess for e (=0) and setup the RHS
   #if 1
    zero_grid(domain,0,e_id);                // ee = 0
   scale_grid(domain,0,R_id,1.0,F_id);       // R_id = F_id
   #else

miniGMG.c

//int log2_coarse_dim = 2; // i.e. coarsen to 4^3
   create_domain(&domain,
                            subdomain_dim_i,subdomain_dim_j,subdomain_dim_k,
                            subdomains_per_rank_in_i,subdomains_per_rank_in_j,subdomains_per_rank_in_k,
                            ranks_in_i,ranks_in_j,ranks_in_k,
                            MPI_Rank,
                            boundary_conditions,
                            __NumGrids,1,levels_in_vcycle);
   double h0=1.0/((double)(domain.dim.i));
-  if(MPI_Rank==0){printf("initializing alpha, beta, RHS for the ``hard problem''...");fflush(stdout);}
+  if(MPI_Rank==0){printf("initializing alpha, beta, RHS for the ``hard problem''.");fflush(stdout);}
   double  a=0.9;       // i.e. good Helmholtz
   double  b=0.9;
   initialize_problem(&domain,0,h0,a,b);
   if(MPI_Rank==0){printf("done\n");fflush(stdout);}
   //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   MGBuild(&domain,a,b,h0); // restrictions, dominant eigenvalue, etc...
   //- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
   int s,sMax=2;

miniGMG.reference_output

 1 MPI Tasks of 1 threads
 truncating the v-cycle at 2^3 subdomains
-creating domain...       done
-  64 x 64 x 64 (per subdomain)
-  128 x 128 x 128 (per process)
-  128 x 128 x 128 (overall)
+creating domain.       done
+  32 x 32 x 32 (per subdomain)
+  64 x 64 x 64 (per process)
+  64 x 64 x 64 (overall)
   1-deep ghost zones
-  allocated 246 MB
-initializing alpha, beta, RHS for the ``hard problem''...
-  average value of f =   0.000000000000e+00
+  allocated 34 MB
+initializing alpha, beta, RHS for the ``hard problem''.
+  average value of f =  -1.110223024625e-15
 done
-MGBuild...                      
-  level= 0, eigenvalue_max ~= 1.000000e+00
-  level= 1, eigenvalue_max ~= -1.000000e+00
-  level= 2, eigenvalue_max ~= -1.000000e+00
-  level= 3, eigenvalue_max ~= -1.000000e+00
-  level= 4, eigenvalue_max ~= -1.000000e+00
-  level= 5, eigenvalue_max ~= -1.000000e+00
+MGBuild.                      
+  level= 0, eigenvalue_max ~= 1.999996e+00
+  level= 1, eigenvalue_max ~= 1.999984e+00
+  level= 2, eigenvalue_max ~= 1.999935e+00
+  level= 3, eigenvalue_max ~= 1.999740e+00
+  level= 4, eigenvalue_max ~= 1.998958e+00
 done
-MGSolve...                      done
-MGSolve...                      done
-Error test: h = 7.812500e-03, max = 0.000000e+00
-Error test: h = 7.812500e-03, L2  = 0.000000e+00
-deallocating domain...   done
+MGSolve.                      
+v-cycle= 1, norm=0.00104306838695688285 (1.043068e-03)
+v-cycle= 2, norm=0.00008508164981969001 (8.508165e-05)
+v-cycle= 3, norm=0.00001867173484947213 (1.867173e-05)
+v-cycle= 4, norm=0.00000204150498158241 (2.041505e-06)
+v-cycle= 5, norm=0.00000036859010963302 (3.685901e-07)
+v-cycle= 6, norm=0.00000004500376124219 (4.500376e-08)
+v-cycle= 7, norm=0.00000000777275418452 (7.772754e-09)
+v-cycle= 8, norm=0.00000000096441180332 (9.644118e-10)
+v-cycle= 9, norm=0.00000000016688265714 (1.668827e-10)
+v-cycle=10, norm=0.00000000002055567632 (2.055568e-11)
+done
+MGSolve.                      
+v-cycle= 1, norm=0.00104306838695688285 (1.043068e-03)
+v-cycle= 2, norm=0.00008508164981969001 (8.508165e-05)
+v-cycle= 3, norm=0.00001867173484947213 (1.867173e-05)
+v-cycle= 4, norm=0.00000204150498158241 (2.041505e-06)
+v-cycle= 5, norm=0.00000036859010963302 (3.685901e-07)
+v-cycle= 6, norm=0.00000004500376124219 (4.500376e-08)
+v-cycle= 7, norm=0.00000000777275418452 (7.772754e-09)
+v-cycle= 8, norm=0.00000000096441180332 (9.644118e-10)
+v-cycle= 9, norm=0.00000000016688265714 (1.668827e-10)
+v-cycle=10, norm=0.00000000002055567632 (2.055568e-11)
+done
+Error test: h = 1.562500e-02, max = 3.529710e-04
+Error test: h = 1.562500e-02, L2  = 5.559935e-05
+deallocating domain.   done
 exit 0

misc.inc

#ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid = domain->subdomains[box].levels[level].grids[grid_id] + ghosts*(1+pencil+plane);
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=-ghosts;k<dim_k+ghosts;k++){
      for(j=-ghosts;j<dim_j+ghosts;j++){
       for(i=-ghosts;i<dim_i+ghosts;i++){
         int ijk = i + j*pencil + k*plane;
         grid[ijk] = 0.0;
     }}}
   }
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
 }
 
 
 
 //------------------------------------------------------------------------------------------------------------------------------
 void initialize_grid_to_scalar(domain_type *domain, int level, int grid_id, double scalar){
   // initializes the grid to a scalar while zero'ing the ghost zones...
   uint64_t _timeStart = CycleTime();
   int CollaborativeThreadingBoxSize = 100000; // i.e. never
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid = domain->subdomains[box].levels[level].grids[grid_id] + ghosts*(1+pencil+plane);
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=-ghosts;k<dim_k+ghosts;k++){
      for(j=-ghosts;j<dim_j+ghosts;j++){
       for(i=-ghosts;i<dim_i+ghosts;i++){
         int ijk = i + j*pencil + k*plane;
         int ghostZone = (i<0) || (j<0) || (k<0) || (i>=dim_i) || (j>=dim_j) || (k>=dim_k);
         grid[ijk] = ghostZone ? 0.0 : scalar;
     }}}
   }
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
 }
 
 
 //------------------------------------------------------------------------------------------------------------------------------
 void add_grids(domain_type *domain, int level, int id_c, double scale_a, int id_a, double scale_b, int id_b){ // c=scale_a*id_a + scale_b*id_b
   uint64_t _timeStart = CycleTime();
 
   int CollaborativeThreadingBoxSize = 100000; // i.e. never
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid_c = domain->subdomains[box].levels[level].grids[id_c] + ghosts*(1+pencil+plane);
     double * __restrict__ grid_a = domain->subdomains[box].levels[level].grids[id_a] + ghosts*(1+pencil+plane);
     double * __restrict__ grid_b = domain->subdomains[box].levels[level].grids[id_b] + ghosts*(1+pencil+plane);
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<dim_k;k++){
      for(j=0;j<dim_j;j++){
       for(i=0;i<dim_i;i++){
         int ijk = i + j*pencil + k*plane;
         grid_c[ijk] = scale_a*grid_a[ijk] + scale_b*grid_b[ijk];
     }}}
   }
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
 }
 
 
 //------------------------------------------------------------------------------------------------------------------------------
 void mul_grids(domain_type *domain, int level, int id_c, double scale, int id_a, int id_b){ // id_c=scale*id_a*id_b
   uint64_t _timeStart = CycleTime();
 
   int CollaborativeThreadingBoxSize = 100000; // i.e. never
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid_c = domain->subdomains[box].levels[level].grids[id_c] + ghosts*(1+pencil+plane);
     double * __restrict__ grid_a = domain->subdomains[box].levels[level].grids[id_a] + ghosts*(1+pencil+plane);
     double * __restrict__ grid_b = domain->subdomains[box].levels[level].grids[id_b] + ghosts*(1+pencil+plane);
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<dim_k;k++){
      for(j=0;j<dim_j;j++){
       for(i=0;i<dim_i;i++){
         int ijk = i + j*pencil + k*plane;
         grid_c[ijk] = scale*grid_a[ijk]*grid_b[ijk];
     }}}
   }
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
 }
 
 
 
 //------------------------------------------------------------------------------------------------------------------------------
 void scale_grid(domain_type *domain, int level, int id_c, double scale_a, int id_a){ // c[]=scale_a*a[]
   uint64_t _timeStart = CycleTime();
 
   int CollaborativeThreadingBoxSize = 100000; // i.e. never
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid_c = domain->subdomains[box].levels[level].grids[id_c] + ghosts*(1+pencil+plane);
     double * __restrict__ grid_a = domain->subdomains[box].levels[level].grids[id_a] + ghosts*(1+pencil+plane);
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<dim_k;k++){
      for(j=0;j<dim_j;j++){
       for(i=0;i<dim_i;i++){
         int ijk = i + j*pencil + k*plane;
         grid_c[ijk] = scale_a*grid_a[ijk];
     }}}
   }
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
 }
 
 
 //------------------------------------------------------------------------------------------------------------------------------
 double dot(domain_type *domain, int level, int id_a, int id_b){
   uint64_t _timeStart = CycleTime();
 
   int CollaborativeThreadingBoxSize = 100000; // i.e. never
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
 
   int box;
   double a_dot_b_domain =  0.0;
 #ifdef OMP
   // FIX, schedule(static) is a stand in to guarantee reproducibility...
   #pragma omp parallel for private(box) if(omp_across_boxes) reduction(+:a_dot_b_domain) schedule(static)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid_a = domain->subdomains[box].levels[level].grids[id_a] + ghosts*(1+pencil+plane); // i.e. [0] = first non ghost zone point
     double * __restrict__ grid_b = domain->subdomains[box].levels[level].grids[id_b] + ghosts*(1+pencil+plane);
     double a_dot_b_box = 0.0;
+#ifdef OMP
     #pragma omp parallel for private(i,j,k) if(omp_within_a_box) collapse(2) reduction(+:a_dot_b_box) schedule(static)
+#endif
     for(k=0;k<dim_k;k++){
     for(j=0;j<dim_j;j++){
     for(i=0;i<dim_i;i++){
       int ijk = i + j*pencil + k*plane;
       a_dot_b_box += grid_a[ijk]*grid_b[ijk];
     }}}
     a_dot_b_domain+=a_dot_b_box;
   }
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
 
   #ifdef __MPI
   uint64_t _timeStartAllReduce = CycleTime();
   double send = a_dot_b_domain;
   MPI_Allreduce(&send,&a_dot_b_domain,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
   uint64_t _timeEndAllReduce = CycleTime();
   domain->cycles.collectives[level]   += (uint64_t)(_timeEndAllReduce-_timeStartAllReduce);
double norm(domain_type *domain, int level, int grid_id){ // implements the max norm
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
 
   int box;
   double max_norm =  0.0;
 #ifdef OMP
   // FIX, schedule(static) is a stand in to guarantee reproducibility...
   #pragma omp parallel for private(box) if(omp_across_boxes) reduction(max:max_norm) schedule(static)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid   = domain->subdomains[box].levels[level].grids[ grid_id] + ghosts*(1+pencil+plane); // i.e. [0] = first non ghost zone point
     double box_norm = 0.0;
+#ifdef OMP
     #pragma omp parallel for private(i,j,k) if(omp_within_a_box) collapse(2) reduction(max:box_norm) schedule(static)
+#endif
     for(k=0;k<dim_k;k++){
     for(j=0;j<dim_j;j++){
     for(i=0;i<dim_i;i++){
       int ijk = i + j*pencil + k*plane;
       double fabs_grid_ijk = fabs(grid[ijk]);
       if(fabs_grid_ijk>box_norm){box_norm=fabs_grid_ijk;} // max norm
     }}}
     if(box_norm>max_norm){max_norm = box_norm;}
   } // box list
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
 
   #ifdef __MPI
   uint64_t _timeStartAllReduce = CycleTime();
   double send = max_norm;
   MPI_Allreduce(&send,&max_norm,1,MPI_DOUBLE,MPI_MAX,MPI_COMM_WORLD);
   uint64_t _timeEndAllReduce = CycleTime();
   domain->cycles.collectives[level]   += (uint64_t)(_timeEndAllReduce-_timeStartAllReduce);
double mean(domain_type *domain, int level, int id_a){
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
 
   int box;
   double sum_domain =  0.0;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes) reduction(+:sum_domain)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid_a = domain->subdomains[box].levels[level].grids[id_a] + ghosts*(1+pencil+plane); // i.e. [0] = first non ghost zone point
     double sum_box = 0.0;
+#ifdef OMP
     #pragma omp parallel for private(i,j,k) if(omp_within_a_box) collapse(2) reduction(+:sum_box)
+#endif
     for(k=0;k<dim_k;k++){
     for(j=0;j<dim_j;j++){
     for(i=0;i<dim_i;i++){
       int ijk = i + j*pencil + k*plane;
       sum_box += grid_a[ijk];
     }}}
     sum_domain+=sum_box;
   }
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
   double ncells_domain = (double)domain->dim.i*(double)domain->dim.j*(double)domain->dim.k;
 
   #ifdef __MPI
   uint64_t _timeStartAllReduce = CycleTime();
   double send = sum_domain;
   MPI_Allreduce(&send,&sum_domain,1,MPI_DOUBLE,MPI_SUM,MPI_COMM_WORLD);
   uint64_t _timeEndAllReduce = CycleTime();
#ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
 
   int box;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid_c = domain->subdomains[box].levels[level].grids[id_c] + ghosts*(1+pencil+plane); // i.e. [0] = first non ghost zone point
     double * __restrict__ grid_a = domain->subdomains[box].levels[level].grids[id_a] + ghosts*(1+pencil+plane); // i.e. [0] = first non ghost zone point
 
+#ifdef OMP
     #pragma omp parallel for private(i,j,k) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<dim_k;k++){
     for(j=0;j<dim_j;j++){
     for(i=0;i<dim_i;i++){
       int ijk = i + j*pencil + k*plane;
       grid_c[ijk] = grid_a[ijk] + shift_a;
     }}}
   }
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
 }
 
 //------------------------------------------------------------------------------------------------------------------------------
 void project_cell_to_face(domain_type *domain, int level, int id_cell, int id_face, int dir){
   uint64_t _timeStart = CycleTime();
   int CollaborativeThreadingBoxSize = 100000; // i.e. never
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double * __restrict__ grid_cell = domain->subdomains[box].levels[level].grids[id_cell] + ghosts*(1+pencil+plane);
     double * __restrict__ grid_face = domain->subdomains[box].levels[level].grids[id_face] + ghosts*(1+pencil+plane);
     int stride;
     switch(dir){
       case 0: stride =      1;break;//i-direction
       case 1: stride = pencil;break;//j-direction
       case 2: stride =  plane;break;//k-direction
     }
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<=dim_k;k++){ // <= to ensure you do low and high faces
      for(j=0;j<=dim_j;j++){
       for(i=0;i<=dim_i;i++){
         int ijk = i + j*pencil + k*plane;
         grid_face[ijk] = 0.5*(grid_cell[ijk-stride] + grid_cell[ijk]); // simple linear interpolation
     }}}
   }
-#endif
   domain->cycles.blas1[level] += (uint64_t)(CycleTime()-_timeStart);
 }

residual.inc

void residual(domain_type * domain, int level,  int res_id, int phi_id, int rhs_id, double a, double b){
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int pencil = domain->subdomains[box].levels[level].pencil;
     int  plane = domain->subdomains[box].levels[level].plane;
     int ghosts = domain->subdomains[box].levels[level].ghosts;
     int  dim_k = domain->subdomains[box].levels[level].dim.k;
     int  dim_j = domain->subdomains[box].levels[level].dim.j;
     int  dim_i = domain->subdomains[box].levels[level].dim.i;
     double h2inv = 1.0/(domain->h[level]*domain->h[level]);
     double * __restrict__ phi    = domain->subdomains[box].levels[level].grids[  phi_id] + ghosts*(1+pencil+plane); // i.e. [0] = first non ghost zone point
     double * __restrict__ rhs    = domain->subdomains[box].levels[level].grids[  rhs_id] + ghosts*(1+pencil+plane);
     double * __restrict__ alpha  = domain->subdomains[box].levels[level].grids[__alpha ] + ghosts*(1+pencil+plane);
     double * __restrict__ beta_i = domain->subdomains[box].levels[level].grids[__beta_i] + ghosts*(1+pencil+plane);
     double * __restrict__ beta_j = domain->subdomains[box].levels[level].grids[__beta_j] + ghosts*(1+pencil+plane);
     double * __restrict__ beta_k = domain->subdomains[box].levels[level].grids[__beta_k] + ghosts*(1+pencil+plane);
     double * __restrict__ res    = domain->subdomains[box].levels[level].grids[  res_id] + ghosts*(1+pencil+plane);
 
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<dim_k;k++){
     for(j=0;j<dim_j;j++){
     for(i=0;i<dim_i;i++){
       int ijk = i + j*pencil + k*plane;
       double helmholtz =  a*alpha[ijk]*phi[ijk]
                          -b*h2inv*(
                             beta_i[ijk+1     ]*( phi[ijk+1     ]-phi[ijk       ] )
                            -beta_i[ijk       ]*( phi[ijk       ]-phi[ijk-1     ] )
                            +beta_j[ijk+pencil]*( phi[ijk+pencil]-phi[ijk       ] )
                            -beta_j[ijk       ]*( phi[ijk       ]-phi[ijk-pencil] )
                            +beta_k[ijk+plane ]*( phi[ijk+plane ]-phi[ijk       ] )
                            -beta_k[ijk       ]*( phi[ijk       ]-phi[ijk-plane ] )
                           );
       res[ijk] = rhs[ijk]-helmholtz;
     }}}
   }
-#endif
   domain->cycles.residual[level] += (uint64_t)(CycleTime()-_timeStart);
 }
 
 #if 1
 //------------------------------------------------------------------------------------------------------------------------------
 // This version maximizes parallelism by parallelizing over the resultant coarse grid.  
 // Thus, 
 //  one parallelizes over the list of 2x2 fine-grid bars,
void residual_and_restriction(domain_type *domain, int level_f, int phi_id, int rhs_id, int level_c, int res_id, double a, double b){
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level_f].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level_f].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int kk,jj;
     int pencil_c = domain->subdomains[box].levels[level_c].pencil;
     int  plane_c = domain->subdomains[box].levels[level_c].plane;
     int ghosts_c = domain->subdomains[box].levels[level_c].ghosts;
     int  dim_k_c = domain->subdomains[box].levels[level_c].dim.k;
     int  dim_j_c = domain->subdomains[box].levels[level_c].dim.j;
     int  dim_i_c = domain->subdomains[box].levels[level_c].dim.i;
for(box=0;box<domain->subdomains_per_rank;box++){
     double * __restrict__ phi    = domain->subdomains[box].levels[level_f].grids[  phi_id] + ghosts_f*(1+pencil_f+plane_f); // i.e. [0] = first non ghost zone point
     double * __restrict__ rhs    = domain->subdomains[box].levels[level_f].grids[  rhs_id] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ alpha  = domain->subdomains[box].levels[level_f].grids[__alpha ] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ beta_i = domain->subdomains[box].levels[level_f].grids[__beta_i] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ beta_j = domain->subdomains[box].levels[level_f].grids[__beta_j] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ beta_k = domain->subdomains[box].levels[level_f].grids[__beta_k] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ res    = domain->subdomains[box].levels[level_c].grids[  res_id] + ghosts_c*(1+pencil_c+plane_c);
 
+#ifdef OMP
     #pragma omp parallel for private(kk,jj) if(omp_within_a_box) collapse(2)
+#endif
     for(kk=0;kk<dim_k_f;kk+=2){
     for(jj=0;jj<dim_j_f;jj+=2){
       int i,j,k;
       for(i=0;i<dim_i_c;i++){
         int ijk_c = (i) + (jj>>1)*pencil_c + (kk>>1)*plane_c;
         res[ijk_c] = 0.0;
       }
       for(k=kk;k<kk+2;k++){
for(jj=0;jj<dim_j_f;jj+=2){
                              -beta_j[ijk_f         ]*( phi[ijk_f         ]-phi[ijk_f-pencil_f] )
                              +beta_k[ijk_f+plane_f ]*( phi[ijk_f+plane_f ]-phi[ijk_f         ] )
                              -beta_k[ijk_f         ]*( phi[ijk_f         ]-phi[ijk_f-plane_f ] )
                             );
         res[ijk_c] += (rhs[ijk_f]-helmholtz)*0.125;
       }
     }}}}
   }
-#endif
   domain->cycles.residual[level_f] += (uint64_t)(CycleTime()-_timeStart);
 }
 #else
 //------------------------------------------------------------------------------------------------------------------------------
 // This version performs a 1D parallelization over the coarse-grid k-dimension (every two fine-grid planes)
 // It first zeros the coarse grid plane, then increments with restrictions from the fine grid
 //------------------------------------------------------------------------------------------------------------------------------
 void residual_and_restriction(domain_type *domain, int level_f, int phi_id, int rhs_id, int level_c, int res_id, double a, double b){
   // exchange the boundary for x in prep for Ax...
   // for 7-point stencil, only needs to be a 1-deep ghost zone & faces only
   exchange_boundary(domain,level_f,phi_id,1,0,0); 
 
   uint64_t _timeStart = CycleTime();
   int CollaborativeThreadingBoxSize = 100000; // i.e. never
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level_f].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level_f].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 
+#ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int pencil_c = domain->subdomains[box].levels[level_c].pencil;
     int  plane_c = domain->subdomains[box].levels[level_c].plane;
     int ghosts_c = domain->subdomains[box].levels[level_c].ghosts;
     int  dim_k_c = domain->subdomains[box].levels[level_c].dim.k;
     int  dim_j_c = domain->subdomains[box].levels[level_c].dim.j;
     int  dim_i_c = domain->subdomains[box].levels[level_c].dim.i;
 
for(box=0;box<domain->subdomains_per_rank;box++){
     double * __restrict__ rhs    = domain->subdomains[box].levels[level_f].grids[  rhs_id] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ alpha  = domain->subdomains[box].levels[level_f].grids[__alpha ] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ beta_i = domain->subdomains[box].levels[level_f].grids[__beta_i] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ beta_j = domain->subdomains[box].levels[level_f].grids[__beta_j] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ beta_k = domain->subdomains[box].levels[level_f].grids[__beta_k] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ res    = domain->subdomains[box].levels[level_c].grids[  res_id] + ghosts_c*(1+pencil_c+plane_c);
 
     int kk;
+#ifdef OMP
     #pragma omp parallel for private(kk) if(omp_within_a_box)
+#endif
+
     for(kk=0;kk<dim_k_f;kk+=2){
       int i,j,k; 
       // zero out the next coarse grid plane
       for(j=0;j<dim_j_c;j++){
       for(i=0;i<dim_i_c;i++){
           int ijk_c = (i) + (j)*pencil_c + (kk>>1)*plane_c;
           res[ijk_c] = 0.0;
       }}

restriction.inc

void restriction(domain_type *domain, int level_f, int id_c, int id_f){
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level_c].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level_c].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int ghosts_c = domain->subdomains[box].levels[level_c].ghosts;
     int pencil_c = domain->subdomains[box].levels[level_c].pencil;
     int  plane_c = domain->subdomains[box].levels[level_c].plane;
     int  dim_i_c = domain->subdomains[box].levels[level_c].dim.i;
     int  dim_j_c = domain->subdomains[box].levels[level_c].dim.j;
     int  dim_k_c = domain->subdomains[box].levels[level_c].dim.k;
   
     int ghosts_f = domain->subdomains[box].levels[level_f].ghosts;
     int pencil_f = domain->subdomains[box].levels[level_f].pencil;
     int  plane_f = domain->subdomains[box].levels[level_f].plane;
   
     double * __restrict__ grid_f = domain->subdomains[box].levels[level_f].grids[id_f] + ghosts_f*(1+pencil_f+plane_f);
     double * __restrict__ grid_c = domain->subdomains[box].levels[level_c].grids[id_c] + ghosts_c*(1+pencil_c+plane_c);
-  
+#ifdef OMP  
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<dim_k_c;k++){
     for(j=0;j<dim_j_c;j++){
     for(i=0;i<dim_i_c;i++){
       int ijk_c = (i   ) + (j   )*pencil_c + (k   )*plane_c;
       int ijk_f = (i<<1) + (j<<1)*pencil_f + (k<<1)*plane_f;
       grid_c[ijk_c] = ( grid_f[ijk_f                   ]+grid_f[ijk_f+1                 ] +
                         grid_f[ijk_f  +pencil_f        ]+grid_f[ijk_f+1+pencil_f        ] +
                         grid_f[ijk_f           +plane_f]+grid_f[ijk_f+1         +plane_f] +
                         grid_f[ijk_f  +pencil_f+plane_f]+grid_f[ijk_f+1+pencil_f+plane_f] ) * 0.125;
     }}}
   }
-#endif
   domain->cycles.restriction[level_f] += (uint64_t)(CycleTime()-_timeStart);
 }
 
 
 
 //------------------------------------------------------------------------------------------------------------------------------
 // restriction of face-centered data...
 //------------------------------------------------------------------------------------------------------------------------------
 void restriction_betas(domain_type * domain, int level_f, int level_c){
   uint64_t _timeStart = CycleTime();
   int CollaborativeThreadingBoxSize = 100000; // i.e. never
   #ifdef __COLLABORATIVE_THREADING
     CollaborativeThreadingBoxSize = 1 << __COLLABORATIVE_THREADING;
   #endif
   int omp_across_boxes = (domain->subdomains[0].levels[level_c].dim.i <  CollaborativeThreadingBoxSize);
   int omp_within_a_box = (domain->subdomains[0].levels[level_c].dim.i >= CollaborativeThreadingBoxSize);
   int box;
 
 #ifdef OMP
   #pragma omp parallel for private(box) if(omp_across_boxes)
+#endif
   for(box=0;box<domain->subdomains_per_rank;box++){
     int i,j,k;
     int ghosts_c = domain->subdomains[box].levels[level_c].ghosts;
     int pencil_c = domain->subdomains[box].levels[level_c].pencil;
     int  plane_c = domain->subdomains[box].levels[level_c].plane;
     int  dim_i_c = domain->subdomains[box].levels[level_c].dim.i;
     int  dim_j_c = domain->subdomains[box].levels[level_c].dim.j;
     int  dim_k_c = domain->subdomains[box].levels[level_c].dim.k;
for(i=0;i<dim_i_c;i++){
       int ijk_f = (i<<1) + (j<<1)*pencil_f + (k<<1)*plane_f;
       beta_c[ijk_c] = ( beta_f[ijk_f        ]+beta_f[ijk_f+pencil_f        ] +
                         beta_f[ijk_f+plane_f]+beta_f[ijk_f+pencil_f+plane_f] ) * 0.25;
     }}}
 
     // restrict beta_j (== face in ik)
     beta_f = domain->subdomains[box].levels[level_f].grids[__beta_j] + ghosts_f*plane_f + ghosts_f*pencil_f + ghosts_f;
     beta_c = domain->subdomains[box].levels[level_c].grids[__beta_j] + ghosts_c*plane_c + ghosts_c*pencil_c + ghosts_c;
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<dim_k_c;k++){
     for(j=0;j<dim_j_c;j++){
     for(i=0;i<dim_i_c;i++){
       int ijk_c = (i   ) + (j   )*pencil_c + (k   )*plane_c;
       int ijk_f = (i<<1) + (j<<1)*pencil_f + (k<<1)*plane_f;
       beta_c[ijk_c] = ( beta_f[ijk_f        ]+beta_f[ijk_f+1        ] +
                         beta_f[ijk_f+plane_f]+beta_f[ijk_f+1+plane_f] ) * 0.25;
     }}}
 
     // restrict beta_k (== face in ij)
     beta_f = domain->subdomains[box].levels[level_f].grids[__beta_k] + ghosts_f*plane_f + ghosts_f*pencil_f + ghosts_f;
     beta_c = domain->subdomains[box].levels[level_c].grids[__beta_k] + ghosts_c*plane_c + ghosts_c*pencil_c + ghosts_c;
+#ifdef OMP
     #pragma omp parallel for private(k,j,i) if(omp_within_a_box) collapse(2)
+#endif
     for(k=0;k<dim_k_c;k++){
     for(j=0;j<dim_j_c;j++){
     for(i=0;i<dim_i_c;i++){
       int ijk_c = (i   ) + (j   )*pencil_c + (k   )*plane_c;
       int ijk_f = (i<<1) + (j<<1)*pencil_f + (k<<1)*plane_f;
       beta_c[ijk_c] = ( beta_f[ijk_f         ]+beta_f[ijk_f+1         ] +
                         beta_f[ijk_f+pencil_f]+beta_f[ijk_f+1+pencil_f] ) * 0.25;
     }}}
   }
-#endif
   domain->cycles.restriction[level_f] += (uint64_t)(CycleTime()-_timeStart);
 }