 
   // Execution configurations for small and large pencil tiles
 
-  grid_sp[0]  = dim3(my / sPencils, mz, 1);
-  block_sp[0] = dim3(mx, sPencils, 1);
-
-  grid_lp[0]  = dim3(my / lPencils, mz, 1);
-  block_lp[0] = dim3(mx, sPencils, 1);
-
-  grid_sp[1]  = dim3(mx / sPencils, mz, 1);
-  block_sp[1] = dim3(sPencils, my, 1);
-
-  grid_lp[1]  = dim3(mx / lPencils, mz, 1);
+  grid[0][0]  = dim3(my / sPencils, mz, 1);
+  block[0][0] = dim3(mx, sPencils, 1);
+
+  grid[0][1]  = dim3(my / lPencils, mz, 1);
+  block[0][1] = dim3(mx, sPencils, 1);
+
+  grid[1][0]  = dim3(mx / sPencils, mz, 1);
+  block[1][0] = dim3(sPencils, my, 1);
+
+  grid[1][1]  = dim3(mx / lPencils, mz, 1);
   // we want to use the same number of threads as above,
   // so when we use lPencils instead of sPencils in one
   // dimension, we multiply the other by sPencils/lPencils
-  block_lp[1] = dim3(lPencils, my * sPencils / lPencils, 1);
-
-  grid_sp[2]  = dim3(mx / sPencils, my, 1);
-  block_sp[2] = dim3(sPencils, mz, 1);
-
-  grid_lp[2]  = dim3(mx / lPencils, my, 1);
-  block_lp[2] = dim3(lPencils, mz * sPencils / lPencils, 1);
+  block[1][1] = dim3(lPencils, my * sPencils / lPencils, 1);
+
+  grid[2][0]  = dim3(mx / sPencils, my, 1);
+  block[2][0] = dim3(sPencils, mz, 1);
+
+  grid[2][1]  = dim3(mx / lPencils, my, 1);
+  block[2][1] = dim3(lPencils, mz * sPencils / lPencils, 1);
 }
 
 void initInput(float *f, int dim)

                              sPencils, mz,
                              lPencils, mz };

  float milliseconds;
  cudaEvent_t startEvent, stopEvent;
  checkCuda( cudaEventCreate(&startEvent) );
  checkCuda( cudaEventCreate(&stopEvent) );

  double error, maxError;

  printf("%c derivatives\n\n", (char)(0x58 + dimension));

  float f[mx*my*mz];
  float df[mx*my*mz];
  float sol[mx*my*mz];
  
  initInput(f, dimension);
  initSol(sol, dimension);

  const int nReps = 20;

  // device arrays
  int bytes = mx*my*mz * sizeof(float);

  checkCuda( cudaMalloc((void**)&d_f, bytes) );
  checkCuda( cudaMalloc((void**)&d_df, bytes) );

  const int nReps = 20;
  float milliseconds;
  cudaEvent_t startEvent, stopEvent;
  checkCuda( cudaEventCreate(&startEvent) );
  checkCuda( cudaEventCreate(&stopEvent) );

  double error, maxError;

  printf("%c derivatives\n\n", (char)(0x58 + dimension));

  for (int fp = 0; fp < 2; fp++) { 
    checkCuda( cudaMemcpy(d_f, f, bytes, cudaMemcpyHostToDevice) );  
    checkCuda( cudaMemset(d_df, 0, bytes) );
    
    fpDeriv[fp]<<<grid[dimension][fp],block[dimension][fp]>>>(d_f, d_df); // warm up
    checkCuda( cudaEventRecord(startEvent, 0) );
    for (int i = 0; i < nReps; i++)
       fpDeriv[fp]<<<grid[dimension][fp],block[dimension][fp]>>>(d_f, d_df);
    
    checkCuda( cudaEventRecord(stopEvent, 0) );
    checkCuda( cudaEventSynchronize(stopEvent) );
