GPGPU 加速波 PDE
4.93/5 (50投票s)
使用 GPGPU 功能的波 PDE 模拟
摘要
本文旨在利用GPGPU（GP2U）的计算能力来改进物理和科学模拟。为了让读者理解所有步骤，我们将逐步解释一个基于著名的波动方程的简单物理模拟。我们将开发一个经典的CPU实现，并最终提出一个使用GPU工作负载的版本。
波动方程PDE
PDE代表“偏微分方程”，它指的是方程的项中包含一个或多个自变量的偏导数。PDE的阶数通常由其中最高的偏导数来定义。以下是一个二阶PDE
通常，一个 n 阶、具有 m 个变量 xi (i=1,2…m) 的 PDE 表示为
一种表示
的紧凑形式是ux，同样适用于
（uxx）和（uxy）。
波动方程是一个二阶偏微分方程，用于描述水波、光波、声波等。它是电磁学、流体动力学和声学等领域的基本方程。其应用涉及模拟弦的振动或飞机运动产生的气流动力学。
通过考虑一根紧绷在x轴上两个固定端点之间的柔性弹性弦，我们可以推导出（即弦的振动）一维波动方程。
吉他弦的恢复力与其拉伸程度成正比。假设在忽略重力的情况下，我们在弦上施加一个y方向的位移（即我们稍微拉动它）。让我们只考虑它在x和x+
之间的一小段 ：
让我们写下上面图示中红色小段的
。我们假设弦的线密度为
（线密度是每单位长度的质量的度量，常用于一维物体）。回想一下，如果一个一维物体由长度为L且总质量为m的均匀物质构成，则线密度为
我们有m= 。
力（如前所述，我们忽略重力、空气阻力和其它微小力）是两端的张力T。由于弦实际上在波动，我们不能假设这两个T向量相互抵消：这就是我们正在寻找的
力。现在让我们做一些假设来简化问题：首先，我们考虑我们正在寻找的净力是垂直的（实际上它不完全是垂直的，但非常接近）。
此外，我们认为波的幅度很小。也就是说，在小弦段末端x+处的张力的垂直分量
是
。
斜率为
，如果我们考虑dx和dy作为上述图像的水平和垂直分量。由于我们认为波的幅度非常小，因此我们可以假设
。这极大地帮助了我们：
，以及两端张力产生的总垂直力变为
在极限
时，等式变得精确。
我们看到y是x的函数，但它也是t的函数：y=y(x,t)。当一个变量保持不变时，相对于另一个变量求导的标准约定（我们正在查看某一时刻的力之和）使我们能够写出
最后一步是使用牛顿第二定律并将所有内容放在一起：所有力的总和、质量（用线密度乘以段长度代替）和加速度（即仅
，因为我们只在垂直方向上移动，请记住小幅度近似）。
我们终于得到了：波动方程。使用空间拉普拉斯算子，将y函数（取决于x和t）表示为u，并代入
（一个固定常数），我们得到常见的紧凑形式
二维波动方程是通过向方程添加第二个空间项获得的，如下所示：
常数c具有距离/时间的单位，因此代表速度。我们在此不证明c实际上是波在弦上或通过表面传播的速度（尽管这肯定值得注意）。这很有意义，因为波速随着介质的张力增加而增加，并随着介质的密度而减小。
在物理模拟中，我们需要考虑除表面张力之外的其他力：真实流体表面波的平均幅度会减小。因此，我们可以通过引入一个与表面上某点的速度方向相反的力来向方程添加粘性阻尼力：
其中非负常数
代表流体的粘度（它控制波在表面平静下来的时间，小的允许波长时间存在，如水一样，而大的会导致波迅速减弱，如浓稠的油）。
使用有限差分法求解波动方程PDE
为了实际实现一个使用波动方程PDE的物理模拟，我们需要找到一种求解方法。让我们用阻尼力来求解我们最后写的方程：
其中t是时间，x和y是二维波的空间坐标，c2是波速，是阻尼因子。u=u(t,x,y)是我们需要的计算的函数：总的来说，它可以代表各种效应，例如池塘水的高度变化、电磁波中的电势等。
解决这个问题的一个常用方法是使用有限差分法。其基本思想是用有限差分代替导数，而这些导数可以在离散算法中轻松计算。如果有一个函数f=f(x)，我们想计算它关于x变量的导数，我们可以写：
如果h是一个小的离散值（非零），那么我们可以用以下方法近似导数：
这种方法的误差可以从泰勒定理导出，但这并非本文的目的。
有限差分方法使用以下三种近似方法之一来表示导数：
- 向前差分
- 向后差分
- 中心差分
让我们坚持使用后者（即中心差分）；这种方法可以泛化，因此如果我们想计算f''(x)的离散版本，我们可以先写：
然后计算f''(x)如下：
同样的方法用于估计以下偏导数：
让我们回到波动方程。我们有以下：
让我们将其应用我们刚才得到的偏导数公式（记住u=u(t,x,y)，即u依赖于t,x和y）：
这是一个相当长的表达式。我们只是用前面得到的公式替换了二阶导数和一阶导数。
如果我们现在考虑
，我们基本上强制导数计算的区间在x和y方向上相同（我们可以大大简化表达式）：
为了提高可读性，让我们用字母替换一些术语：
我们得到：
如果我们查看t变量，我们有类似的东西：
这告诉我们，新的波的高度（在t+1时）将是当前波的高度（t）加上之前的波的高度（t-1）加上现在（在t时）的值，该值仅取决于我们正在考虑的波点周围的值。
这可以被可视化为一系列时间帧，一个接一个，其中我们考虑的表面的一个点演变：
对象
有一个中心点，代表时间t时的表面上的一个点（t,x,y）。由于我们之前称为something_at_t(x,y)的项实际上是：
中心点的值受五个项的影响，最后一个项是其自身值（-4ut,x,y）乘以-4。
创建算法
如前所述，波动方程PDE可以有效地用有限差分方法求解。然而，我们仍然需要编写一些分辨率代码，才能设置一个真正的物理模拟。在最后一段中，我们最终得到了：
然后我们用以下方法简化了这个表达式：
这确实是一种递归形式，可以用以下伪代码进行建模：
对于 t 中的 时间
{
ut+1 <- ut+ut-1+something_at_t(x,y)
// 对新的 ut+1 进行一些操作，例如，您可以在此处绘制波
ut-1<- ut
ut <- ut+1
}
上面的伪代码应该以保留模式运行，因此应用程序可以简单地通过调用draw_wave(ut+1)来绘制我们正在考虑的表面的每个点上的波。
假设我们正在模拟水波的演变，因此让u函数代表波高（0 - 水平线）：我们现在可以写下伪代码的开头。
// 表面数据
高度 = 500;
宽度 = 500;
// 波参数
c = 1; // 波速
dx = 1; // 空间步长
dt = 0.05; // 时间步长
k=0.002 // 衰减因子
kdt=k*dt; // 每个时间步的衰减因子，请记住 q = 2 - kdt, r = -1 +kdt
c1 = dt^2*c^2/dx^2; // b 因子
// 初始条件
u_current_t = zeros(height,width); // 创建一个高度 x 宽度零矩阵
u_previous_t = u_current_t;
我们基本上定义了一个将绘制波演变效果的表面（一个500x500的区域），并初始化了我们在前几段中看到的波参数（请记住我们之前进行的q、r和b替换）。初始条件是零矩阵（u_current_t），因此整个表面是平静的，没有波在演变。
考虑到我们正在处理矩阵表面（位于(x;y)坐标的点由一个u值描述，表示波的高度），我们需要编写代码来实现
for循环中的那一行。实际上，上述表达式是以下表达式的简化形式：
我们需要实现最后一个。我们可以写如下：
for(t=0; t<A_LOT_OF_TIME; t += dt)
{
u_next_t = (2-kdt)*u_current_t+(kdt-1)*u_previous_t+c1*something_at_t(x,y)
u_previous_t = u_current_t; // 当前变为旧
u_current_t = u_next_t; // 新变为当前
// 绘制波
draw_wave(u_current_t)
}
也就是说，一个for循环，其索引变量t在每一步增加dt。现在一切都应该很熟悉了，因为（2-kdt）、（kdt-1）和c1项是通常的q、r和b替换。最后需要实现的是something_at_t(x,y)项，也称为：
我们开始的波动方程PDE是：
而我们现在感兴趣的项是这个：
在我们的例子中，它是：
由于我们有一个代表我们表面的点矩阵，我们完全处于离散域，并且由于我们需要对离散点矩阵执行二阶导数，因此我们的问题相当于有一个像素强度值为I(x,y)的图像，并且需要计算其拉普拉斯算子。
这是一个常见的图像处理任务，通常通过将卷积滤波器应用于我们感兴趣的图像（在本例中：代表表面的矩阵）来解决。用于近似拉普拉斯算子定义的二阶导数的一个常见小核是：
因此，为了实现
项，我们需要将D拉普拉斯核作为滤波器应用于我们的 图像（即u_current_t）：
u_next_t=(2-kdt)*u_current_t+(kdt-1)*u_previous_t+c1* convolution(u_current_t,D);
事实上，在前面我们看到的元素中 
红色点的元素被加权，并用D核的二维卷积计算。
在用D核与我们的u_current_t矩阵进行卷积时，一个需要记住的重要元素是，外部光晕中的每个值（在卷积中涉及但超出u_current_t矩阵边界的每个值）都应为零，如下图所示：
在上图中，红色边框是u_current_t矩阵的周长，而蓝色3x3矩阵是D拉普拉斯核，红色边框之外的所有值都为零。这一点很重要，因为我们希望我们的表面像容器中的水一样，水中的波撞击容器边缘时会“弹回”。通过将表面矩阵之外的所有值置零，我们不会从我们边界外部接收波的贡献，也不会以任何方式影响其外部。此外，波的“能量”不会扩散，并且会被方程部分地“弹回”。
现在算法几乎完成了：我们的PDE确保我们表面上的每个波峰都会被正确地“传输”成真实的波。唯一的问题是：从一个处处为零的矩阵开始并让它演变只会产生什么都没有。永远。
我们将向我们的表面添加一小滴水来扰动它并启动引擎。
为了尽可能逼真地模拟液滴落入表面的效果，我们将在表面矩阵中引入一个“包”。这意味着我们将添加一个代表离散高斯函数（类似于高斯核）的矩阵到我们的表面矩阵。请注意，我们将“添加”，而不是“卷积”。
鉴于二维高斯公式： 
我们有A是高斯幅度（因此将是我们的液滴幅度），x0 和y0 是中心的坐标，而
和
是斑点（blob）的展开宽度。将负幅度置于如下图像所示的位置，我们可以模拟一个刚刚落入表面的液滴。
生成液滴矩阵，我们可以使用以下伪代码：
// 生成液滴矩阵
dsz = 3; // 液滴大小
da=0.07; // 液滴幅度
[X,Y] = generate_XY_planes(dsz,da);
DropletMatrix = -da * exp( -(X/dsz)^2 -(Y/dsz)^2);
da是液滴幅度（并且如我们刚才所说，它是负的），而dsz是液滴大小，也就是说，高斯“钟形”的半径。X和Y是两个矩阵，表示X和Y平面的离散值。
因此，上面图像的X和Y矩阵是：
最终的DropletMatrix 类似于以下：
其中中心值是-0.0700。如果你绘制上述矩阵，你将得到一个三维高斯函数，现在它可以模拟我们的液滴。
我们波动算法的最终伪代码如下：
// 表面数据
高度 = 500;
宽度 = 500;
// 波参数
c = 1; // 波速
dx = 1; // 空间步长
dt = 0.05; // 时间步长
k=0.002 // 衰减因子
kdt=k*dt; // 每个时间步的衰减因子，请记住 q = 2 - kdt, r = -1 +kdt
c1 = dt^2*c^2/dx^2; // b 因子
// 初始条件
u_current_t=zeros(height,width); // 创建一个高度 x 宽度零矩阵
u_previous_t = u_current_t;
// 生成液滴矩阵
dsz = 3; // 液滴大小
da=0.07; // 液滴幅度
[X,Y] = generate_XY_planes(dsz,da);
DropletMatrix = -da * exp( -(X/dsz)^2 -(Y/dsz)^2);
// 此变量用于只添加一个液滴
One_single_droplet_added = 0;
for(t=0; t<A_LOT_OF_TIME; t = t + dt)
{
u_next_t = (2-kdt)*u_current_t+(kdt-1)*u_previous_t+c1*convolution(u_current_t,D);
u_previous_t = u_current_t; // 当前变为旧
u_current_t = u_next_t; // 新变为当前
// 绘制波
draw_wave(u_current_t)
if(One_single_droplet_added == 0)
{
One_single_droplet_added = 1; // 不再添加液滴
addMatrix(u_current_t, DropletMatrix, [+100;+100]);
}
}
变量One_single_droplet_added用于检查液滴是否已插入（我们只需要一个液滴）。addMatrix函数将DropletMatrix的值添加到u_current_t表面矩阵值中，将DropletMatrix置于点(100;100)的中心。请记住，DropletMatrix小于（或等于）表面矩阵，因此我们只需将DropletMatrix的值添加到落在DropletMatrix尺寸内的u_current_t的值中。
现在算法已经完成，尽管它仍然是一个理论模拟。我们将很快用实际代码实现它。
实现波模拟
我们现在将讨论上述算法如何在实际C++项目中实现，以创建完全功能的OpenGL物理模拟。
下面的序列图显示了程序的骨架，该程序基本上由三个主要部分组成：包含启动函数和创建整个窗口的矩阵图像的内核函数的主模块；一个封装GLUT库函数和回调处理程序的OpenGL渲染器包装器；以及一个手动编写的矩阵类，用于简化矩阵数据访问和操作。尽管序列图需要遵循标准的软件工程方法，并且其使用受到预定规则的严格执行，但我们仍将使用它作为抽象来建模程序的控制流。
程序在main()开始，并创建一个openglrenderer对象，该对象将处理所有图形通信与GLUT库和回调事件（如鼠标移动、键盘按键事件等）。OpenGL不提供用于与窗口系统交互的例程，因此在此项目中我们将依赖GLUT库，该库提供了一个平台无关的接口来管理窗口和输入事件。为了尽可能快地创建动画，我们将使用glutIdleFunc()函数设置一个空闲回调。我们将在稍后详细解释这一点。
最初，算法设置其初始化变量（时间步长、空间步长、液滴幅度、拉普拉斯2D核等……基本上是我们在理论部分看到的所有内容），并且与要渲染的图像对应的每个矩阵都被归零。高斯液滴矩阵也被预处理。OpenGL渲染器头文件中定义的一个结构包含了应在图像渲染之间保留的所有数据。
typedef struct kernelData
{
float a1; // 2-kdt
float a2; // kdt-1
float c1; // decayment factor
sMatrix* u;
sMatrix* u0;
sMatrix* D;
int Ddim; // Droplet matrix width/height
int dsz; // Gaussian radius
sMatrix* Zd;
} kernelData;
该结构在执行每个时间步时都会更新，因为它包含描述波随时间演变的当前和先前矩阵。由于此结构同时被OpenGL渲染器和主模块用于初始化，因此该变量被声明为external，并在OpenGL渲染器cpp文件中定义（因此其范围超出了单个翻译单元）。在所有设置完成后，将调用openglRenderer类的方法startRendering()，并且OpenGL主循环开始获取事件。我们在主模块的kernel()函数中看到了算法的核心，该函数在每次调度OpenGL空闲事件时被调用（也就是说，屏幕被更新，并且更改只有在空闲回调完成后才会显示，因此此处执行的渲染工作量应最小化以避免性能损失）。
kernel函数的代码如下：
// This kernel is called at each iteration
// It implements the main loop algorithm and someway "rasterize" the matrix data
// to pass to the openGL renderer. It also adds droplets in the waiting queue
//
void kernel(unsigned char *ptr, kernelData& data)
{
// Wave evolution
sMatrix un(DIM,DIM);
// The iterative discrete update (see documentation)
un = (*data.u)*data.a1 + (*data.u0)*data.a2 + convolutionCPU((*data.u),(*data.D))*data.c1;
// Step forward in time
(*data.u0) = (*data.u);
(*data.u) = un;
// Prepare matrix data for rendering
matrix2Bitmap( (*data.u), ptr );
if(first_droplet == 1) // By default there's just one initial droplet
{
first_droplet = 0;
int x0d= DIM / 2; // Default droplet center
int y0d= DIM / 2;
// Place the (x0d;y0d) centered Zd droplet on the wave data (it will be added at the next iteration)
for(int Zdi=0; Zdi < data.Ddim; Zdi++)
{
for(int Zdj=0; Zdj < data.Ddim; Zdj++)
{
(*data.u)(y0d-2*data.dsz+Zdi,x0d-2*data.dsz+Zdj) += (*data.Zd)(Zdi, Zdj);
}
}
}
// Render the droplet queue
m_renderer->renderWaitingDroplets();
}
我们在波动PDE演化中遵循的模式很容易在计算密集型代码行中识别：
un = (*data.u)*data.a1 + (*data.u0)*data.a2 + convolutionCPU((*data.u),(*data.D))*data.c1;
这基本上执行了著名的迭代步：
为提高性能，所有常数都已预处理。
需要注意的是，这一行将大的矩阵相加，在代码中称为sMatrix对象。sMatrix类是一个手动编写的简单类，它公开简单的运算符重载，以方便使用矩阵。除了应该记住大型矩阵操作应避免按值传递参数，并且要创建一个新矩阵并将其复制到目标需要复制构造函数（以避免获得浅拷贝而没有实际数据）之外，代码相当直接，因此不再赘述。
// This class handles matrix objects
class sMatrix
{
public:
int rows, columns;
float *values;
sMatrix(int height, int width)
{
if (height == 0 || width == 0)
throw "Matrix constructor has 0 size";
rows = height;
columns = width;
values = new float[rows*columns];
}
// Copy constructor, this is needed to perform a deep copy (not a shallow one)
sMatrix(const sMatrix& mt)
{
this->rows = mt.rows;
this->columns = mt.columns;
this->values = new float[rows*columns];
// Copy the values
memcpy(this->values, mt.values, this->rows*this->columns*sizeof(float));
}
~sMatrix()
{
delete [] values;
}
// Allows matrix1 = matrix2
sMatrix& operator= (sMatrix const& m)
{
if(m.rows != this->rows || m.columns != this->columns)
{
throw "Size mismatch";
}
memcpy(this->values,m.values,this->rows*this->columns*sizeof(float));
return *this; // Since "this" continues to exist after the function call, it is perfectly legit to return a reference
}
// Allows both matrix(3,3) = value and value = matrix(3,3)
float& operator() (const int row, const int column)
{
// May be suppressed to slightly increase performances
if (row < 0 || column < 0 || row > this->rows || column > this->columns)
throw "Size mismatch";
return values[row*columns+column]; // Since the float value continues to exist after the function call, it is perfectly legit to return a reference
}
// Allows scalar*matrix (e.g. 3*matrix) for each element
sMatrix operator* (const float scalar)
{
sMatrix result(this->rows, this->columns);
// Multiply each value by the scalar
for(int i=0; i<rows*columns; i++)
{
result.values[i] = this->values[i]*scalar;
}
return result; // Copy constructor
}
// Allows matrix+matrix (if same size)
sMatrix operator+ (const sMatrix& mt)
{
if (this->rows != mt.rows || this->columns != mt.columns)
throw "Size mismatch";
sMatrix result(this->rows, this->columns);
// Sum each couple of values
for(int i=0; i<rows; i++)
{
for(int j=0; j<columns; j++)
result.values[i*columns+j] = this->values[i*columns+j] + mt.values[i*columns+j];
}
return result; // Copy constructor
}
};
卷积是用以下代码执行的（一种经典方法）：
// Returns the convolution between the matrix A and the kernel matrix B,
// A's size is preserved
//
sMatrix convolutionCPU(sMatrix& A, sMatrix& B)
{
sMatrix result(A.rows, A.columns);
int kernelHradius = (B.rows - 1) / 2;
int kernelWradius = (B.columns - 1) / 2;
float convSum, convProd, Avalue;
for(int i=0; i<A.rows; i++)
{
for(int j=0; j<A.columns; j++)
{
//
// --------j--------->
// _ _ _ _ _ _ _
// | | |
// | | |
// | A | i
// | | |
// | | |
// _ _ _ _ _ _ _| v
//
convSum = 0;
for(int bi=0; bi<B.rows; bi++)
{
for(int bj=0; bj<B.columns; bj++)
{
// A's value respect to the kernel center
int relpointI = i-kernelHradius+bi;
int relpointJ = j-kernelWradius+bj;
if(relpointI < 0 || relpointJ < 0 || relpointI >= A.rows || relpointJ >= A.columns)
Avalue = 0;
else
Avalue = A(i-kernelHradius+bi,j-kernelWradius+bj);
convProd = Avalue*B(bi,bj);
convSum += convProd;
}
}
// Store the convolution result
result(i,j) = convSum;
}
}
return result;
} 在计算完系统演化后，时间经过是通过将新矩阵与旧矩阵交换并丢弃先前状态来模拟的，如我们之前所述。
然后调用matrix2Bitmap()执行矩阵到位图的转换，正如其名称所示，更准确地说，整个模拟区域由一个大的sMatrix对象描述，该对象显然包含浮点值。为了将这些值实际渲染为像素单位，我们需要将每个值转换为其相应的RGBA值，并将其传递给openglRenderer类（该类反过来将整个位图缓冲区传递给GLUT库）。简而言之：我们需要执行浮点到RGB颜色的映射。
由于在物理模拟中我们假设静止水面高度为0，任何扰动都会升高或降低该值（特别是液滴高斯矩阵将其降低最多-0.07），我们正在寻找一个[-1;1]到颜色的映射。HSV颜色模型可以更好地模拟我们用肉眼实际体验到的渐变色过渡，但这需要将其转换回RGB值，以便设置一个位图图以传回GLUT包装器。出于性能原因，我们选择为每个值分配一种颜色（颜色映射）。一个初步的解决方案是实现一个完整的[-1;1] -> [0x0;0xFFFFFF]映射，以涵盖RGB格式的所有可能颜色。
// Returns a RGB color from a value in the interval between [-1;+1]
RGB getColorFromValue(float value)
{
RGB result;
if(value <= -1.0f)
{
result.r = 0x00;
result.g = 0x00;
result.b = 0x00;
}
else if(value >= 1.0f)
{
result.r = 0xFF;
result.g = 0xFF;
result.b = 0xFF;
}
else
{
float step = 2.0f/0xFFFFFF;
unsigned int cvalue = (unsigned int)((value + 1.0f)/step);
if(cvalue < 0)
cvalue = 0;
else if(cvalue > 0xFFFFFF)
cvalue = 0xFFFFFF;
result.r = cvalue & 0xFF;
result.g = (cvalue & 0xFF00) >> 8;
result.b = (cvalue & 0xFF0000) >> 16;
}
return result;
}
然而，上述方法既耗费性能，又不能很好地渲染平滑的颜色过渡：让我们看看以这种方式映射的液滴。
看起来更像分形而不是液滴，所以上述解决方案行不通。一种改进性能（和图像外观）的更好方法是将颜色映射硬编码在一个数组中，并在需要时使用它。
float pp_step = 2.0f / (float)COLOR_NUM;
// The custom colormap, you can customize it to whatever you like
unsigned char m_colorMap[COLOR_NUM][3] =
{
{0,0,143},{0,0,159},{0,0,175},{0,0,191},{0,0,207},{0,0,223},{0,0,239},{0,0,255},
{0,16,255},{0,32,255},{0,48,255},{0,64,255},{0,80,255},{0,96,255},{0,111,255},{0,128,255},
{0,143,255},{0,159,255},{0,175,255},{0,191,255},{0,207,255},{0,223,255},{0,239,255},{0,255,255},
{16,255,239},{32,255,223},{48,255,207},{64,255,191},{80,255,175},{96,255,159},{111,255,143},{128,255,128},
{143,255,111},{159,255,96},{175,255,80},{191,255,64},{207,255,48},{223,255,32},{239,255,16},{255,255,0},
{255,239,0},{255,223,0},{255,207,0},{255,191,0},{255,175,0},{255,159,0},{255,143,0},{255,128,0},
{255,111,0},{255,96,0},{255,80,0},{255,64,0},{255,48,0},{255,32,0},{255,16,0},{255,0,0},
{239,0,0},{223,0,0},{207,0,0},{191,0,0},{175,0,0},{159,0,0},{143,0,0},{128,0,0},
{100,0,20},{80,0,40},{60,0,60},{40,0,80},{20,0,100},{0,0,120}
};
// Returns a RGB color from a value in the interval between [-1;+1] using the above colormap
RGB getColorFromValue(float value)
{
RGB result;
unsigned int cvalue = (unsigned int)((value + 1.0f)/pp_step);
if(cvalue < 0)
cvalue = 0;
else if(cvalue >= COLOR_NUM)
cvalue = COLOR_NUM-1;
result.r = m_colorMap[cvalue][0];
result.g = m_colorMap[cvalue][1];
result.b = m_colorMap[cvalue][2];
return result;
}
创建颜色映射并不难，并且网络上有各种不同的颜色映射可供选择，它们会产生不同的过渡效果。这次结果要好得多（稍后见截图），性能（尽管每次值转换总是一项繁重的任务）也大大提高。
涉及屏幕渲染的最后一部分是在用户用光标点击窗口的任何地方添加液滴。一个液滴会自动添加到表面的中心（您可以在kernel()函数中找到代码，它由first_droplet变量控制），但用户可以点击表面上的任何地方（几乎任何地方）添加另一个液滴。为了实现这一点，实现了一个队列，最多可以包含60个液滴中心，高斯矩阵将放置在这些中心（请注意，矩阵将被添加到表面值中，而不是简单地替换它们）。
#define MAX_DROPLET 60
typedef struct Droplet
{
int mx;
int my;
} Droplet;
Droplet dropletQueue[MAX_DROPLET];
int dropletQueueCount = 0; 队列系统之所以实现，是有原因的：与我们之前编写的伪代码算法不同，使用OpenGL渲染场景需要程序控制立即模式下要显示的物体：这意味着程序需要在实际渲染执行之前处理要绘制的内容，它不能简单地将液滴放入要绘制的数据中，因为它可能正在使用中（您可以在保留模式下执行此操作）。此外，我们不知道液滴何时会被添加，因为它完全取决于用户。因此，每次kernel()完成时，液滴队列都会被清空，液滴高斯矩阵会被添加到要渲染的数据（表面数据）中。执行此操作的代码如下：
void openGLRenderer::renderWaitingDroplets()
{
// If there are droplets waiting to be rendered, schedule for the next rendering
while(dropletQueueCount > 0)
{
dropletQueueCount--;
addDroplet(dropletQueue[dropletQueueCount].mx,dropletQueue[dropletQueueCount].my);
}
}
void addDroplet( int x0d, int y0d )
{
y0d = DIM - y0d;
// Place the (x0d;y0d) centered Zd droplet on the wave data (it will be added at the next iteration)
for(int Zdi=0; Zdi< m_simulationData.Ddim; Zdi++)
{
for(int Zdj=0; Zdj< m_simulationData.Ddim; Zdj++)
{
(*m_simulationData.u)(y0d-2*m_simulationData.dsz+Zdi,x0d-2*m_simulationData.dsz+Zdj)
+= (*m_simulationData.Zd)(Zdi, Zdj);
}
}
}
代码现在应该很熟悉了：addDroplet函数只是将Zd二维高斯矩阵添加到模拟数据（内核数据）中，在“当前”时间，即代表表面的u矩阵。
代码循环直到键盘回调处理程序（由openglrenderer定义）检测到Esc键按下，之后应用程序收到终止信号，循环结束。程序释放分配的资源，然后发出退出信号，但由于程序终止会释放其先前分配的所有资源，因此这可能不是必需的。
添加液滴功能和所有处理程序设置后，代码就可以运行了。这次结果比以前好得多，因为我们使用了更平滑的颜色映射（请看下面的图像）。请注意，卷积项如何产生“波传播”和“反弹”效果，当计算与表面矩阵之外的填充零数据的值时（即当波撞击窗口边缘并被反射回来时）。第一张图像是模拟的早期阶段，即一些液滴刚刚添加时，第二张图像代表后期阶段，此时表面即将平静下来（在我们的颜色映射中，蓝色值高于红色值）。
由于我们引入了阻尼因子（请回顾理论部分），波最终会停止，表面迟早会再次平静下来。整个模拟（除了颜色）相当逼真，但也很慢。这是因为整个表面矩阵正在被系统彻底更新和计算。kernel()函数连续运行，更新渲染缓冲区。对于512x512的图像，CPU必须处理大量数据，并且还必须执行512x512的浮点卷积。使用性能分析器（如VS的集成性能分析器）显示，程序大部分时间花费在kernel()调用中（正如我们所预期的），并且卷积操作是最耗费CPU资源的操作。
同样有趣的是，添加大量液滴时，模拟速度会显著下降。
在真实的科学模拟环境中，需要在大约短时间间隔内处理海量数据。这时GPGPU计算就派上用场了。我们将简要总结这个缩写词的含义，然后介绍一个GPU加速的波动方程PDE模拟版本。
GPGPU架构
GPGPU代表通用图形处理器计算，它指的是使用图形处理器和设备执行高并行化的计算，而这些计算通常由CPU设备处理。利用图形设备帮助CPU处理工作负载的想法并不新鲜，但直到最近的架构和框架（如CUDA（©NVIDIA供应商特定）或OpenCL）出现之前，程序员不得不依赖一系列的变通方法和技巧，以不方便且不直观的方法和数据结构进行工作。开发人员应该考虑将CPU原生代码移植到新GPU版本的原因在于CPU和GPU之间的架构设计差异。虽然CPU（多核）通过顺序执行（流水线、缓存、控制流等）来获得性能优势，但GPU则以多核方式发展：它们倾向于在更高的数据带宽下运行，并选择大幅增加执行线程数量。近年来，GPGPU已被视为一种利用图形处理单元作为代数协处理器的新机会，能够处理大规模并行化和精确浮点计算。GPGPU架构背后的思想是让CPU处理程序中的顺序部分，而GPU处理可并行化的部分。事实上，许多科学应用和系统通过这种方法提高了性能，GPGPU现在是许多领域（如医学成像、物理模拟、信号处理、密码学、入侵检测、环境科学等）的基本技术。
使用CUDA实现大规模并行化
我们选择使用CUD-Architecture来并行化我们代码的某些部分。使用GPGPU技术并行化意味着从顺序设计的代码转变为并行设计的代码，这通常也意味着不得不重写代码的大部分。整个模拟算法中最明显可以从并行方法中受益的部分是表面矩阵与二维拉普拉斯核的卷积。
注意：为简洁起见，我们将假设读者已经熟悉CUDA C/C++。
CUDA SDK附带了大量的示例，说明如何实现高效的矩阵卷积以及如何将核作为图像滤波器应用。然而，我们决定实现自己的版本，而不是依赖标准技术。原因有很多：
- 我们想展示如何设计和实现大规模并行算法。
- 由于我们的核非常小（3x3二维拉普拉斯核，基本上是9个浮点值），使用Victor Podlozhnyuk描述的FFT方法效率不高。
- 二维高斯核是唯一径向对称且可分离的函数，我们的核不是，因此我们无法使用可分离卷积方法。
- 这样一个小的核似乎非常适合在卷积运算中“积极缓存”。一旦描述了设计的CUDA核，我们将对此进行详细说明。
执行卷积的最明显方法（尽管效率极低）是委托每个GPU线程对整个图像的每个元素执行多个卷积。
请看下图，我们将使用一个9x9的线程网格（为简化起见，只有一个块）来执行卷积。紫色方块是我们的9x9网格，而红色网格对应于9x9的核。每个线程在其元素内执行卷积，然后X轴移动，整个块“虚拟地”向右移动。当X坐标完成（即覆盖了整个图像的水平区域）时，Y坐标增加，然后重复此过程直到完成。在边界区域，图像外的每个值都将设置为零。
这种简单方法的代码如下，其中A是图像矩阵，B是核。
__global__ void convolutionGPU(float *A, int A_cols, int A_rows, float *B, int B_Wradius, int B_Hradius,
int B_cols, int B_rows, float *result)
{
// Initial position
int threadXabs = blockIdx.x*blockDim.x + threadIdx.x;
int threadYabs = blockIdx.y*blockDim.y + threadIdx.y;
int threadXabsInitialPos = threadXabs;
float convSum;
while(threadYabs < A_rows)
{
while(threadXabs < A_cols)
{
// If we are still in the image, start the convolution
convSum = 0.0f;
// relative x coord to the absolute thread
#pragma unroll
for(int xrel=-B_Wradius;xrel<(B_Wradius+1); xrel++)
{
#pragma unroll
for(int yrel=-B_Hradius;yrel<(B_Hradius+1); yrel++)
{
// Check the borders, 0 if outside
float Avalue;
if(threadXabs + xrel < 0 || threadYabs + yrel <0 || threadXabs + xrel >= A_cols || threadYabs + yrel >= A_rows)
Avalue = 0;
else
Avalue = A[ (threadYabs+yrel)*A_cols + (threadXabs + xrel) ];
// yrel+b_hradius makes the interval positive
float Bvalue = B[ (yrel+B_Hradius)*B_cols + (xrel+B_Wradius) ];
convSum += Avalue * Bvalue;
}
}
// Store the result and proceed ahead in the grid
result[threadYabs*A_cols + threadXabs ] = convSum;
threadXabs += blockDim.x * gridDim.x;
}
// reset X pos and forward the Y pos
threadXabs = threadXabsInitialPos;
threadYabs += blockDim.y * gridDim.y;
}
// convolution finished
} 如前所述，这种简单方法有几个缺点：
- 核非常小，将其保留在全局内存中并为每次卷积访问它效率极低。
- 虽然矩阵读取部分合并（coalesced），但在边界区域活动的线程和不活动的线程之间可能会出现显著的线程发散。
- 线程之间没有协作行为，尽管它们基本上使用相同的核并共享大部分围边区域。
因此，设计了一种更好的GPU卷积方法，同时考虑到上述几点。
思路很简单：让每个线程将围边区域和数据区域的一部分加载到共享内存中，从而最大化读取合并并减少发散。
在GPU版本模拟中执行卷积的代码如下：
// For a 512x512 image the grid is 170x170 blocks 3x3 threads each one
__global__ void convolutionGPU(float *A, float *result)
{
__shared__ float data[laplacianD*2][laplacianD*2];
// Absolute position into the image
const int gLoc = threadIdx.x + IMUL(blockIdx.x,blockDim.x) + IMUL(threadIdx.y,DIM) + IMUL(blockIdx.y,blockDim.y)*DIM;
// Image-relative position
const int x0 = threadIdx.x + IMUL(blockIdx.x,blockDim.x);
const int y0 = threadIdx.y + IMUL(blockIdx.y,blockDim.y);
// Load the apron and data regions
int x,y;
// Upper left square
x = x0 - kernelRadius;
y = y0 - kernelRadius;
if(x < 0 || y < 0)
data[threadIdx.x][threadIdx.y] = 0.0f;
else
data[threadIdx.x][threadIdx.y] = A[ gLoc - kernelRadius - IMUL(DIM,kernelRadius)];
// Upper right square
x = x0 + kernelRadius + 1;
y = y0 - kernelRadius;
if(x >= DIM || y < 0)
data[threadIdx.x + blockDim.x][threadIdx.y] = 0.0f;
else
data[threadIdx.x + blockDim.x][threadIdx.y] = A[ gLoc + kernelRadius+1 - IMUL(DIM,kernelRadius)];
// Lower left square
x = x0 - kernelRadius;
y = y0 + kernelRadius+1;
if(x < 0 || y >= DIM)
data[threadIdx.x][threadIdx.y + blockDim.y] = 0.0f;
else
data[threadIdx.x][threadIdx.y + blockDim.y] = A[ gLoc - kernelRadius + IMUL(DIM,(kernelRadius+1))];
// Lower right square
x = x0 + kernelRadius+1;
y = y0 + kernelRadius+1;
if(x >= DIM || y >= DIM)
data[threadIdx.x + blockDim.x][threadIdx.y + blockDim.y] = 0.0f;
else
data[threadIdx.x + blockDim.x][threadIdx.y + blockDim.y] = A[ gLoc + kernelRadius+1 + IMUL(DIM,(kernelRadius+1))];
__syncthreads();
float sum = 0;
x = kernelRadius + threadIdx.x;
y = kernelRadius + threadIdx.y;
// Execute the convolution in the shared memory (kernel is in constant memory)
#pragma unroll
for(int i = -kernelRadius; i<=kernelRadius; i++)
for(int j=-kernelRadius; j<=kernelRadius; j++)
sum += data[x+i][y+j] * gpu_D[i+kernelRadius][j+kernelRadius];
// Transfer the risult to global memory
result[gLoc] = sum;
}
核只接收表面矩阵和存储卷积图像的结果。没有提供核，因为它被放入了一个称为“constant内存”的特殊内存中，该内存是核只读的，预取并且高度优化，让所有线程以最小的延迟从特定位置读取。缺点是这种内存非常有限（在64Kb的量级），所以应该明智地使用。将我们的3x3核声明为常量内存可以给我们带来显著的速度优势。
__device__ __constant__ float gpu_D[laplacianD][laplacianD]; // 拉普拉斯2D核
下图有助于确定线程如何从全局内存中的表面矩阵加载数据，并将它们存储到更快的片上共享内存中，然后再实际使用它们进行卷积运算。紫色3x3方块是核窗口，中心元素是我们实际旋转的值。网格是172x172个块，每个块有3x3个线程；每个3x3线程块在进入卷积循环之前需要完成四个阶段：将上半部分围边和图像数据加载到共享内存（从核紫色窗口的上半部分红色方块），加载上半部分区域（红色方块），加载下半部分区域（红色方块），以及加载下半部分区域（同上）。由于共享内存只对块中的线程可用，每个块加载自己的共享区域。请注意，我们选择让每个线程从全局内存读取一些内容以最大化合并，但我们实际上并不使用所有元素。图像显示了一个黄色区域和一个灰色区域：黄色数据将用于我们正在考虑的块执行的每个元素紫色核方块的卷积运算（它包括围边和数据），而灰色区域不会被该块执行的任何卷积使用。
在填充每个块的共享内存数组后，CUDA线程会进行同步以最小化它们的分歧。然后执行卷积算法：共享数据与恒定核数据相乘，从而产生高度优化的操作。
`#pragma unroll`指令指示编译器（在可能的情况下）展开循环，以减少周期控制开销并提高性能。循环展开的一个小例子：以下循环
for(int i=0;i<1000;i++)
a[i] = b[i] + c[i];
可以通过展开进行优化
for(int i=0;i<1000;i+=2)
{
a[i] = b[i] + c[i];
a[i+1] = b[i+1] + c[i+1];
}
以便执行的控制指令更少，并且整个循环提高了性能。需要注意的是，CUDA代码中的几乎所有优化都需要经过仔细而彻底的测试，因为不同的架构和不同的程序控制流可能会产生不同的结果（以及不同的编译器优化，不幸的是，这些优化并不总是可信的）。
另请注意，代码中使用了IMUL宏，它定义为
#define IMUL(a,b) __mul24(a,b)
在CUDA计算能力为1.x的设备上，32位整数乘法使用多个指令实现，因为它不是本机支持的。24位整数乘法通过__ [u] mul24内在函数本机支持。然而，在计算能力为2.0的设备上，32位整数乘法是本机支持的，但24位整数乘法不是。因此，__ [u] mul24使用多个指令实现，不应使用。所以，如果你打算在2.x设备上使用该代码，请确保重新定义宏指令。
调用我们刚刚编写的核的典型代码可能是：
sMatrix convolutionGPU_i(sMatrix& A, sMatrix& B)
{
unsigned int A_bytes = A.rows*A.columns*sizeof(float);
sMatrix result(A.rows, A.columns);
float *cpu_result = (float*)malloc(A_bytes);
// Copy A data to the GPU global memory (B aka the kernel is already there)
cudaError_t chk;
chk = cudaMemcpy(m_gpuData->gpu_matrixA, A.values, A_bytes, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT TRANSFER MEMORY TO GPU");
return result;
}
// Call the convolution kernel
dim3 blocks(172,172);
dim3 threads(3,3);
convolutionGPU<<<blocks,threads>>>(m_gpuData->gpu_matrixA, m_gpuData->gpu_matrixResult);
// Copy back the result
chk = cudaMemcpy(cpu_result, m_gpuData->gpu_matrixResult, A_bytes, cudaMemcpyDeviceToHost);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT TRANSFER MEMORY FROM GPU");
return result;
}
// Allocate a sMatrix and return it with the GPU data
free(result.values);
result.values = cpu_result;
return result;
}
显然，CUDA内存应该在程序开始时使用cudaMalloc分配，并在GPU工作完成后才释放。
然而，正如我们之前所说，将顺序设计的程序转换为并行程序并非易事，通常需要比简单的函数到函数转换更多的东西（这取决于应用程序）。在我们的例子中，仅仅用GPU卷积函数替换CPU卷积函数是行不通的。事实上，即使我们从CPU版本中更好地分配了工作负载（请参阅下面的CPU-GPU独占时间百分比图），**我们实际上减慢了整个模拟速度**。
原因很简单：我们的kernel()函数在每次调度绘图事件时都会被调用，所以需要非常频繁地调用它。尽管CUDA核比CPU卷积函数快，并且尽管GPU内存带宽高于CPU，但在（可能已分页）主机内存和全局设备内存之间来回传输会扼杀我们的模拟性能。从CUDA方法中受益更多的应用程序通常执行一次计算量大的计算工作负载，然后将结果传回。实时应用程序可能从并发核方法中受益，但这需要2.x能力的设备。
为了真正加速我们的模拟，需要对代码进行更大的修订。
另一个需要考虑的更微妙的事情是GPU代码的CPU优化：看一下下面CPU版本的代码行
un = (*data.u)*data.a1 + (*data.u0)*data.a2 + convolutionCPU((*data.u),(*data.D))*data.c1;
000000013F2321EF mov r8,qword ptr [m_simulationData+20h (13F234920h)]
000000013F2321F6 mov rdx,qword ptr [m_simulationData+10h (13F234910h)]
000000013F2321FD lea rcx,[rbp+2Fh]
000000013F232201 call convolutionCPU (13F231EC0h)
000000013F232206 nop
000000013F232207 movss xmm2,dword ptr [m_simulationData+8 (13F234908h)]
000000013F23220F lea rdx,[rbp+1Fh]
000000013F232213 mov rcx,rax
000000013F232216 call sMatrix::operator* (13F2314E0h)
000000013F23221B mov rdi,rax
000000013F23221E movss xmm2,dword ptr [m_simulationData+4 (13F234904h)]
000000013F232226 lea rdx,[rbp+0Fh]
000000013F23222A mov rcx,qword ptr [m_simulationData+18h (13F234918h)]
000000013F232231 call sMatrix::operator* (13F2314E0h)
000000013F232236 mov rbx,rax
000000013F232239 movss xmm2,dword ptr [m_simulationData (13F234900h)]
000000013F232241 lea rdx,[rbp-1]
000000013F232245 mov rcx,qword ptr [m_simulationData+10h (13F234910h)]
000000013F23224C call sMatrix::operator* (13F2314E0h)
000000013F232251 nop
000000013F232252 mov r8,rbx
000000013F232255 lea rdx,[rbp-11h]
000000013F232259 mov rcx,rax
000000013F23225C call sMatrix::operator+ (13F2315B0h)
000000013F232261 nop
000000013F232262 mov r8,rdi
000000013F232265 lea rdx,[rbp-21h]
000000013F232269 mov rcx,rax
000000013F23226C call sMatrix::operator+ (13F2315B0h)
000000013F232271 nop
000000013F232272 cmp dword ptr [rax],1F4h
000000013F232278 jne kernel+33Fh (13F2324CFh)
000000013F23227E cmp dword ptr [rax+4],1F4h
000000013F232285 jne kernel+33Fh (13F2324CFh)
000000013F23228B mov r8d,0F4240h
000000013F232291 mov rdx,qword ptr [rax+8]
000000013F232295 mov rcx,r12
000000013F232298 call memcpy (13F232DDEh)
000000013F23229D nop
000000013F23229E mov rcx,qword ptr [rbp-19h]
000000013F2322A2 call qword ptr [__imp_operator delete (13F233090h)]
000000013F2322A8 nop
000000013F2322A9 mov rcx,qword ptr [rbp-9]
000000013F2322AD call qword ptr [__imp_operator delete (13F233090h)]
000000013F2322B3 nop
000000013F2322B4 mov rcx,qword ptr [rbp+7]
000000013F2322B8 call qword ptr [__imp_operator delete (13F233090h)]
000000013F2322BE nop
000000013F2322BF mov rcx,qword ptr [rbp+17h]
000000013F2322C3 call qword ptr [__imp_operator delete (13F233090h)]
000000013F2322C9 nop
000000013F2322CA mov rcx,qword ptr [rbp+27h]
000000013F2322CE call qword ptr [__imp_operator delete (13F233090h)]
000000013F2322D4 nop
000000013F2322D5 mov rcx,qword ptr [rbp+37h]
000000013F2322D9 call qword ptr [__imp_operator delete (13F233090h)]
现在看看GPU版本的代码行
un = (*data.u)*data.a1 + (*data.u0)*data.a2 + convolutionGPU_i ((*data.u),(*data.D))*data.c1;
000000013F7E23A3 mov rax,qword ptr [data]
000000013F7E23AB movss xmm0,dword ptr [rax+8]
000000013F7E23B0 movss dword ptr [rsp+0A8h],xmm0
000000013F7E23B9 mov rax,qword ptr [data]
000000013F7E23C1 mov r8,qword ptr [rax+20h]
000000013F7E23C5 mov rax,qword ptr [data]
000000013F7E23CD mov rdx,qword ptr [rax+10h]
000000013F7E23D1 lea rcx,[rsp+70h]
000000013F7E23D6 call convolutionGPU_i (13F7E1F20h)
000000013F7E23DB mov qword ptr [rsp+0B0h],rax
000000013F7E23E3 mov rax,qword ptr [rsp+0B0h]
000000013F7E23EB mov qword ptr [rsp+0B8h],rax
000000013F7E23F3 movss xmm0,dword ptr [rsp+0A8h]
000000013F7E23FC movaps xmm2,xmm0
000000013F7E23FF lea rdx,[rsp+80h]
000000013F7E2407 mov rcx,qword ptr [rsp+0B8h]
000000013F7E240F call sMatrix::operator* (13F7E2B20h)
000000013F7E2414 mov qword ptr [rsp+0C0h],rax
000000013F7E241C mov rax,qword ptr [rsp+0C0h]
000000013F7E2424 mov qword ptr [rsp+0C8h],rax
000000013F7E242C mov rax,qword ptr [data]
000000013F7E2434 movss xmm0,dword ptr [rax+4]
000000013F7E2439 movaps xmm2,xmm0
000000013F7E243C lea rdx,[rsp+50h]
000000013F7E2441 mov rax,qword ptr [data]
000000013F7E2449 mov rcx,qword ptr [rax+18h]
000000013F7E244D call sMatrix::operator* (13F7E2B20h)
000000013F7E2452 mov qword ptr [rsp+0D0h],rax
000000013F7E245A mov rax,qword ptr [rsp+0D0h]
000000013F7E2462 mov qword ptr [rsp+0D8h],rax
000000013F7E246A mov rax,qword ptr [data]
000000013F7E2472 movss xmm2,dword ptr [rax]
000000013F7E2476 lea rdx,[rsp+40h]
000000013F7E247B mov rax,qword ptr [data]
000000013F7E2483 mov rcx,qword ptr [rax+10h]
000000013F7E2487 call sMatrix::operator* (13F7E2B20h)
000000013F7E248C mov qword ptr [rsp+0E0h],rax
000000013F7E2494 mov rax,qword ptr [rsp+0E0h]
000000013F7E249C mov qword ptr [rsp+0E8h],rax
000000013F7E24A4 mov r8,qword ptr [rsp+0D8h]
000000013F7E24AC lea rdx,[rsp+60h]
000000013F7E24B1 mov rcx,qword ptr [rsp+0E8h]
000000013F7E24B9 call sMatrix::operator+ (13F7E2BF0h)
000000013F7E24BE mov qword ptr [rsp+0F0h],rax
000000013F7E24C6 mov rax,qword ptr [rsp+0F0h]
000000013F7E24CE mov qword ptr [rsp+0F8h],rax
000000013F7E24D6 mov r8,qword ptr [rsp+0C8h]
000000013F7E24DE lea rdx,[rsp+90h]
000000013F7E24E6 mov rcx,qword ptr [rsp+0F8h]
000000013F7E24EE call sMatrix::operator+ (13F7E2BF0h)
000000013F7E24F3 mov qword ptr [rsp+100h],rax
000000013F7E24FB mov rax,qword ptr [rsp+100h]
000000013F7E2503 mov qword ptr [rsp+108h],rax
000000013F7E250B mov rdx,qword ptr [rsp+108h]
000000013F7E2513 lea rcx,[un]
000000013F7E2518 call sMatrix::operator= (13F7E2A90h)
000000013F7E251D nop
000000013F7E251E lea rcx,[rsp+90h]
000000013F7E2526 call sMatrix::~sMatrix (13F7E2970h)
000000013F7E252B nop
000000013F7E252C lea rcx,[rsp+60h]
000000013F7E2531 call sMatrix::~sMatrix (13F7E2970h)
000000013F7E2536 nop
000000013F7E2537 lea rcx,[rsp+40h]
000000013F7E253C call sMatrix::~sMatrix (13F7E2970h)
000000013F7E2541 nop
000000013F7E2542 lea rcx,[rsp+50h]
000000013F7E2547 call sMatrix::~sMatrix (13F7E2970h)
000000013F7E254C nop
000000013F7E254D lea rcx,[rsp+80h]
000000013F7E2555 call sMatrix::~sMatrix (13F7E2970h)
000000013F7E255A nop
000000013F7E255B lea rcx,[rsp+70h]
000000013F7E2560 call sMatrix::~sMatrix(13F7E2970h)
代码虽然涉及的数据几乎相同，但看起来要臃肿得多（甚至有一些无意义的操作，请看地址000000013F7E23DB）。可能让CPU在GPU完成工作后继续计算不是个好主意。
由于还有其他函数可以并行化（如matrix2Bitmap()函数），我们需要将尽可能多的工作负载转移到设备上。
首先，我们需要在程序开始时在设备上分配内存，并在完成后取消分配。像我们算法参数这样的小数据可以存储在常量内存中以提高性能，而矩阵大数据的存储更适合全局内存（常量内存大小非常有限）。
// Initialize all data used by the device
// and the rendering simulation data as well
void initializeGPUData()
{
/* Algorithm parameters */
// Time step
float dt = (float)0.05;
// Speed of the wave
float c = 1;
// Space step
float dx = 1;
// Decay factor
float k = (float)0.002;
// Droplet amplitude (Gaussian amplitude)
float da = (float)0.07;
// Initialize u0
sMatrix u0(DIM,DIM);
for(int i=0; i<DIM; i++)
{
for(int j=0; j<DIM; j++)
{
u0(i,j) = 0.0f; // The corresponding color in the colormap for 0 is green
}
}
// Initialize the rendering img to the u0 matrix
CPUsMatrix2Bitmap(u0, renderImg);
// Decayment per timestep
float kdt=k*dt;
// c1 constant
float c1=pow(dt,2)*pow(c,2)/pow(dx,2);
// Droplet as gaussian
// This code creates a gaussian discrete droplet, see the documentation for more information
const int dim = 4*dropletRadius+1;
sMatrix xd(dim, dim);
sMatrix yd(dim, dim);
for(int i=0; i<dim; i++)
{
for(int j=-2*dropletRadius; j<=2*dropletRadius; j++)
{
xd(i,j+2*dropletRadius) = j;
yd(j+2*dropletRadius,i) = j;
}
}
float m_Zd[dim][dim];
for(int i=0; i<dim; i++)
{
for(int j=0; j<dim; j++)
{
// Calculate Gaussian centered on zero
m_Zd[i][j] = -da*exp(-pow(xd(i,j)/dropletRadius,2)-pow(yd(i,j)/dropletRadius,2));
}
}
/* GPU data initialization */
// Allocate memory on the GPU for u and u0 matrices
unsigned int UU0_bytes = DIM*DIM*sizeof(float);
cudaError_t chk;
chk = cudaMalloc((void**)&m_gpuData.gpu_u, UU0_bytes);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT ALLOCATE GPU MEMORY");
return;
}
chk = cudaMalloc((void**)&m_gpuData.gpu_u0, UU0_bytes);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT ALLOCATE GPU MEMORY");
return;
}
// Allocate memory for ris0, ris1, ris2 and ptr matrices
chk = cudaMalloc((void**)&m_gpuData.ris0, UU0_bytes);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT ALLOCATE GPU MEMORY");
return;
}
chk = cudaMalloc((void**)&m_gpuData.ris1, UU0_bytes);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT ALLOCATE GPU MEMORY");
return;
}
chk = cudaMalloc((void**)&m_gpuData.ris2, UU0_bytes);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT ALLOCATE GPU MEMORY");
return;
}
chk = cudaMalloc((void**)&m_gpuData.gpu_ptr, DIM*DIM*4);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT ALLOCATE GPU MEMORY");
return;
}
// Initialize to zero both u and u0
chk = cudaMemcpy(m_gpuData.gpu_u0, u0.values, UU0_bytes, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT TRANSFER MEMORY TO GPU");
return;
}
chk = cudaMemcpy(m_gpuData.gpu_u, u0.values, UU0_bytes, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT TRANSFER MEMORY TO GPU");
return;
}
// Preload Laplacian kernel
float m_D[3][3];
m_D[0][0] = 0.0f; m_D[1][0] = 1.0f; m_D[2][0]=0.0f;
m_D[0][1] = 1.0f; m_D[1][1] = -4.0f; m_D[2][1]=1.0f;
m_D[0][2] = 0.0f; m_D[1][2] = 1.0f; m_D[2][2]=0.0f;
// Copy Laplacian to constant memory
chk = cudaMemcpyToSymbol((const char*)gpu_D, m_D, 9*sizeof(float), 0, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCONSTANT MEMORY TRANSFER FAILED");
return;
}
// Store all static algorithm parameters in constant memory
const float a1 = (2-kdt);
chk = cudaMemcpyToSymbol((const char*)&gpu_a1, &a1, sizeof(float), 0, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCONSTANT MEMORY TRANSFER FAILED");
return;
}
const float a2 = (kdt-1);
chk = cudaMemcpyToSymbol((const char*)&gpu_a2, &a2, sizeof(float), 0, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCONSTANT MEMORY TRANSFER FAILED");
return;
}
chk = cudaMemcpyToSymbol((const char*)&gpu_c1, &c1, sizeof(float), 0, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCONSTANT MEMORY TRANSFER FAILED");
return;
}
const int ddim = dim;
chk = cudaMemcpyToSymbol((const char*)&gpu_Ddim, &ddim, sizeof(int), 0, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCONSTANT MEMORY TRANSFER FAILED");
return;
}
const int droplet_dsz = dropletRadius;
chk = cudaMemcpyToSymbol((const char*)&gpu_dsz, &droplet_dsz, sizeof(int), 0, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCONSTANT MEMORY TRANSFER FAILED");
return;
}
chk = cudaMemcpyToSymbol((constchar*)&gpu_Zd, &m_Zd, sizeof(float)*dim*dim, 0, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCONSTANT MEMORY TRANSFER FAILED");
return;
}
//
// Initialize colormap and ppstep in constant memory
chk = cudaMemcpyToSymbol((const char*)&gpu_pp_step, &pp_step, sizeof(float), 0, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCONSTANT MEMORY TRANSFER FAILED");
return;
}
chk = cudaMemcpyToSymbol((const char*)&gpu_m_colorMap, &m_colorMap, sizeof(unsigned char)*COLOR_NUM*3, 0, cudaMemcpyHostToDevice);
if(chk != cudaSuccess)
{
printf("\nCONSTANT MEMORY TRANSFER FAILED");
return;
}
}
void deinitializeGPUData()
{
// Free everything from device memory
cudaFree(m_gpuData.gpu_u);
cudaFree(m_gpuData.gpu_u0);
cudaFree(m_gpuData.gpu_ptr);
cudaFree(m_gpuData.ris0);
cudaFree(m_gpuData.ris1);
cudaFree(m_gpuData.ris2);
}
在初始化GPU内存后，可以像往常一样启动openglRenderer来调用kernel()函数，以获得一个有效的可渲染表面图像矩阵。但现在有一个区别，就在openglRenderer构造函数中：
openGLRenderer::openGLRenderer(void)
{
//. . .
// Sets up the CUBLAS
cublasStatus_t status = cublasInit();
if (status != CUBLAS_STATUS_SUCCESS)
{
// CUBLAS initialization error
printf("\nCRITICAL: CUBLAS LIBRARY FAILED TO LOAD");
return;
}
// Set up the bitmap data with page-locked memory for fast access
//no more : renderImg = new unsigned char[DIM*DIM*4];
cudaError_t chk = cudaHostAlloc((void**)&renderImg, DIM*DIM*4*sizeof(char), cudaHostAllocDefault);
if(chk != cudaSuccess)
{
printf("\nCRITICAL: CANNOT ALLOCATE PAGELOCKED MEMORY");
return ;
}
} 首先，我们决定使用CUBLAS库执行矩阵加法，原因有两个：
- 我们在设备上的行主序数据已准备好供CUBLAS函数使用（cublasMalloc只是cudaMalloc的包装器）。
- CUBLAS库对于大型矩阵运算经过高度优化；我们的矩阵不是那么大，但这有助于为未来版本扩展架构。
使用我们的sMatrix包装器不再是有效的选择，在设备上工作时我们需要摆脱它，尽管我们仍然可以在初始化阶段使用它。
我们需要注意的第二个基本事情是，在openglRenderer构造函数中，我们使用cudaHostAlloc而不是经典的malloc分配了主机端内存（将包含要渲染的数据的内存）。正如文档所述，使用此函数分配的内存可确保CUDA驱动程序跟踪使用此函数分配的虚拟内存范围，并加速对cudaMemCpy等函数的调用。使用cudaHostAlloc分配的主机内存通常被称为“固定内存”，不能被分页（因此，分配过多的内存可能会降低系统性能，因为它减少了系统可用于分页的内存量）。这一技巧将为设备和主机之间的内存传输带来额外的速度。
我们现在准备好一窥修订后的kernel()函数：
// This kernel is called at each iteration
// It implements the main loop algorithm and someway "rasterize" the matrix data
// to be passed to the openGL renderer. It also adds droplets in the waiting queue
//
void kernel(unsigned char *ptr)
{
// Set up the grid
dim3 blocks(172,172);
dim3 threads(3,3); // 516x516 img is 172x172 (3x3 thread) blocks
// Implements the un = (*data.u)*data.a1 + (*data.u0)*data.a2 + convolution((*data.u),(*data.D))*data.c1;
// line by means of several kernel calls
convolutionGPU<<<blocks,threads>>>(m_gpuData.gpu_u, m_gpuData.ris0);
// Now multiply everything by c1 constant
multiplyEachElementby_c1<<<blocks,threads>>>(m_gpuData.ris0, m_gpuData.ris1);
// First term is ready, now u*a1
multiplyEachElementby_a1<<<blocks,threads>>>(m_gpuData.gpu_u, m_gpuData.ris0);
// u0*a2
multiplyEachElementby_a2<<<blocks,threads>>>(m_gpuData.gpu_u0, m_gpuData.ris2);
// Perform the matrix addition with the CUBLAS library
// un = ris0 + ris2 + ris1
// Since everything is already stored as row-major device vectors, we don't need to do anything to pass it to the CUBLAS
cublasSaxpy(DIM*DIM, 1.0f, m_gpuData.ris0, 1, m_gpuData.ris2, 1);
cublasSaxpy(DIM*DIM, 1.0f, m_gpuData.ris2, 1, m_gpuData.ris1, 1);
// Result is not in m_gpuData.ris1
// Step forward in time
cudaMemcpy(m_gpuData.gpu_u0, m_gpuData.gpu_u, DIM*DIM*sizeof(float), cudaMemcpyDeviceToDevice);
cudaMemcpy(m_gpuData.gpu_u, m_gpuData.ris1, DIM*DIM*sizeof(float), cudaMemcpyDeviceToDevice);
// Draw the u surface matrix and "rasterize" it into gpu_ptr
gpuMatrix2Bitmap<<<blocks,threads>>>(m_gpuData.gpu_u, m_gpuData.gpu_ptr);
// Back on the pagelocked host memory
cudaMemcpy(ptr, m_gpuData.gpu_ptr, DIM*DIM*4, cudaMemcpyDeviceToHost);
if(first_droplet == 1) // By default there's just one initial droplet
{
first_droplet = 0;
int x0d= DIM / 2; // Default droplet center
int y0d= DIM / 2;
cudaMemcpy(m_gpuData.ris0, m_gpuData.gpu_u, DIM*DIM*sizeof(float), cudaMemcpyDeviceToDevice);
addDropletToU<<<blocks,threads>>>(m_gpuData.ris0, x0d,y0d, m_gpuData.gpu_u);
}
// Add all the remaining droplets in the queue
while(dropletQueueCount >0)
{
dropletQueueCount--;
int y0d = DIM - dropletQueue[dropletQueueCount].my;
// Copy from u to one of our buffers
cudaMemcpy(m_gpuData.ris0, m_gpuData.gpu_u, DIM*DIM*sizeof(float), cudaMemcpyDeviceToDevice);
addDropletToU<<<blocks,threads>>>(m_gpuData.ris0, dropletQueue[dropletQueueCount].mx,y0d, m_gpuData.gpu_u);
}
// Synchronize to make sure all kernels executions are done
cudaThreadSynchronize();
}
行：
un = (*data.u)*data.a1 + (*data.u0)*data.a2 + convolution((*data.u),(*data.D))*data.c1;
已被多个内核调用完全取代，这些内核分别执行卷积运算，将矩阵数据乘以算法常数，并通过CUBLAS执行矩阵-矩阵加法。所有操作都在设备上执行，包括点到RGB值映射（这是一个高度可并行化的操作，因为必须对表面图像矩阵中的每个值执行）。通过设备方法实现时间向前推进。最终，数据被复制回分页固定的主机内存，并且队列中等待的液滴将在下一次迭代中添加到u表面模拟数据矩阵中。
kernel()函数调用的CUDA内核如下：
/******************************************************************************
CUDA KERNELS
/******************************************************************************/
// For a 512x512 image the grid is 170x170 blocks 3x3 threads each one
__global__ void convolutionGPU(float *A, float *result)
{
__shared__ float data[laplacianD*2][laplacianD*2];
// Absolute position into the image
const int gLoc = threadIdx.x + IMUL(blockIdx.x,blockDim.x) + IMUL(threadIdx.y,DIM) + IMUL(blockIdx.y,blockDim.y)*DIM;
// Image-relative position
const int x0 = threadIdx.x + IMUL(blockIdx.x,blockDim.x);
const int y0 = threadIdx.y + IMUL(blockIdx.y,blockDim.y);
// Load the apron and data regions
int x,y;
// Upper left square
x = x0 - kernelRadius;
y = y0 - kernelRadius;
if(x < 0 || y < 0)
data[threadIdx.x][threadIdx.y] = 0.0f;
else
data[threadIdx.x][threadIdx.y] = A[ gLoc - kernelRadius - IMUL(DIM,kernelRadius)];
// Upper right square
x = x0 + kernelRadius + 1;
y = y0 - kernelRadius;
if(x >= DIM || y < 0)
data[threadIdx.x + blockDim.x][threadIdx.y] = 0.0f;
else
data[threadIdx.x + blockDim.x][threadIdx.y] = A[ gLoc + kernelRadius+1 - IMUL(DIM,kernelRadius)];
// Lower left square
x = x0 - kernelRadius;
y = y0 + kernelRadius+1;
if(x < 0 || y >= DIM)
data[threadIdx.x][threadIdx.y + blockDim.y] = 0.0f;
else
data[threadIdx.x][threadIdx.y + blockDim.y] = A[ gLoc - kernelRadius + IMUL(DIM,(kernelRadius+1))];
// Lower right square
x = x0 + kernelRadius+1;
y = y0 + kernelRadius+1;
if(x >= DIM || y >= DIM)
data[threadIdx.x + blockDim.x][threadIdx.y + blockDim.y] = 0.0f;
else
data[threadIdx.x + blockDim.x][threadIdx.y + blockDim.y] = A[ gLoc + kernelRadius+1 + IMUL(DIM,(kernelRadius+1))];
__syncthreads();
float sum = 0;
x = kernelRadius + threadIdx.x;
y = kernelRadius + threadIdx.y;
// Execute the convolution in the shared memory (kernel is in constant memory)
#pragma unroll
for(int i = -kernelRadius; i<=kernelRadius; i++)
for(int j=-kernelRadius; j<=kernelRadius; j++)
sum += data[x+i][y+j] * gpu_D[i+kernelRadius][j+kernelRadius];
// Transfer the risult to global memory
result[gLoc] = sum;
}
__global__ void multiplyEachElementby_c1(float *matrix, float *result)
{
// Absolute position into the image
const int gLoc = threadIdx.x + IMUL(blockIdx.x,blockDim.x) + IMUL(threadIdx.y,DIM) + IMUL(blockIdx.y,blockDim.y)*DIM;
// Multiply by c1 each matrix's element
result[gLoc] = matrix[gLoc]*gpu_c1;
}
__global__ void multiplyEachElementby_a1(float *matrix, float *result)
{
// Absolute position into the image
const int gLoc = threadIdx.x + IMUL(blockIdx.x,blockDim.x) + IMUL(threadIdx.y,DIM) + IMUL(blockIdx.y,blockDim.y)*DIM;
// Multiply by c1 each matrix's element
result[gLoc] = matrix[gLoc]*gpu_a1;
}
__global__ void multiplyEachElementby_a2(float *matrix, float *result)
{
// Absolute position into the image
const int gLoc = threadIdx.x + IMUL(blockIdx.x,blockDim.x) + IMUL(threadIdx.y,DIM) + IMUL(blockIdx.y,blockDim.y)*DIM;
// Multiply by c1 each matrix'selement
result[gLoc] = matrix[gLoc]*gpu_a2;
}
// Associate a colormap RGB value to each point
__global__ void gpuMatrix2Bitmap(float *matrix, BYTE *bitmap)
{
// Absolute position into the image
const int gLoc = threadIdx.x + IMUL(blockIdx.x,blockDim.x) + IMUL(threadIdx.y,DIM) + IMUL(blockIdx.y,blockDim.y)*DIM;
int cvalue = (int)((matrix[gLoc] + 1.0f)/gpu_pp_step);
if(cvalue < 0)
cvalue = 0;
else if(cvalue >= COLOR_NUM)
cvalue = COLOR_NUM-1;
bitmap[gLoc*4] = gpu_m_colorMap[cvalue][0];
bitmap[gLoc*4 + 1] = gpu_m_colorMap[cvalue][1];
bitmap[gLoc*4 + 2] = gpu_m_colorMap[cvalue][2];
bitmap[gLoc*4 + 3] = 0xFF; // Alpha
}
// Add a gaussian 2D droplet matrix to the surface data
// Warning: this kernel has a high divergence factor, it is meant to be seldom called
__global__ void addDropletToU(float *matrix, int x0d, int y0d, float *result)
{
// Absolute position into the image
const int gLoc = threadIdx.x + IMUL(blockIdx.x,blockDim.x) + IMUL(threadIdx.y,DIM) + IMUL(blockIdx.y,blockDim.y)*DIM;
// Image relative position
const int x0 = threadIdx.x + IMUL(blockIdx.x,blockDim.x);
const int y0 = threadIdx.y + IMUL(blockIdx.y,blockDim.y);
// Place the (x0d;y0d) centered Zd droplet on the wave data (it will be added at the next iteration)
if (x0 >= x0d-gpu_dsz*2 && y0 >= y0d-gpu_dsz*2 && x0 <= x0d+gpu_dsz*2 && y0 <= y0d+gpu_dsz*2)
{
// Add to result the matrix value plus the Zd corresponding value
result[gLoc] = matrix[gLoc] + gpu_Zd[x0 -(x0d-gpu_dsz*2)][y0 - (y0d-gpu_dsz*2)];
}
else
result[gLoc] = matrix[gLoc]; // This value shouln't be changed
}
请注意，我们宁愿“硬编码”常量值的用法与不同的内核，而不是引入条件分支导致发散。唯一增加线程发散的内核是addDropletToU，因为只有少数线程实际执行高斯包-启动例程（参见前面几段描述的理论算法），但这不成问题，因为其调用频率很低。
性能对比
计时测量和性能比较已在以下系统上进行：
Intel 2 quad cpu Q9650 @ 3.00 Ghz
6 GB RAM
64 位操作系统
NVIDIA GeForce GTX 285 (1GB DDR3 @ 1476 Mhz, 240 CUDA 核)
我们用于编译项目的CUDA版本是4.2，如果您遇到问题，请确保安装正确的版本或按照README文件中的说明进行更改。
为了对CUDA核执行进行基准测试，我们使用了CUDA版本附带的cudaEventCreate / cudaEventRecord / cudaEventSynchronize / cudaEventElapsedTime函数；而对于CPU版本，我们使用了两个Windows平台相关的API：QueryPerformanceFrequency和QueryPerformanceCounter。
我们将基准测试分为四个阶段：启动阶段，此时图像中唯一的液滴是默认液滴；第二阶段，CPU和GPU版本都稳定下来；第三阶段，我们向渲染队列添加60-70个液滴；第四阶段，应用程序运行15-20分钟。我们发现，在所有测试中，CPU的性能都比GPU版本差，GPU版本可以依赖一个大的线程网格，准备好拆分即将到来的显著工作负载并提供固定的渲染时间。另一方面，在长期运行中，虽然GPU仍然表现更好，但CPU版本显示出微小的性能提升，也许得益于缓存机制。
请注意，处理更大数据的应用程序肯定会从大规模并行化方法中获得更大的优势。我们的波动方程PDE模拟确实非常简单，并且不需要显著的工作负载，因此限制了可以实现的性能增益。
一劳永逸：将顺序设计的算法转换为并行算法没有通用规则，每个情况都必须在自己的上下文和架构中进行评估。此外，使用CUDA可以在科学模拟中提供巨大优势，但不应将GPU视为CPU的替代品，而应将其视为可以依赖大规模数据并行化的代数协处理器。结合CPU的顺序代码部分和GPU的并行代码部分是成功的关键。
CUDA核函数最佳实践
本文的最后（但并非最不重要）部分提供了一个清单，列出了在编写CUDA内核时应遵循的最佳实践和要避免的错误，以最大限度地发挥GPU加速应用程序的性能。
- 最小化主机<->设备传输，特别是设备->主机传输（后者较慢），即使这意味着在设备上运行的内核在CPU上不会更慢。
- 在主机端使用固定内存（分页锁定）以利用带宽优势。但要注意不要滥用它，否则会降低整个系统的性能。
- cudaMemcpy是一个阻塞函数，如果可能，请异步使用它，并结合固定内存+CUDA流，以便将传输与内核执行重叠。
- 如果您的显卡是集成显卡，零拷贝内存（即使用cudaHostAllocMapped标志分配的固定内存）始终能带来优势。如果不是，则不确定，因为GPU不会缓存该内存。
- 始终检查设备的计算能力（使用cudaGetDeviceProperties），如果< 2.0，则每个块不能使用超过512个线程（最多65535个块）。
- 查看显卡规格：当启动一个AxB块网格内核时，每个SM（流多处理器）可以服务一部分。如果您的卡最多可以容纳1024个线程/SM，您应该调整块的大小，以便尽可能多地填充它们，但不要太多（否则会产生调度延迟）。每个warp通常是32个线程（尽管这不是固定值，并且取决于体系结构），并且是每个SM的最小调度单元（在任何时候只有一个warp执行，并且warp中的所有线程执行相同的指令-SIMD），所以您应该考虑以下示例：在一块GT200卡上，您需要执行矩阵乘法。您应该使用8x8、16x16还是32x32线程/块？
对于8x8块，我们有64个线程/块。由于每个SM最多可以容纳1024个线程，因此每个SM可以容纳16个块（1024/64）。然而，每个SM最多只能容纳8个块。因此，每个SM只有512个线程（64*8）进入每个SM -> SM执行资源未充分利用；用于调度长延迟操作的wrap数量较少。
对于16x16块，我们有256个线程/块。由于每个SM最多可以容纳1024个线程，因此它可以容纳多达4个块（1024/256），并且不会达到8个块的限制 -> 每个SM的线程容量已满，并且用于调度长延迟操作的wrap数量最大（1024/32=32个wrap）。
对于32x32块，我们有1024个线程/块 -> 甚至一个也无法放入SM！（每个块有512个线程的限制）。
- 检查每个块的SM寄存器限制，并除以线程数：您将得到一个线程中可以使用的最大寄存器数。仅一个线程超过此限制将导致更少的warp被调度，性能下降。
- 检查每个块的共享内存，并确保不超过该值。超过此值将导致更少的warp被调度，性能下降。
- 始终检查线程数是否小于设备支持的最大线程数。
- 在可能的情况下，使用归约（reduction）和局部性（locality）技术。
- 如果可能，切勿将内核代码分割成条件分支（if-then-else），不同的路径将导致同一个warp执行更多的次数，从而产生以下开销。
- 在可能的情况下，使用具有发散最小化的归约技术（即，尝试使用warp执行每次周期合并读取的归约，如Kirk和Hwu书籍的第6章所述）。
- 合并是通过强制硬件读取连续数据来实现的。如果warp中的每个线程访问连续内存，则合并会显著增加（warp中一半的线程应同时访问全局内存），这就是为什么对于大型矩阵（行主序），按列读取比按行读取更好。通过局部性（即线程协作加载其他线程所需的数据）可以增加合并。内核应执行具有局部性目的的合并读取，以最大化内存吞吐量。将值存储到共享内存（如果可能）也是一个好习惯。
- 确保您没有未使用的线程/块，无论是由于设计原因还是代码原因。如前所述，显卡有SM上的最大线程数和每个SM的最大块数等限制。设计网格时不咨询卡规格强烈不推荐。
- 如前所述，添加比卡最大寄存器限制更多的寄存器是性能损失的根源。不过，添加寄存器也可能导致指令的添加，即：更多的时间来并行化传输或调度warp，以及更好的性能。再说一遍：没有通用规则，您应该遵循设计代码时的最佳实践，然后自己进行实验。
- 数据预取是指预加载您当前不需要的数据，以在未来操作中获得性能（与局部性密切相关）。将数据预取与矩阵瓦片（tiles）结合起来可以解决许多长延迟内存访问问题。
- 如果适用（例如，循环很小），则最好展开循环。理想情况下，循环应由编译器自动展开，但始终检查以确保这一点是更好的选择。
- 在处理矩阵时，使用矩形瓦片（Kirk和Hwu书的第6章）减小线程粒度，以避免不同块从全局内存中多次读取行/列。
- 纹理是缓存内存，如果使用得当，它们比全局内存快得多，也就是说：纹理更适合二维空间局部性访问（例如，多维数组），并在某些特定情况下表现更好（同样，这是一个逐例规则）。
- 尽量使至少25%的总寄存器被占用，否则无法通过其他warp的计算来隐藏访问延迟。
- 每个块的线程数应始终是warp大小的倍数，以有利于合并读取。
- 通常，如果一个SM支持多个块，则应使用超过64个线程/块。
- 开始尝试内核网格的起始值是每个块128到256个线程。
- 高延迟问题可以通过使用更多块、每个SM（流多处理器）上的线程更少来解决，而不是只使用一个块、每个SM上的线程很多。这对于经常调用__syncthreads()的内核尤其重要。
- 如果内核失败，请使用cudaGetLastError()检查错误值，您可能使用了过多的寄存器/过多的共享或常量内存/过多的线程。
- CUDA函数和硬件在float数据类型下表现最佳。强烈建议使用它。
- 整数除法和模运算符(%)是昂贵的运算符，尽可能用位运算替换它们。如果n是2的幂，
- 如果可能，避免从double到float的自动数据转换。
- 前面带有__的数学函数是硬件实现的，它们更快但精度较低。如果精度不是关键目标，可以使用它们（例如，__sinf(x)而不是sinf(x)）。
- 在循环中使用有符号整数而不是无符号整数，因为某些编译器对有符号整数的优化更好（有符号整数溢出会导致未定义行为，编译器可能会对此进行激进优化）。
- 请记住，由于舍入误差，浮点数学不是可结合的。因此，大规模并行计算的结果可能会有所不同。
- 1) 如果您的设备支持并发内核执行（请参阅 concurrentKernels 设备属性），您可能可以通过并行运行内核来获得额外的性能。请检查您的设备规格和CUDA编程指南。
本文到此结束。目的是探讨GPGPU应用程序在改进科学模拟和旨在处理大数据操作的程序方面的能力。
参考文献
- Mathematics for 3D Game Programming and Computer Graphics Third Edition - Eric Lengyel
- Programming Massively Parallel Processors – David B. Kirk, Wen-mei W. Hwu
- CUDA by Example - Jason Sanders, Edward Kandrot
- Analyzing Waves on a String - Michael Fowler (http://galileo.phys.virginia.edu/classes/152.mf1i.spring02/)
- Frank Wattenberg - Estimating Second Derivatives Numerically
- Maxim Vedenev (http://simulations.narod.ru/) - solve 2d wave equation with FDM
- 2008 Fall CS 525 CUDA Convolution - Sangyoon (James) Lee, University of Illinois (http://www.evl.uic.edu/sjames/cs525/final.html)