发表于

前言

最近发现了coursera上比较好的几门课，是来自北京大学的程序设计与算法系列课程，于是想重新回顾一下C++这门及其注重效率的语言的一些特性。同时对自己的相关知识查漏补缺，以及做一些记录。course

目录

1， 内联成员函数和缺省参数的函数
2， 构造函数

内联成员函数和缺省参数的函数

之前我们知道，内联函数存在的目的是，当一些比较简单又调用次数特别多的函数的时候，为了减少调用函数入栈时候的额外开销，我们设置一种内联函数，编译器在碰到内联函数的时候，直接把函数的代码粘过来运行。这是C++这样一种追求极致效率的语言的一种特殊机制。
那么内联成员函数的声明方法：
1， inline + 成员函数
2，整个函数体出现在类定义内部，也就是说当你把函数体写在类定义里面时，这些默认为内联函数。

1
2
3
4
5
6
class B{
    inline void func1();
    void func2(){};
    void init(int x = 1, int y = 2);
    private: int x, int y;
}

上面的func1()和func2()都是内联函数。

而缺省参数的路由是在函数的参数列表中已经给一些参数一个默认值，如之上的init()函数。那么当你调用函数的参数不够时，剩下的（靠近右边的）用默认参数。那么这种设定有什么用呢？是当你在开发较为简单的功能的时候，考虑到程序的可扩展性，即之后可能还会将某些现在是固定值的变量变成一个可控变量。则可以用到这种技术。

用重载函数和缺省参数的函数时，一定要避免函数的二义性。

构造函数

1， 一个类可以有多个构造函数。
2， 当定义了构造函数，则编译器不生成默认的无参数的构造函数（肯定会有一个构造函数）
3， 不分配内存，是在new一个对象时对这个对象进行作用。

复制构造函数

只有一个参数，即对同类对象的引用。形如X::X( X& )或者X::X(const X& )。
如果没有定义复制构造函数，那么编译器生成默认的复制构造函数。（所以当我们不写构造函数的时候，编译器会至少为我们写两个构造函数）而当你自己写了一个复制构造函数的时候，编译器就不会帮你生成复制构造函数了。

注意：
1，当初始化的时候调用 A a1 = a2;这一句是初始化语句。a1 = a2这一句是赋值语句，体会这两者的不同。
2， 当一个函数中有一个参数是类A的对象，那么当该函数被调用时，类A的复制构造函数将会被调用（对形参）。
这个时候如果你自己写了这个类的复制构造函数，那么可能这个时候实参和形参没有经历过复制工作。
3，如果函数的返回值是类A的对象时，则函数返回式，A的复制构造函数被调用。

在什么时候要自己写复制构造函数。（to be continued..）
引用作为参数时是没有生成对象的，也即不会有复制的过程，也不会调用复制构造函数。

类型转换构造函数

目的：实现类型的自动转换。
特点：只有一个参数，不是复制构造函数。
编译器自动调用->转换构造函数->建立一个临时对象。

析构函数

名字与类名相同。
前面加上~
没有返回值，没有参数，只有一个。

在对象消亡时自动被调用。用于做“善后”工作，释放内存等。

注意：对象数组的生命期结束时，每个元素的析构函数都会被调用。而对象的生命作用域是离他最近的一对花括号。当跳出作用域时，就会调用析构函数。当然，static静态对象和全局对象则在最后依次被消亡。

静态成员变量和静态成员函数

1，静态成员对象和普通成员对象的本质区别是，普通成员变量每个对象有格子的一份，而静态成员变量一共就一份，为所有对象共享。这也是sizeof不会算static变量的原因。
2，普通成员函数必须作用于某个对象，而静态成员函数并不具体作用于某个对象。
3，因此静态成员不需要通过对象就能访问。

如果访问静态成员？

不管怎么访问，但它都不属于某个特定的对象。

静态成员变量本质上是全局变量，哪怕一个对象都不存在，类的静态成员变量也存在；
静态成员函数本质上是全局函数；
设置静态成员这种机制的目的是将和某些类紧密相关的全局变量和函数写到类里面，看上去像一个整体，易于维护和理解。（换句话说是对全局变量分类）

注意：
1，在C++中，必须在定义类的文件中对静态成员变量进行一次说明或初始化，否则编译能通过链接不能通过。
2，在静态成员函数中，不能访问非静态成员变量，也不能调用非静态成员函数。

