《GPU并行计算与CUDA编程》学习笔记

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16

cmake_minimum_required(VERSION 3.8)
project(CUDA_MAT_MUL LANGUAGES CXX CUDA)

aux_source_directory(cuda_src CUDA_SRC_LIST)
aux_source_directory(cuda_head CUDA_INCLUDE_LIST)
include_directories(cuda_head)
add_library(cudaPractice ${CUDA_SRC_LIST} ${CUDA_INCLUDE_LIST})
target_compile_features(cudaPractice PUBLIC cxx_std_11)

#opencv
find_package(OpenCV 3.4.3 REQUIRED)
include_directories( ${OpenCV_INCLUDE_DIRS} )

add_executable(main main.cc)
target_link_libraries(main cudaPractice)
target_link_libraries(main ${OpenCV_LIBS})

1
2
3
4
5
6
7
8

__global__ void add_kernel(double *a, double *b, double *c) {
	//block id
	int tid = blockIdx.x;
	if (tid < N)
	{
		c[tid] = a[tid] + b[tid];
	}	
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

err1 = cudaMalloc((void**)&dev_a, N * sizeof(double));
err2 = cudaMalloc((void**)&dev_b, N * sizeof(double));
err3 = cudaMalloc((void**)&dev_c, N * sizeof(double));

//表示 N 个block， 每个block分配 1个 thread
add_kernel << <N, 1 >> > (dev_a, dev_b, dev_c);////在GPU上相加操作

////用完设备指针要释放
cudaFree(dev_a);
cudaFree(dev_b);
cudaFree(dev_c);

1
2
3
4
5
6

__global__ void kernel(unsigned char *ptr) {
	int x = blockIdx.x;
	int y = blockIdx.y;
	int offset = x + y * gridDim.x;
	//...
}

1
2
3
4
5
6
7
8

unsigned char *dev_bitmap;

HANDLE_ERROR(cudaMalloc((void**)&dev_bitmap, bitmap.image_size()));

dim3 grid(DIM1, DIM2); ////实际上是DIM1*DIM2*1的三维线程格
//三维grid, 1个thread
kernel << <grid, 1 >> > (dev_bitmap);
HANDLE_ERROR(cudaFree(dev_bitmap));

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116

//thread 1D
__global__ void testThread1(int *c, const int *a, const int *b)
{
    int i = threadIdx.x;
    c[i] = b[i] - a[i];
}

//thread 2D
__global__ void testThread2(int *c, const int *a, const int *b)
{
    int i = threadIdx.x + threadIdx.y*blockDim.x;
    c[i] = b[i] - a[i];
}

//thread 3D
__global__ void testThread3(int *c, const int *a, const int *b)
{
    int i = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.x*blockDim.y;
    c[i] = b[i] - a[i];
}

//block 1D
__global__ void testBlock1(int *c, const int *a, const int *b)
{
    int i = blockIdx.x;
    c[i] = b[i] - a[i];
}

//block 2D
__global__ void testBlock2(int *c, const int *a, const int *b)
{
    int i = blockIdx.x + blockIdx.y*gridDim.x;
    c[i] = b[i] - a[i];
}

//block 3D
__global__ void testBlock3(int *c, const int *a, const int *b)
{
    int i = blockIdx.x + blockIdx.y*gridDim.x + blockIdx.z*gridDim.x*gridDim.y;
    c[i] = b[i] - a[i];
}

//block-thread 1D-1D
__global__ void testBlockThread1(int *c, const int *a, const int *b)
{
    int i = threadIdx.x + blockDim.x*blockIdx.x;
    c[i] = b[i] - a[i];
}

//block-thread 1D-2D
__global__ void testBlockThread2(int *c, const int *a, const int *b)
{
    int threadId_2D = threadIdx.x + threadIdx.y*blockDim.x;
    int i = threadId_2D+ (blockDim.x*blockDim.y)*blockIdx.x;
    c[i] = b[i] - a[i];
}

//block-thread 1D-3D
__global__ void testBlockThread3(int *c, const int *a, const int *b)
{
    int threadId_3D = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.x*blockDim.y;
    int i = threadId_3D + (blockDim.x*blockDim.y*blockDim.z)*blockIdx.x;
    c[i] = b[i] - a[i];
}

//block-thread 2D-1D
__global__ void testBlockThread4(int *c, const int *a, const int *b)
{
    int blockId_2D = blockIdx.x + blockIdx.y*gridDim.x;
    int i = threadIdx.x + blockDim.x*blockId_2D;
    c[i] = b[i] - a[i];
}

//block-thread 3D-1D
__global__ void testBlockThread5(int *c, const int *a, const int *b)
{
    int blockId_3D = blockIdx.x + blockIdx.y*gridDim.x + blockIdx.z*gridDim.x*gridDim.y;
    int i = threadIdx.x + blockDim.x*blockId_3D;
    c[i] = b[i] - a[i];
}

//block-thread 2D-2D
__global__ void testBlockThread6(int *c, const int *a, const int *b)
{
    int threadId_2D = threadIdx.x + threadIdx.y*blockDim.x;
    int blockId_2D = blockIdx.x + blockIdx.y*gridDim.x;
    int i = threadId_2D + (blockDim.x*blockDim.y)*blockId_2D;
    c[i] = b[i] - a[i];
}

//block-thread 2D-3D
__global__ void testBlockThread7(int *c, const int *a, const int *b)
{
    int threadId_3D = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.x*blockDim.y;
    int blockId_2D = blockIdx.x + blockIdx.y*gridDim.x;
    int i = threadId_3D + (blockDim.x*blockDim.y*blockDim.z)*blockId_2D;
    c[i] = b[i] - a[i];
}

//block-thread 3D-2D
__global__ void testBlockThread8(int *c, const int *a, const int *b)
{
    int threadId_2D = threadIdx.x + threadIdx.y*blockDim.x;
    int blockId_3D = blockIdx.x + blockIdx.y*gridDim.x + blockIdx.z*gridDim.x*gridDim.y;
    int i = threadId_2D + (blockDim.x*blockDim.y)*blockId_3D;
    c[i] = b[i] - a[i];
}

//block-thread 3D-3D
__global__ void testBlockThread9(int *c, const int *a, const int *b)
{
    int threadId_3D = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.x*blockDim.y;
    int blockId_3D = blockIdx.x + blockIdx.y*gridDim.x + blockIdx.z*gridDim.x*gridDim.y;
    int i = threadId_3D + (blockDim.x*blockDim.y*blockDim.z)*blockId_3D;
    c[i] = b[i] - a[i];
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54

	//testThread1<<<1, size>>>(dev_c, dev_a, dev_b);

    //uint3 s;s.x = size/5;s.y = 5;s.z = 1;
    //testThread2 <<<1,s>>>(dev_c, dev_a, dev_b);

    //uint3 s; s.x = size / 10; s.y = 5; s.z = 2;
    //testThread3<<<1, s >>>(dev_c, dev_a, dev_b);

    //testBlock1<<<size,1 >>>(dev_c, dev_a, dev_b);

    //uint3 s; s.x = size / 5; s.y = 5; s.z = 1;
    //testBlock2<<<s, 1 >>>(dev_c, dev_a, dev_b);

    //uint3 s; s.x = size / 10; s.y = 5; s.z = 2;
    //testBlock3<<<s, 1 >>>(dev_c, dev_a, dev_b);

    //testBlockThread1<<<size/10, 10>>>(dev_c, dev_a, dev_b);

    //uint3 s1; s1.x = size / 100; s1.y = 1; s1.z = 1;
    //uint3 s2; s2.x = 10; s2.y = 10; s2.z = 1;
    //testBlockThread2 << <s1, s2 >> >(dev_c, dev_a, dev_b);

    //uint3 s1; s1.x = size / 100; s1.y = 1; s1.z = 1;
    //uint3 s2; s2.x = 10; s2.y = 5; s2.z = 2;
    //testBlockThread3 << <s1, s2 >> >(dev_c, dev_a, dev_b);

    //uint3 s1; s1.x = 10; s1.y = 10; s1.z = 1;
    //uint3 s2; s2.x = size / 100; s2.y = 1; s2.z = 1;
    //testBlockThread4 << <s1, s2 >> >(dev_c, dev_a, dev_b);

    //uint3 s1; s1.x = 10; s1.y = 5; s1.z = 2;
    //uint3 s2; s2.x = size / 100; s2.y = 1; s2.z = 1;
    //testBlockThread5 << <s1, s2 >> >(dev_c, dev_a, dev_b);

    //uint3 s1; s1.x = size / 100; s1.y = 10; s1.z = 1;
    //uint3 s2; s2.x = 5; s2.y = 2; s2.z = 1;
    //testBlockThread6 << <s1, s2 >> >(dev_c, dev_a, dev_b);

    //uint3 s1; s1.x = size / 100; s1.y = 5; s1.z = 1;
    //uint3 s2; s2.x = 5; s2.y = 2; s2.z = 2;
    //testBlockThread7 << <s1, s2 >> >(dev_c, dev_a, dev_b);

    //uint3 s1; s1.x = 5; s1.y = 2; s1.z = 2;
    //uint3 s2; s2.x = size / 100; s2.y = 5; s2.z = 1;
    //testBlockThread8 <<<s1, s2 >>>(dev_c, dev_a, dev_b);

    uint3 s1; s1.x = 5; s1.y = 2; s1.z = 2;
    uint3 s2; s2.x = size / 200; s2.y = 5; s2.z = 2;
    testBlockThread9<<<s1, s2 >>>(dev_c, dev_a, dev_b);
    
    //或者：
    dim3    blocks(DIM/16,DIM/16，1); ////二维线程块
    dim3    threads(16,16,1); ////二维线程
    func_kernel<<<blocks,threads>>>（参数）；

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

__global__ void add_kernel(double *a, double *b, double *c) {
	int tid = threadIdx.x + blockIdx.x * blockDim.x;
	if (tid < N)
	{
		c[tid] = a[tid] + b[tid];
		tid += blockDim.x * gridDim.x;
 	}

}

add_kernel << < 128, 128 >> > (dev_a, dev_b, dev_c);////在GPU上相加操作

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

struct cuComplex {
	float r;
	float i;
	__device__ cuComplex(float a, float b) :r(a), i(b) {}
	__device__ float magnitude2(void) {
		return r*r + i*i;
	} ////返回复数的模的平方
	__device__ cuComplex operator*(const cuComplex& a) {
		return cuComplex(r*a.r - i*a.i, i*a.r + r*a.i);
	}
	__device__ cuComplex operator+(const cuComplex& a) {
		return cuComplex(r + a.r, i + a.i);
	}
};

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

__global__ void dot(float *a, float *b, float *c) {
	__shared__ float cache[threadsPerBlock];  ////共享内存缓存区（驻留在物理GPU上），编译器为每一个线程块生成一个共享变量的副本。同一线程块中的每个线程共享这块内存。
	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
	////归约，将共享内存区域每一个储存的值相加起来，由于规约每次迭代数量减半，要求 threadsPerBlock 是2 的指数倍。
	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];  ////把一个线程块中的最后计算得到的相加值返还给全局变量。
}

1
2

dot << <blocksPerGrid, threadsPerBlock >> >(dev_a, dev_b,
		dev_partial_c);

1

__constant__ Sphere s[num]

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21

#define SPHERES 20
__constant__ Sphere s[SPHERES];

int main(void) {
// allocate temp memory, initialize it, copy to constant
	// memory on the GPU, then free our temp memory
	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拷贝到常量内存
	HANDLE_ERROR(cudaMemcpyToSymbol(s, temp_s,
		sizeof(Sphere) * SPHERES));
	free(temp_s);
}

1
2
3
4
5
6

while (i != 0) {
		if (cacheIndex < i)
			cache[cacheIndex] += cache[cacheIndex + i];
		__syncthreads(); 
		i /= 2;
	}

1
2
3
4
5
6
7
8

while (i != 0) {
		if (cacheIndex < i)
		{
			cache[cacheIndex] += cache[cacheIndex + i];
		    __syncthreads();
		}
		i /= 2;
	}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15

cudaEvent_t     start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
*****************
在GPU上执行一些工作（包括前后的设备内存复制）
*****************
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);
float   elapsedTime;
cudaEventElapsedTime(&elapsedTime,start, stop);
printf("Time to generate:  %3.1f ms\n", elapsedTime);

