CUDA并行计算加速体积计算

Question

• 设备：Tesla C2050
• 操作系统：Windows 7 企业版
• IDE：VS 2010
• CUDA：5.0（最新）

第一次在这里提问。我在 CUDA 程序中遇到了一些问题。

我有数百万个四面体，一点在 (0,0,0)，所以我可以使用公式：

得到四面体的体积。

所以，这里是代码：

struct Triangle
{
    double x1;
    double y1;
    double z1;
    double x2;
    double y2;
    double z2;
    double x3;
    double y3;
    double z3;
};

还有 CUDA 代码：

__global__ void getResult(double *d_volume ,Triangle *d_triangles, Origin *d_point)
{
    extern __shared__ Triangle s_data[];
    int tid = threadIdx.x;
    int i =  blockDim.x * blockIdx.x + threadIdx.x;
    s_data[tid] = d_triangles[i];
    __syncthreads();
    d_volume[i] =s_data[tid].x1 * s_data[tid].y2 * s_data[tid].z3 + \
                s_data[tid].y1 * s_data[tid].z2 * s_data[tid].x3 + \
                s_data[tid].x2 * s_data[tid].y3 * s_data[tid].z1 - \
                s_data[tid].x3 * s_data[tid].y2 * s_data[tid].z1 - \
                s_data[tid].x2 * s_data[tid].y1 * s_data[tid].z3 - \
                s_data[tid].y3 * s_data[tid].z2 * s_data[tid].x1;
}

我从其他函数中获得了数百万个四面体作为数组。

// Host
Triangle *h_triangles = triangles;
double *h_volume;
// Device
Triangle *d_triangles;
double *d_volume;

// define grid and block size
int numThreadsPerBlock = numTriangles;
int numBlocks = numTrianges / 512;

// Shard memory size
int sharedMemSize = numThreadsPerBlock * sizeof(Triangle);

// allocate host and device memory
size_t memSize_triangles = numBlocks * numThreadsPerBlock * sizeof(Triangle);
size_t memSize_volume = numBlocks * numThreadsPerBlock * sizeof(double);

cudaMalloc( (void **) &d_triangles, memSize_triangles );
cudaMalloc( (void **) &d_volume, memSize_volume );

// Copy host array to device array
cudaMemcpy( d_triangles, h_triangles, memSize_triangles, cudaMemcpyHostToDevice );
cudaMemcpy( d_point, h_point, memSize_point, cudaMemcpyHostToDevice );

// launch kernel
dim3 dimGrid(numBlocks);
dim3 dimBlock(numThreadsPerBlock);

getResult<<< dimGrid, dimBlock, sharedMemSize >>>( d_volume, d_triangles);

// block until the device has completed
cudaThreadSynchronize();

// device to host copy
cudaMemcpy( h_volume, d_volume, memSize_volume, cudaMemcpyDeviceToHost );

// free device memory
cudaFree(d_triangles);
cudaFree(d_volume);

// free host memory
free(h_triangles); 
free(h_volume);

到目前为止，一切正常。但是我花费了比我想象的更多的时间来获得这本书。 我的设备是 Tesla C2050（515Gflops），比我的 CPU（单核，20.25Gflops）快 20 倍。 但只提速10倍左右（不包括设备和主机之间复制内存的时间。）

我想知道如何使它比 CPU 代码快 20 倍（for 循环获取音量。）。

谢谢！

PS：也许 cudaMallocPitch() 会帮助我，但是三角形不是矩阵，我不能使用 cudaMemcpy2D() 而不是 cudaMemcpy() 来复制内存。谁能帮我解决这个问题？

Answer 1

与 CPU 相比，GPU 的峰值性能通常更难获得。原因之一是许多内核受带宽限制而不是计算限制。

因为你的内核的计算复杂度是 O(n)。您可能应该使用带宽指标来计算理论峰值性能，如下所示

1024*1024*64 * sizeof(double) * (9  +   1)     / (144e9    *    8/9)     = 42 ms
#tetrahedron                     #input #output   peak mem bw   ECC cost

另一方面，您的内核可以进一步优化。

• 请谨慎选择 blockDim/gridDim，错误的数字有时会导致 20% 的性能损失。
• 您可以为每个线程计算 多个 卷，而不是每个线程计算 一个 卷，这将减少线程启动开销。
• 由于你不在线程之间共享数据，__syncthreads() 或许可以消除。
• 由于未合并的内存访问，结构数组 (AoS) 通常比 GPU 上的数组结构 (SoA) 慢。你也可以尝试改变你的数据结构。

