cuda学习笔记2

本文主要介绍如何用cuda进行并行计算,还有一个有趣的实验:生成julia分形图。

本节要点

  • CUDA实现并行性
  • 一种错误处理的宏定义

基于GPU的矢量求和

看代码:

 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
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <math.h>
#include <helper_cuda.h>
using namespace std;

#define N 100


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

__global__ void value_init_kernel(double *a, double *b) {
	int tid = blockIdx.x;
	if (tid < N)
	{
		a[tid] = 1.0*tid;
		b[tid] = (1.0*tid*tid);
	}

}

int main(void)
{
	cudaError_t err1 = cudaSuccess, err2 = cudaSuccess, err3 = cudaSuccess;
	double a[N], b[N], c[N];
	double *dev_a, *dev_b, *dev_c;

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

	
	if (err1 != cudaSuccess || err2 != cudaSuccess || err3 != cudaSuccess)
	{
		fprintf(stderr, "Failed to allocate device value  (error code (%s,%s,%s))!\n", cudaGetErrorString(err1), cudaGetErrorString(err2), cudaGetErrorString(err3));
		exit(EXIT_FAILURE);
	}
	
	value_init_kernel <<<N, 1 >>> (dev_a, dev_b);////在GPU上赋值操作
	add_kernel <<<N, 1 >>> (dev_a, dev_b, dev_c);////在GPU上相加操作

	err1 = cudaMemcpy(a, dev_a, N * sizeof(double), cudaMemcpyDeviceToHost);
	err2 = cudaMemcpy(b, dev_b, N * sizeof(double), cudaMemcpyDeviceToHost);
	err3 = cudaMemcpy(c, dev_c, N * sizeof(double), cudaMemcpyDeviceToHost);
	
	if (err1 != cudaSuccess || err2 != cudaSuccess || err3 != cudaSuccess)
	{
		fprintf(stderr, "Failed to copy device value to host value (error code (%s,%s,%s))!\n", cudaGetErrorString(err1), cudaGetErrorString(err2), cudaGetErrorString(err3));
		exit(EXIT_FAILURE);
	}
	
	
	for (int i = 0; i < N; i++)
	{
		printf("%f + %f = %f\n", a[i], b[i], c[i]);
	}
	
	////释放GPU内存
	cudaFree(dev_a);
	cudaFree(dev_b);
	cudaFree(dev_c);


	return 0;	
	
}

代码实现在GPU上对变量赋值,然后相加返回给主机函数,上面每次对错误处理的代码太冗长了,可以用一个宏定义来简化:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
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__ ))

julia分形图案

最后是一个有趣的例子:Julia分形图案:

  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
117
118
119
#include <stdio.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include "../common/book.h" ////GPU高性能编程CUDA实战代码
#include "../common/image.h"


#include "book.h"
#include "image.h"

#define DIM1 5760
#define DIM2 5760 //每一维度的长度
#define iter_N 200

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);
	}
};

__device__ int julia(int x, int y) {
	const float scale = 1.5;
	float jx = scale * (float)(DIM1 / 2 - x) / (DIM1/ 2);
	float jy = scale * (float)(DIM2 / 2 - y) / (DIM2 / 2);

	cuComplex c(-0.8, 0.156); //-0.8,0.156;
	cuComplex a(jx, jy);

	
	for (int i = 1; i < iter_N; i++) {
		a = a * a + c;
		if (a.magnitude2() > 1000)
			return i;
	}
	return 0;
}

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

	int juliaValue = julia(x, y);
	////美工部分。。。。。。
	if (juliaValue ==0)
	{
		ptr[offset * 4 + 0] = 0;
		ptr[offset * 4 + 1] = 0;
		ptr[offset * 4 + 2] = 0;
		ptr[offset * 4 + 3] = 255;
	}

	if (juliaValue < 90 && juliaValue >= 1)
	{
		ptr[offset * 4 + 0] = (int)(255 * juliaValue / (2.0 * iter_N));
		ptr[offset * 4 + 1] = 0;
		ptr[offset * 4 + 2] = 0;
		ptr[offset * 4 + 3] = 255;
	}

	if (juliaValue < 120 && juliaValue >=90)
	{
		ptr[offset * 4 + 0] = 255;
		ptr[offset * 4 + 1] = 255 - (int)(255 * juliaValue  / (5.0 * iter_N));
		ptr[offset * 4 + 2] = 255 - (int)(255 * juliaValue / (5.0 * iter_N));
		ptr[offset * 4 + 3] = 255;
	}

	if (juliaValue < 180 && juliaValue >=120)
	{
		ptr[offset * 4 + 0] = 10;
		ptr[offset * 4 + 1] = 215;
		ptr[offset * 4 + 2] = 200;
		ptr[offset * 4 + 3] = 255;
	}

	if (juliaValue <= 255 && juliaValue >=180)
	{
		ptr[offset * 4 + 0] = (int)(255 * juliaValue / (1.0 * iter_N));
		ptr[offset * 4 + 1] = 0;
		ptr[offset * 4 + 2] = 0;
		ptr[offset * 4 + 3] = 255;
	}
	

}

struct DataBlock {
	unsigned char   *dev_bitmap;
};

int main(void) {
	DataBlock data;
	IMAGE bitmap(DIM1, DIM2);
	unsigned char *dev_bitmap;

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

	dim3 grid(DIM1, DIM2); ////实际上是DIM1*DIM2*1的三维线程格
	kernel << <grid, 1 >> > (dev_bitmap);
	HANDLE_ERROR(cudaMemcpy(bitmap.get_ptr(), dev_bitmap,
		bitmap.image_size(),
		cudaMemcpyDeviceToHost));
	HANDLE_ERROR(cudaFree(dev_bitmap));
	imwrite("C:/Users/Lenovo/Pictures/image/julia.png", bitmap.image);
	bitmap.show_image();
	
}
updatedupdated2019-12-282019-12-28