成员对象和封闭类

成员对象：一个类的成员变量是另一个类的对象。
那么包含成员对象的类叫做封闭类（Enclosing）.
封闭类必须要自己定义构造函数，否则就会编译报错。所以要通过初始化列表来实现。

初始化列表
类名::构造函数（参数表）:成员变量1(参数表)，成员变量2(参数表),…
{
}
而成员对象初始化列表中的参数.1可以是任意复杂的表达式.2函数/变量等。

当封闭类对象生成时，先执行所有成员对象的构造函数，再执行封闭类的构造函数。（顺序和成员对象在类中的说明顺序一致，和初始化列表无关。）
当封闭类对象消亡时，先执行封闭类的析构函数，再执行所有成员对象的析构函数。
析构函数顺序和构造函数的调用顺序相反

this指针

作用：指向成员函数所作用的对象。是一种特殊的指针，相当于指向一个对象入口的指针。
注意，在静态成员函数中是不能使用this指针的。

常量对象，常量成员函数

如果不希望某个对象的值被被改变，则定义该对象的时候可以再前面加const关键字。
而常量成员函数在执行期间不应该修改其所作用的对象。因此常量成员函数不能修改成员变量的值（静态成员变量除外），也不能调用同类的非常量成员函数（静态成员函数除外）。

注意：如果两个名字一样的函数一个定义了const一个没有定义，那么算是重载而不是重复定义。分别使用在const对象和非const对象上。
在引用作为参数的是很常用的，而在有些情况下，为了确保对象不会被改变可以使用常引用。

发表于

枚举，基于已有的知识进行答案猜测的一种问题求解策略。
从可能的集合中一一列举各元素，根据已经知道的知识，给一个猜测的答案。

Ex.以找到比N小的最大素数为例

枚举算法

对问题可能解集的每一项根据问题给定的检验条件判定哪些是成立的，使条件成立的即是问题的解。

枚举过程

1, 判断猜测的答案是否正确。（例如，2是小于N的最大素数吗）
2, 进行新的猜测：注意要保证，（1）猜测的结果必须是前面的猜测中没有出现过的。（每次猜测的素数要比之前的要大）（2）猜测的过程中要尽早的排除错误答案。（除2之外，只有奇数才可能是素数。）

所以要考虑到三个关键问题。

1， 给出解空间，建立简洁的数学模型。
2， 减少搜索的空间。 利用知识缩小模型中各个变量的范围。
3， 采用合适的搜索顺序。

先附上之前问题的解法

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
//找到小于N的最大素数
//用这种枚举的方法，要比去判断n-1, n-2...是否是质数要快。
int findNum(const int N)
{
	int ret;
	std::vector<int> prim;
	prim.push_back(2);
	bool flag = true; //标识当前的数是否是质数
	for (int i = 3; i < N; i = i + 2/* 排除偶数*/)
	{
		for (int j = 0; j < prim.size() - 1; ++j)
		{
			if (i % prim[j] == 0)
			{
				flag = false;
				break;
			}
		}
		if (flag) prim.push_back(i);//如果当前i是质数，把它Push到质数数组里
		flag = true;
	}
	return prim[prim.size() - 1];
}

例一 熄灯问题

问题描述

有一个由按钮组成的矩阵，其中每行有6个按钮，共5行，每个按钮上有一盏灯。当你按下一个按钮后，该按钮以及其周围上下左右四盏灯都会改变，亮-->暗或者暗-->亮。
要求对矩阵中的每盏灯设置一个初始状态，请你写一个程序，确定需要按下哪些按钮，恰好使得所有的灯都熄灭。

