CUDA By Example（三）

        GPU性能的瓶颈往往不在于芯片的数学计算吞吐量，而在于芯片的内存宽带。GPU有非常多的数字逻辑单元（ALU），因此有时输入数据的速率无法维持如此高的计算速率。
         CUDA C除了可以使用全局内存和共享内存，还支持常量内存。常量内存用于保持核函数执行期间不会发生变化的数据。NVIDA提供了64KB的常量内存。在某些情况下，用常量内存替代全局内存能有效减少内存宽带。

 

         下面的例子演示光线跟踪。

         光线跟踪（Ray Tracing）是从三维对象场景中生成二维图像的一种方式。（OpenGL、DirectX中有一种相同目的的技术：光栅化，Rasterization）。光线跟踪的原理是：在场景中选择一个位置放一台假象的相机，这台相机用光传感器来生成图像，图像的每个像素与命中传感器的光线有着相同的颜色和强度。命中的光线可能来自场景中的任意位置，因此采用逆向计算跟容易一些。从像素中投射出的光线穿过场景，直到光像命中某个物体，然后计算这个像素的颜色，这里称像素将“看到”这个物体。光线跟踪中的大部分计算都是光线与场景中物体的相交运算。
         更复杂的光线跟踪模型中，场景中会有反射光线和折射光线出现，这将生成二次射线、三次射线等。

#include "cuda_runtime.h"  
#include "device_launch_parameters.h"  
#include "device_functions.h" 
#include <stdio.h>
#include <math.h>
#include "bitmap.h"

#define INF 2e10f
#define SPHERES 20
#define DIM 1024
#define rnd( x )  (x*rand()/RAND_MAX)

struct Sphere{
	float r,b,g;
	float radius;
	float x,y,z;
	__device__ float hit(float ox,float oy,float *n){  
	    	float dx=ox-x;
	      	float dy=oy-y;
	    	if(dx*dx+dy*dy<radius*radius){
	    	   float dz=sqrtf(radius*radius-dx*dx-dy*dy);
	    	   *n=dz/sqrtf(radius*radius);
	    	   return dz+z;
    	    }
    	    return -INF;
    	 }
    };

__global__ void kernel(unsigned char*ptr ,Sphere *s);

int main(){
	Sphere *s;
	cudaEvent_t start,stop;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);
	cudaEventRecord(start,0);

	Bitmap bitmap(DIM ,DIM);
	unsigned char *dev_bitmap;
	cudaMalloc((void**)&dev_bitmap,bitmap.image_size());  //给dev_bitmap分配内存
	cudaMalloc((void **)&s,sizeof(Sphere)*SPHERES);  //给s分配内存

	Sphere *temp_s=(Sphere*)malloc(sizeof(Sphere)*SPHERES);//在CPU上先给s一个拷贝
	for(int i=0;i<SPHERES;i++){  //生成20个随机数组给temp_s赋值
		temp_s[i].r=rnd(1.0f);
		temp_s[i].g=rnd(1.0f);
		temp_s[i].b=rnd(1.0f);
		temp_s[i].x=rnd(1000.0f)-500;
		temp_s[i].y=rnd(1000.0f)-500;
		temp_s[i].z=rnd(1000.0f)-500;
		temp_s[i].radius=rnd(100.0f)+20;
	}

	cudaMemcpy(s,temp_s,sizeof(Sphere)*SPHERES,cudaMemcpyHostToDevice);
	free(temp_s);

	//数据都分配到GPU上了，接下来启动核函数
	dim3 grids(DIM/16,DIM/16);
	dim3 threads(16,16);
	kernel<<<grids,threads>>>(dev_bitmap,s);

	cudaMemcpy(bitmap.get_ptr(),dev_bitmap,bitmap.image_size(),cudaMemcpyDeviceToHost);

	cudaEventRecord(stop,0);
	cudaEventSynchronize(stop);
	float elapsedTime;
	cudaEventElapsedTime(&elapsedTime,start,stop);
	printf("time to generate:%3.1f 毫秒",elapsedTime);
	cudaEventDestroy(start);
	cudaEventDestroy(stop);

	bitmap.display_and_exit();
	cudaFree(dev_bitmap);
	cudaFree(s);
}


__global__ void kernel(unsigned char*ptr,Sphere *s){
	int x=threadIdx.x+blockDim.x*blockIdx.x;
	int y=threadIdx.y+blockDim.y*blockIdx.y;
	int offset=x+y*blockDim.x*gridDim.x;
	float ox=(x-DIM/2);
	float oy=(y-DIM/2);

	float r=0,g=0,b=0;
	float maxz=-INF;
	for(int i=0;i<SPHERES;i++){
		float n;
		float t=s[i].hit(ox,oy,&n);
		if(t>maxz){
			float fscale=n;
			r=s[i].r*fscale;
			g=s[i].g*fscale;
			b=s[i].b*fscale;
			maxz=t;
		}
	}

	ptr[offset*4+0]=(int)(r*255);
	ptr[offset*4+1]=(int)(g*255);
	ptr[offset*4+2]=(int)(b*255);
	ptr[offset*4+3]=255;
}

 

        程序思路大概如下：数据结构Sphere用来保存球的信息（颜色，位置，半径），Sphere中的hit用来计算（ox,oy）的像素的光线是否与这个球面相交。如果相交，那么将这个方法计算从相机到光线命中球面处的距离。然后随机产生一些Sphere，拷贝到GPU上。kernel是一个简单的光线跟踪模型，hit来判断像素点是否看到Sphere，如果比相机更近，那么作为新的最接近球面。

        运行结果：

 

 

 

