Accelerated Ray Tracing （十）

https://developer.nvidia.com/blog/finite-difference-methods-cuda-cc-part-1/

有限差分方法

通过使用共享内存来合并全局内存访问，在有效带宽上获得了大约一个数量级的改进。本片如何使用共享内存以有限差分代码增强数据重用。除shared memory,，我们还将讨论constant memory，这是一种缓存的read-only memory，已优化用于跨块（或扭曲）中的线程进行统一访问。

Problem Statement: 3D Finite Difference

为简单起见，我们假定周期边界条件并且仅考虑一阶导数，尽管扩展代码以使用其他类型的边界条件计算高阶导数很简单。

有限差分法实质上是在相邻点使用函数值的加权求和来近似特定点的导数。对于在x方向上具有均匀间距Δx的（2N + 1）点模版，以下方程式给出了x导数的中心有限差分方案。其他方向的方程式相似。

系数C_i通常是从泰勒级数展开中生成的，可以选择以获得具有所需特性（例如精度）的方案，并且可以在偏微分方程，色散和耗散的情况下使用。对于诸如上述类型的显式有限差分方案，较大的模板通常具有较高的精度。在本文中，我们使用具有八阶精度的九点模板。我们还选择了对称的模具，如以下等式所示。

在这里，我们使用网格索引i，j，k而不是物理坐标x，y，z在计算网格上指定函数的值。这里的系数是ax = 4 /（5 ∆x），bx = −1 /（5 ∆x），cx = 4 /（105 ∆x）和dx = − 1 /（280 ∆x），这是典型的八阶方案。对于y和z方向的导数，上式中的索引偏移量仅适用于j和k索引，并且系数相同，只是使用∆y和∆z代替∆x。

因为我们计算了64^3周期网格上每个点的导数的近似值，所以每个点的f值被使用了八次：上面表达式中的每个右侧项一次。在设计派生内核时，我们的目标是通过将数据块存储在共享内存中来利用这种数据重用，以确保我们尽可能少地从全局内存中获取f的值。

Data Reuse and Shared Memory

我们采用tiling方法，其中每个线程块将多维网格中的数据块加载到共享内存中，以便该块中的每个线程都可以根据需要访问共享内存块的所有元素。我们如何选择最佳的tile形状和尺寸？需要进行一些实验，但是有限差分模具的特征和网格尺寸提供了方向。

选择用于模板计算的图块形状时，图块通常重叠模板大小的一半，如下图左图所示。

 

 在这里，为了计算16×16 tile（黄色）中的导数，每个线程块必须不仅加载该tile，而且还要加载两个附加的4×16部分（橙色）中的f值。总体而言，橙色部分中的f值被加载两次，一次由计算该位置导数的线程块加载，一次由每个相邻线程块加载。结果，从全局内存中两次加载16×16中的8×16值或一半的值。此外，在计算能力为2.0或更高的设备上合并时，使用16×16 tile将是次优的，因为在此类设备上进行完美的合并需要每个负载访问全局内存中32个连续元素内的数据。

如上图右图所示，针对我们的问题和数据大小的瓦片（和线程块）大小的更好选择是64×4tile。当为64^3网格上的一维模板计算x导数时，此图块完全避免重叠，因为该图块包含导数方向上的所有点。最小的图块将只有一根pencil，或一个方向上所有点的一维数组。这将仅对应于64个线程的线程块，因此，从occupancy 的角度来看，每个tile使用多支pencils 是有益的。我们将铅笔的数量参数化以进行实验。除了仅将f的每个值加载一次之外，每个线程warp都使用此tile从全局内存加载连续数据，因此可以实现对全局内存的完美合并访问。

X-Derivative Implementation

在这里，mx，my和mz是设置为64的数组维度，而sPencils是4，这是用于制作共享内存图块的铅笔数。

Y and Z Derivatives

我们可以轻松地修改x导数代码以在其他方向操作。在x导数中，每个线程块都会在64×sPencils元素的x，y切片中计算导数。对于y导数，我们可以使用一个线程块来计算x，y中sPencils×64个元素的tile上的导数，如下图左图所示。

#include <stdio.h>
#include <assert.h>