ji
如果穷举所有的状态，即为2^(56)种，肯定是不现实的。
*那么此时穷举的技巧是，如果存在某个局部，一旦这个局部的状态被确定，那么剩余其他部分的状态只能是确定的一种，或者不多的n种，那么就只需要枚举这个局部的状态就可以。
本题是有这样的“局部”的，当第一行的开关状态确定后，如果要熄灭第一行的所有灯，只能依靠第二行在第一行亮灯的对应位置下Press,故在第一行状态后，第二行的状态实际上已经确定好了；依次类推，第三行第四行…所有有一个思路是枚举第一行的所有的状态，然后判断第五行是否恰好都熄灭。以下是代码：

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
//-------------------------------------------------
//**熄灯问题**
//枚举第一行的所有情况 2^6 = 64种。
//对于每一种情况，下面的1-4行结果唯一，最后判断第五行是否正好都别熄灭，若是，则得到正确解。
int press[6][8], puzzle[6][8];
bool guess()
{
	for (int i = 1; i < 7; ++i)
	{
		for (int j = 0; j < 5; ++j)
		{
			//根据Press第一行和puzzle数组，计算press其它行的值
			press[i + 1][j] = (puzzle[i][j] + press[i - 1][j] + press[i][j] + press[i][j + 1] + press[i][j - 1]) % 2;
		}
	}
	for (int i = 1; i < 7; ++i)
	{
		//判断所计算的press数组能否熄灭第五行的所有灯
		if ((puzzle[5][i] + press[5][i] + press[5][i - 1] + press[5][i + 1]) % 2 == 1)
			return false;
	}
	return true;
}
void enumerate()
{
	int c;
	for (c = 1; c < 7; ++c)
	{
		press[1][c] = 0;
	}
	while (!guess())
	{
		//对press第一行的元素的各种取值进行枚举，依次考虑
		c = 1;
		press[1][c]++;
		//累加进位
		while (press[1][c] > 1)
		{
			press[1][c] = 0;
			c++;
			press[1][c] += 1;
		}
	}
	
	return;
}
int main(int argc, char** argv)
{
	int cases, i, r, c;
	scanf("%d", &cases); 
	//将扩大的矩阵的边缘用0填充
	for (r = 0; r < 6; r++)
	{
		press[r][0] = 0;
		press[r][7] = 0;
	}
	for (c = 0; c < 8; c++)
	{
		press[0][c] = 0;
	}
	for (i = 0; i < cases; ++i)
	{
		for (r = 1; r < 6; ++r)
		{
			for (c = 1; c < 7; ++c)
			{
				scanf("%d", &puzzle[r][c]);
			}
		}
		enumerate();
		
		printf("PUZZLE# %d\n", i + 1);
		for (r = 1; r < 6; ++r)
		{
			for (c = 1; c < 7; ++c)
			{
				printf("%d ", press[r][c]);
			}
			printf("\n");
		}
	}
	return 0;
	
	
}

讨厌的青蛙

问题描述
有一种青蛙，会跳跃稻田，从而踩踏稻子。青蛙总是沿着一条直线跳跃稻田，并且每次跳跃的距离都相同。（可竖直也可倾斜）。来看看标准输入输出。



那么，我们枚举的是什么？肯定不能枚举所有可能出现的情况。
我们知道，两点确定一条直线，所以基本思路就是枚举第一个点和第二个点。
先直接粘贴代码，注释中有比较详细的思路。总结在之后。

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
#include <stdio.h>
#include <stdlib.h>
#include <algorithm>
//-------------------------------------------------
//**讨厌的青蛙**
//定义了一个结构体用于描述踩踏水稻的横纵坐标
struct Plant
{
	int x;
	int y;
};
Plant plants[5000];
int r, c, n;
//申请一个动态的空间来存放稻田是否被踩踏
int **field;
int searchPath(int dx, int dy, Plant p)
{
	int step = 1;
	int x = p.x;
	int y = p.y;
	while (x >= 1 && x <= c && y >= 1 && y <= r)
	{
		step++;
		//这一步也可以用二叉查找来找plants里面有没有x, y相对应的点
		if(field[x][y] != 1 /* std::!binary_search(plants, plants + n, plants)*/) return 0;
		x += dx;
		y += dy;
	}
	return step;
}
void input()
{
	scanf("%d %d %d", &r, &c, &n);
	field = new int*[r + 1];
	for (int i = 0; i < r + 1; ++i)
	{
		field[i] = new int[c + 1];
	}
	//先初始化其为0
	for (int i = 1; i < r + 1; ++i)
	{
		for (int j = 1; j < c + 1; ++j)
		{
			field[i][j] = 0;
		}
	}
	for (int i = 0; i < n; ++i)
	{
		scanf("%d %d", &plants[i].x, &plants[i].y);
		field[plants[i].x][plants[i].y] = 1;
	}
	for (int i = 1; i < r + 1; ++i)
	{
		for (int j = 1; j < c + 1; ++j)
		{
			printf("%d ", field[i][j]);
		}
		printf("\n");
	}
}
int hatingFlog()
{
	int dx, dy;
	int max = 2;
	int tmpStep;
	input();
	std::sort(plants, plants + n);
	for (int i = 0; i < n - 1; ++i)
	{
		for (int j = i + 1; j < n; ++j)
		{
			dx = plants[j].x - plants[i].x;
			dy = plants[j].y - plants[i].y;
			//如果第一个点的前一个点还在field内，则肯定不是最优答案的起点，换第二个点。
			if (plants[i].x - dx >= 1 && plants[i].x - dx <= c && plants[i].y - dy >= 1 && plants[i].y <= r)
				continue;
			//如果第一个点经过比max还小的step就越界了，肯定也不可能满足条件。
			//由于是按x排好序的，再往后换第二个点必然也不满足条件，所以换第一个点。
			if (plants[i].x + dx * (max - 1) > r)
				break;
			//如果y方向过早越界
			if (plants[i].y + dy * (max - 1) > c || plants[i].y + dy * (max - 1) < 1)
				continue;
			//如果三种提前的排除无法排除该店，那么算这条路径的step，看其是否为最大。
			int tmpStep = searchPath(dx, dy, plants[i]);
			if (tmpStep > max) max = tmpStep;
		}
	}
	if (max <= 2) max = 0;
	return max;
}
bool operator < (const Plant& p1, const Plant& p2)
{
	if (p1.x == p2.x) return p1.y < p2.y;
	return p1.x < p2.x;
}