cudaEventDestroy(start);
cudaEventDestroy(stop);

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

extern "C" __global__ void kernel_func(int * counter, char * data_0)
{
    // 计算线程号
    unsigned int block_index = blockIdx.x + blockIdx.y * gridDim.x + blockIdx.z * gridDim.x * gridDim.y;
    unsigned int thread_index = block_index * blockDim.x * blockDim.y * blockDim.z + \
        threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y;

    // 统计结果
    int value = data_0[thread_index];
    counter[value] ++;  
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11

extern "C" __global__ void kernel_func(int * counter, char * data_0)
{
    // 计算线程号
    unsigned int block_index = blockIdx.x + blockIdx.y * gridDim.x + blockIdx.z * gridDim.x * gridDim.y;
    unsigned int thread_index = block_index * blockDim.x * blockDim.y * blockDim.z + \
        threadIdx.x + threadIdx.y * blockDim.x + threadIdx.z * blockDim.x * blockDim.y;

    // 统计结果
    int value = data_0[thread_index];
    atomicAdd(&counter[value], 1);
}

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67

uint8_t* bgraBuffer;
uint8_t* yuvBuffer;
uint8_t* deviceBgraBuffer;
uint8_t* deviceYuvBuffer;

const int dataSizeBgra = 7680 * 4320 * 4;
const int dataSizeYuv = 7680 * 4320 * 3;

cudaMallocHost(&bgraBuffer, dataSizeBgra);
cudaMallocHost(&yuvBuffer, dataSizeYuv);
cudaMalloc(&deviceBgraBuffer, dataSizeBgra);
cudaMalloc(&deviceYuvBuffer, dataSizeYuv);

//随机生成8K的BGRA图像
GenerateBgra8K(bgraBuffer, dataSizeBgra);

//Stream的数量，这里用8个
const int nStreams = 8;

//Stream的初始化
cudaStream_t streams[nStreams];
for (int i = 0; i < nStreams; i++) {
  cudaStreamCreate(&streams[i]);
}

//计算每个Stream处理的数据量。这里只是简单将数据分成8等分
//这里不会出现不能整除的情况，但实际中要小心
int brgaOffset = 0;
int yuvOffset = 0;
const int brgaChunkSize = dataSizeBgra / nStreams;
const int yuvChunkSize = dataSizeYuv / nStreams;

//这个循环依次启动 nStreams 个 Stream
for(int i=0; i<nStreams; i++)
{
  
  brgaOffset = brgaChunkSize*i;
  yuvOffset = yuvChunkSize*i;

  //CPU到GPU的数据拷贝（原始数据），Stream i
  cudaMemcpyAsync(  deviceBgraBuffer+brgaOffset,
                    bgraBuffer+brgaOffset,
                    brgaChunkSize,
                    cudaMemcpyHostToDevice,
                    streams[i] );

  //数值计算，Stream i
  convertPixelFormat<<<4096, 1024, 0, streams[i]>>>(
                                 deviceBgraBuffer+brgaOffset, 
                                 deviceYuvBuffer+yuvOffset, 
                                 brgaChunkSize/4 );

  //GPU到CPU的数据拷贝（计算结果），Stream i
  cudaMemcpyAsync(  yuvBuffer+yuvOffset,
                    deviceYuvBuffer+yuvOffset,
                    yuvChunkSize,
                    cudaMemcpyDeviceToHost,
                    streams[i] );
}

//等待所有操作完成
cudaDeviceSynchronize();

cudaFreeHost(bgraBuffer);
cudaFreeHost(yuvBuffer);
cudaFree(deviceBgraBuffer);
cudaFree(deviceYuvBuffer);

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31

/****************************1. 并行编程的基础知识***********************************/
	//1. hello cuda
	showhello();

	//2. 矩阵相乘
	showMatMul();

	//3. 求平方
	showsquare();

	//4. 多数之和 输入参数0：使用全局内存 输入参数1：使用共享内存
	showreduce(0);

	//5. 数组扫描 更新数组，每一个数组的值是前面所有遍历过的数组的值之和 输入参数0：使用全局内存 输入参数1：使用共享内存
	showscan(1);
	//上面两个例子的global memory 都比 share memory 的时间小，不太懂为什么，是跑的数据量比较小吗

	//6. 分段直方图统计 输入参数0：错误示范 输入参数1：不分段，适用于bin数比较多的情况 输入参数2：分段，适用于bin数比较少的情况
	//分段直方图思想：假设开的线程数是THREADS_COUNT，直方图分组的数是BIN_COUNT；
	//每个线程不只处理一个数据，而是平均地处理ARRAY_SIZE/THREADS_COUNT个数据。
	//这时需要创建一个THREADS_COUNT×BIN_COUNT×sizeof(int)的共享内存，每个线程处理ARRAY_SIZE/THREADS_COUNT个数据，然后把结果放在
	//BIN_COUNT×sizeof(int)的共享内存里面；将所有数据处理完后再做一次reduce即可。
	//速度比不分段的慢，但开的线程很少，节省资源
	showhisto(1);
	/*****************************2. 图像处理中的并行编程基础**********************************/
        //7. 彩色图转灰度图
	showrgb2grey();
	//8. 均值滤波
	showblur();
	//9. 图像翻转
	showreverse();