更新

获得了具有大型 L1 缓存设置和最佳 blockDim/gridDim 选择的新内核。它快 15%。这是代码和配置文件结果。我的设备是 M2090。

#include <stdlib.h>
#include <thrust/transform.h>
#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <iterator>
#include <thrust/inner_product.h>

using namespace thrust::placeholders;

struct Triangle
{
    double x1;
    double y1;
    double z1;
    double x2;
    double y2;
    double z2;
    double x3;
    double y3;
    double z3;
};

__global__ void getResultNoSMem(double *d_volume, Triangle *d_triangles)
{
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    d_volume[i] = d_triangles[i].x1 * d_triangles[i].y2 * d_triangles[i].z3 +
            d_triangles[i].y1 * d_triangles[i].z2 * d_triangles[i].x3 +
            d_triangles[i].x2 * d_triangles[i].y3 * d_triangles[i].z1 -
            d_triangles[i].x3 * d_triangles[i].y2 * d_triangles[i].z1 -
            d_triangles[i].x2 * d_triangles[i].y1 * d_triangles[i].z3 -
            d_triangles[i].y3 * d_triangles[i].z2 * d_triangles[i].x1;
}

__global__ void getResult(double *d_volume, Triangle *d_triangles)
{
    extern __shared__ Triangle s_data[];
    int tid = threadIdx.x;
    int i = blockDim.x * blockIdx.x + threadIdx.x;
    s_data[tid] = d_triangles[i];
//  __syncthreads();
    d_volume[i] = s_data[tid].x1 * s_data[tid].y2 * s_data[tid].z3 +
            s_data[tid].y1 * s_data[tid].z2 * s_data[tid].x3 +
            s_data[tid].x2 * s_data[tid].y3 * s_data[tid].z1 -
            s_data[tid].x3 * s_data[tid].y2 * s_data[tid].z1 -
            s_data[tid].x2 * s_data[tid].y1 * s_data[tid].z3 -
            s_data[tid].y3 * s_data[tid].z2 * s_data[tid].x1;
}

__global__ void getResultOpt(double *d_volume, Triangle *d_triangles, int len)
{
    const int gridSize = blockDim.x * gridDim.x;
    int i = blockDim.x * blockIdx.x + threadIdx.x;

    while (i < len)
    {
        d_volume[i] = d_triangles[i].x1 * d_triangles[i].y2 * d_triangles[i].z3 +
                d_triangles[i].y1 * d_triangles[i].z2 * d_triangles[i].x3 +
                d_triangles[i].x2 * d_triangles[i].y3 * d_triangles[i].z1 -
                d_triangles[i].x3 * d_triangles[i].y2 * d_triangles[i].z1 -
                d_triangles[i].x2 * d_triangles[i].y1 * d_triangles[i].z3 -
                d_triangles[i].y3 * d_triangles[i].z2 * d_triangles[i].x1;
        i += gridSize;
    }
}