思路值得我们继续巩固的点有：
1，通过思考得到“局部”枚举的方式。
2，对标准库里面函数的运用，sort, 包括binary_search的运用；当然，有很多时候我们对一些复杂的数据结构进行排序时，需要我们自己进行操作符重载。

发表于

##前言

由于一直在学习图形学，很多时候，图形学中的计算，如最近碰到的问题，计算随机点的中垂面，每个点的计算K近邻树，都是独立的。如果用GPU并行加速，可能会有好的提速效果。再加上导师最想希望我和实验室学长开始做一个基于Cuda的项目，故想开始学CUDA。

结构目录

  • Cuda是什么
  • Cuda程序的编译过程
  • 我的第一个Cuda程序
  • 利用GPU加速Matrix Multiplication
  • 优化策略1 : shared memory
  • 优化策略2 : texture memroy
  • 总结

Cuda是什么

CUDA(Compute Unified Device Architecture)，是显卡厂商NVIDIA推出的运算平台。 CUDA™是一种由NVIDIA推出的通用并行计算架构，该架构使GPU能够解决复杂的计算问题。而我们关心的事情有2点：

  1. Cuda是对GPU运算能力的并行利用；
  2. Cuda很好的兼容了C/C++，在此之上做了扩展；

听到过一个不错的比喻：　
CPU可以类比于一个强壮的男人，随机现在硬件能力的不断提升，单个CPU已经拥有了很强大的能力，甚至CPU也有了多核，同样也有很多基于CPU的多线程编程。
而GPU则是一群小孩，每一个能力不大，但是由于GPU可以同时开启成千上万个线程，<线程格，线程块，线程>，因此可以同时做很多事情，极大的提升程序的运行效率，这就是并行的好处。

tips : GPU是不适合做太多的逻辑判断的。它更希望能做一些简单的计算工作。

Cuda程序的编译过程

Cuda程序一般分为两个部分，其中一部分是用Nvidia的编译器进行编译，在GPU上运行；另一部分，用原本的VS（我使用的编译器是Visual Studio）编译，在CPU上跑。
Nvcc的编译过程
Cuda C是支持C/C++语言的。它只对C语言做了一个很小的扩展并且提供了一个C runtime library.
想要知道Cuda是怎么运行的，我们首先要知道Cuda程序的编译过程。

NVCC 的工作流主要分下面几步：

1，将程序中的host code 和 device code 区别开来。
2，将device code进行转化可装配形式（assembly form (PTX code)），进而转化成2进制流，用于交给GPU处理。
3，将Host code中不符合C语言标准的代码进行替换，然后按照正常的编译过程进行编译链接，在CPU中处理。

Kernel 核函数