接下来是使用常量的版本。这里讲球Sphere数组设置为常量，常量需要使用需要注意的问题在程序中有说明。

#include "cuda_runtime.h"  
#include "device_launch_parameters.h"  
#include "device_functions.h" 
#include <stdio.h>
#include <math.h>
#include "bitmap.h"

#define INF 2e10f
#define SPHERES 200
#define DIM 1024
#define rnd( x )  (x*rand()/RAND_MAX)

struct Sphere{
	float r,b,g;
	float radius;
	float x,y,z;
	__device__ float hit(float ox,float oy,float *n){
	    	float dx=ox-x;
	      	float dy=oy-y;
	    	if(dx*dx+dy*dy<radius*radius){
	    	   float dz=sqrtf(radius*radius-dx*dx-dy*dy);
	    	   *n=dz/sqrtf(radius*radius);
	    	   return dz+z;
    	    }
    	    return -INF;
    	 }
    };

__global__ void kernel(unsigned char*ptr); //常量这里不需要将指针传递到函数里面去

__constant__ Sphere s[SPHERES];  //常量定义在外面

int main(){
	//Sphere *s;
	cudaEvent_t start,stop;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);
	cudaEventRecord(start,0);

	Bitmap bitmap(DIM ,DIM);
	unsigned char *dev_bitmap;
	cudaMalloc((void**)&dev_bitmap,bitmap.image_size());
	//常量不需要分配内存，变量需要分配内存

	Sphere *temp_s=(Sphere*)malloc(sizeof(Sphere)*SPHERES);
	for(int i=0;i<SPHERES;i++){
		temp_s[i].r=rnd(1.0f);
		temp_s[i].g=rnd(1.0f);
		temp_s[i].b=rnd(1.0f);
		temp_s[i].x=rnd(1000.0f)-500;
		temp_s[i].y=rnd(1000.0f)-500;
		temp_s[i].z=rnd(1000.0f)-500;
		temp_s[i].radius=rnd(100.0f)+20;
	}

	cudaMemcpyToSymbol(s,temp_s,sizeof(Sphere)*SPHERES);  //常量拷贝过程用的函数不一样
	free(temp_s);

	dim3 grids(DIM/16,DIM/16);
	dim3 threads(16,16);
	kernel<<<grids,threads>>>(dev_bitmap);

	cudaMemcpy(bitmap.get_ptr(),dev_bitmap,bitmap.image_size(),cudaMemcpyDeviceToHost);

	cudaEventRecord(stop,0);
	cudaEventSynchronize(stop);
	float elapsedTime;
	cudaEventElapsedTime(&elapsedTime,start,stop);
	printf("time to generate:%3.1f 毫秒",elapsedTime);
	cudaEventDestroy(start);
	cudaEventDestroy(stop);

	bitmap.display_and_exit();
	cudaFree(dev_bitmap);
	//常量不需要free
}

__global__ void kernel(unsigned char*ptr){
	int x=threadIdx.x+blockDim.x*blockIdx.x;
	int y=threadIdx.y+blockDim.y*blockIdx.y;
	int offset=x+y*blockDim.x*gridDim.x;
	float ox=(x-DIM/2);
	float oy=(y-DIM/2);

	float r=0,g=0,b=0;
	float maxz=-INF;
	for(int i=0;i<SPHERES;i++){
		float n;
		float t=s[i].hit(ox,oy,&n);
		if(t>maxz){
			float fscale=n;
			r=s[i].r*fscale;
			g=s[i].g*fscale;
			b=s[i].b*fscale;
		}
	}

	ptr[offset*4+0]=(int)(r*255);
	ptr[offset*4+1]=(int)(g*255);
	ptr[offset*4+2]=(int)(b*255);
	ptr[offset*4+3]=255;
}

 

         __constant__把变量的访问限制为只读。与全局内存中读取数据相比，从常量内存中读取相同的数据可以节约内存带宽，因为：对常量内存的单次读取可以广播到其他的“近邻”线程；常量内存的数据将缓存起来，因此对相同地址的连续读操作将不会产生额外的内存通信量。
        这里所谓的“近邻”涉及到“线程束（wrap）”的概念。线程束可以看成是一组线程通过交织而形成的一个整体。CUDA中，线程束是指一个包含32个线程的集合，这个线程集合交织到一起并以步调一致的形式执行。在程序中的每一行，线程束中的每一个线程都将在不同的数据上执行相同的命令。
         当处理常量时，硬件将单次内存读取广播到每半个线程束（包含16个线程），这样，需要的内存流量大概是原始方式的1/16。硬件主动将这个常量数据存在GPU上，在第一次读取常量后，当其他半线程束请求同一数据的地址时，就会命中缓存。

        如果16个线程读取相同地址时候能够极大提升性能，如果16个线程分别读取不同地址可能会降低性能。

 

这里的cudaEvent用来统计执行时间。cudaEventRecord()视为一条记录当前时间的语句。当cudaEventSynchronize()返回时，stop之前的工作就完成了。这时就可以计算时间戳。

CUDA中的event是直接在GPU上实现的，他们不适合计算host与device的混合计时。