// Convenience function for checking CUDA runtime API results
// can be wrapped around any runtime API call. No-op in release builds.
inline
cudaError_t checkCuda(cudaError_t result)
{
#if defined(DEBUG) || defined(_DEBUG)
    if (result != cudaSuccess) {
        fprintf(stderr, "CUDA Runtime Error: %s\n", cudaGetErrorString(result));
        assert(result == cudaSuccess);
    }
#endif
    return result;
}

float fx = 1.0f, fy = 1.0f, fz = 1.0f;
const int mx = 64, my = 64, mz = 64;

// shared memory tiles will be m*-by-*Pencils
// sPencils is used when each thread calculates the derivative at one point
// lPencils is used for coalescing in y and z where each thread has to 
//     calculate the derivative at mutiple points
const int sPencils = 4;  // small # pencils
const int lPencils = 32; // large # pencils

dim3 grid[3][2], block[3][2];

/*
* 常量存储器是一种专用的存储空间，只能从设备代码读取，但可以由主机读取和写入。恒定内存驻留在设备DRAM中，并缓存在芯片上。
* 常量存储器高速缓存只有一个读取端口，但可以通过warp从该端口广播数据。
* 这意味着，当warp中的所有线程都读取相同的地址时，对常量的访问是有效的，但是当warp中的线程读取了不同的地址时，则将读取序列化。
* 常量存储器非常适合系数和其他在线程间统一使用的数据，例如系数c_ax，c_bx等。
* 在CUDA C / C ++中，常量数据必须在全局范围内声明，并且可以（仅）从设备代码读取，也可以由主机代码读取或写入。
* 设备代码中使用常量内存的方式与使用任何CUDA C变量或数组/指针的方式相同，但必须使用cudaMemcpyToSymbol或其变体之一从主机代码中对其进行初始化。
* 在我们的有限差分代码中，我们有以下常量声明，它们使用__constant__声明说明符。
*/
// stencil coefficients
__constant__ float c_ax, c_bx, c_cx, c_dx;
__constant__ float c_ay, c_by, c_cy, c_dy;
__constant__ float c_az, c_bz, c_cz, c_dz;