核函数是作用在GPU上，用____global____定义，相当于是cuda提供给C/C++的一个接口，每次可以同时调用多个核函数，每个线程都会执行这个核函数，线程块和线程的数量可以自己根据计算的规模和硬件条件来设定。核函数是cuda编程的核心，在之后的例子中我们也能看到。

Initialization 初始化过程：

一般来说，Cuda并没有一个明确的开始标识，当第一个runtime function被调用的时候，GPU section就被初始化了。在初始化的过程中，会创建一个cuda contexet，初始创建的这个cuda环境是主要环境，被所有的host所共享。*cudaDeviceReset()*可以destory当前的上下文，直至下一个runtime function被调用时，将重新创建primary context。

Device Memory 设备内存：

设备内存相当于是GPU中的内存，一般的分配方法有两种，线性内存或者cuda arrays，其中cuda arrays是和纹理内存相关的，我们后续再谈。
线性内存通常使用*cudaMalloc()*分配，*cudaFree()*释放，*cudaMemcpy()*在主机和设备之间传递数据。

我的第一个Cuda程序

给出一个最基本的例子，也是在Visual Studio上创建一个.cu文件时，会自动生成的一个示例代码。

输入：两个长度一样的数组，假设长度为5；
输出：一个结果数组，数组里面的每一个数是2个输入数组里对应位置的和。

CPU思路，这是一道简单得不能再简单的题，通常CPU中的算法必然是定义一个循环，遍历2个输入数组并将他们对应元素的和相加，得到结果数组。

1
2
3
4
for(int i = 0; i < 5; ++i)
{
    _c[i] = a[i] + b[i];
}

可以看到，CPU中进行的是五次串行运算。

那么回到我们想要实现的GPU并行计算。之前也提到过，GPU是通过调用核函数来调用进程的。

1
2
3
4
5
6
7
8
9
 // Launch a kernel on the GPU with one thread for each element.
  addKernel<<<1, 5>>>(dev_c, dev_a, dev_b);
    
    __global__ void addKernel(int *c, const int *a, const int *b)
{
    int i = threadIdx.x;
    c[i] = a[i] + b[i];
}

可以大致想象，我们开启了5个线程，然后通过thread.x来获取所在线程的信息，同时计算结果数组C，以此来实现并行的效果。当然，在GPU中需要单独分配和释放内存，同时数据也会在CPU和GPU中来回传递，这个Cuda都已经为我们提供好了相应的接口，后续我们会继续介绍。

利用GPU加速Matrix Multiplication

1，矩阵乘法的CPU思路

在这里我是默认读者已经有基本的矩阵运算基础的。所以不难想到，对于CPU的矩阵乘法算法如下：
由于，对于C中的每一个数，都要经过n次乘法运算，因此总的时间复杂度应该是O(n^3)。
对于计算两个1024*1024的输入方阵，CPU串行计算时间大概在9000ms左右。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
//CPU version 其中C是结果矩阵，AB是输入矩阵。
void MatrixMulCPU(float *_C, const float* _A, const float* _B, int WA, int HA, int WB, int HB)
{
	if (WA != HB)
	{
		printf("the matrix A and B cannot be multipled!");
			exit(0);
	}
	for (int i = 0; i < HA; ++i)
	{
		for (int j = 0; j < WB; ++j)
		{
			for (int k = 0; k < WA; ++k)
			{
				_C[i * WA + j] += _A[i * WA + k] * _B[k * WB + j];
			}
		}
	}
}

2，GPU优化思路。

借矩阵乘法的GPU算法，我在这里完整地说一下GPU调用核函数的方法。
首先申请的是CPU上的内存，也即是host memory。

1
2
3
4
5
6
7
8
   const int width_A = 1024;
const int height_A = 1024;
const int width_B = 1024;
const int height_B = 1024;
float *B = (float *)malloc(sizeof(float) * height_B * width_B);
float *A = (float *)malloc(sizeof(float) * height_A * width_A);
float *C = (float *)malloc(sizeof(float) * height_A * width_B);

我是给了A,B矩阵0-100的随机数。

1
2
3
4
5
//产生随机数生成器
	srand((unsigned)time(0));
	randomInit(B, height_B * width_B);
	randomInit(A, height_A * width_A);

之后申请device memory，并且把A, B所在内存拷贝过GPU上。

1
2
3
4
5
6
7
8
9
10
11
12
   float *dev_a = 0;