int main(void)
{
    const int m = 1024 * 1024;
    thrust::host_vector<Triangle> data(m);
    for (int i = 0; i < m; i++)
    {
        data[i].x1 = (double) rand() / RAND_MAX;
        data[i].y1 = (double) rand() / RAND_MAX;
        data[i].z1 = (double) rand() / RAND_MAX;
        data[i].x2 = (double) rand() / RAND_MAX;
        data[i].y2 = (double) rand() / RAND_MAX;
        data[i].z2 = (double) rand() / RAND_MAX;
        data[i].x3 = (double) rand() / RAND_MAX;
        data[i].y3 = (double) rand() / RAND_MAX;
        data[i].z3 = (double) rand() / RAND_MAX;
    }

    thrust::device_vector<Triangle> triangles = data;
    thrust::device_vector<double> volume(m);
    thrust::device_vector<double> volumeOpt(m);

    Triangle* dTriangles = thrust::raw_pointer_cast(&triangles[0]);
    double* dVolume = thrust::raw_pointer_cast(&volume[0]);
    double* dVolumeOpt = thrust::raw_pointer_cast(&volumeOpt[0]);

    int g;
    int b;

    int threadUpperLimit = 48 * 1024 / sizeof(Triangle);

    //for (b = 32; b <= 1024; b += 32)
    {
        b = 64;
        int gridDim = (m + b - 1) / b;
        getResultNoSMem<<<gridDim, b, 0, 0>>>(dVolume, dTriangles);
    }

    //  for (b = 32; b <= threadUpperLimit; b += 32)
    {
        b = 64;
        int gridDim = (m + b - 1) / b;
        getResult<<<gridDim, b, b * sizeof(Triangle), 0>>>(dVolume, dTriangles);
    }

    //for (g = 32; g <= 512; g += 32)
    //  for (b = 32; b <= 1024; b += 32)
    {
        b = 64;
        g = 64;
        getResultOpt<<<g, b, 0, 0>>>(dVolumeOpt, dTriangles, m);
    }

    //for (g = 32; g <= 512; g += 32)
    //  for (b = 32; b <= 1024; b += 32)
    {
        b = 64;
        g = 512;
        cudaFuncSetCacheConfig(getResultOpt, cudaFuncCachePreferL1);
        getResultOpt<<<g, b, 0, 0>>>(dVolumeOpt, dTriangles, m);
    }

    thrust::device_vector<double> X = volume;
    thrust::device_vector<double> Y = volumeOpt;
    thrust::transform(X.begin(), X.end(), Y.begin(), X.begin(), _1 - _2);
    double result = thrust::inner_product(X.begin(), X.end(), X.begin(), 0.0);

    std::cout << "difference: " << result << std::endl;

    return 0;
}

Answer 2

正如 Eric 的回答所暗示的，您的内核需要 9 个 64 位加载和每个线程 64 位存储，但每个线程只执行 17 次 FLOP。这可能意味着您的代码受内存带宽限制，而不是计算限制，并且您不应该期望这种类型的代码能够达到峰值 FLOP/s 吞吐量。

因此，优化性能的关键很可能在于内存带宽优化。目前，您的内核有几个明显的问题，我在评论中谈到了其中一个问题。你真的不需要内核中的共享内存 - 它比寄存器慢，并且使用它没有内存带宽改进。使用__syncthreads() 也会给内核增加不必要的延迟。你的代码可以写成这样：

__global__ void getresult2(double *d_volume, Triangle *d_triangles)
{
    int i =  blockDim.x * blockIdx.x + threadIdx.x;
    Triangle t = d_triangles[i];
    d_volume[i] = t.x1 * t.y2 * t.z3 + 
                  t.y1 * t.z2 * t.x3 +
                  t.x2 * t.y3 * t.z1 -
                  t.x3 * t.y2 * t.z1 -
                  t.x2 * t.y1 * t.z3 - 
                  t.y3 * t.z2 * t.x1;
}

[免责声明，从未编译或运行，使用风险自负]

我希望它的性能比共享内存版本更好。

第二个问题是内存合并问题。您拥有的结构相当大，并且每个加载Triangle 完整实例的线程对内存合并或缓存重用都不是很友好。您可以尝试使用共享内存来提高内存加载性能并通过执行以下操作来合并写入：

__global__ 
void __launch_bounds__(288, 3)getresult3(double *d_volume, Triangle *d_triangles)
{
    __shared__ double s_data[9*32];

    int i =  blockDim.x * blockIdx.x + threadIdx.x;
    double * t_data = reinterpret_cast<double *>(d_triangles);
    s_data[threadIdx.x] = t_data[i];
    __syncthreads();

    if (threadIdx.x < 32) {
        Triangle * t = reinterpret_cast<Triangle *>(&s_data[9*threadIdx.x]);
        d_volume[i] = t->x1 * t->y2 * t->z3 + 
            t->y1 * t->z2 * t->x3 +
            t->x2 * t->y3 * t->z1 -
            t->x3 * t->y2 * t->z1 -
            t->x2 * t->y1 * t->z3 - 
            t->y3 * t->z2 * t->x1;
    }
}

[免责声明，从未编译或运行，使用风险自负]

在这里，288 个线程的块在合并负载中将 32 个Triangle 实例提取到共享内存，然后块中的前 32 个线程执行计算并存储 32 个结果。如果内核确实没有实现很高比例的全局内存带宽吞吐量，那么像这样的方案可能会更快。 CUDA 工具包分析工具有一些相当有用的性能指标分析，可以帮助查明代码中的性能瓶颈。像所有优化练习一样，关键是仔细的分析和基准测试，这是只有您可以做的事情。