为什么我的 CUDA 实现与我的 CPU 实现一样快
Why is my CUDA implementation equally fast as my CPU implementation
我在标准 C++ 和 CUDA 中创建了一些代码来对 1300x1300 灰度图像和 15x15 内核进行二维卷积。两个版本：
CPU:
#include <iostream>
#include <exception>
#define N 1300
#define K 15
#define K2 ((K - 1) / 2)
template<int mx, int my>
inline int index(int x, int y)
{
return x*my + y;
}
int main() {
double *image = new double[N * N];
double *kernel = new double[K * K];
double *result = new double[N * N];
for (int x=0; x<N; ++x)
for (int y=0; y<N; ++y)
{
double r = 0;
for(int i=0; i<K; ++i)
for(int j=0; j<K; ++j)
{
if (x + i - K2 >= 0 and
x + i - K2 < N and
y + j - K2 >= 0 and
y + j - K2 < N)
{
r += kernel[index<K,K>(i,j)] * image[index<N,N>(x+i-K2, y+j-K2)];
}
}
result[index<N,N>(x, y)] = r;
}
delete[] image;
delete[] kernel;
delete[] result;
}
显卡：
#include <iostream>
#include <exception>
// ignore, just for error handling
struct ErrorHandler {
int d_line;
char const *d_file;
ErrorHandler(int line, char const *file) : d_line(line), d_file(file) {};
};
#define EH ErrorHandler(__LINE__, __FILE__)
ErrorHandler operator<<(ErrorHandler eh, cudaError_t err)
{
if (err != cudaSuccess)
{
std::cerr << cudaGetErrorString( err ) << " in " << eh.d_file << " at line " << eh.d_line << '\n';
throw std::exception();
}
return eh;
}
// end.
#define N 1300
#define K 15
#define K2 ((K - 1) / 2)
template<int mx, int my>
__device__ inline int index(int x, int y)
{
return x*my + y;
}
__global__ void kernelkernel(double *image, double *kernel, double *result)
{
int x = blockIdx.x;
int y = blockIdx.y; // becomes: int y = threadIdx.x;
double r = 0;
for(int i=0; i<K; ++i)
for(int j=0; j<K; ++j)
{
if (x + i - K2 >= 0 and
x + i - K2 < N and
y + j - K2 >= 0 and
y + j - K2 < N)
{
r += kernel[index<K,K>(i,j)] * image[index<N,N>(x+i-K2, y+j-K2)];
}
}
result[index<N,N>(x, y)] = r;
}
int main() {
double *image = new double[N * N];
double *kernel = new double[K * K];
double *result = new double[N * N];
double *image_cuda;
double *kernel_cuda;
double *result_cuda;
EH << cudaMalloc((void **) &image_cuda, N*N*sizeof(double));
EH << cudaMalloc((void **) &kernel_cuda, K*K*sizeof(double));
EH << cudaMalloc((void **) &result_cuda, N*N*sizeof(double));
EH << cudaMemcpy(image_cuda, image, N*N*sizeof(double), cudaMemcpyHostToDevice);
EH << cudaMemcpy(kernel_cuda, kernel, K*K*sizeof(double), cudaMemcpyHostToDevice);
dim3 grid ( N, N );
kernelkernel<<<grid, 1>>>(image_cuda, kernel_cuda, result_cuda);
// replace previous 2 statements with:
// kernelkernel<<<N, N>>>(image_cuda, kernel_cuda, result_cuda);
EH << cudaMemcpy(result, result_cuda, N*N*sizeof(double), cudaMemcpyDeviceToHost);
cudaFree( image_cuda );
cudaFree( kernel_cuda );
cudaFree( result_cuda );
delete[] image;
delete[] kernel;
delete[] result;
}
我希望 cuda 代码会快很多，但是：
$ nvprof ./gpuversion
==17806== NVPROF is profiling process 17806, command: ./gpuversion
==17806== Profiling application: ./gpuversion
==17806== Profiling result:
Time(%) Time Calls Avg Min Max Name
99.89% 3.83149s 1 3.83149s 3.83149s 3.83149s kernelkernel(double*, double*, double*)
0.07% 2.6420ms 1 2.6420ms 2.6420ms 2.6420ms [CUDA memcpy DtoH]
0.04% 1.5111ms 2 755.54us 736ns 1.5103ms [CUDA memcpy HtoD]
并且：
$ time ./cpuversion
real 0m3.382s
user 0m3.371s
sys 0m0.012s
它们的差异在统计上不显着。 CUDA 内核大约需要 3-4 秒，为什么它没有快很多？我的代码 运行 是并行的吗？
PS：我是 CUDA 的新手，所以我可能会遗漏一些微不足道的东西。
解决方案
我发现，CUDA 不允许您随意从块访问内存。我猜CUDA编程的大体策略是：
- 使用 cudaMalloc 和 cudaMemCpy 将内存从 RAM 分配和复制到 cuda
- 以不同块访问的内存不会重叠太多的方式在块和线程之间划分工作负载。
- 如果块使用的内存之间存在重叠，则通过在 shared 数组中复制内存来启动每个块。请注意：
- 编译时必须知道这个数组的大小
- 大小有限
- 此内存由一个块中的每个线程共享，因此 __shareddouble foo[10] 为每个块分配 10 个双精度数。
- 将一个block所需要的内存复制到内核中的shared变量中。当然，你使用不同的线程来做这个'efficiently'
- 同步线程，这样所有数据在使用前都已存在。
- 处理数据，写入结果。它到内核的输出数组
- 再次同步，我不知道为什么，但互联网上的每个人都在这样做 :S
- 将 GPU 内存复制回 RAM
- 清理 GPU 内存。
这给出了以下代码。它是 mex 代码，用于 Matlab 的结构相似性，它也通过滑动内核工作，但超过 2 个图像并且与点积具有不同的聚合。
// author: Herbert Kruitbosch, CC: be nice, include my name in documentation/papers/publications when used
#include <matrix.h>
#include <mex.h>
#include <cmath>
#include <iostream>
#include <fstream>
#include <iostream>
#include <stdio.h>
static void HandleError(
cudaError_t err,
const char *file,
int line )
{
if (err != cudaSuccess)
{
printf( "%s in %s at line %d\n", cudaGetErrorString( err ), file, line );
exit( EXIT_FAILURE );
}
}
#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))
#define TILE_WIDTH 31
__device__ inline double sim(double v0, double v1, double c)
{
return (c + 2*v0*v1) / (c + v1*v1 + v0*v0);
}
__device__ inline int index(int rows, int cols, int row, int col)
{
return row + col*rows;
}
__global__ void ssimkernel(double *test, double *reference, const double * __restrict__ kernel, double *ssim, int k, int rows, int cols, int tile_batches_needed)
{
int radius = k / 2;
int block_width = TILE_WIDTH - k + 1;
__shared__ double tile_test [TILE_WIDTH][TILE_WIDTH];
__shared__ double tile_reference[TILE_WIDTH][TILE_WIDTH];
for(int offset=0; offset < tile_batches_needed; ++offset)
{
int dest = block_width*block_width*offset + threadIdx.y * block_width + threadIdx.x;
int destRow = dest / TILE_WIDTH;
int destCol = dest % TILE_WIDTH;
int srcRow = blockIdx.y * block_width + destRow - radius;
int srcCol = blockIdx.x * block_width + destCol - radius;
int src = srcCol * rows + srcRow;
if (destRow < TILE_WIDTH)
{
if (srcRow >= 0 and srcRow < rows and
srcCol >= 0 and srcCol < cols)
{
tile_test [destRow][destCol] = test [src];
tile_reference[destRow][destCol] = reference[src];
}
else
{
tile_test [destRow][destCol] = 0;
tile_reference[destRow][destCol] = 0;
}
}
}
__syncthreads();
double mean_test = 0;
double mean_reference = 0;
for(int i=0; i<k; ++i)
for(int j=0; j<k; ++j)
{
double w = kernel[i * k + j];
mean_test += w * tile_test [threadIdx.y+i][threadIdx.x+j];
mean_reference += w * tile_reference[threadIdx.y+i][threadIdx.x+j];
}
double var_test = 0;
double var_reference = 0;
double correlation = 0;
for(int i=0; i<k; ++i)
for(int j=0; j<k; ++j)
{
double w = kernel[i * k + j];
double a = (tile_test [threadIdx.y+i][threadIdx.x+j] - mean_test );
double b = (tile_reference[threadIdx.y+i][threadIdx.x+j] - mean_reference);
var_test += w * a * a;
var_reference += w * b * b;
correlation += w * a * b;
}
int destRow = blockIdx.y * block_width + threadIdx.y;
int destCol = blockIdx.x * block_width + threadIdx.x;
if (destRow < rows and destCol < cols)
ssim[destCol * rows + destRow] = sim(mean_test, mean_reference, 0.01) * (0.03 + 2*correlation) / (0.03 + var_test + var_reference);
__syncthreads();
}
template<typename T>
inline T sim(T v0, T v1, T c)
{
return (c + 2*v0*v1) / (c + v1*v1 + v0*v0);
}
inline int upperdiv(int a, int b) {
return (a + b - 1) / b;
}
void mexFunction(int nargout, mxArray *argout[], int nargin, const mxArray *argin[])
{
mwSize rows = mxGetDimensions(argin[0])[0];
mwSize cols = mxGetDimensions(argin[0])[1];
mwSize k = mxGetDimensions(argin[2])[0];
mwSize channels = mxGetNumberOfDimensions(argin[0]) <= 2 ? 1 : mxGetDimensions(argin[0])[2];
int dims[] = {rows, cols, channels};
argout[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
double *test = (double *)mxGetData(argin[0]);
double *reference = (double *)mxGetData(argin[1]);
double *gaussian = (double *)mxGetData(argin[2]);
double *ssim = (double *)mxGetData(argout[0]);
double *test_cuda;
double *reference_cuda;
double *gaussian_cuda;
double *ssim_cuda;
HANDLE_ERROR( cudaMalloc((void **) &test_cuda, rows*cols*sizeof(double)) );
HANDLE_ERROR( cudaMalloc((void **) &reference_cuda, rows*cols*sizeof(double)) );
HANDLE_ERROR( cudaMalloc((void **) &gaussian_cuda, k*k*sizeof(double)) );
HANDLE_ERROR( cudaMalloc((void **) &ssim_cuda, rows*cols*sizeof(double)) );
HANDLE_ERROR( cudaMemcpy(gaussian_cuda, gaussian, k*k*sizeof(double), cudaMemcpyHostToDevice) );
int block_width = TILE_WIDTH - k + 1;
int tile_batches_needed = upperdiv(TILE_WIDTH*TILE_WIDTH, block_width*block_width);
for(int c=0; c<channels; ++c)
{
HANDLE_ERROR( cudaMemcpy(test_cuda, test + rows*cols*c, rows*cols*sizeof(double), cudaMemcpyHostToDevice) );
HANDLE_ERROR( cudaMemcpy(reference_cuda, reference + rows*cols*c, rows*cols*sizeof(double), cudaMemcpyHostToDevice) );
dim3 dimGrid(upperdiv(cols, block_width), upperdiv(rows, block_width), 1);
dim3 dimBlock(block_width, block_width, 1);
ssimkernel<<<dimGrid, dimBlock>>>(test_cuda, reference_cuda, gaussian_cuda, ssim_cuda, k, rows, cols, tile_batches_needed);
HANDLE_ERROR( cudaMemcpy(ssim + rows*cols*c, ssim_cuda, rows*cols*sizeof(double), cudaMemcpyDeviceToHost) );
}
cudaFree( test_cuda );
cudaFree( reference_cuda );
cudaFree( gaussian_cuda );
cudaFree( ssim_cuda );
}
如果您在 CUDA 中使用全局内存，所有数据访问都将在队列之类的东西中同步，您将获得几乎线性的解决方案，而不是并行的。
此外，将大型数据集从 RAM 内存传输到 GPU 内存也需要很多时间（总线速度有限）。
所以，我认为你必须以某种方式在你的 GPU 中跨计算单元并行你的数据（将它们分成共享内存）。
检查 this 以查看在与您的情况类似的情况下如何提高 GPU 内存使用率的解决方案。
kernelkernel<<<grid, 1>>>
这是一个重大问题； nVidia GPU 上的线程以 32 个线程的 warp 工作。但是，您只为每个块分配了一个线程，这意味着其中 31 个线程将处于空闲状态，而只有一个线程在工作。通常，对于具有灵活性的内核，您通常希望每个块有多个 warp 而不是一个。
通过使用 N 个块和每个块的 N 个线程，而不是使用 N^2 个块，您可以立即获得加速。
实际上，N 可能太大了，因为每个块的线程数有上限。尽管您可以选择合适的 M，以便每个块使用 N/M 个线程，并且 N * M 个块。
事实上，您可能会在这方面获得最佳结果，方法是选择一些 M（我猜 256 可能接近最佳）并启动 L=ceiling(N*N/M) 个块和每个线程 M 个块。然后每个线程数字根据其块和线程 ID 在 [0, M*L) 中重建一个索引，然后那些索引在 [0,N*N) 中的线程将继续将该索引拆分为 x 和 y 坐标并进行工作。
由于延迟，访问内核中的全局内存成本很高。全局内存请求（读取和写入）需要数百个时钟周期才能完成。您希望最大程度地减少访问全局内存的次数，并在连续的块中访问它。
如果每条数据只被访问一次，则延迟无关紧要，但这种情况很少见。在您的代码中绝对不是这种情况，其中 kernel 数组被所有线程以相同的模式访问，并且很多 image 也被多个线程访问。
解决方案是通过从高延迟 全局内存 获取数据到低延迟 共享内存来启动内核。 共享内存是多处理器上的一块内存，其延迟与寄存器相当。所以大多数简单的内核都遵循这样的结构：
每个线程从全局内存中获取数据到共享内存中。如果可能，您希望以连续的顺序获取数据，因为全局内存是通过事务访问的。如果没有足够的数据供所有线程获取，请让其中一些线程空闲。
线程对共享内存中的数据进行操作。
数据从共享内存写回全局内存，其模式与在步骤 1 中获取的模式相同。
共享内存由线程块内的所有线程共享。这引出了您代码中的第二个大问题：您根本没有使用线程块。在一个多处理器上一个块 运行 中的线程，共享共享内存，可以相互同步等。您需要将线程很好地组织到块中以充分利用它们。
块网格只是一种能够在一次调用中 运行 更多块的机制。并行指令执行和共享内存访问的所有好处都在 一个块内。块的网格只是 "yeah, sorry, my data's so big a single block won't do, just run many of them."
你做的恰恰相反：你的块每个都有一个线程，这意味着在每个步骤中，多处理器上的每个 warp 运行s 只有一个线程（基于你设备的计算能力和可用的 warp 调度器的数量，这意味着一个多处理器上最多有 2-4 个线程。
您必须重新构建线程以反映数据访问模式，并将数据预取到共享内存中。这将为您带来预期的性能提升。
以上只是一个简短的总结。有关块组织、共享内存和全局内存事务的详细信息，请参阅 CUDA 编程指南。
我在标准 C++ 和 CUDA 中创建了一些代码来对 1300x1300 灰度图像和 15x15 内核进行二维卷积。两个版本：
CPU:
#include <iostream>
#include <exception>
#define N 1300
#define K 15
#define K2 ((K - 1) / 2)
template<int mx, int my>
inline int index(int x, int y)
{
return x*my + y;
}
int main() {
double *image = new double[N * N];
double *kernel = new double[K * K];
double *result = new double[N * N];
for (int x=0; x<N; ++x)
for (int y=0; y<N; ++y)
{
double r = 0;
for(int i=0; i<K; ++i)
for(int j=0; j<K; ++j)
{
if (x + i - K2 >= 0 and
x + i - K2 < N and
y + j - K2 >= 0 and
y + j - K2 < N)
{
r += kernel[index<K,K>(i,j)] * image[index<N,N>(x+i-K2, y+j-K2)];
}
}
result[index<N,N>(x, y)] = r;
}
delete[] image;
delete[] kernel;
delete[] result;
}
显卡：
#include <iostream>
#include <exception>
// ignore, just for error handling
struct ErrorHandler {
int d_line;
char const *d_file;
ErrorHandler(int line, char const *file) : d_line(line), d_file(file) {};
};
#define EH ErrorHandler(__LINE__, __FILE__)
ErrorHandler operator<<(ErrorHandler eh, cudaError_t err)
{
if (err != cudaSuccess)
{
std::cerr << cudaGetErrorString( err ) << " in " << eh.d_file << " at line " << eh.d_line << '\n';
throw std::exception();
}
return eh;
}
// end.
#define N 1300
#define K 15
#define K2 ((K - 1) / 2)
template<int mx, int my>
__device__ inline int index(int x, int y)
{
return x*my + y;
}
__global__ void kernelkernel(double *image, double *kernel, double *result)
{
int x = blockIdx.x;
int y = blockIdx.y; // becomes: int y = threadIdx.x;
double r = 0;
for(int i=0; i<K; ++i)
for(int j=0; j<K; ++j)
{
if (x + i - K2 >= 0 and
x + i - K2 < N and
y + j - K2 >= 0 and
y + j - K2 < N)
{
r += kernel[index<K,K>(i,j)] * image[index<N,N>(x+i-K2, y+j-K2)];
}
}
result[index<N,N>(x, y)] = r;
}
int main() {
double *image = new double[N * N];
double *kernel = new double[K * K];
double *result = new double[N * N];
double *image_cuda;
double *kernel_cuda;
double *result_cuda;
EH << cudaMalloc((void **) &image_cuda, N*N*sizeof(double));
EH << cudaMalloc((void **) &kernel_cuda, K*K*sizeof(double));
EH << cudaMalloc((void **) &result_cuda, N*N*sizeof(double));
EH << cudaMemcpy(image_cuda, image, N*N*sizeof(double), cudaMemcpyHostToDevice);
EH << cudaMemcpy(kernel_cuda, kernel, K*K*sizeof(double), cudaMemcpyHostToDevice);
dim3 grid ( N, N );
kernelkernel<<<grid, 1>>>(image_cuda, kernel_cuda, result_cuda);
// replace previous 2 statements with:
// kernelkernel<<<N, N>>>(image_cuda, kernel_cuda, result_cuda);
EH << cudaMemcpy(result, result_cuda, N*N*sizeof(double), cudaMemcpyDeviceToHost);
cudaFree( image_cuda );
cudaFree( kernel_cuda );
cudaFree( result_cuda );
delete[] image;
delete[] kernel;
delete[] result;
}
我希望 cuda 代码会快很多，但是：
$ nvprof ./gpuversion
==17806== NVPROF is profiling process 17806, command: ./gpuversion
==17806== Profiling application: ./gpuversion
==17806== Profiling result:
Time(%) Time Calls Avg Min Max Name
99.89% 3.83149s 1 3.83149s 3.83149s 3.83149s kernelkernel(double*, double*, double*)
0.07% 2.6420ms 1 2.6420ms 2.6420ms 2.6420ms [CUDA memcpy DtoH]
0.04% 1.5111ms 2 755.54us 736ns 1.5103ms [CUDA memcpy HtoD]
并且：
$ time ./cpuversion
real 0m3.382s
user 0m3.371s
sys 0m0.012s
它们的差异在统计上不显着。 CUDA 内核大约需要 3-4 秒，为什么它没有快很多？我的代码 运行 是并行的吗？
PS：我是 CUDA 的新手，所以我可能会遗漏一些微不足道的东西。
解决方案
我发现，CUDA 不允许您随意从块访问内存。我猜CUDA编程的大体策略是：
- 使用 cudaMalloc 和 cudaMemCpy 将内存从 RAM 分配和复制到 cuda
- 以不同块访问的内存不会重叠太多的方式在块和线程之间划分工作负载。
- 如果块使用的内存之间存在重叠，则通过在 shared 数组中复制内存来启动每个块。请注意：
- 编译时必须知道这个数组的大小
- 大小有限
- 此内存由一个块中的每个线程共享，因此 __shareddouble foo[10] 为每个块分配 10 个双精度数。
- 将一个block所需要的内存复制到内核中的shared变量中。当然，你使用不同的线程来做这个'efficiently'
- 同步线程，这样所有数据在使用前都已存在。
- 处理数据，写入结果。它到内核的输出数组
- 再次同步，我不知道为什么，但互联网上的每个人都在这样做 :S
- 将 GPU 内存复制回 RAM
- 清理 GPU 内存。
这给出了以下代码。它是 mex 代码，用于 Matlab 的结构相似性，它也通过滑动内核工作，但超过 2 个图像并且与点积具有不同的聚合。
// author: Herbert Kruitbosch, CC: be nice, include my name in documentation/papers/publications when used
#include <matrix.h>
#include <mex.h>
#include <cmath>
#include <iostream>
#include <fstream>
#include <iostream>
#include <stdio.h>
static void HandleError(
cudaError_t err,
const char *file,
int line )
{
if (err != cudaSuccess)
{
printf( "%s in %s at line %d\n", cudaGetErrorString( err ), file, line );
exit( EXIT_FAILURE );
}
}
#define HANDLE_ERROR( err ) (HandleError( err, __FILE__, __LINE__ ))
#define TILE_WIDTH 31
__device__ inline double sim(double v0, double v1, double c)
{
return (c + 2*v0*v1) / (c + v1*v1 + v0*v0);
}
__device__ inline int index(int rows, int cols, int row, int col)
{
return row + col*rows;
}
__global__ void ssimkernel(double *test, double *reference, const double * __restrict__ kernel, double *ssim, int k, int rows, int cols, int tile_batches_needed)
{
int radius = k / 2;
int block_width = TILE_WIDTH - k + 1;
__shared__ double tile_test [TILE_WIDTH][TILE_WIDTH];
__shared__ double tile_reference[TILE_WIDTH][TILE_WIDTH];
for(int offset=0; offset < tile_batches_needed; ++offset)
{
int dest = block_width*block_width*offset + threadIdx.y * block_width + threadIdx.x;
int destRow = dest / TILE_WIDTH;
int destCol = dest % TILE_WIDTH;
int srcRow = blockIdx.y * block_width + destRow - radius;
int srcCol = blockIdx.x * block_width + destCol - radius;
int src = srcCol * rows + srcRow;
if (destRow < TILE_WIDTH)
{
if (srcRow >= 0 and srcRow < rows and
srcCol >= 0 and srcCol < cols)
{
tile_test [destRow][destCol] = test [src];
tile_reference[destRow][destCol] = reference[src];
}
else
{
tile_test [destRow][destCol] = 0;
tile_reference[destRow][destCol] = 0;
}
}
}
__syncthreads();
double mean_test = 0;
double mean_reference = 0;
for(int i=0; i<k; ++i)
for(int j=0; j<k; ++j)
{
double w = kernel[i * k + j];
mean_test += w * tile_test [threadIdx.y+i][threadIdx.x+j];
mean_reference += w * tile_reference[threadIdx.y+i][threadIdx.x+j];
}
double var_test = 0;
double var_reference = 0;
double correlation = 0;
for(int i=0; i<k; ++i)
for(int j=0; j<k; ++j)
{
double w = kernel[i * k + j];
double a = (tile_test [threadIdx.y+i][threadIdx.x+j] - mean_test );
double b = (tile_reference[threadIdx.y+i][threadIdx.x+j] - mean_reference);
var_test += w * a * a;
var_reference += w * b * b;
correlation += w * a * b;
}
int destRow = blockIdx.y * block_width + threadIdx.y;
int destCol = blockIdx.x * block_width + threadIdx.x;
if (destRow < rows and destCol < cols)
ssim[destCol * rows + destRow] = sim(mean_test, mean_reference, 0.01) * (0.03 + 2*correlation) / (0.03 + var_test + var_reference);
__syncthreads();
}
template<typename T>
inline T sim(T v0, T v1, T c)
{
return (c + 2*v0*v1) / (c + v1*v1 + v0*v0);
}
inline int upperdiv(int a, int b) {
return (a + b - 1) / b;
}
void mexFunction(int nargout, mxArray *argout[], int nargin, const mxArray *argin[])
{
mwSize rows = mxGetDimensions(argin[0])[0];
mwSize cols = mxGetDimensions(argin[0])[1];
mwSize k = mxGetDimensions(argin[2])[0];
mwSize channels = mxGetNumberOfDimensions(argin[0]) <= 2 ? 1 : mxGetDimensions(argin[0])[2];
int dims[] = {rows, cols, channels};
argout[0] = mxCreateNumericArray(3, dims, mxDOUBLE_CLASS, mxREAL);
double *test = (double *)mxGetData(argin[0]);
double *reference = (double *)mxGetData(argin[1]);
double *gaussian = (double *)mxGetData(argin[2]);
double *ssim = (double *)mxGetData(argout[0]);
double *test_cuda;
double *reference_cuda;
double *gaussian_cuda;
double *ssim_cuda;
HANDLE_ERROR( cudaMalloc((void **) &test_cuda, rows*cols*sizeof(double)) );
HANDLE_ERROR( cudaMalloc((void **) &reference_cuda, rows*cols*sizeof(double)) );
HANDLE_ERROR( cudaMalloc((void **) &gaussian_cuda, k*k*sizeof(double)) );
HANDLE_ERROR( cudaMalloc((void **) &ssim_cuda, rows*cols*sizeof(double)) );
HANDLE_ERROR( cudaMemcpy(gaussian_cuda, gaussian, k*k*sizeof(double), cudaMemcpyHostToDevice) );
int block_width = TILE_WIDTH - k + 1;
int tile_batches_needed = upperdiv(TILE_WIDTH*TILE_WIDTH, block_width*block_width);
for(int c=0; c<channels; ++c)
{
HANDLE_ERROR( cudaMemcpy(test_cuda, test + rows*cols*c, rows*cols*sizeof(double), cudaMemcpyHostToDevice) );
HANDLE_ERROR( cudaMemcpy(reference_cuda, reference + rows*cols*c, rows*cols*sizeof(double), cudaMemcpyHostToDevice) );
dim3 dimGrid(upperdiv(cols, block_width), upperdiv(rows, block_width), 1);
dim3 dimBlock(block_width, block_width, 1);
ssimkernel<<<dimGrid, dimBlock>>>(test_cuda, reference_cuda, gaussian_cuda, ssim_cuda, k, rows, cols, tile_batches_needed);
HANDLE_ERROR( cudaMemcpy(ssim + rows*cols*c, ssim_cuda, rows*cols*sizeof(double), cudaMemcpyDeviceToHost) );
}
cudaFree( test_cuda );
cudaFree( reference_cuda );
cudaFree( gaussian_cuda );
cudaFree( ssim_cuda );
}
如果您在 CUDA 中使用全局内存，所有数据访问都将在队列之类的东西中同步，您将获得几乎线性的解决方案，而不是并行的。
此外，将大型数据集从 RAM 内存传输到 GPU 内存也需要很多时间（总线速度有限）。
所以，我认为你必须以某种方式在你的 GPU 中跨计算单元并行你的数据（将它们分成共享内存）。
检查 this 以查看在与您的情况类似的情况下如何提高 GPU 内存使用率的解决方案。
kernelkernel<<<grid, 1>>>
这是一个重大问题； nVidia GPU 上的线程以 32 个线程的 warp 工作。但是，您只为每个块分配了一个线程，这意味着其中 31 个线程将处于空闲状态，而只有一个线程在工作。通常，对于具有灵活性的内核，您通常希望每个块有多个 warp 而不是一个。
通过使用 N 个块和每个块的 N 个线程，而不是使用 N^2 个块，您可以立即获得加速。
实际上，N 可能太大了，因为每个块的线程数有上限。尽管您可以选择合适的 M，以便每个块使用 N/M 个线程，并且 N * M 个块。
事实上，您可能会在这方面获得最佳结果，方法是选择一些 M（我猜 256 可能接近最佳）并启动 L=ceiling(N*N/M) 个块和每个线程 M 个块。然后每个线程数字根据其块和线程 ID 在 [0, M*L) 中重建一个索引，然后那些索引在 [0,N*N) 中的线程将继续将该索引拆分为 x 和 y 坐标并进行工作。
由于延迟，访问内核中的全局内存成本很高。全局内存请求（读取和写入）需要数百个时钟周期才能完成。您希望最大程度地减少访问全局内存的次数，并在连续的块中访问它。
如果每条数据只被访问一次，则延迟无关紧要，但这种情况很少见。在您的代码中绝对不是这种情况，其中 kernel 数组被所有线程以相同的模式访问，并且很多 image 也被多个线程访问。
解决方案是通过从高延迟 全局内存 获取数据到低延迟 共享内存来启动内核。 共享内存是多处理器上的一块内存，其延迟与寄存器相当。所以大多数简单的内核都遵循这样的结构：
每个线程从全局内存中获取数据到共享内存中。如果可能，您希望以连续的顺序获取数据，因为全局内存是通过事务访问的。如果没有足够的数据供所有线程获取，请让其中一些线程空闲。
线程对共享内存中的数据进行操作。
数据从共享内存写回全局内存，其模式与在步骤 1 中获取的模式相同。
共享内存由线程块内的所有线程共享。这引出了您代码中的第二个大问题：您根本没有使用线程块。在一个多处理器上一个块 运行 中的线程，共享共享内存，可以相互同步等。您需要将线程很好地组织到块中以充分利用它们。
块网格只是一种能够在一次调用中 运行 更多块的机制。并行指令执行和共享内存访问的所有好处都在 一个块内。块的网格只是 "yeah, sorry, my data's so big a single block won't do, just run many of them."
你做的恰恰相反：你的块每个都有一个线程，这意味着在每个步骤中，多处理器上的每个 warp 运行s 只有一个线程（基于你设备的计算能力和可用的 warp 调度器的数量，这意味着一个多处理器上最多有 2-4 个线程。
您必须重新构建线程以反映数据访问模式，并将数据预取到共享内存中。这将为您带来预期的性能提升。
以上只是一个简短的总结。有关块组织、共享内存和全局内存事务的详细信息，请参阅 CUDA 编程指南。