// host routine to set constant data
void setDerivativeParameters()
{
    // check to make sure dimensions are integral multiples of sPencils
    if ((mx % sPencils != 0) || (my % sPencils != 0) || (mz % sPencils != 0)) {
        printf("'mx', 'my', and 'mz' must be integral multiples of sPencils\n");
        exit(1);
    }

    if ((mx % lPencils != 0) || (my % lPencils != 0)) {
        printf("'mx' and 'my' must be multiples of lPencils\n");
        exit(1);
    }

    // stencil weights (for unit length problem)
    float dsinv = mx - 1.f;

    float ax = 4.f / 5.f * dsinv;//我们使用以下代码初始化x系数，并且类似地初始化y和z系数。
    float bx = -1.f / 5.f * dsinv;
    float cx = 4.f / 105.f * dsinv;
    float dx = -1.f / 280.f * dsinv;
    checkCuda(cudaMemcpyToSymbol(c_ax, &ax, sizeof(float), 0, cudaMemcpyHostToDevice));
    checkCuda(cudaMemcpyToSymbol(c_bx, &bx, sizeof(float), 0, cudaMemcpyHostToDevice));
    checkCuda(cudaMemcpyToSymbol(c_cx, &cx, sizeof(float), 0, cudaMemcpyHostToDevice));
    checkCuda(cudaMemcpyToSymbol(c_dx, &dx, sizeof(float), 0, cudaMemcpyHostToDevice));

    dsinv = my - 1.f;

    float ay = 4.f / 5.f * dsinv;
    float by = -1.f / 5.f * dsinv;
    float cy = 4.f / 105.f * dsinv;
    float dy = -1.f / 280.f * dsinv;
    checkCuda(cudaMemcpyToSymbol(c_ay, &ay, sizeof(float), 0, cudaMemcpyHostToDevice));
    checkCuda(cudaMemcpyToSymbol(c_by, &by, sizeof(float), 0, cudaMemcpyHostToDevice));
    checkCuda(cudaMemcpyToSymbol(c_cy, &cy, sizeof(float), 0, cudaMemcpyHostToDevice));
    checkCuda(cudaMemcpyToSymbol(c_dy, &dy, sizeof(float), 0, cudaMemcpyHostToDevice));

    dsinv = mz - 1.f;

    float az = 4.f / 5.f * dsinv;
    float bz = -1.f / 5.f * dsinv;
    float cz = 4.f / 105.f * dsinv;
    float dz = -1.f / 280.f * dsinv;
    checkCuda(cudaMemcpyToSymbol(c_az, &az, sizeof(float), 0, cudaMemcpyHostToDevice));
    checkCuda(cudaMemcpyToSymbol(c_bz, &bz, sizeof(float), 0, cudaMemcpyHostToDevice));
    checkCuda(cudaMemcpyToSymbol(c_cz, &cz, sizeof(float), 0, cudaMemcpyHostToDevice));
    checkCuda(cudaMemcpyToSymbol(c_dz, &dz, sizeof(float), 0, cudaMemcpyHostToDevice));

    // Execution configurations for small and large pencil tiles

    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[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)
{
    const float twopi = 8.f * (float)atan(1.0);

    for (int k = 0; k < mz; k++) {
        for (int j = 0; j < my; j++) {
            for (int i = 0; i < mx; i++) {
                switch (dim) {
                case 0:
                    f[k * mx * my + j * mx + i] = cos(fx * twopi * (i - 1.f) / (mx - 1.f));
                    break;
                case 1:
                    f[k * mx * my + j * mx + i] = cos(fy * twopi * (j - 1.f) / (my - 1.f));
                    break;
                case 2:
                    f[k * mx * my + j * mx + i] = cos(fz * twopi * (k - 1.f) / (mz - 1.f));
                    break;
                }
            }
        }
    }
}

void initSol(float* sol, int dim)
{
    const float twopi = 8.f * (float)atan(1.0);

    for (int k = 0; k < mz; k++) {
        for (int j = 0; j < my; j++) {
            for (int i = 0; i < mx; i++) {
                switch (dim) {
                case 0:
                    sol[k * mx * my + j * mx + i] = -fx * twopi * sin(fx * twopi * (i - 1.f) / (mx - 1.f));
                    break;
                case 1:
                    sol[k * mx * my + j * mx + i] = -fy * twopi * sin(fy * twopi * (j - 1.f) / (my - 1.f));
                    break;
                case 2:
                    sol[k * mx * my + j * mx + i] = -fz * twopi * sin(fz * twopi * (k - 1.f) / (mz - 1.f));
                    break;
                }
            }
        }
    }
}

void checkResults(double& error, double& maxError, float* sol, float* df)
{
    // error = sqrt(sum((sol-df)**2)/(mx*my*mz))
    // maxError = maxval(abs(sol-df))
    maxError = 0;
    error = 0;
    for (int k = 0; k < mz; k++) {
        for (int j = 0; j < my; j++) {
            for (int i = 0; i < mx; i++) {
                float s = sol[k * mx * my + j * mx + i];
                float f = df[k * mx * my + j * mx + i];
                //printf("%d %d %d: %f %f\n", i, j, k, s, f);
                error += (s - f) * (s - f);
                if (fabs(s - f) > maxError) maxError = fabs(s - f);
            }
        }
    }
    error = sqrt(error / (mx * my * mz));
}


// -------------
// x derivatives
// -------------

__global__ void derivative_x(float* f, float* df)//x导数内核的实现如下。
{
    __shared__ float s_f[sPencils][mx + 8]; // 4-wide halo;在x维度的每一端声明共享内存tile，并在其中填充4个元素，以适应在端点处计算导数所需的周期性图像。

    //该内核以64×sPencils线程块启动，该线程在64×sPencils元素的x×y切片上计算导数，因此每个线程在一个点计算导数。
    int i = threadIdx.x;
    int j = blockIdx.x * blockDim.y + threadIdx.y;
    int k = blockIdx.y;//索引i，j和k对应于64^3网格中的坐标。
    int si = i + 4;       // local i for shared memory access + halo offset;索引si和sj是共享内存tile的本地（行，列）坐标。
    int sj = threadIdx.y; // local j for shared memory access

    //我们计算线性输入数组的索引，然后从全局内存读取到共享内存磁贴。
    int globalIdx = k * mx * my + j * mx + i;//在这里，mx，my和mz是设置为64的数组维度，而sPencils是4，这是用于制作共享内存图块的Pencils数。

    s_f[sj][si] = f[globalIdx];

    __syncthreads();
    
    //来自全局内存的这些读取已完美合并。通过以下代码将数据复制到共享内存中，以在x方向上填充周期性图像。每个元素只能由前一条语句从全局存储器读取一次。
    // fill in periodic images in shared memory array 
    if (i < 4) {
        s_f[sj][si - 4] = s_f[sj][si + mx - 5];
        s_f[sj][si + mx] = s_f[sj][si + 1];
    }
    /*
    * 请注意，在此操作中，我们从共享内存中读取，而不是从全局内存中读取。
    * 由于线程从另一个线程写入的共享内存中读取数据，因此在写入周期图像之前，我们需要__syncthreads（）调用。
    * 同样，在写入周期图像之后，我们需要一个__syncthreads（）调用，因为它们是由不同的线程访问的。
    * 在第二个__syncthreads（）调用之后，使用以下代码计算我们的有限差分近似值。
    */

    __syncthreads();

    df[globalIdx] =
        (c_ax * (s_f[sj][si + 1] - s_f[sj][si - 1])//常数c_ax，c_bx，c_cx和c_dx在此内核的外部常量存储器中声明。
            + c_bx * (s_f[sj][si + 2] - s_f[sj][si - 2])
            + c_cx * (s_f[sj][si + 3] - s_f[sj][si - 3])
            + c_dx * (s_f[sj][si + 4] - s_f[sj][si - 4]));
}


// this version uses a 64x32 shared memory tile, 
// still with 64*sPencils threads

__global__ void derivative_x_lPencils(float* f, float* df)
{/*
 在准确性方面，我们获得了x导数的相同最大误差，但对于基本相同的函数，则获得了不同的（实际上并不差）RMS误差。
 这种差异是由于在主机上完成累加的顺序所致，简单地交换主机代码错误计算中的循环顺序将产生相同的结果。
 恢复完美合并的一种方法是扩展共享内存块中的Pencil数量。对于Tesla K20（Kepler），这将需要32支铅笔。
 使用这种方法将需要9216字节的共享内存块，这不是问题。但是，如果将线程与要计算导数的元素进行一对一映射，则需要2048个线程的线程块。
 K20的限制是每个线程块1024个线程，每个多处理器2048个线程。我们可以通过让每个线程计算多个点的导数来解决这些限制。
 如果我们使用一个线程的32×8×1块,每个线程计算导数在8点,而不是一个线程块4×64×1,并且每个线程计算导数在只有一个点,我们推出一个内核与相同数量的块和线程/块,但恢复完美的合并。
 */

    __shared__ float s_f[lPencils][mx + 8]; // 4-wide halo

    int i = threadIdx.x;
    int jBase = blockIdx.x * lPencils;
    int k = blockIdx.y;
    int si = i + 4; // local i for shared memory access + halo offset

    for (int sj = threadIdx.y; sj < lPencils; sj += blockDim.y) {
        int globalIdx = k * mx * my + (jBase + sj) * mx + i;
        s_f[sj][si] = f[globalIdx];
    }

    __syncthreads();

    // fill in periodic images in shared memory array 
    if (i < 4) {
        for (int sj = threadIdx.y; sj < lPencils; sj += blockDim.y) {
            s_f[sj][si - 4] = s_f[sj][si + mx - 5];
            s_f[sj][si + mx] = s_f[sj][si + 1];
        }
    }

    __syncthreads();

    for (int sj = threadIdx.y; sj < lPencils; sj += blockDim.y) {
        int globalIdx = k * mx * my + (jBase + sj) * mx + i;
        df[globalIdx] =
            (c_ax * (s_f[sj][si + 1] - s_f[sj][si - 1])
                + c_bx * (s_f[sj][si + 2] - s_f[sj][si - 2])
                + c_cx * (s_f[sj][si + 3] - s_f[sj][si - 3])
                + c_dx * (s_f[sj][si + 4] - s_f[sj][si - 4]));
    }
}
/*
* 这里我们让l铅笔等于32。这段代码与上面的y导数代码类似，除了变量j是一个循环索引，而不是计算一次。
当每个线程块使用大量共享内存时，我们应该检查共享内存使用是否限制了驻留在多处理器上的线程数量，以及在多大程度上限制了这些线程的数量。
每个线程块需要9216字节的共享内存，并且每个多处理器的共享内存限制为48 KB，这意味着最多可以在一个多处理器上同时驻留5个线程块。3
(使用——ptxas-options=-v进行编译表明每个线程使用10个寄存器，这不会限制每个多处理器的线程块数量。)
对于一个包含32个8个线程的线程块，这个内核的多处理器上可以驻留的线程数是1280，而K20上的最大线程数是2048。
占用率为1280/2048或0.63，这通常不是问题，特别是因为我们通过j上的循环使用了8倍指令级并行性。
*/
// -------------
// y derivatives
// -------------

__global__ void derivative_y(float* f, float* df)
{//同样，对于z导数，线程块可以计算sPencils x 64个元素的x，z切片中的导数。下面的内核显示了使用这种方法的y导数内核。
    __shared__ float s_f[my + 8][sPencils];//通过使用共享内存tile，我们可以确保全局内存中的每个元素仅被读取一次。
    /*
    * 这种方法的缺点是，对于这些tile，使用sPencils或x中的4个点，我们将不再具有完美的合并效果。性能结果证明了这一点。我们的性能大约是x派生内核的一半。
    */
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = threadIdx.y;
    int k = blockIdx.y;
    int si = threadIdx.x;
    int sj = j + 4;

    int globalIdx = k * mx * my + j * mx + i;

    s_f[sj][si] = f[globalIdx];

    __syncthreads();

    if (j < 4) {
        s_f[sj - 4][si] = s_f[sj + my - 5][si];
        s_f[sj + my][si] = s_f[sj + 1][si];
    }

    __syncthreads();

    df[globalIdx] =
        (c_ay * (s_f[sj + 1][si] - s_f[sj - 1][si])
            + c_by * (s_f[sj + 2][si] - s_f[sj - 2][si])
            + c_cy * (s_f[sj + 3][si] - s_f[sj - 3][si])
            + c_dy * (s_f[sj + 4][si] - s_f[sj - 4][si]));
}

// y derivative using a tile of 32x64,
// launch with thread block of 32x8
__global__ void derivative_y_lPencils(float* f, float* df)
{
    __shared__ float s_f[my + 8][lPencils];

    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int k = blockIdx.y;
    int si = threadIdx.x;

    for (int j = threadIdx.y; j < my; j += blockDim.y) {
        int globalIdx = k * mx * my + j * mx + i;
        int sj = j + 4;
        s_f[sj][si] = f[globalIdx];
    }

    __syncthreads();

    int sj = threadIdx.y + 4;
    if (sj < 8) {
        s_f[sj - 4][si] = s_f[sj + my - 5][si];
        s_f[sj + my][si] = s_f[sj + 1][si];
    }

    __syncthreads();

    for (int j = threadIdx.y; j < my; j += blockDim.y) {
        int globalIdx = k * mx * my + j * mx + i;
        int sj = j + 4;
        df[globalIdx] =
            (c_ay * (s_f[sj + 1][si] - s_f[sj - 1][si])
                + c_by * (s_f[sj + 2][si] - s_f[sj - 2][si])
                + c_cy * (s_f[sj + 3][si] - s_f[sj - 3][si])
                + c_dy * (s_f[sj + 4][si] - s_f[sj - 4][si]));
    }
}


// ------------
// z derivative
// ------------

__global__ void derivative_z(float* f, float* df)
{
    __shared__ float s_f[mz + 8][sPencils];

    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = blockIdx.y;
    int k = threadIdx.y;
    int si = threadIdx.x;
    int sk = k + 4; // halo offset

    int globalIdx = k * mx * my + j * mx + i;

    s_f[sk][si] = f[globalIdx];

    __syncthreads();

    if (k < 4) {
        s_f[sk - 4][si] = s_f[sk + mz - 5][si];
        s_f[sk + mz][si] = s_f[sk + 1][si];
    }

    __syncthreads();

    df[globalIdx] =
        (c_az * (s_f[sk + 1][si] - s_f[sk - 1][si])
            + c_bz * (s_f[sk + 2][si] - s_f[sk - 2][si])
            + c_cz * (s_f[sk + 3][si] - s_f[sk - 3][si])
            + c_dz * (s_f[sk + 4][si] - s_f[sk - 4][si]));
}

__global__ void derivative_z_lPencils(float* f, float* df)
{
    __shared__ float s_f[mz + 8][lPencils];

    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = blockIdx.y;
    int si = threadIdx.x;

    for (int k = threadIdx.y; k < mz; k += blockDim.y) {
        int globalIdx = k * mx * my + j * mx + i;
        int sk = k + 4;
        s_f[sk][si] = f[globalIdx];
    }

    __syncthreads();

    int k = threadIdx.y + 4;
    if (k < 8) {
        s_f[k - 4][si] = s_f[k + mz - 5][si];
        s_f[k + mz][si] = s_f[k + 1][si];
    }

    __syncthreads();

    for (int k = threadIdx.y; k < mz; k += blockDim.y) {
        int globalIdx = k * mx * my + j * mx + i;
        int sk = k + 4;
        df[globalIdx] =
            (c_az * (s_f[sk + 1][si] - s_f[sk - 1][si])
                + c_bz * (s_f[sk + 2][si] - s_f[sk - 2][si])
                + c_cz * (s_f[sk + 3][si] - s_f[sk - 3][si])
                + c_dz * (s_f[sk + 4][si] - s_f[sk - 4][si]));
    }
}

// Run the kernels for a given dimension. One for sPencils, one for lPencils
void runTest(int dimension)
{
    void (*fpDeriv[2])(float*, float*);

    switch (dimension) {
    case 0:
        fpDeriv[0] = derivative_x;
        fpDeriv[1] = derivative_x_lPencils;
        break;
    case 1:
        fpDeriv[0] = derivative_y;
        fpDeriv[1] = derivative_y_lPencils;
        break;
    case 2:
        fpDeriv[0] = derivative_z;
        fpDeriv[1] = derivative_z_lPencils;
        break;
    }

    int sharedDims[3][2][2] = { mx, sPencils,
                                mx, lPencils,
                                sPencils, my,
                                lPencils, my,
                                sPencils, mz,
                                lPencils, mz };

    float* f = new float[mx * my * mz];
    float* df = new float[mx * my * mz];
    float* sol = new float[mx * my * mz];

    initInput(f, dimension);
    initSol(sol, dimension);

    // device arrays
    int bytes = mx * my * mz * sizeof(float);
    float* d_f, * d_df;
    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));
        checkCuda(cudaEventElapsedTime(&milliseconds, startEvent, stopEvent));

        checkCuda(cudaMemcpy(df, d_df, bytes, cudaMemcpyDeviceToHost));

        checkResults(error, maxError, sol, df);

        printf("  Using shared memory tile of %d x %d\n",
            sharedDims[dimension][fp][0], sharedDims[dimension][fp][1]);
        printf("   RMS error: %e\n", error);
        printf("   MAX error: %e\n", maxError);
        printf("   Average time (ms): %f\n", milliseconds / nReps);
        printf("   Average Bandwidth (GB/s): %f\n\n",
            2.f * 1e-6 * mx * my * mz * nReps * sizeof(float) / milliseconds);
    }

    checkCuda(cudaEventDestroy(startEvent));
    checkCuda(cudaEventDestroy(stopEvent));

    checkCuda(cudaFree(d_f));
    checkCuda(cudaFree(d_df));

    delete[] f;
    delete[] df;
    delete[] sol;
}


// This the main host code for the finite difference 
// example.  The kernels are contained in the derivative_m module

int main(void)
{
    // Print device and precision
    cudaDeviceProp prop;
    checkCuda(cudaGetDeviceProperties(&prop, 0));
    printf("\nDevice Name: %s\n", prop.name);
    printf("Compute Capability: %d.%d\n\n", prop.major, prop.minor);

    setDerivativeParameters(); // initialize 

    runTest(0); // x derivative
    runTest(1); // y derivative
    runTest(2); // z derivative

    return 0;
}