float *dev_b = 0;
float *dev_c = 0;
cudaError_t cudaStatus;
// Choose which GPU to run on, change this on a multi-GPU system.
   cudaStatus = cudaSetDevice(0);
   
   cudaStatus = cudaMalloc((void**)&dev_c, HA * WB * sizeof(float));
   
   cudaStatus = cudaMemcpy(dev_a, a, HA * WA * sizeof(float), cudaMemcpyHostToDevice);
   cudaStatus = cudaMemcpy(dev_b, b, HB * WB * sizeof(float), cudaMemcpyHostToDevice);

到这一步后，我们的准备工作就算是做好了，之后的工作是自己来设定，我们需要在GPU上开启多少个线程，这个你可以把它也视为核函数的一个参数，并且把这些device上的内存地址当作参数传递给核函数。

在这个基本的GPU算法中，我们为C的每一行每一列都设置一个线程。

1
2
3
4
5
6
7
8
//为每一个C[i][j]设置一个线程进行计算
int block_size = 16;
dim3 Threads(block_size, block_size);
dim3 Blocks(WB / block_size, HA / block_size);
//用这个方式来调用核函数
MatrixMulGPU_1 << <Blocks, Threads >>>(dev_c, dev_a, dev_b, WA, WB);

核函数，我们通过blockIdx和threadIdx索引获得当前是哪个线程，来计算它所负责的部分。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
__global__ void MatrixMulGPU_1(float *c, const float *a, const float *b, unsigned int WA, unsigned int WB)
{
	float sum = 0;
	//找出该线程所在的行和列
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	//线程Thread(row, col)负责计算C(row, col)
	for (int i = 0; i < WB; ++i)
	{
		sum += a[row * WA + i] * b[i * WB + col];
	}
	c[row * WB + col] = sum;
}

可以看到，每个线程只进行了N次运算，而每个线程之前是并行的，所以时间理论上只是O(n)，虽然在拷贝内存上还需要花费一些时间，但是从O(n^3)->O(n)这个提升也是巨大的，对时间的统计数据显示也是如此。从9000ms提升到50ms左右。

优化政策1，shared memory

那么在实现了基本的矩阵GPU算法后，我们还有没有什么方法继续优化我们的算法，使得速度进一步提高呢？结果是肯定的。在之前的算法中，我们的每一个C[i][j]都算了n次，而在这n次取数的过程中，还有我们可以做的很多事情。首先我们可以对每一行每一列进行细分，划分成更细的block_size大小（可以自己来设定），然后再计算每个block_size*block_size时，我们设计一个共享区，然后让数的计算存这个共享区中取数，就可以减少取数的时间。下面是相应的代码，可以在英伟达的标准实例代码中找到：

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
template<int BLOCK_SIZE> __global__ void MatrixMulGPU_2(float *c, const float *a, const float *b, unsigned int WA, unsigned int WB)
{
	// Block index
	int bx = blockIdx.x;
	int by = blockIdx.y;
	// Thread index
	int tx = threadIdx.x;
	int ty = threadIdx.y;
	// Index of the first sub-matrix of A processed by the block
	int aBegin = WA * BLOCK_SIZE * by;
	// Index of the last sub-matrix of A processed by the block
	int aEnd = aBegin + WA - 1;
	// Step size used to iterate through the sub-matrices of A
	int aStep = BLOCK_SIZE;
	// Index of the first sub-matrix of B processed by the block
	int bBegin = BLOCK_SIZE * bx;
	// Step size used to iterate through the sub-matrices of B 
	int bStep = BLOCK_SIZE * WB;
	// Csub is used to store the element of the block sub-matrix
	// that is computed by the thread
	float Csub = 0;
	// Loop over all the sub-matrices of A and B
	// required to compute the block sub-matrix
	for (int i = aBegin, j = bBegin;
		i <= aEnd;
		i += aStep, j += bStep)
	{
		// Declaration of the shared memory array As used to
		// store the sub-matrix of A
		__shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
		// Declaration of the shared memory array Bs used to
		// store the sub-matrix of B
		__shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];
		// Load the matrices from device memory
		// to shared memory; each thread loads
		// one element of each matrix
		As[ty][tx] = a[i + WA * ty + tx];
		Bs[ty][tx] = b[j + WB * ty + tx];
		// Synchronize to make sure the matrices are loaded
		__syncthreads();
		// Multiply the two matrices together;
		// each thread computes one element
		// of the block sub-matrix
#pragma unroll
		for (int k = 0; k < BLOCK_SIZE; ++k)
		{
			Csub += As[ty][k] * Bs[k][tx];
		}
		// Synchronize to make sure that the preceding
		// computation is done before loading two new
		// sub-matrices of A and B in the next iteration
		__syncthreads();
	}
	// Write the block sub-matrix to device memory;
	// each thread writes one element
	int k = WB * BLOCK_SIZE * by + BLOCK_SIZE * bx;
	c[k + WB * ty + tx] = Csub;
}

