SoCVista深入浅出谈CUDA
深入浅出谈CUDA
Hotball
CUDA 是 NVIDIA 的 GPGPU 模型，它使用 C 语言为基础，可以直接以大多数人熟悉的 C
”
CUDA是什么？能吃吗？
编者注：NVIDIA的GeFoce 8800GTX发布后，它的通用计算架构CUDA经过一年多的推广后，
易了解CUDA，我们征得Hotball的本人同意，发表他最近亲自撰写的本文。这篇文章的特点是深入浅出，也包含了hotball本人编写一些简单CUDA程序的亲身体验，对于希望了解CUDA的读者来说是非常不错的入门文章，PCINLIFE对本文的发表没有作任何的删减，主要是把一些台湾的词汇转换成大陆的词汇以及作了若干"编者注"的注释。
现代的显示芯片已经具有高度的可程序化能力，由于显示芯片通常具有相当高的内存带宽，以及大量的执行单元，因此开始有利用显示芯片来帮助进行一些计算工作的想法，即 GPGPU。CUDA 即是NVIDIA 的 GPGP U 模型。
NVIDIA 的新一代显示芯片，包括 GeForce 8 系列及更新的显示芯片都支持 CUDA。NVIDIA 免费提供 CUDA 的开发工具（包括 Wind ows 版本和 Li nux 版本）、程序范例、文件等等，可以在CUDA Zone 下载。
GPGPU 的优缺点
使用显示芯片来进行运算工作，和使用 CPU 相比，主要有几个好处：
1.显示芯片通常具有更大的内存带宽。例如，NVIDIA 的 GeF orce 8800GTX 具有超过
50GB/s 的内存带宽，而目前高阶 CPU 的内存带宽则在 10GB/s 左右。.
2.显示芯片具有更大量的执行单元。例如 G eForce 8800GTX 具有 128 个 "strea m
processors"，频率为 1.35GHz。CPU 频率通常较高，但是执行单元的数目则要少得多。
3.和高阶 CPU相比，显卡的价格较为低廉。例如目前一张 G eForce 8800GT 包括
512MB 内存的价格，和一颗 2.4GHz 四核心 CP U 的价格相若。
当然，使用显示芯片也有它的一些缺点：
1.显示芯片的运算单元数量很多，因此对于不能高度并行化的工作，所能带来的帮助就
不大。
2.显示芯片目前通常只支持 32 bits 浮点数，且多半不能完全支持 IEEE 754 规格，有
些运算的精确度可能较低。目前许多显示芯片并没有分开的整数运算单元，因此整数运算的效率较差。
3.显示芯片通常不具有分支预测等复杂的流程控制单元，因此对于具有高度分支的程序，
效率会比较差。
4.目前 GPGP U 的程序模型仍不成熟，也还没有公认的标准。例如 NVIDIA和
AMD/ATI 就有各自不同的程序模型。
整体来说，显示芯片的性质类似 str eam processor，适合一次进行大量相同的工作。CPU 则比较有弹性，能同时进行变化较多的工作。
CUDA 架构
CUDA 是 NVIDI A 的 GPGPU 模型，它使用 C 语言为基础，可以直接以大多数人熟悉的 C 语言，写出在显示芯片上执行的程序，而不需要去学习特定的显示芯片的指令或是特殊的结构。
在 CUDA 的架构下，一个程序分为两个部份：host 端和 devi ce 端。Host 端是指在 CPU 上执行的部份，而 de vice 端则是在显示芯片上执行的部份。Device 端的程序又称为 "kern el"。通常 hos t 端程序会将数据准备好后，复制到显卡的内存中，再由显示芯片执行 de vice 端程序，完成后再由 ho st 端程序将结果从显卡的内存中取回。
由于 CPU 存取显卡内存时只能透过 PCI Express 接口，因此速度较慢（PCI Express x16 的理论带宽是双向各 4G B/s），因此不能太常进行这类动作，以免降低效率。
在 CUDA 架构下，显示芯片执行时的最小单位是thread。数个 threa d 可以组成一个block。一个 bl ock 中的 thre ad 能存取同一块共享的内存，而且可以快速进行同步的动作。
每一个 bloc k 所能包含的 thre ad 数目是有限的。不过，执行相同程序的 blo ck，可以组成grid。不同 b lock 中的 t hread 无法存取同一个共享的内存，因此无法直接互通或进行同步。因此，不同 b lock 中的 t hread 能合作的程度是比较低的。不过，利用这个模式，可以让程序不用担心显示芯片实际上能同时执行的 thre ad 数目限制。例如，一个具有很少量执行单元的显示芯片，可能会把各个 bl ock 中的 th read 顺序执行，而非同时执行。不同的 gri d 则可以执行不同的程序（即 ker nel）。
Grid、block 和 thre ad 的关系，如下图所示：
每个 thre ad 都有自己的一份 register和 lo cal me mory 的空间。同一个 bl ock 中的每个thread 则有共享的一份 share memory。此外，所有的 thre ad（包括不同 b lock 的 thre ad）都共享一份 g lobal memory、constant memory、和 text ure memory。不同的 gri d 则有各自的 gl obal memory、constant memory 和 te xture memory。这些不同的内存的差别，会在之后讨论。
执行模式
由于显示芯片大量并行计算的特性，它处理一些问题的方式，和一般 CPU 是不同的。主要的特点包括：
1.内存存取 la tency 的问题：CPU 通常使用 ca che 来减少存取主内存的次数，以避免
内存 la tency 影响到执行效率。显示芯片则多半没有 ca che（或很小），而利用并行化执行的方式来隐藏内存的 la tency（即，当第一个 thre ad 需要等待内存读取结果时，则开始执行第二个 t hread，依此类推）。
2.分支指令的问题：CPU 通常利用分支预测等方式来减少分支指令造成的 p ipeline
bubble。显示芯片则多半使用类似处理内存 la tency 的方式。不过，通常显示芯片处理分支的效率会比较差。
因此，最适合利用 CUDA 处理的问题，是可以大量并行化的问题，才能有效隐藏内存的latency，并有效利用显示芯片上的大量执行单元。使用 CUDA 时，同时有上千个 t hread 在执行是很正常的。因此，如果不能大量并行化的问题，使用 CUDA 就没办法达到最好的效率了。
CUDA Toolkit的安装
目前 NVIDI A 提供的 CUDA Toolkit（可从这里下载）支持 Wind ows （32 bits 及 64 bits 版本）及许多不同的 Li nux 版本。
CUDA Toolkit 需要配合 C/C++ compiler。在 Wi ndows 下，目前只支持 V isual Studio 7.x 及Visual St udio 8（包括免费的Visual Studio C++ 2005 Ex press）。Visual S tudio 6 和 gcc 在Windows 下是不支援的。在 Li nux 下则只支援 gcc。
这里简单介绍一下在 Wi ndows 下设定并使用 CUDA 的方式。
下载及安装
在 Win dows 下，CUDA Toolkit 和 CU DA SDK 都是由安装程序的形式安装的。CUDA Toolkit 包括 CUDA 的基本工具，而 CUDA SDK 则包括许多范例程序以及链接库。基本上要写CUDA 的程序，只需要安装 CUDA Toolkit 即可。不过 CUDA SDK 仍值得安装，因为里面的许多范例程序和链接库都相当有用。
CUDA Toolkit 安装完后，预设会安装在 C:\CUDA 目录里。其中包括几个目录：
?bin -- 工具程序及动态链接库
?doc -- 文件
?include -- header 檔
?lib -- 链接库档案
?open64 -- 基于 Op en64 的 CUDA compiler
?src -- 一些原始码
安装程序也会设定一些环境变量，包括：
?CUDA_BIN_PATH -- 工具程序的目录，默认为 C:\CUDA\b in
?CUDA_INC_PATH -- header 文件的目录，默认为 C:\CUDA\in c
?CUDA_LIB_PATH -- 链接库文件的目录，默认为 C:\CUDA\li b
在 Visual Studio 中使用 CUDA
CUDA 的主要工具是 nvcc，它会执行所需要的程序，将 CUDA 程序代码编译成执行档(或object 檔) 。在 Visu al Studio 下，我们透过设定custom build tool 的方式，让 Visu al Studio 会自动执行 nvcc。
这里以Visual Studio 2005 为例：
1.首先，建立一个 W in32 Console 模式的 pro ject（在 A pplication Settings 中记得勾选
Empty project），并新增一个档案，例如 ma in.cu。
2.在 ma in.cu 上右键单击，并选择Properties。点选General，确定Tool的部份是选
择Custom Build Tool。
3.选择Custom Build Step，在 C ommand Line 使用以下设定：
o Release 模式："$(CUDA_BIN_PATH)\nvcc.exe" -ccbin "$(VCIn stallDir)bin" -c
-DWIN32 -D_CONSOLE -D_MBCS -Xcompiler /EHsc,/W3,/nologo,/Wp64,/O2,/Zi,/MT
-I"$(CUDA_INC_PATH)" -o $(ConfigurationName)\$(InputName).obj $(InputFileName) o Debug 模式："$(CUDA_BIN_PATH)\nvcc.exe" -ccb in "$(VCIn stallDir)bin" -c
-D_DEBUG -D WIN32 -D_CONSOLE -D_MBCS -Xco mpiler
/EHsc,/W3,/nologo,/Wp64,/Od,/Zi,/RTC1,/MTd -I"$(CUDA_INC_PATH)" -o
$(ConfigurationName)\$(InputName).obj $(InputFileName)
4.如果想要使用软件仿真的模式，可以新增两个额外的设定：
o EmuRelease 模式："$(CUDA_BIN_PATH)\nvcc.exe" -ccb in "$(VCInstallDir)bin"
-deviceemu -c -DWIN32 -D_C ONSOLE -D_MBCS -Xcompiler
/EHsc,/W3,/nologo,/Wp64,/O2,/Zi,/MT -I"$(CUDA_INC_PATH)" -o
$(ConfigurationName)\$(InputName).obj $(InputFileName)
o EmuDebug 模式："$(CUDA_BIN_PATH)\nvcc.exe" -ccb in "$(VCInstallDir)bin"
-deviceemu -c -D_DEBUG -DWIN32 -D_CONSOLE -D_MBCS -Xc ompiler
/EHsc,/W3,/nologo,/Wp64,/Od,/Zi,/RTC1,/MTd -I"$(CUDA_INC_PATH)" -o
$(ConfigurationName)\$(InputName).obj $(InputFileName)
5.对所有的配置文件，在Custom Build Step的Outputs中加入
$(ConfigurationName)\$(InputName).obj。
6.选择 project，右键单击选择Properties，再点选Linker。对所有的配置文件修改以
下设定：
o General/Enable Incremental Linking：No
o General/Additional Library Directories：$(CUDA_LIB_PATH)
o Input/Additional Dependencies：cudart.lib
这样应该就可以直接在 Visu al Studio 的 IDE 中，编辑 CUDA 程序后，直接 bui ld 以及执行程序了。
第一个CUDA程序
CUDA 目前有两种不同的 API：Runtime API 和 Driv er API，两种 API 各有其适用的范围。由于 runtim e API 较容易使用，一开始我们会以 rune time API 为主。
CUDA 的初始化
首先，先建立一个档案 fir st_cuda.cu。如果是使用 V isual St udio 的话，则请先按照这里的设定方式设定 proje ct。
要使用 runtime API 的时候，需要 in clude cuda_runtime.h。所以，在程序的最前面，加上
#include
#include
接下来是一个 In itCUDA 函式，会呼叫 runtim e API 中，有关初始化 CUDA 的功能：
bool InitCUDA()
{
int count;
cudaGetDeviceCount(&count);
if(count == 0) {
fprintf(stderr, "There is no device.\n");
return false;
}
int i;
for(i = 0; i < count; i++) {
cudaDeviceProp prop;
if(cudaGetDeviceProperties(&prop, i) == cudaSuccess) {
if(prop.major >= 1) {
break;
}
}
}
if(i == count) {
fprintf(stderr, "There is no device supporting CUDA 1.x.\n");
return false;
}
cudaSetDevice(i);
return true;
}
这个函式会先呼叫 cuda GetDeviceCount 函式，取得支持 CUDA 的装置的数目。如果系统上没有支持 CUDA 的装置，则它会传回1，而 d evice 0 会是一个仿真的装置，但不支持 CUDA 1.0 以上的功能。所以，要确定系统上是否有支持 CUDA 的装置，需要对每个 dev ice 呼叫cudaGetDeviceProperties 函式，取得装置的各项数据，并判断装置支持的 CUDA 版本（prop.major 和 prop.minor 分别代表装置支持的版本号码，例如 1.0 则 prop.m ajor 为 1 而prop.minor 为0）。
透过 cudaGe tDeviceProperties 函式可以取得许多数据，除了装置支持的 CUDA 版本之外，还有装置的名称、内存的大小、最大的 thre ad 数目、执行单元的频率等等。详情可参考NVIDIA 的CUDA Programming Guide。
在找到支持 CUDA 1.0 以上的装置之后，就可以呼叫 cuda SetDevice 函式，把它设为目前要使用的装置。
最后是 ma in 函式。在 main 函式中我们直接呼叫刚才的 In itCUDA 函式，并显示适当的讯息：
int main()
{
if(!InitCUDA()) {
return 0;
}
printf("CUDA initialized.\n");
return 0;
}
这样就可以利用 nv cc 来 com pile 这个程序了。使用 Visu al Studio 的话，若按照先前的设定方式，可以直接 Bu ild Project 并执行。
nvcc 是 CU DA 的 com pile 工具，它会将 .c u 檔拆解出在 GPU 上执行的部份，及在 ho st 上执行的部份，并呼叫适当的程序进行 com pile 动作。在 GP U 执行的部份会透过 NVIDI A 提供的 com piler 编译成中介码，而 ho st 执行的部份则会透过系统上的 C++ compiler 编译（在Windows 上使用 Visua l C++ 而在 Li nux 上使用 gcc）。
编译后的程序，执行时如果系统上有支持 CUDA的装置，应该会显示 CUDA initialized. 的讯息，否则会显示相关的错误讯息。
利用 CUDA 进行运算
到目前为止，我们的程序并没有做什么有用的工作。所以，现在我们加入一个简单的动作，就是把一大堆数字，计算出它的平方和。
首先，把程序最前面的 i nclude 部份改成：
#include
#include
#include
#define DATA_SIZE 1048576
int data[DATA_SIZE];
并加入一个新函式 Gen erateNumbers：
void GenerateNumbers(int *number, int size)
{
for(int i = 0; i < size; i++) {
number[i] = rand() % 10;
}
}
这个函式会产生一大堆 0 ~ 9 之间的随机数。
要利用 CUDA 进行计算之前，要先把数据复制到显卡内存中，才能让显示芯片使用。因此，需要取得一块适当大小的显卡内存，再把产生好的数据复制进去。在 ma in 函式中加入：GenerateNumbers(data, DATA_SIZE);
int* gpudata, *result;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
上面这段程序会先呼叫 Gen erateNumbers 产生随机数，并呼叫 cu daMalloc 取得一块显卡内存（result 则是用来存取计算结果，在稍后会用到），并透过 cudaMem cpy 将产生的随机数复制到显卡内存中。cudaMalloc 和 cu daMemcpy 的用法和一般的 m alloc 及 me mcpy 类似，不过 cudaMem cpy 则多出一个参数，指示复制内存的方向。在这里因为是从主内存复制到显卡内存，所以使用 cudaM emcpyHostToDevice。如果是从显卡内存到主内存，则使用cudaMemcpyDeviceToHost。这在之后会用到。
接下来是要写在显示芯片上执行的程序。在 CUDA 中，在函式前面加上 __global__ 表示这个函式是要在显示芯片上执行的。因此，加入以下的函式：
__global__ static void sumOfSquares(int *num, int* result)
{
int sum = 0;
int i;
for(i = 0; i < DATA_SIZE; i++) {
sum += num[i] * num[i];
}
*result = sum;
}
在显示芯片上执行的程序有一些限制，例如它不能有传回值。其它的限制会在之后提到。
接下来是要让 CUDA 执行这个函式。在 CUDA 中，要执行一个函式，使用以下的语法：函式名称<<
呼叫完后，还要把结果从显示芯片复制回主内存上。在 ma in 函式中加入以下的程序：sumOfSquares<<<1, 1, 0>>>(gpudata, result);
int sum;
cudaMemcpy(&sum, result, sizeof(int), cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
printf("sum: %d\n", sum);
因为这个程序只使用一个 thread，所以 bloc k 数目、thread 数目都是1。我们也没有使用到任何 shared memory，所以设为0。编译后执行，应该可以看到执行的结果。
为了确定执行的结果正确，我们可以加上一段以 CPU 执行的程序代码，来验证结果：sum = 0;
for(int i = 0; i < DATA_SIZE; i++) {
sum += data[i] * data[i];
}
printf("sum (CPU): %d\n", sum);
编译后执行，确认两个结果相同。
计算运行时间
CUDA 提供了一个 cl ock 函式，可以取得目前的 ti mestamp，很适合用来判断一段程序执行所花费的时间（单位为 GPU 执行单元的频率）。这对程序的优化也相当有用。要在我们的程序中记录时间，把 sumO fSquares 函式改成：
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
int sum = 0;
int i;
clock_t start = clock();
for(i = 0; i < DATA_SIZE; i++) {
sum += num[i] * num[i];
}
*result = sum;
*time = clock() - start;
}
把 ma in 函式中间部份改成：
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int));
cudaMalloc((void**) &time, sizeof(clock_t));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<1, 1, 0>>>(gpudata, result, time);
int sum;
clock_t time_used;
cudaMemcpy(&sum, result, sizeof(int), cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t),
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
printf("sum: %d time: %d\n", sum, time_used);
编译后执行，就可以看到执行所花费的时间了。
如果计算实际运行时间的话，可能会注意到它的执行效率并不好。这是因为我们的程序并没有利用到 CUDA 的主要的优势，即并行化执行。在下一段文章中，会讨论如何进行优化的动作。
改良第一个 CUDA程序
在上一篇文章中，我们做了一个计算一大堆数字的平方和的程序。不过，我们也提到这个程序的执行效率并不理想。当然，实际上来说，如果只是要做计算平方和的动作，用 CPU 做会比用 GPU 快得多。这是因为平方和的计算并不需要太多运算能力，所以几乎都是被内存带宽所限制。因此，光是把数据复制到显卡内存上的这个动作，所需要的时间，可能已经和直接在 CPU 上进行计算差不多了。
不过，如果进行平方和的计算，只是一个更复杂的计算过程的一部份的话，那么当然在 GPU 上计算还是有它的好处的。而且，如果数据已经在显卡内存上（例如在 GPU 上透过某种算法产生），那么，使用 GPU 进行这样的运算，还是会比较快的。
刚才也提到了，由于这个计算的主要瓶颈是内存带宽，所以，理论上显卡的内存带宽是相当大的。这里我们就来看看，倒底我们的第一个程序，能利用到多少内存带宽。
程序的并行化
我们的第一个程序，并没有利用到任何并行化的功能。整个程序只有一个 thre ad。在 GeFo rce 8800GT 上面，在 GPU 上执行的部份（称为"kernel"）大约花费 640M 个频率。GeForce 8800GT 的执行单元的频率是 1.5GHz，因此这表示它花费了约 0.43 秒的时间。1M 个 32 bits 数字的数据量是 4MB，因此，这个程序实际上使用的内存带宽，只有 9.3MB/s 左右！这是非常糟糕的表现。
为什么会有这样差的表现呢？这是因为 GPU 的架构特性所造成的。在 CUDA 中，一般的数据复制到的显卡内存的部份，称为global memory。这些内存是没有 cache 的，而且，存取global memory 所需要的时间（即 lat ency）是非常长的，通常是数百个 cycles。由于我们的程序只有一个 th read，所以每次它读取 g lobal memory 的内容，就要等到实际读取到数据、累加到 sum之后，才能进行下一步。这就是为什么它的表现会这么的差。
由于 gl obal memory 并没有 cac he，所以要避开巨大的 la tency 的方法，就是要利用大量的threads。假设现在有大量的 threa ds 在同时执行，那么当一个 thre ad 读取内存，开始等待结果的时候，GPU 就可以立刻切换到下一个 thre ad，并读取下一个内存位置。因此，理想上当thread 的数目够多的时候，就可以完全把 gl obal memory 的巨大 la tency 隐藏起来了。
要怎么把计算平方和的程序并行化呢？最简单的方法，似乎就是把数字分成若干组，把各组数字分别计算平方和后，最后再把每组的和加总起来就可以了。一开始，我们可以把最后加总的动作，由 CPU 来进行。
首先，在 first_cuda.cu 中，在 #d efine DATA_SIZE 的后面增加一个 #d efine，设定 th read 的数目：
#define DATA_SIZE 1048576
#define THREAD_NUM 256
接着，把 ker nel 程序改成：
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
const int tid = threadIdx.x;
const int size = DATA_SIZE / THREAD_NUM;
int sum = 0;
int i;
clock_t start;
if(tid == 0) start = clock();
for(i = tid * size; i < (tid + 1) * size; i++) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0) *time = clock() - start;
}
程序里的 thr eadIdx 是 CUDA 的一个内建的变量，表示目前的 t hread 是第几个 thre ad（由 0 开始计算）。以我们的例子来说，会有 256 个 thre ads，所以同时会有 256 个 sumOfSquares 函式在执行，但每一个的 thre adIdx.x 则分别会是0 ~ 255。利用这个变量，我们就可以让每
一份函式执行时，对整个数据不同的部份计算平方和。另外，我们也让计算时间的动作，只在 thre ad 0（即 thre adIdx.x = 0 的时候）进行。
同样的，由于会有 256 个计算结果，所以原来存放 result 的内存位置也要扩大。把 ma in 函式中的中间部份改成：
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int) * THREAD_NUM);
cudaMalloc((void**) &time, sizeof(clock_t));
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<<1, THREAD_NUM, 0>>>(gpudata, result, time);
int sum[THREAD_NUM];
clock_t time_used;
cudaMemcpy(&sum, result, sizeof(int) * THREAD_NUM,
cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t),
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
可以注意到我们在呼叫 su mOfSquares 函式时，指定 THRE AD_NUM 为 threa d 的数目。最后，在 CPU 端把计算好的各组数据的平方和进行加总：
int final_sum = 0;
for(int i = 0; i < THREAD_NUM; i++) {
final_sum += sum[i];
}
printf("sum: %d time: %d\n", final_sum, time_used);
final_sum = 0;
for(int i = 0; i < DATA_SIZE; i++) {
sum += data[i] * data[i];
}
printf("sum (CPU): %d\n", final_sum);
编译后执行，确认结果和原来相同。
这个版本的程序，在 GeFo rce 8800GT 上执行，只需要约 8.3M cycles，比前一版程序快了 77倍！这就是透过大量 thr ead 来隐藏 la tency 所带来的效果。
不过，如果计算一下它使用的内存带宽，就会发现其实仍不是很理想，大约只有 723MB/s 而已。这和 GeFo rce 8800GT 所具有的内存带宽是很大的差距。为什么会这样呢？
内存的存取模式
显卡上的内存是 DRAM，因此最有效率的存取方式，是以连续的方式存取。前面的程序，虽然看起来是连续存取内存位置（每个 thre ad 对一块连续的数字计算平方和），但是我们要考
虑到实际上 t hread 的执行方式。前面提过，当一个 thre ad 在等待内存的数据时，GPU 会切换到下一个 t hread。也就是说，实际上执行的顺序是类似
thread 0 -> thread 1 -> thread 2 -> ...
因此，在同一个 thre ad 中连续存取内存，在实际执行时反而不是连续了。要让实际执行结果是连续的存取，我们应该要让 t hread 0 读取第一个数字，thread 1 读取第二个数字…依此类推。所以，我们可以把 kerne l 程序改成如下：
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
const int tid = threadIdx.x;
int sum = 0;
int i;
clock_t start;
if(tid == 0) start = clock();
for(i = tid; i < DATA_SIZE; i += THREAD_NUM) {
sum += num[i] * num[i];
}
result[tid] = sum;
if(tid == 0) *time = clock() - start;
}
编译后执行，确认结果相同。
仅仅是这样简单的修改，实际执行的效率就有很大的差别。在 GeFo rce 8800GT 上，上面的程序执行需要的频率是 2.6M cycles，又比前一版程序快了三倍。不过，这样仍只有 2.3GB/s 的带宽而已。
这是因为我们使用的 thr ead 数目还是不够多的原因。理论上 256 个 t hreads 最多只能隐藏256 cycles 的 la tency。但是 GPU 存取 gl obal memory 时的 laten cy 可能高达 500 cycles 以上。如果增加 thre ad 数目，就可以看到更好的效率。例如，可以把 THREAD_NUM 改成 512。在 GeFo rce 8800GT 上，这可以让执行花费的时间减少到 1.95M cycles。有些改进，但是仍不够大。不幸的是，目前 GeFo rce 8800GT 一个 b lock 最多只能有 512 个 threa ds，所以不能再增加了，而且，如果 t hread 数目增加太多，那么在 CPU 端要做的最后加总工作也会变多。更多的并行化
前面提到了 blo ck。在之前介绍呼叫 CUDA 函式时，也有提到 "b lock 数目" 这个参数。到目前为止，我们都只使用一个 bl ock。究竟 bl ock 是什么呢？
在 CUDA 中，thread 是可以分组的，也就是 bl ock。一个 bl ock 中的 th read，具有一个共享的 shared memory，也可以进行同步工作。不同 b lock 之间的 thre ad 则不行。在我们的程序中，其实不太需要进行 threa d 的同步动作，因此我们可以使用多个 bl ock 来进一步增加thread 的数目。
首先，在 #defin e DATA_SIZE 的地方，改成如下：
#define DATA_SIZE 1048576
#define BLOCK_NUM 32
#define THREAD_NUM 256
这表示我们会建立 32个 bl ocks，每个 bl ocks 有 256 个 t hreads，总共有 32*256 = 8192 个threads。
接着，我们把 ker nel 部份改成：
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int sum = 0;
int i;
if(tid == 0) time[bid] = clock();
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
sum += num[i] * num[i];
}
result[bid * THREAD_NUM + tid] = sum;
if(tid == 0) time[bid + BLOCK_NUM] = clock();
}
blockIdx.x 和 thre adIdx.x 一样是 CUDA 内建的变量，它表示的是目前的 bl ock 编号。另外，注意到我们把计算时间的方式改成每个 bl ock 都会记录开始时间及结束时间。
main 函式部份，修改成：
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result,
sizeof(int) * THREAD_NUM * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<<
int sum[THREAD_NUM * BLOCK_NUM];
clock_t time_used[BLOCK_NUM * 2];
cudaMemcpy(&sum, result, sizeof(int) * THREAD_NUM * BLOCK_NUM, cudaMemcpyDeviceToHost);
cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2,
cudaMemcpyDeviceToHost);
cudaFree(gpudata);
cudaFree(result);
cudaFree(time);
int final_sum = 0;
for(int i = 0; i < THREAD_NUM * BLOCK_NUM; i++) {
final_sum += sum[i];
}
clock_t min_start, max_end;
min_start = time_used[0];
max_end = time_used[BLOCK_NUM];
for(int i = 1; i < BLOCK_NUM; i++) {
if(min_start > time_used[i])
min_start = time_used[i];
if(max_end < time_used[i + BLOCK_NUM])
max_end = time_used[i + BLOCK_NUM];
}
printf("sum: %d time: %d\n", final_sum, max_end - min_start);
基本上我们只是把 result 的大小变大，并修改计算时间的方式，把每个 blo ck 最早的开始时间，和最晚的结束时间相减，取得总运行时间。
这个版本的程序，执行的时间减少很多，在 GeFo rce 8800GT 上只需要约 150K cycles，相当于 40GB/s 左右的带宽。不过，它在 CPU 上执行的部份，需要的时间加长了（因为 CPU 现在需要加总 8192 个数字）。为了避免这个问题，我们可以让每个 bl ock 把自己的每个 t hread 的计算结果进行加总。
Thread 的同步
前面提过，一个 bl ock 内的 thre ad 可以有共享的内存，也可以进行同步。我们可以利用这一点，让每个 b lock 内的所有 thre ad 把自己计算的结果加总起来。把 ker nel 改成如下：
__global__ static void sumOfSquares(int *num, int* result,
clock_t* time)
{
extern __shared__ int shared[];
const int tid = threadIdx.x;
const int bid = blockIdx.x;
int i;
if(tid == 0) time[bid] = clock();
shared[tid] = 0;
for(i = bid * THREAD_NUM + tid; i < DATA_SIZE;
i += BLOCK_NUM * THREAD_NUM) {
shared[tid] += num[i] * num[i];
}
__syncthreads();
if(tid == 0) {
for(i = 1; i < THREAD_NUM; i++) {
shared[0] += shared[i];
}
result[bid] = shared[0];
}
if(tid == 0) time[bid + BLOCK_NUM] = clock();
}
利用 __share d__ 声明的变量表示这是 shared m emory，是一个 blo ck 中每个 thread 都共享的内存。它会使用在 GPU 上的内存，所以存取的速度相当快，不需要担心 la tency 的问题。__syncthreads() 是一个 CUDA 的内部函数，表示 bl ock 中所有的 thre ad 都要同步到这个点，才能继续执行。在我们的例子中，由于之后要把所有 t hread 计算的结果进行加总，所以我们需要确定每个 thre ad 都已经把结果写到 shared[tid] 里面了。
接下来，把 ma in 函式的一部份改成：
int* gpudata, *result;
clock_t* time;
cudaMalloc((void**) &gpudata, sizeof(int) * DATA_SIZE);
cudaMalloc((void**) &result, sizeof(int) * BLOCK_NUM);
cudaMalloc((void**) &time, sizeof(clock_t) * BLOCK_NUM * 2);
cudaMemcpy(gpudata, data, sizeof(int) * DATA_SIZE,
cudaMemcpyHostToDevice);
sumOfSquares<< THREAD_NUM * sizeof(int)>>>(gpudata, result, time); int sum[BLOCK_NUM]; clock_t time_used[BLOCK_NUM * 2]; cudaMemcpy(&sum, result, sizeof(int) * BLOCK_NUM, cudaMemcpyDeviceToHost); cudaMemcpy(&time_used, time, sizeof(clock_t) * BLOCK_NUM * 2, cudaMemcpyDeviceToHost); cudaFree(gpudata); cudaFree(result); cudaFree(time); int final_sum = 0; for(int i = 0; i < BLOCK_NUM; i++) { final_sum += sum[i]; } 可以注意到，现在 CPU 只需要加总 BLOCK_NUM 也就是 32个数字就可以了。 不过，这个程序由于在 GPU 上多做了一些动作，所以它的效率会比较差一些。在 GeFo rce 8800GT 上，它需要约 164K cycles。 当然，效率会变差的一个原因是，在这一版的程序中，最后加总的工作，只由每个 blo ck 的thread 0来进行，但这并不是最有效率的方法。理论上，把 256 个数字加总的动作，是可以并行化的。最常见的方法，是透过树状的加法： 把 ker nel 改成如下： __global__ static void sumOfSquares(int *num, int* result, clock_t* time) { extern __shared__ int shared[]; const int tid = threadIdx.x; const int bid = blockIdx.x; int i; int offset = 1, mask = 1; if(tid == 0) time[bid] = clock(); shared[tid] = 0; for(i = bid * THREAD_NUM + tid; i < DATA_SIZE; i += BLOCK_NUM * THREAD_NUM) { shared[tid] += num[i] * num[i]; } __syncthreads(); while(offset < THREAD_NUM) { if((tid & mask) == 0) { shared[tid] += shared[tid + offset]; } offset += offset; mask = offset + mask; __syncthreads(); } if(tid == 0) { result[bid] = shared[0]; time[bid + BLOCK_NUM] = clock(); } } 后面的 wh ile 循环就是进行树状加法。main 函式则不需要修改。 这一版的程序，在 GeFo rce 8800GT 上执行需要的时间，大约是 140K cycles（相当于约43GB/s），比完全不在 GPU 上进行加总的版本还快！这是因为，在完全不在 GPU上进行 加总的版本，写入到 g lobal memory 的数据数量很大（8192 个数字），也对效率会有影响。所以，这一版程序不但在 CPU 上的运算需求降低，在 GPU 上也能跑的更快。 进一步改善 上一个版本的树状加法是一般的写法，但是它在 GPU 上执行的时候，会有 share memory 的bank conflict 的问题（详情在后面介绍 GPU 架构时会提到）。采用下面的方法，可以避免这个问题： offset = THREAD_NUM / 2; while(offset > 0) { if(tid < offset) { shared[tid] += shared[tid + offset]; } offset >>= 1; __syncthreads(); } 这样同时也省去了 ma sk 变数。因此，这个版本的执行的效率就可以再提高一些。在 GeForce 8800GT 上，这个版本执行的时间是约 137K cycles。当然，这时差别已经很小了。如果还要再提高效率，可以把树状加法整个展开： if(tid < 128) { shared[tid] += shared[tid + 128]; } __syncthreads(); if(tid < 64) { shared[tid] += shared[tid + 64]; } __syncthreads(); if(tid < 32) { shared[tid] += shared[tid + 32]; } __syncthreads(); if(tid < 16) { shared[tid] += shared[tid + 16]; } __syncthreads(); if(tid < 8) { shared[tid] += shared[tid + 8]; } __syncthreads(); if(tid < 4) { shared[tid] += shared[tid + 4]; } __syncthreads(); if(tid < 2) { shared[tid] += shared[tid + 2]; } __syncthreads(); if(tid < 1) { shared[tid] += shared[tid + 1]; } __syncthreads(); 当然这只适用于 THREAD_NUM 是 256 的情形。这样可以再省下约 1000 cycl es 左右（约44GB/s）。最后完整的程序文件可以从这里下载。 第二个 CUDA程序 前面介绍的计算平方和的程序，似乎没有什么实用价值。所以我们的第二个 CUDA 程序，要做一个确实有（某些）实用价值的程序，也就是进行矩阵乘法。而且，这次我们会使用浮点数。 虽然矩阵乘法有点老套，不过因为它相当简单，而且也可以用来介绍一些有关 CUDA 的有趣性质。 矩阵乘法 为了单纯起见，我们这里以方形的矩阵为例子。基本上，假设有两个矩阵 A 和B，则计算 AB = C 的方法如下： for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { C[i][j] = 0; for(k = 0; k < n; k++) { C[i][j] += A[i][k] * B[k][j]; } } } 一开始，我们先准备好产生数据、设定 CUDA 等等的工作： int main() { float *a, *b, *c, *d; int n = 1000; if(!InitCUDA()) return 0; a = (float*) malloc(sizeof(float) * n * n); b = (float*) malloc(sizeof(float) * n * n); c = (float*) malloc(sizeof(float) * n * n); d = (float*) malloc(sizeof(float) * n * n); srand(0); matgen(a, n, n); matgen(b, n, n); clock_t time = matmultCUDA(a, n, b, n, c, n, n); matmult(a, n, b, n, d, n, n); compare_mat(c, n, d, n, n); double sec = (double) time / CLOCKS_PER_SEC; printf("Time used: %.2f (%.2lf GFLOPS)\n", sec, 2.0 * n * n * n / (sec * 1E9)); return 0; } InitCUDA 函式和第一个 CUDA 程序一样，可以直接参考前面的文章。以下是上面用到的一些其它的函式： 产生矩阵： void matgen(float* a, int lda, int n) { int i, j; for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { a[i * lda + j] = (float) rand() / RAND_MAX + (float) rand() / (RAND_MAX * RAND_MAX); } } } 这个函式只是利用随机数生成器把矩阵填满 0 ~ 1 之间的数字。特别注意到因为 C 语言中无法声明变动大小的二维矩阵，所以我们使用i * lda + j 的方式。 进行矩阵乘法： void matmult(const float* a, int lda, const float* b, int ldb, float* c, int ldc, int n) { int i, j, k; for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { double t = 0; for(k = 0; k < n; k++) { t += a[i * lda + k] * b[k * ldb + j]; } c[i * ldc + j] = t; } } } 这是以 CPU 进行矩阵乘法、用来进行验证答案正确与否的程序。特别注意到它用 do uble 来储存暂时的计算结果，以提高精确度。 验证结果： void compare_mat(const float* a, int lda, const float* b, int ldb, int n) { float max_err = 0; float average_err = 0; int i, j; for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { if(b[i * ldb + j] != 0) { float err = fabs((a[i * lda + j] - b[i * ldb + j]) / b[i * ldb + j]); if(max_err < err) max_err = err; average_err += err; } } } printf("Max error: %g Average error: %g\n", max_err, average_err / (n * n)); } 这个函式计算两个矩阵的最大相对误差和平均相对误差，并把结果印出来。 最后是 CUDA 的矩阵乘法的部份： #define NUM_THREADS 256 clock_t matmultCUDA(const float* a, int lda, const float* b, int ldb, float* c, int ldc, int n) { float *ac, *bc, *cc; clock_t start, end; start = clock(); cudaMalloc((void**) &ac, sizeof(float) * n * n); cudaMalloc((void**) &bc, sizeof(float) * n * n); cudaMalloc((void**) &cc, sizeof(float) * n * n); cudaMemcpy2D(ac, sizeof(float) * n, a, sizeof(float) * lda, sizeof(float) * n, n, cudaMemcpyHostToDevice); cudaMemcpy2D(bc, sizeof(float) * n, b, sizeof(float) * ldb, sizeof(float) * n, n, cudaMemcpyHostToDevice); int blocks = (n + NUM_THREADS - 1) / NUM_THREADS; matMultCUDA<< (ac, n, bc, n, cc, n, n); cudaMemcpy2D(c, sizeof(float) * ldc, cc, sizeof(float) * n, sizeof(float) * n, n, cudaMemcpyDeviceToHost); cudaFree(ac); cudaFree(bc); cudaFree(cc); end = clock(); return end - start; } 这个函式相当单纯，就是在显卡内存中配置存放矩阵的内存，然后把主内存中的矩阵数据复制到显卡内存上。不过，因为我们的矩阵乘法函式可以指定 pit ch（即 lda、ldb、和 ldc），所以如果用一般的 cu daMemcpy 函式来复制内存的话，会需要每个 row 都分开复制，那会需要呼叫很多次 cudaM emcpy 函式，会使效率变得很差。因此，在这里我们用了一个新的cudaMemcpy2D 函式，它是用来复制二维数组，可以指定数组的 pi tch。这样就可以透过一次函数调用就可以了。 进行计算的 kerne l 如下： __global__ static void matMultCUDA(const float* a, size_t lda, const float* b, size_t ldb, float* c, size_t ldc, int n) { const int tid = threadIdx.x; const int bid = blockIdx.x; const int idx = bid * blockDim.x + tid; const int row = idx / n; const int column = idx % n; int i; if(row < n && column < n) { float t = 0; for(i = 0; i < n; i++) { t += a[row * lda + i] * b[i * ldb + column]; } c[row * ldc + column] = t; } } 这个函式一开始先从 bi d 和 ti d 计算出这个 t hread 应该计算的 row 和 co lumn，在判断 row 和 colum n 在范围内之后，就直接进行计算，并把结果写到 c 矩阵中，是非常单纯的函式。在 GeFo rce 8800GT 上实际执行的结果如下： Max error: 2.01484e-006 Average error: 3.36637e-007 Time used: 1.1560 (1.73 GFLOPS) 可以看到两个问题： 1.很明显的，执行效率相当低落。 2.最大相对误差偏高。理想上应该要低于 1e-6。 计算结果的误差偏高的原因是，在 CPU 上进行计算时，我们使用 do uble（即 64 bits 浮点数）来累进计算过程，而在 GPU 上则只能用 float（32 bits 浮点数）。在累加大量数字的时候，由于累加结果很快会变大，因此后面的数字很容易被舍去过多的位数。 由于 CUDA 的浮点数运算，在进行加、减、乘法时是符合 IEEE 754 规定的精确度的，因此，我们可以利用Kahan's Summation Formula 来提高精确度。把程序改成： if(row < n && column < n) { float t = 0; float y = 0; for(i = 0; i < n; i++) { float r; y -= a[row * lda + i] * b[i * ldb + column]; r = t - y; y = (r - t) + y; t = r; } } 修改后的程序的执行结果是： Max error: 1.19209e-007 Average error: 4.22751e-008 Time used: 1.1560 (1.73 GFLOPS) 可以看到相对误差有很大的改善，效率则没什么变化。 高性能计算实验报告 学生姓名：X X 学号：XXXXXXXXXX 班号：116122 指导教师：郭明强 中国地质大学（武汉）信息工程学院 第一题 1.编写console程序 2.由下图看出，电脑是双核CPU 3.多线程程序，利用windowsAPI函数创建线程 代码 #include"stdafx.h" #include } 第二题多线程实现计算e和π的乘积 代码 #include"stdafx.h" #include"windows.h" #define num_steps 2000000 #include ANSYS高性能并行计算 作者：安世亚太雷先华 高性能并行计算主要概念 ·高性能并行计算机分类 并行计算机主要可以分为如下四类：对称多处理共享存储并行机（SMP，Symmetric Multi-Processor）、分布式共享存储多处理机（DSM，Distributied Shared Memory）、大规模并行处理机（MPP，Massively Parallel Processor）和计算机集群系统（Cluster）。 这四类并行计算机也正好反映了高性能计算机系统的发展历程，前三类系统由于或多或少需要在CPU、内存、封装、互联、操作系统等方面进行定制，因而成本非常昂贵。最后一类，即计算机集群系统，由于几乎全采用商业化的非定制系统，具有极高的性能价格比，因而成为现代高性能并行计算的主流系统。它通过各种互联技术将多个计算机系统连接在一起，利用所有被连接系统的综合计算能力来处理大型计算问题，所以又通常被称为高性能计算集群。高性能并行计算的基本原理就是将问题分为若干部分，而相连的每台计算机（称为节点）均可同时参与问题的解决，从而显著缩短解决整个问题所需的计算时间。 ·集群互联网络 计算机集群系统的互联网络大体上经历了从Ethernet到Giganet、Myrinet、Infiniband、SCI、Quadrics(Q-net)等发展历程，在“延时”和“带宽”两个最主要指标上有了非常大的改善，下表即是常用的互联方式： ANSYS主要求解器的高性能并行计算特性 ANSYS系列CAE软件体系以功能齐全、多物理场耦合求解、以及协同仿真而著称于世。其核心是一系列面向各个方向应用的高级求解器，并行计算也主要是针对这些求解器而言。 ANSYS的主要求解器包括： Mechanical：隐式有限元方法结构力学求解器； CFX ：全隐式耦合多重网格计算流体力学求解器； AUTODYN：显式有限元混合方法流固耦合高度非线性动力学求解器； LS-DYNA：显式有限元方法非线性结构动力学求解器； FEKO：有限元法、矩量法、高频近似方法相互混合的计算电磁学求解器； ·高性能并行计算的典型应用 现代CAE计算的发展方向主要有两个：系统级多体耦合计算和多物理场耦合计算，前者摒弃了以往只注重零部件级CAE仿真的传统，将整个对象的完整系统（如整机、整车）一次性纳入计算范畴；后者在以往只注重单一物理场分析（如结构力学、流体力学）的基础上，将影响系统性能的所有物理因素一次性纳入计算范畴，考虑各物理因素综合起来对分析对象的影响。因此，可以说，高性能并行计算也是CAE的发展方向，因为它是大规模CAE 应用的基石。例如，在航空航天领域，需要高性能并行计算的典型CAE应用有： –飞机/火箭/导弹等大型对象整体结构静力、动力响应、碰撞、安全性分析，整体外流场分析，多天线系统电磁兼容性及高频波段RCS分析，全模型流体－结构－电磁耦合分析；–航空发动机多级转子/静子联合瞬态流动分析，流体－结构－热耦合分析； –大型运载火箭/导弹发射过程及弹道分析…… · ANSYS求解器对高性能并行计算的支持 作为大型商用CAE软件的领头雁，ANSYS在对高性能并行计算的支持方面也走在所有CAE软件的前列，其各个求解器对高性能并行系统的支持可用下表描述： 高性能计算和网格技术 实验报告 实验题目OpenMP和MPI编程姓名 学号 专业计算机系统结构 指导教师 助教 所在学院计算机科学与工程学院论文提交日期 一、实验目的 本实验的目的是通过练习掌握OpenMP 和MPI 并行编程的知识和技巧。 1、熟悉OpenMP 和MPI 编程环境和工具的使用； 2、掌握并行程序编写的基本步骤； 3、了解并行程序调试和调优的技巧。 二、实验要求 1、独立完成实验内容； 2、了解并行算法的设计基础； 3、熟悉OpenMP和MPI的编程环境以及运行环境； 4、理解不同线程数，进程数对于加速比的影响。 三、实验内容 3.1、矩阵LU分解算法的设计： 参考文档sy6.doc所使用的并行算法： 在LU分解的过程中，主要的计算是利用主行i对其余各行j，(j>i)作初等行变换，各行计算之间没有数据相关关系，因此可以对矩阵A 按行划分来实现并行计算。考虑到在计算过程中处理器之间的负载均衡，对A采用行交叉划分：设处理器个数为p，矩阵A的阶数为n，??p =，对矩阵A行交叉划分后，编号为i(i=0,1,…,p-1)的处理器存有m/ n A的第i, i+p,…, i+(m-1)p行。然后依次以第0,1,…,n-1行作为主行，将 其广播给所有处理器，各处理器利用主行对其部分行向量做行变换，这实际上是各处理器轮流选出主行并广播。若以编号为my_rank的处理器的第i行元素作为主行，并将它广播给所有处理器，则编号大于等于my_rank的处理器利用主行元素对其第i+1,…,m-1行数据做行变换，其它处理器利用主行元素对其第i,…,m-1行数据做行变换。 根据上述算法原理用代码表示如下（关键代码）： for(k = 0;k 第十章高性能计算和并行算法 §10．1 引言 计算机的运算速度在日新月异地增长，计算机的市场价格却不断地下降。 当前的计算机技术仍然远远不能满足物理问题计算的需要。 高性能计算机是一个所有最先进的硬件，软件，网络和算法的综合概念，“高性能”的标准是随着技术的发展而发展的。 高性能计算系统中最为关键的要素是单处理器的最大计算速度，存贮器访问速度和内部处理器通讯速度，多处理器系统稳定性，计算能力与价格比，以及整机性能等。 传统的计算机是冯.纽曼(Von Newmann)计算机，它是由中央处理器、内存器和输入/输出设备构成。 为了要超越这个冯.纽曼“瓶颈”，人们发展了两种计算机体系结构和相关软件技术的应用原则。一个是并行算法(parallelism)，另一个是流水线技术(pipelining)。 由于高性能计算机与当前能够应用的新计算技术相关联，因而它与并行算法和流水线技术有着密切的联系。 §10. 2并行计算机和并行算法 并行计算机是由多个处理器组成，并能够高速、高效率地进行复杂问题计算的计算机系统。 串行计算机是指只有单个处理器，顺序执行计算程序的计算机，也称为顺序计算机。 并行计算作为计算机技术，该技术的应用已经带来单机计算能力的巨大改进。 并行计算就是在同一时间内执行多条指令，或处理多个数据的计算。并行计算机是并行计算的载体。 为什么要采用并行计算呢？ z并行计算可以大大加快运算速度，即在更短的时间内完成相同的计算量，或解决原来根本不能计算的非常复杂的问题。 z提高传统的计算机的计算速度一方面受到物理上光速极限和量子效应的限制，另一方面计算机器件产品和材料的生产受到加工工艺的限制，其尺寸不可能做得无限小。因此我们只能转向并行算法。 File Name: hello.f program hello print *, 'hello series word!' !$OMP PARALLEL print *,'hello parallel world!' !$OMP END PARALLEL print *, 'hello series word!' stop end program hello ------------------------------------------- ! File Name: hp1.f program hello_parallel1 !$OMP PARALLEL print *,'hello world!' !$OMP END PARALLEL stop end program hello_parallel1 ------------------------------------------- ! File Name: hp2.f program hello_parallel_2 implicit none include 'omp_lib.h' integer :: idcpu,mcpu call OMP_SET_NUM_THREADS(3) idcpu=OMP_GET_THREAD_NUM() mcpu=OMP_GET_NUM_THREADS() print *,'------before parallel' print '(a,i4,a,i4,a)','Hello from thread',idcpu,' in',mcpu,' CPUs' print * !$OMP PARALLEL DEFAULT(NONE) PRIV ATE(IDCPU,MCPU) 高性能计算的应用 随着高性能计算技术的发展，高性能计算开始广泛应用于各个领域。在核电，气象，工业工程，水下工程，建筑，生物医学，社会科学等方面均有重要的应用。 1、核电工程领域 在核电工程领域中，核电压力容器分析，开孔安全壳环向应力分析，核电厂房抗震分析，核反应堆压力容器与管道温度分析，核电流固耦合分析，核安全防护分析等方面均需要大规模的计算[1]。通过高性能计算，对工业仿真流程进行分析，直接减少了计算时间，降低了成本，提高了企业的竞争力。 2、气象 在气象领域中，数值天气预报模式的科学研究和业务运行需要高性能计算。目前，数值预报模式的水平分辨率已达到了15~20公里，而未来的3-5年内几乎世界各国的全球数值预报模式的水平分辨率都将要提高10~20公里[2]，为适应其快速发展，气象部门需要引进和更新高性能计算机系统用以支持气象应用。 3、工业工程 对于工业和工程领域来说，使用高性能计算对于计算数学特别是用力学计算仿真手段来模拟实际产品制造、产品运行环境和工程建设环境具有不可代替的作用[3]。高性能计算降低了物理原型和实验的数量，提高了设计质量和效率，提升了企业解决复杂技术难题的手段和能力。 在石油勘测方面，由于地震波法勘测收集的数据通常都以TB计，在海洋勘测过程中的数据容量更是达到了PB级别量[4]，面对这些海量的数据，只有借助性能出色的高性能计算机系统，才能缩短时间，以实现最佳的勘测效益。 在高光谱遥感数据处理方面，高光谱遥感数据的海量特性严重制约了应用的拓展和实际工程应用效率的提高，大量数据操作和处理的复杂性决定了高光谱遥感图像处理具有很强的计算性[5]，普通计算机远远无法满足遥感数据处理的增长需求，因此高性能计算是解决海量数据处理效率低的有效方法。 在飞机设计方面，首先，飞机设计需要做大量的气动力预测工作，采用高性能计算比采用传统的风洞试验成本要低得多，而且在提升飞机性能时，常规基于雷诺平均方程的CFD技术并不能有效处理，因为它需要的计算网格约10亿量级，需求的计算能力比常规计算高出2个量级以上。其次，精确噪声预测，螺旋桨滑流研究，需要的计算能力比常规计算高出2个量级以上。而现代军用飞机对雷达散射截面计算的要求十分严格，只有基于高性能计算的电磁数值仿真技术有望解决RCS预测难题[6]。 在岩石力学课程教学方面，由于岩石力学需要将工程实例，实验模型和理论模型相结合，才能增强教学效果。尤其在涉及破裂问题上，从变形，损伤演化到最终失稳的过程对数值模拟而言需要网格重划分、单元消去与再生、节点释放和数据存储管理，串行CPU和内存无力做到精细表征[7]。因此，需要通过高性能计算，充分利用数值模拟辅助岩石力学教学 在电力系统工程方面，现代的电力系统分析需要越来越多的计算，包括仿真，优化，控制和分析。人们需要寻求新的方式追求计算更快的方法来进行效率性更高的计算及解决问题，以确保电力网格系统的安全性和可靠性[8]。一个很明显 切实实现高性能并行计算应用分析 高性能并行计算的应用软件位于高性能计算生态系统的最上层，针对不同的行业有专业的产品，针对各个领域的科学与工程计算应用，直接为用户创造价值。这些软件原来大多运行在大型主机上，是面向多个处理器、多进程、多任务的单节点软件，进程之间的通信通过大型主机操作系统的消息机制进行，消息机制的启动通过函数进行调用。 本系统中，应用软件面向教学和科研应用领域的多个方面基于多节点IA架构系统，进程或任务之间的通信，基于多节点集群的中间件提供的并行通信库MPI，物理层是基于标准互联以太网系统。并行库的启动，通过特定的程序语句进行调用。 高性能应用软件总体概括分类: √多媒体运算 主要使用整型和双精度运算。包括图形图像处理和三维图像生成的高性能计算系统，强调计算节点的多媒体计算功能。计算科学院的大气科学和流体力学应用中需要的许多模拟仿真计算都属于这类计算。 √科学计算 主要使用浮点运算功能，这也是目前高性能计算系统的最主要应用领域。比如:高分子运动分析、石油勘测分析等。计算科学院的大气科学、固体力学、分子力学、流体力学、有限元分析等的主体计算都属于这类计算，这类计算需要系统具有强大的浮点运算能力。本项目的计算属于此类应用。 √数据库应用 主要使用逻辑计算和I/O操作。包括数据库集群系统和网格数据库系统的应用。强调计算节点有很强的I/O处理能力，同时，整个高性能计算系统具有足够的外接存储空间。本系统结合此类应用，奠定未来网格计算的基础。 INTEL和宝德技术人员针对华南理工的项目特点和目标，投入极高的专注和热情，在华南理工项目前期进行了详细的测试分析，提出系统优化和移植的策略，帮助客户将微分方程数值计算并行模拟器勘测系统移植到IA平台上。 Intel还提供了系列的优化工具、编译工具、集群工具等众多高性能计算组件和虚拟技术，为IA架构、标准互联的高性能计算系统应用提供高效率的保证，成为本次HPC项目成功实施的关键。 解决方案高性能计算报告
ANSYS高性能并行计算
并行计算实验报告(高性能计算与网格技术)
高性能计算和并行算法-计算物理课件
多核并行高性能计算OpenMP第二章源程序
高性能计算的应用
切实实现高性能并行计算应用分析