10.3 典型并行编程环境
本节主要介绍数据并行SIMD编程、共享存储编程模型Pthreads和OpenMP以及消息传递编程模型MPI等。
10.3.1 数据并行SIMD编程
工业界广泛应用的单指令流多数据流（Single Instruction Multiple Data，简称SIMD）并行就是典型的数据并行技术。相比于传统的标量处理器上的单指令流单数据流（Single Instruction Single Data，简称SISD）指令，一条SIMD指令可以同时对一组数据进行相同的计算。比如将两个数组SRC0[8]和SRC1[8]中的每个对应元素求和，将结果放入数组RESULT中，对于传统的标量处理器平台，C语言实现如下：
for（i = 0; i < 8; i++）
RESULT[i] = SRC0[i] + SRC1[i];
也就是通过循环遍历需要求和的8组对应数据，对SRC0和SRC1的各对应项求和，将结果存入RESULT数组的对应项中。在龙芯处理器平台上，用机器指令（汇编代码）实现该运算的代码如下（这里假设$src0、$src1、$result分别为存储了SRC0、SRC1和RESULT数组起始地址的通用寄存器）：
li $r4, 0x0 li $r5, 0x8 1: addi.d $src0,$src0,$r4 addi.d $src1,$src1,$4 addi.d $result,$result,$4 ld.b $r6,$src0, 0x0 ld.b $r7,$src1, 0x0 add.d $r6,$r6,$r7 st.b $r6,$result, 0x0 addi.d $r4,$r4, 0x1 blt $r4,$r5, 1b
如果采用龙芯处理器的SIMD指令编写程序的话，上述两个数组的求和只需要将上述两个源操作数数组SRC0[8]和SRC1[8]一次性加载到龙芯处理器的向量寄存器（龙芯向量寄存器复用了浮点寄存器）中，然后只需要一条paddb指令就可以完成上述8个对应数组项的求和，最后只需要一条store指令就可以将结果存回RESULT[8]数组所在的内存空间中。该实现的机器指令序列如下：
vld $vr0,$src0,0 vld $vr1,$src1,0 vadd.b $vr0,$vr0,$vr1 vst $vr0,$result,0
图10.6简要示意了采用传统SISD指令和SIMD指令实现上述8个对应数组项求和的执行控制流程。
图10.6 SISD和SIMD执行控制流示意图
10.3.2 POSIX编程标准
POSIX（Portable Operating System Interface）属于早期的操作系统接口标准。Pthreads代表官方IEEE POSIX1003.1C_1995线程标准，是由IEEE标准委员会所建立的，主要包含线程管理、线程调度、同步等原语定义，体现为C语言的一套函数库。下面只简介其公共性质。
1.线程管理
线程库用于管理线程，Pthreads中基本线程管理原语如表10.1所示。其中pthread_create()在进程内生成新线程，新线程执行带有变元arg的myroutine，如果pthread_create()生成，则返回0并将新线程之ID置入thread_id，否则返回指明错误类型的错误代码；pthread_exit()结束调用线程并执行清场处理；pthread_self()返回调用线程的ID；pthread_join()等待其他线程结束。
表10.1 Pthreads中基本线程管理一览表
2.线程调度
pthread_yield()的功能是使调用者将处理器让位于其他线程；pthread_cancel()的功能是中止指定的线程。
3.线程同步
Pthreads中的同步原语见表10.2。重点讨论互斥变量mutex（Mutual Exclusion）和条件变量cond（Conditional）。前者类似于信号灯结构，后者类似于事件结构。注意，使用同步变量之前需被初始化（生成），用后应销毁。
表10.2 Pthreads中线程相互作用原语
如果mutex未被上锁，pthread_mutex_lock()将锁住mutex；如果mutex已被上锁，调用线程一直被阻塞到mutex变成有效。pthead_mutex_trylock()的功能是尝试对mutex上锁。pthread_mutex_lock()和pthead_mutex_trylock()的区别是：前者会阻塞等待mutex被解锁；后者尝试去加锁，如果不成功就返回非0，如果成功返回0，不会产生阻塞。pthead_mutex_unlock()解锁先前上锁的mutex，当mutex被解锁，它就能由别的线程获取。
pthread_cond_wait()自动阻塞等待条件满足的现行线程，并开锁mutex。pthread_cond_timedwait()与pthread_cond_wait()类似，除了当等待时间达到时限它将解除阻塞外。pthread_cond_signal()解除一个等待条件满足的已被阻塞的线程的阻塞。pthread_cond_broadcast()将所有等待条件满足的已被阻塞的线程解除阻塞。
4.示例
以下程序示例用数值积分法求π的近似值。
我们使用梯形规则来求解这个积分。其基本思想是用一系列矩形填充一条曲线下的区域，也就是求出在区间[0,1]内函数曲线4/（1+x2）下的面积，此面积就是π的近似值。为此先将区间[0,1]划分成N个等间隔的子区间，每个子区间的宽度为1.0/N；然后计算出各子区间中点处的函数值；再将各子区间面积相加就可得出π的近似值。N的值越大，π值的误差越小。图10.7为进行π值计算的C语言描述的串行代码。为简化起见，将积分中的迭代次数固定为1000000。
图10.7 利用梯形规则计算π的C语言串行代码
图10.8所示为采用Pthreads的并行代码。
图10.8 利用梯形规则计算π的Pthreads并行代码
下面举一个矩阵乘的Pthreads并行代码例子，如图10.9所示。该例子将两个n阶的方阵A和B相乘，结果存放在方阵C中。
图10.9 矩阵乘的Pthreads并行代码
10.3.3 OpenMP标准
OpenMP是由OpenMP Architecture Review Board（ARB，结构审议委员会）牵头提出的，是一种用于共享存储并行系统的编程标准。最初的OpenMP标准形成于1997年，2002年发布了OpenMP 2.0标准，2008年发布了OpenMP 3.0标准，2013年发布了OpenMP 4.0标准。实际上，OpenMP不是一种新语言，是对基本编程语言进行编译制导（Compiler Directive）扩展，支持C/C++和Fortran。由于OpenMP制导嵌入到C/C++、Fortran语言中，所以具体语言不同会有所区别，本书介绍主要参考支持C/C++的OpenMP 4.0标准。
OpenMP标准中定义了制导指令、运行库和环境变量，使得用户可以按照标准逐步将已有串行程序并行化。制导语句是对程序设计语言的扩展，提供了对并行区域、工作共享、同步构造的支持；运行库和环境变量使用户可以调整并行程序的执行环境。程序员通过在程序源代码中加入专用pragma制导语句（以“#pragmaomp”字符串开头）来指明自己的意图，支持OpenMP标准的编译器可以自动将程序进行并行化，并在必要之处加入同步互斥以及通信。当选择忽略这些pragma，或者编译器不支持OpenMP时，程序又可退化为普通程序（一般为串行），代码仍然可以正常运行，只是不能利用多线程来加速程序执行。
由于OpenMP标准具有简单、移植性好和可扩展等优点，目前已被广泛接受，主流处理器平台均支持OpenMP编译器，如Intel、AMD、IBM、龙芯等。开源编译器GCC也支持OpenMP标准。
1.OpenMP的并行执行模型
OpenMP是一个基于线程的并行编程模型，一个OpenMP进程由多个线程组成，使用fork-join并行执行模型。OpenMP程序开始于一个单独的主线程（Master Thread），主线程串行执行，遇到一个并行域（Parallel Region）开始并行执行。接下来的过程如下：
1）fork（分叉）。主线程派生出一队并行的线程，并行域的代码在主线程和派生出的线程间并行执行。
2）join（合并）。当派生线程在并行域中执行完后，它们或被阻塞或被中断，所计算的结果会被主线程收集，最后只有主线程在执行。
实际上，OpenMP的并行化都是使用嵌入到C/C++或者Fortran语言的制导语句来实现的。图10.10为OpenMP程序的并行结构。
图10.10 OpenMP程序的并行结构
2.编译制导语句
下面介绍编译制导语句的格式。参看前面的OpenMP程序并行结构的例子，在并行开始部分需要语句“#pragma omp parallel private（var1,var2）shared（var3）”。表10.3是编译制导语句的格式及解释。
表10.3 编译制导语句的解释
3.并行域结构
一个并行域就是一个能被多个线程执行的程序块，它是最基本的OpenMP并行结构。并行域的具体格式为：
#pragma omp parallel [if（scalar_expression）| num_threads（integer-expression）|default （shared|none）|private（list）|firstprivate（list）|shared（list）| copyin（list）|reduction（operator: list）| proc_bind（master|close|spread）] newline
当一个线程执行到parallel这个指令时，线程就会生成一列线程，线程号依次从0到n-1，而它自己会成为主线程（线程号为0）。当并行域开始时，程序代码就会被复制，每个线程都会执行该代码。这就意味着，到了并行域结束就会有一个栅障，且只有主线程能够通过这个栅障。
4.共享任务结构
共享任务结构将其内封闭的代码段划分给线程队列中的各线程执行。它不产生新的线程，在进入共享任务结构时不存在栅障，但是在共享任务结构结束时存在一个隐含的栅障。图10.11显示了3种典型的共享任务结构。其中：do/for将循环分布到线程列中执行，可看作是一种表达数据并行的类型；sections把任务分割成多个各个部分（section），每个线程执行一个section，可很好地表达任务并行；single由线程队列中的一个线程串行执行。
图10.11 共享任务类型
下面具体来看一下。
1）for编译制导语句。for语句（即C/C++中的for语句），表明若并行域已经初始化了，后面的循环就在线程队列中并行执行，否则就会顺序执行。语句格式如下：
#pragma omp for [private（list）|firstprivate（list）|lastprivate（list）|reduction（reduction- identifier:list）|schedule（kind[,chunk_size]）|collapse（n）|ordered|nowait] newline
其中，schedule子句描述如何在线程队列中划分循环。kind为static时，将循环划分为chunk_size大小的循环块，静态分配给线程执行，若chunk_size没有声明，则尽量将循环在线程队列中平分；kind为dynamic时，线程会动态请求循环块来执行，执行完一个循环块后就申请下一个循环块，直到所有循环都执行完，循环块的大小为chunk_size，若chunk_size没有声明，则默认的块长度为1；kind为guide时，线程会动态请求循环块来执行，循环块的大小为未调度的循环数除以线程数，但循环块大小不能小于chunk_size（除了最后一块），若chunk_size没有声明，则默认为1。
2）sections编译制导语句。该语句是非循环的共享任务结构，它表明内部的代码是被线程队列分割的。语句格式如下：
#pragmaompsections[private（list）|firstprivate（list）|lastprivate（list）|reduction（re-
duction-identifier:list）|nowait] newline
{
[#pragma omp section newline]
Structured_block
[#pragma omp section newline
Structured_block]
…
}
值得注意的是，在没有nowait子句时，sections后面有栅障。
3）single编译制导语句。该语句表明内部的代码只由一个线程执行。语句格式如下：
#pragma omp single [private（list）|firstprivate（list）| copyprivate（list）|nowait] newline
Structured_block
若没有nowait子句，线程列中没有执行single语句的线程，会一直等到代码栅障同步才会继续往下执行。
5.组合的并行共享任务结构
下面介绍两种将并行域制导和共享任务制导组合在一起的编译制导语句。
1）parallel for编译制导语句。该语句表明一个并行域包含一个单独的for语句。语句格式如下：
#pragma omp parallel for [if（scalar_expression）|num_threads（integer-expression |default （shared|none）|private（list）|firstprivate（list）|lastprivate（list）|shared（list）|copyin （list）|reduction（Structured_block:list）|proc_bind（master |close |spread）|schedule（kind[, chunk_size]）|collapse（n）|ordered] newline for_loop
该语句的子句可以是parallel和for语句的任意子句组合，除了nowait子句。
2）parallel sections编译制导语句。该语句表明一个并行域包含单独的一个sections语句。语句格式如下：
#pragmaompparallelsections
[if（scalar_expression）|num_threads（integer-expression）|default（shared|none）|private
（list）|firstprivate（list）|lastprivate（list）|shared（list）|copyin（list）|reduction（Struc-
tured_block:list）|proc_bind（master|close|spread）]
{
[#progma omp section newline]
Structured_block
[#progma omp section newline
Structured_block]
…
}
同样，该语句的子句可以是parallel和for语句的任意子句组合，除了nowait子句。
6.同步结构
OpenMP提供了多种同步结构来控制与其他线程相关的线程的执行。下面列出几种常用的同步编译制导语句。
1）master编译制导语句。该语句表明一个只能被主线程执行的域。线程队列中所有其他线程必须跳过这部分代码的执行，语句中没有栅障。语句格式如下：
#pragma omp master newline
2）critical编译制导语句。该语句表明域中的代码一次只能由一个线程执行。语句格式如下：
#pragma omp critical[name] newline
3）barrier编译指导语句。该语句同步线程队列中的所有线程。当有一个barrier语句时，线程必须要等到所有的其他线程也到达这个栅障时才能继续执行。然后所有线程并行执行栅障之后的代码。语句格式如下：
#pragma omp barrier newline
4）atomic编译制导语句。该语句表明一个特别的存储单元只能原子地更新，而不允许让多个线程同时去写。语句格式如下：
#pragma omp atomic newline
另外，还有flush、order等语句。
7.数据环境
OpenMP中提供了用来控制并行域在多线程队列中执行时的数据环境的制导语句和子句。下面选择主要的进行简介。
1）threadprivate编译制导语句。该语句表明变量是复制的，每个线程都有自己私有的备份。这条语句必须出现在变量序列定义之后。每个线程都复制这个变量块，所以一个线程的写数据对其他线程是不可见的。语句格式如下：
#pragma omp threadprivate（list）
2）数据域属性子句。OpenMP的数据域属性子句用来定义变量的范围，它包括private、firstprivate、lastprivate、shared、default、reduction和copyin等。数据域变量与编译制导语句parallel、for、sections等配合使用，可控制变量的范围。它们在并行结构执行过程中控制数据环境。例如：哪些串行部分的数据变量被传到程序的并行部分以及如何传送；哪些变量对所有的并行部分是可见的；哪些变量是线程私有的；等等。具体说明如下。
·private子句：表示它列出的变量对于每个线程是局部的，即线程私有的。其格式为：
private（list）
·shared子句：表示它列出的变量被线程队列中的所有线程共享，程序员可以使多线程对其进行读写（例如通过critical语句）。其格式为：
shared（list）
·default子句：该子句让用户可以规定在并行域的词法范围内所有变量的一个默认属性（如可以是private、shared、none）。其格式为：
default（shared|none）
·firstprivate子句：该子句包含private子句的操作，并将其列出的变量的值初始化为并行域外同名变量的值。其格式为：
firstprivate（list）
·lastprivate子句：该子句包含private子句的操作，并将值复制给并行域外的同名变量。其格式为：
lastprivate（list）
·copyin子句：该子句赋予线程中变量与主线程中threadprivate同名变量的值。其格式为：
copyin（list）
·reduction子句：该子句用来归约其列表中出现的变量。归约操作可以是加、减、乘、与（and）、或（or）、相等（eqv）、不相等（neqv）、最大（max）、最小（min）等。其格式为：
reduction（reduction-identifier:list）
图10.12给出了一个计算π的采用OpenMP并行的C语言代码示例。
图10.12 利用梯形规则计算π的OpenMP并行代码
图10.13给出了矩阵乘的OpenMP并行代码例子，将两个n阶的方阵A和B相乘，结果存放在方阵C中。
10.3.4 MPI消息传递编程接口
MPI（Message Passing Interface）定义了一组消息传递函数库的编程接口标准。1994年发布了MPI第1版MPI-1,1997年发布了扩充版MPI-2，2012年发布了MPI-3标准。有多种支持MPI标准的函数库实现，开源实现有MPICH（由Argonne National Laboratory（ANL）和Mississippi State University开发）、Open MPI和LAM/MPI（由Ohio超算中心开发）等；商业实现来自Intel、Microsoft、HP公司等。MPI编译器用于编译和链接MPI程序，支持C、C++、Fortran语言，如mpicc支持C语言、mpic++支持C++语言、mpif90支持Fortran90。MPI具有高可移植性和易用性，对运行的硬件要求简单，是目前国际上最流行的并行编程环境之一。
图10.13 矩阵乘的OpenMP并行代码
在MPI编程模型中，计算由一个或多个通过调用库函数进行消息收/发通信的进程所组成。在绝大部分MPI实现中，一组固定的进程在程序初始化时生成，在一个处理器核上通常只生成一个进程。这些进程可以执行相同或不同的程序（相应地称为单程序多数据（SPMD）或多程序多数据（MPMD）模式）。进程间的通信可以是点到点的或者集合（Collective）的。MPI只是为程序员提供了一个并行环境库，程序员用标准串行语言编写代码，并在其中调用MPI的库函数来实现消息通信，进行并行处理。
1.最基本的MPI
MPI是个复杂的系统，包括129个函数（根据1994年发布的MPI标准）。事实上，1997年修订的MPI-2标准中函数已超过200个，其中最常用的有约30个，但只需要6个最基本的函数就能编写MPI程序求解许多问题，如表10.4所示。
表10.4 MPI的6个最基本的函数
图10.14显示了这6个基本函数的功能及参数情况。其中，标号IN表明函数使用但是不能修改参数；OUT表明函数不使用但是可以修改参数；INOUT表明函数既可以使用也可以修改参数。
图10.14 基本MPI函数说明
图10.15是一个简单C语言的MPI程序的例子，其中MPI_COMM_WORLD是一个缺省的进程组，它指明所有的进程都参与计算。
图10.15 一个简单的MPI程序的例子
2.集体通信
并行程序中经常需要一些进程组间的集体通信（Collective Communication），包括：①栅障（MPI_Barrier），同步所有进程；②广播（MPI_Bcast），从一个进程发送一条数据给所有进程；③收集（MPI_Gather），从所有进程收集数据到一个进程；④散播（MPI_Scatter），从一个进程散发多条数据给所有进程；⑤归约（MPI_Reduce、MPI_Allreduce），包括求和、求积等。这些函数的功能及参数描述参见MPI 3.0标准。不同于点对点通信，所有的进程都必须执行集体通信函数。集体通信函数不需要同步操作就能使所有进程同步，因此可能造成死锁。这意味着集体通信函数必须在所有进程上以相同的顺序执行。
3.通信域
通信域（Communicator）提供了MPI中独立的安全的消息传递。MPI通信域包含进程组（Process Group）和通信上下文（Context）。其中进程组是参加通信的一个有限并有序的进程集合，如果一共有N个进程参加通信，则进程的编号从0到N-1。通信上下文提供一个相对独立的通信区域，消息总是在其被发送的上下文内被接收，不同上下文的消息互不干涉。通信上下文可以将不同的通信区别开来。MPI提供了一个预定义的通信域MPI_COMM_WORLD，MPI初始化后就会产生，它包含了初始化时可得的全部进程，进程由它们在MPI_COMM_WORLD组中的进程号所标识。
用户可以在原有通信域的基础上定义新的通信域。MPI提供的通信域函数概括：①MPI_Comm_dup，它生成一个新的通信域，具有相同的进程组和新的上下文，这可确保不同目的通信不会混淆；②MPI_Comm_split，它生成一个新的通信域，但只是给定进程组的子集，这些进程可相互通信而不必担心与其他并发计算相冲突；③MPI_Intercomm_create，它构造一个进程组之间的通信域，该通信域链接两组内的进程；④MPI_Comm_free，它用来释放上述三个函数所生成的通信域。
4.MPI点对点通信
点到点通信（Point-to-Point Communication）是MPI中较复杂的部分，其数据传送有阻塞（Blocking）和非阻塞（Non_Blocking）两种机制。在阻塞方式中，它必须等到消息从本地送出之后才可以执行后续的语句，保证了缓冲区等资源可再用；对于非阻塞方式，它无须等到消息从本地送出就可执行后续的语句，从而允许通信和计算的重叠，但非阻塞调用的返回并不保证资源的可再用性。
阻塞和非阻塞有四种通信模式：①标准模式，包括阻塞发送MPI_Send、阻塞接收MPI_Recv、非阻塞发送MPI_Isend和非阻塞接收MPI_Irecv；②缓冲模式，包括阻塞缓冲发送MPI_Bsend和非阻塞缓冲发送MPI_Ibsend；③同步模式，包括阻塞同步发送MPI_Ssend非阻塞同步发送MPI_Issend；④就绪模式，包括阻塞就绪发送MPI_Rsend和非阻塞就绪发送MPI_Irsend。在标准通信模式中，MPI根据当前的状况选取其他三种模式或用户定义的其他模式；缓冲模式在相匹配的接收未开始的情况下，将送出的消息放在缓冲区内，这样发送者可以很快地继续计算，然后由系统处理放在缓冲区中的消息，但这占用内存且多了一次内存拷贝；在同步模式中，只有相匹配的接收操作开始后，发送才能返回；在就绪模式下，只有相匹配的接收操作启动后，发送操作才能开始。
在点到点通信中，发送和接收语句必须是匹配的。为了区分不同进程或同一进程发送来的不同消息，在这些语句中采用了通信域Comm和标志位tag来实现成对语句的匹配。
上述函数中，关于MPI_Send和MPI_Recv的功能和定义可以参考图10.14，其他函数的描述可参考MPI 3.0标准。
图10.16是计算π的C语言MPI程序的例子。
图10.16 利用梯形规则计算π的MPI并行代码
图10.17是进行矩阵乘的C语言MPI程序的例子。该例子将两个n阶的方阵A和B相乘，结果存放在方阵C中，A、B、C都在节点0上，采用主从进程的计算方法，主进程将数据发送给从进程，从进程将计算结果返回给主进程。
图10.17 矩阵乘的MPI并行代码
图10.17 （续）