1 共享内存（一个block 48KB=49152B ）
常用于向量点积1
■ 在考虑线程块大小的时候经常用到向上取整，这里使用了技巧 ceil( a / b ) == floor( (a-1) / b) + 1
■ 算法总体想法是在GPU中将很长的向量分段放入GPU的各线程块中，每个线程块利用共享内存和多线程分别计算乘法和加法。结果整理为每个线程块输出一个浮点数，置于全局内存中，这样就将待计算的元素数量降到了 gridDim.x 的水平，再返回CPU中完成剩下的加法。
■ 算法预先规定了每个线程块使用256个线程（blockDim.x == 256），那么使用的线程块数量应该满足 gridDim.x * blockDim.x ≥ N（待计算的向量长度），另外代码中规定线程块数量至少为32（“选择其他的只可能产生更高或更差的性能，这取决于CPU和GPU的相对速度”）
■ 在核函数中使用了既定大小的共享内存 shared float cache[threadsPerBlock]; ，并采用 __syncthreads(); 函数进行线程同步（因为接下来要进行规约运算，前提就是该线程块内所有的线程已经独立计算完毕）。
/*
* Copyright 1993-2010 NVIDIA Corporation. All rights reserved.
*
* NVIDIA Corporation and its licensors retain all intellectual property and
* proprietary rights in and to this software and related documentation.
* Any use, reproduction, disclosure, or distribution of this software
* and related documentation without an express license agreement from
* NVIDIA Corporation is strictly prohibited.
*
* Please refer to the applicable NVIDIA end user license agreement (EULA)
* associated with this source code for terms and conditions that govern
* your use of this NVIDIA software.
*
*/
#include "../common/book.h"
#define imin(a,b) (a<b?a:b)
const int N = 33 * 1024;
const int threadsPerBlock = 256;
const int blocksPerGrid =
imin( 32, (N+threadsPerBlock-1) / threadsPerBlock );
__global__ void dot( float *a, float *b, float *c ) {
__shared__ float cache[threadsPerBlock]; //共享内存申请方式：__shared__ 变量类型 变量名；
int tid = threadIdx.x + blockIdx.x * blockDim.x;
int cacheIndex = threadIdx.x;
float temp = 0;
while (tid < N) {
temp += a[tid] * b[tid];
tid += blockDim.x * gridDim.x;
}
// set the cache values
cache[cacheIndex] = temp;
// synchronize threads in this block
__syncthreads();
// for reductions, threadsPerBlock must be a power of 2
// because of the following code
int i = blockDim.x/2;
while (i != 0) {
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads();
i /= 2;
}
if (cacheIndex == 0)
c[blockIdx.x] = cache[0];
}
int main( void ) {
float *a, *b, c, *partial_c;
float *dev_a, *dev_b, *dev_partial_c;
// allocate memory on the cpu side
a = (float*)malloc( N*sizeof(float) );
b = (float*)malloc( N*sizeof(float) );
partial_c = (float*)malloc( blocksPerGrid*sizeof(float) );
// allocate the memory on the GPU
HANDLE_ERROR( cudaMalloc( (void**)&dev_a,
N*sizeof(float) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_b,
N*sizeof(float) ) );
HANDLE_ERROR( cudaMalloc( (void**)&dev_partial_c,
blocksPerGrid*sizeof(float) ) );
// fill in the host memory with data
for (int i=0; i<N; i++) {
a[i] = i;
b[i] = i*2;
}
// copy the arrays 'a' and 'b' to the GPU
HANDLE_ERROR( cudaMemcpy( dev_a, a, N*sizeof(float),
cudaMemcpyHostToDevice ) );
HANDLE_ERROR( cudaMemcpy( dev_b, b, N*sizeof(float),
cudaMemcpyHostToDevice ) );
dot<<<blocksPerGrid,threadsPerBlock>>>( dev_a, dev_b,
dev_partial_c );
// copy the array 'c' back from the GPU to the CPU
HANDLE_ERROR( cudaMemcpy( partial_c, dev_partial_c,
blocksPerGrid*sizeof(float),
cudaMemcpyDeviceToHost ) );
// finish up on the CPU side
c = 0;
for (int i=0; i<blocksPerGrid; i++) {
c += partial_c[i];
}
#define sum_squares(x) (x*(x+1)*(2*x+1)/6)
printf( "Does GPU value %.6g = %.6g?\n", c,
2 * sum_squares( (float)(N - 1) ) );
// free memory on the gpu side
HANDLE_ERROR( cudaFree( dev_a ) );
HANDLE_ERROR( cudaFree( dev_b ) );
HANDLE_ERROR( cudaFree( dev_partial_c ) );
// free memory on the cpu side
free( a );
free( b );
free( partial_c );
}
错误的同步方式：
while (i != 0) {
if (cacheIndex < i)
cache[cacheIndex] += cache[cacheIndex + i];
__syncthreads(); //
i /= 2;
}
由于if的存在，只有部分线程包含同步操作，部分线程永远都不可能执行到cache[cacheIndex] += cache[cacheIndex + i];这一步，因此就要永远等待下去，因而程序无法执行。
2 常量内存
用来保存在核函数执行期间不会发生变化的数据。
NVIDIA硬件提供了64KB的常量内存，并且常量内存采用了不同于标准全局内存的处理方式。在某些情况下，用常量内存来替换全局内存可以有效地减少内存带宽。
2.1 声明方式
__constant__ int s[10000] //变量修饰符 __constant__ 将变量的访问限制为只读。
2.2 传递参数
我们为变量分配内存时是先声明一个指针，然后通过cudaMalloc()来为指针分配GPU内存。而当我们将其改为常量内存时，则要将这个声明修改为在常量内存中静态地分配空间。我们不再需要对变量指针调用cudaMalloc()或者cudaFree()，而是在编译时为这个变量（如数组s）提交固定的大小。
当从主机内存复制到GPU上的常量内存时，我们需要使用一个特殊版本的cudaMemcpy()
cudaMemcpyToSymbol()
2.3 特点
对常量内存的单次读操作可以广播到其他的“邻近(nearby)”线程，这将节约15次读取操作；
常量内存的数据将缓存起来，因此对于相同地址的连续操作将不会产生额外的内存通信量。
当处理常量内存时，NVIDIA硬件将把单次内存读取操作广播到每个半线程束(Half-Warp)。在半线程束中包含16个线程，即线程束中线程数量的一半。如果在半线程束中的每个线程从常量内存的相同地址上读取数据，那么GPU只会产生一次读取请求并在随后将数据广播到每个线程。如果从常量内存中读取大量数据，那么这种方式产生的内存流量只是使用全局内存时的1/16。
但在读取常量内存时，所节约的并不仅限于减少94%的带宽。由于这块内存的内容是不发生变化的，因此硬件将主动把这个常量数据缓存在GPU上。在第一次从常量内存的某个地址上读取后，当其他半线程束请求同一个地址时，那么将命中缓存，这同样减少了额外的内存流量。
注意事项
半线程束广播功能实际上是一把双刃剑。虽然当所有16个线程都读取相同地址时，这个功能可以极大提升性能，但当所有16个线程分别读取不同的地址时，它实际上会降低性能。因为这16次不同的读取操作会被串行化，从而需要16倍的时间来发出请求。但如果从全局内存中读取，那么这些请求会同时发出。
3 纹理内存
和常量内存一样，纹理内存是另一种类型的只读内存，在特定的访问模式中，纹理内存同样能够提升性能并减少内存流量。
纹理内存缓存在芯片上，因此在某些情况中，它能够减少对内存的请求并提供更高效的内存带宽。纹理缓存是专门为那些在内存访问模式中存在大量空间局部性(Spatial Locality)的图形应用程序而设计的。在某个计算应用程序中，这意味着一个线程读取的位置可能与邻近线程的读取位置“非常接近”，如下图所示。
从数学的角度，上图中的4个地址并非连续的，在一般的CPU缓存中，这些地址将不会缓存。但由于GPU纹理缓存是专门为了加速这种访问模式而设计的，因此如果在这种情况中使用纹理内存而不是全局内存，那么将会获得性能的提升。
3.1 一维纹理内存
以一维纹理内存为例，我们来看看如何使用纹理内存。
3.1.1 声明方式
首先，需要将输入的数据声明为texture类型的引用。比如：
// 变量将位于GPU上
texture<float> texConstSrc;
接下来分配内存：
HANDLE_ERROR( cudaMalloc( (void**)&data.dev_inSrc, imageSize ) );
3.1.2 变量绑定
需要通过cudaBindTexture()将这些变量绑定到内存缓存区。这相当于告诉CUDA运行时两件事：
HANDLE_ERROR( cudaBindTexture( NULL, texConstSrc, data.dev_inSrc, imageSize ) )
我们希望将指定的缓冲区作为纹理来使用；
我们希望将纹理引用作为纹理的“名字”。
然而，当读取核函数中的纹理时，需要通过特殊的函数来告诉GPU将读取请求转发到纹理内存而不是标准的全局内存。因此，当读取内存时，需要使用特殊的方式。
从线性内存中读取（拾取），使用的函数是tex1Dfetch()：
template<class Type>
Type tex1Dfetch(
texture<Type, 1, cudaReadModeElementType> texRef,
int x);
float tex1Dfetch(
texture<unsigned char, 1, cudaReadModeNormalizedFloat> texRef,
int x);
float tex1Dfetch(
texture<signed char, 1, cudaReadModeNormalizedFloat> texRef,
int x);
float tex1Dfetch(
texture<unsigned short, 1, cudaReadModeNormalizedFloat> texRef,
int x);
float tex1Dfetch(
texture<signed short, 1, cudaReadModeNormalizedFloat> texRef,
int x);
这些函数用纹理坐标x拾取绑定到纹理参考texRef 的线性内存的区域。不支持纹理过滤和寻址模式。对于整数型，这些函数将会将整数型转化为单精度浮点型。除了这些函数，还支持2-和4-分量向量的拾取。
3.1.3 清除绑定
最后，当应用程序运行结束后，还要清除纹理的绑定。
cudaUnbindTexture(texConst);
3.2 二维纹理内存
3.2.1 关键步骤
// 热传导
texture<float, 2> texSrc;// 在全局位置上声明纹理引用
float *dev_Src;
cudaMalloc((void**)&dev_Src, DIM*DIM);// 申请和绑定纹理内存
cudaChannelFormatDesc desc = cudaCreateChannelDesc<float>();
cudaBindTexture2D(NULL, texSrc, dev_Src, desc, DIM, DIM, sizeof(float) * DIM*DIM);
float *temp = (float*)malloc(sizeof(float)*DIM*DIM);// 初始化该内存中的内容
//Initalize data in temp and then free(temp)
cudaMemcpy(dev_Src, temp, sizeof(float)*DIM*DIM, cudaMemcpyHostToDevice);
//Do something
cudaUnbindTexture(texSrc);// 解绑和释放内存
cudaFree(dev_Src);
3.2.2 声明方式
texture<float, 2> texSrc;
增加了代表维数的参数2
3.2.3 访问纹理
（对于热传导模型）
- 可以直接通过x和y来访问纹理
- 使用tex2D()时，不需要担心溢出问题。如果x或y小于0，那么tex2D()将返回0处的值。同理，如果某个值大于宽度，tex2D()将返回位于宽度处的值。
int x = threadIdx.x + blockIdx.x * blockDim.x;
int y = threadIdx.y + blockIdx.y * blockDim.y;
float c = tex2D(texSrc, x, y);
3.2.4 完整代码
// 热传导演示
#include "../common/book.h"
#include "../common/cpu_anim.h"
#define DIM 1024
#define PI 3.1415926535897932f
#define MAX_TEMP 1.0f // 最高温度
#define MIN_TEMP 0.0001f // 最低温度
#define SPEED 0.25f // 热量流动速度
__global__ void blend_kernel(float *outSrc, const float *inSrc) {
// map from threadIdx/blockIdx to pixel position
//每个线程都负责计算一个单元，每个线程都将读取对应单元及其邻接单元的温度值
//执行前面给出的更新运算，即用上，下，左，右相邻的单元温度计算当前单元的温度
int x = threadIdx.x + blockIdx.x * blockDim.x;
int y = threadIdx.y + blockIdx.y * blockDim.y;
int offset = x + y * blockDim.x * gridDim.x;
int left = offset - 1;
int right = offset + 1;
if (x == 0) left++;
if (x == DIM - 1) right--;
int top = offset - DIM;
int bottom = offset + DIM;
if (y == 0) top += DIM;
if (y == DIM - 1) bottom -= DIM;
outSrc[offset] = inSrc[offset] + SPEED * (inSrc[top] + inSrc[bottom] + inSrc[left] + inSrc[right] - inSrc[offset] * 4);
}
__global__ void copy_const_kernel(float *iptr, const float *cptr) {
// map from threadIdx/blockIdx to pixel position
int x = threadIdx.x + blockIdx.x * blockDim.x;
int y = threadIdx.y + blockIdx.y * blockDim.y;
int offset = x + y * blockDim.x * gridDim.x;
if (cptr[offset] != 0) iptr[offset] = cptr[offset];//把cptr的热源温度复制到iptr的输入单元中
}
// globals needed by the update routine
struct DataBlock {
unsigned char *output_bitmap;
float *dev_inSrc;
float *dev_outSrc;
float *dev_constSrc;
CPUAnimBitmap *bitmap;// 输出位图的定义
cudaEvent_t start, stop;// 定义记录GPU处理时间的事件
float totalTime;
float frames;
};
void anim_gpu(DataBlock *d, int ticks) {
HANDLE_ERROR(cudaEventRecord(d->start, 0));
dim3 blocks(DIM / 16, DIM / 16);
dim3 threads(16, 16);
CPUAnimBitmap *bitmap = d->bitmap;
for (int i = 0; i<90; i++) {
// 计算模拟过程的一个时间步
copy_const_kernel << <blocks, threads >> >(d->dev_inSrc, d->dev_constSrc);// 算法第一步：复制温度
blend_kernel << <blocks, threads >> >(d->dev_outSrc, d->dev_inSrc);// 算法第二步：计算输出温度
swap(d->dev_inSrc, d->dev_outSrc);// 算法第三步：新的输出作为下次循环的输入
}
float_to_color << <blocks, threads >> >(d->output_bitmap, d->dev_inSrc);//float_to_color在本书带的库book.h中有相应的定义：功能是将温度与颜色建立映射关系
HANDLE_ERROR(cudaMemcpy(bitmap->get_ptr(), d->output_bitmap, bitmap->image_size(), cudaMemcpyDeviceToHost));
HANDLE_ERROR(cudaEventRecord(d->stop, 0));
HANDLE_ERROR(cudaEventSynchronize(d->stop));
float elapsedTime;
HANDLE_ERROR(cudaEventElapsedTime(&elapsedTime, d->start, d->stop));
d->totalTime += elapsedTime;
++d->frames;// 动画帧增加一
printf("Average time per frame: %3.1f msn", d->totalTime / d->frames);
}
void anim_exit(DataBlock *d) {
cudaFree(d->dev_inSrc);
cudaFree(d->dev_outSrc);
cudaFree(d->dev_constSrc);
HANDLE_ERROR(cudaEventDestroy(d->start));
HANDLE_ERROR(cudaEventDestroy(d->stop));
}
int main(void) {
DataBlock data;
CPUAnimBitmap bitmap(DIM, DIM, &data);
data.bitmap = &bitmap;
data.totalTime = 0;
data.frames = 0;
HANDLE_ERROR(cudaEventCreate(&data.start));
HANDLE_ERROR(cudaEventCreate(&data.stop));
HANDLE_ERROR(cudaMalloc((void**)&data.output_bitmap, bitmap.image_size()));
// assume float == 4 chars in size (ie., rgba)
HANDLE_ERROR(cudaMalloc((void**)&data.dev_inSrc, bitmap.image_size()));
HANDLE_ERROR(cudaMalloc((void**)&data.dev_outSrc, bitmap.image_size()));
HANDLE_ERROR(cudaMalloc((void**)&data.dev_constSrc, bitmap.image_size()));
float *temp = (float*)malloc(bitmap.image_size());
for (int i = 0; i<DIM*DIM; i++) {
temp[i] = 0;
int x = i % DIM;
int y = i / DIM;
if ((x>300) && (x<600) && (y>310) && (y<601))
temp[i] = MAX_TEMP;
}
// 初如化热源单元
temp[DIM * 100 + 100] = (MAX_TEMP + MIN_TEMP) / 2;
temp[DIM * 700 + 100] = MIN_TEMP;
temp[DIM * 300 + 300] = MIN_TEMP;
temp[DIM * 200 + 700] = MIN_TEMP;
for (int y = 800; y<900; y++) {
for (int x = 400; x<500; x++) {
temp[x + y*DIM] = MIN_TEMP;
}
}
HANDLE_ERROR(cudaMemcpy(data.dev_constSrc, temp, bitmap.image_size(), cudaMemcpyHostToDevice));
for (int y = 800; y<DIM; y++) {
for (int x = 0; x<200; x++) {
temp[x + y*DIM] = MAX_TEMP;
}
}
HANDLE_ERROR(cudaMemcpy(data.dev_inSrc, temp, bitmap.image_size(), cudaMemcpyHostToDevice));
free(temp);
bitmap.anim_and_exit((void(*)(void*, int))anim_gpu, (void(*)(void*))anim_exit);// 传递函数指针调用核函数处理动画数据
}
两个向量a = [a1, a2,…, an]和b = [b1, b2,…, bn]的点积定义为：
a·b=a1b1+a2b2+……+anbn。 ↩︎
7969
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
