有限元并行计算简介
1. 并行计算基础
1.1 并行计算概念
图2 并行计算把待解问题分解为互相独立的子问题，然后把在不同的计算硬件上计算这些子问题。要实现并行计算，问题的可分解是关键，如何分解是学问。
1.2 并行计算的硬件基础
1.2.1 PC
现在只怕再廉价的CPU也具有多核, 包含GPU的也常见。下图的主卡包含4核的CPU和GPU
图3 Core-i7-4770K
图4 CPU-GPU构造上图图示了这种主板的构造。虽然同为可编程计算硬件，与CPU(central processing unit)相比，GPU(Graphics Processor Unit)拥有大量绿色ALU但是较少的缓存和逻辑控制，适宜于小粒度的大量并行计算(一个ALU 缓存和逻辑控制等组成一个相对独立的计算单元，称为一个核core)。相应地。 CPU的核数要少得多，但是逻辑控制能力强，适宜于完成更灵活多变的计算任务，是计算机的大脑。很多人称多核CPU为multicore，称GPU为manycore processor。
当然，实际上硬件的构成是多种多样的。如图5
图5 各种构成方式在windows command环境下可以通过键入msinfo32确认你的CPU配置，或键入wmic cpu get caption, deviceid, name, numberofcores, maxclockspeed, status。如要确认GPU配置，键入dxdiag命令，在显示器窗口下即可确认，如
键入wmic path win32_VideoController可以得到更加详细的信息。
在linux下可以用 sudo lshw -c cpu或cat /proc/cpuinfo来确认CPU配置， 用lspci | grep VGA来确认GPU配置。
1.2.2 超级并行计算机
超级并行计算机当然不是只由一块主板组成。在大多数情况下，几块Socket装入一个相对独立的方盒子中构成一个node, 然后把成百上千的这样的方盒子插入数十上百的大盒子中，用网线相连构成一个拥有成千上万CPU的超级计算机。当然，这些方盒子也可以理解为在万维网中用网线连接起来的任意一台计算机，由此构成一台地球级巨型计算机。刘慈欣在他的《球状闪电》中就是用的这台计算机^O^。
图
如此，我们可以根据使用内存的共有/非共有把计算机分为两类: shared memory computer 和distributed memory computer。hybrid parallel computer兼具两种的特性。
图6 计算机硬件构造分类GPU是典型的shared memory computer, 多核的CPU也是shared memory computer。shared memory processing（SMP）中的任意的计算核可以通过内存的全局地址访问所有内存。而在distributed memory computer中内存分离设置，一个计算核（core）可以直接访问其本地内存，而访问远端内存则需要使用约定的语言通过网络间通信来实现。
下图是各国的超算力统计。中国的硬件是很强了。具体的应用水平怎么样？
图7 国家算力（From https://www.top500.org/）1.3 并行计算的软件基础
假设我们要实现标准库blas内的一个函数SAXPY
void saxpy(int n, float a, float *x, float *y)
{
for(int i=0; i<n; i )
y[i] = a*x[i]
}
void main{
...
saxpy(n, 3.0, x, y)
...
}这个简单的实现只能在单核上运行。 下面我们来看看并行计算是怎么实现的。
1.3.1 用Distributed memory computer实现并行计算
分散内存计算机间的标准通信规格为MPI（Message Passing Interface）。现在最为常见的实装有OpenMPI， mpich， Intel MPI， Microsoft MPI， 其中Intel MPI和Microsoft MPI都是基于开源的mpich。Microsoft MPI应该是专门为windows实现的。如果你的系统中不包含上述实装，你需要先下载安装一个。
void saxpy(int n, float a, float *x, float *y)
{
for(int i=0; i<n; i )
y[i] = a*x[i]
}
...
#include "mpi.h"
void main(int argc, char **argv)
{
float x[100000],y[100000];
float px[25000],py[25000]
int myrank, nprocs, i;
MPI_Init(&argc,&argv);
MPI_Comm_size(MPI_COMM_WORLD,&nprocs);
MPI_Comm_rank(MPI COMM WORLD,&myrank);
/* 在CPU 0处给数组x，y赋值 */
if (myrank == 0) {
for (i = 0; i < 100000; i )
x[i] = rand();
for (i = 0; i < 100000; i )
y[i] = rand();
}
/* 把CPU 0处的x，y值按顺序传输到CPU0，CPU1，CPU2和CPU3处的px，py数组中 */
MPI_Scatter(x, 25000, MPI_REAL, px, 25000, MPI_REAL, 0, MPI_COMM_WORLD);
MPI_Scatter(y, 25000, MPI_REAL, py, 25000, MPI_REAL, 0, MPI_COMM_WORLD);
saxpy(25000, 3.0, px, py)
/* 把CPU0，CPU1，CPU2和CPU3处的py值合并到CPU0处的y*/
MPI_Gather(py, 25000, MPI_REAL, y, 25000, MPI_REAL, 0, MPI_COMM_WORLD);
...
MPI_Finalize();
}（上面的程序写的很笨拙，仅由于示例MPI的用法）
上面的程序把两个维数为100000的数组的加法分解到4个CPU（CPU0～CPU3）处计算，然后再把结果合并到CPU0。在上面的程序中，nprocs为CPU个数（在此固定为4），myrank是各个CPU的ID或编号（在此为0~3)。使用MPI，往往需要在各CPU间进行通信（如上述MPI_Scatter，MPI_Gather）
1.3.2 多核CPU的并行编程
常用的多核CPU的并行编程方法有Pthreads和OpenMP。绝大多数编译器如gcc, intel icc都支持这两种方法。Pthreads提供了底层的API来控制各个计算核。用OpenMP来实现多核运行则相对简单，但是弱于控制计算的细节。一般而言，优先使用OpenMP。
void saxpy(int n, float a, float *x, float *restrict y)
{
#progma omp parallel for
for(int i=0; i<n; i )
y[i] = a*x[i]
}
void main{
...
saxpy(n, 3.0, x, y)
...
}如上所示，通过增加一指示行#progma omp parallel for即可实现并行。
实现多核CPU的并行方法还有一些，如std::thread, 最新的MPI也增加了多核并行计算的规格。其实最简单的实现多核运算的方法是调用自带多核运算的库函数，如blas，lapack库都是可以多核运算的
1.3.3 GPU的并行编程
常用的GPU的并行编程方法有CUDA和OpenAcc。CUDA编译器由NAVIA提供，是c语言的扩张，它提供了提供了底层的API来控制GPU运算，这一点有Pthreads相似。相应地，OpenAcc与OpenMP相似，易用，通用性较强（也可用于其他厂家的GPU）。目前gcc和PGI编译器支持OpenAcc。
OpenAcc
void saxpy(int n, float a, float *x, float *restrict y)
{
#progma acc parallel copy(y[:n]) copyin(x[:n])
for(int i=0; i<n; i )
y[i] = a*x[i]
}
void main{
...
saxpy(n, 3.0, x, y)
...
}CUDA
__global__ void saxpy(int n, float a, float *x, float *y)
{
int i = threadldx.x blockDim.x*blockLdx;
if(i<n)
y[i] = a*x[i]
}
...
cudaMemcpy(d_x, x, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_y, y, size, cudaMemcpyHostToDevice);
saxpy<<< N/128, 128 >>>(N, 3.0, d_x, d_y);
cudaMemcpy(y, d_y, size, cudaMemcpyDeviceToHost);
...GPU的编程方法还有一些，比如OpenCL, C AMP, 最新的OpenMP 4.0也增加了GPU对应规格。最后，调用自带GPU运算的库函数，如cuBLAS，cuSOLVER是实现GPU并行运算的最简单的方法。
OpenMP 4.0
void saxpy(int n, float a, float * restrict x, float * restrict y)
{
#pragma omp target teams distribute
parallel for is_device_ptr(x,y)
for(int i=0; i<n; i )
y[i] = a*x[i];
}1.3.4 混合并行编程
如果同时用到OpenMP， OpenAcc， MPI，就是混合编程了。如果要搞大型计算，这个是必须的手艺。
图8 混合并行2. 有限元并列计算
有限元并行计算一般采用SIMD(Single Instruction Multiple Data)方式进行。在这里，待解析空间先被分解投入到分布式计算机中，然后在每台分布式计算机中的各有限元单元再次被分解为数个小组分别由不同的计算核来担当计算。由于有限元计算比较繁杂，上述任务一般不适合于由GPU来完成。GPU一般用于大规模的代数运算，就是图9中显示的最下一层。
图9 有限元并行计算的常见方式如上所示，大型有限元并行计算的第一步是把有限元网格分割配分到分布式CPU中。如下图10，我们把16个单元分割为红，蓝，黄3个部分。由于处在各分割区域边界的节点为至少两个区域共有，为了计算方便，很多软件设计将不属于本分割区但是与边界节点相连的单元也记录在本区域的内存中（参见图10）。不属于本区域但被本区域记录的节点称为ghost或halo节点。在有限元计算中，各分割区域间需要进行频繁的通信以保持这些ghost节点在各个区域同步。由于通过网线的通信速度很慢，如何控制上述通信是有限元并行计算的最大的难点。其注意点如下：
各分割区的计算量应大致相同，避免出现通信待机状态。
通信量最小，也就是说ghost节点数要尽量少。
我们一般使用Graph数据处理软件，如metis，zoltan，来优化区域分割。
避免碎片化的通信。力争把大量的通信数据集中在连续的内存中排成一列一次完成。
图10 4*4单元分解到3个Processor中多核或GPU的并行计算不存在通信问题，但要注意同期处理和数据冲突。
参考文献
1 Introduction to Parallel Computing
2. Introduction to parallel computers and parallel programming.