其中syncthreads()是用于保证同步的，每一个线程在运行到这一步时，都会停下来等待，直至所有进程都进行到这才会继续下去，以此来保证同步。

优化政策2， Texture Memory

纹理这一块，我的理解是，一般的内存都是线性存放的，不论是二维或者是三维，但存放的本质就是串行线性。故在访问例如 a[1][2]和a[1][3]时，比访问a[1][2]和a[2][2]时，后者的耗时要多一些。但对于纹理，就不会出现这样的问题，以此来提高访问效率。

首先，先声明纹理

1
2
texture<float, 2, cudaReadModeElementType> texA;
texture<float, 2, cudaReadModeElementType> texB;

然后以cudaArray为介质，将内存中的数放到纹理中

1
2
3
4
5
6
7
8
9
10
11
12
13
//GPU mode3 with texture memory
		cudaChannelFormatDesc channelDescA = cudaCreateChannelDesc((int)sizeof(float) * 8, 0, 0, 0, cudaChannelFormatKindFloat);
		cudaChannelFormatDesc channelDescB = cudaCreateChannelDesc((int)sizeof(float) * 8, 0, 0, 0, cudaChannelFormatKindFloat);
		cudaArray* mat_A;
		cudaArray* mat_B;
		cudaMallocArray(&mat_A, &channelDescA, width_A, height_A);
		cudaMallocArray(&mat_B, &channelDescB, width_B, height_B);
		cudaMemcpyToArray(mat_A, 0, 0, A, sizeof(float) * height_A * width_A, cudaMemcpyHostToDevice);
		cudaMemcpyToArray(mat_B, 0, 0, B, sizeof(float) * height_B * width_B, cudaMemcpyHostToDevice);

接着设置纹理的相关属性

1
2
3
4
5
6
7
8
texA.addressMode[0] = cudaAddressModeWrap;
texA.addressMode[1] = cudaAddressModeWrap;
texA.filterMode = cudaFilterModePoint;
texA.normalized = false;
texB.addressMode[0] = cudaAddressModeWrap;
   texB.addressMode[1] = cudaAddressModeWrap;
texB.filterMode = cudaFilterModePoint;
texB.normalized = false;

最后纹理的相关核函数如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
__global__ void MatrixMul(float *c, unsigned int w, unsigned int h)
{
	float sum = 0;
	//找出该线程所在的行和列
	int row = blockIdx.y * blockDim.y + threadIdx.y;
	int col = blockIdx.x * blockDim.x + threadIdx.x;
	/*
	//计算纹理坐标
	float u = row / (float)w;
	float v = col / (float)h;
	//线程Thread(row, col)负责计算C(row, col)
	u -= 0.5f;
	v -= 0.5f;
	*/
	for (int i = 0; i < w; ++i)
	{
		sum += tex2D(texA, i, row) * tex2D(texB, col, i);
	}
	c[row * w + col] = sum;
}

可以看到用2维纹理来访问我们相关的矩阵，是非常方便的，只需要用tex2D()就行。

最后是这些方法的截图。
结果截图

注：在release版本中，使用texture优化的时间比正常的GPU算法要慢，这个现象由于矩阵乘法本身对取数是 有一定规律性的，而texture优化对一些随机访问内存算法比较有效，在此我们不深究，只为学习纹理内存用。

发表于

Welcome to Hexo! This is your very first post. Check documentation for more info. If you get any problems when using Hexo, you can find the answer in troubleshooting or you can ask me on GitHub.

Quick Start

Create a new post

1
$ hexo new "My New Post"

More info: Writing

Run server

1
$ hexo server

More info: Server

Generate static files

1
$ hexo generate

More info: Generating

Deploy to remote sites

1
$ hexo deploy

More info: Deployment