并行计算及并行算法

一、并行计算

　　简单地说，并行计算就是在并行计算机上所做的计算。从普通意义上讲，它和常说的高性能计算、超级计算等是同义词。并行计算的初衷是为了努力仿真自然世界中一个序列中含有众多同时发生的、复杂且相关事件的事务状态。
　　为了利用并行计算求解一个计算问题，通常基于以下考虑：1.将计算任务分解成多个子任务，有助于同时解决；2.在同一时间，由不同的执行部件可同时执行多个子任务；3.多计算资源下解决问题的耗时要少于单个计算资源下的耗时。
　　并行计算可分为：1.计算密集型：如大型科学工程计算与数值模拟等；2.数据密集型：如数字图书馆、数据仓库、数据挖掘和计算可视化等；3.网络密集型：如协同计算和远程诊断等。

二、并行计算机体系结构

1.并行计算机结构模型

（1）结构类型

• SISD：单指令流单数据流计算机（冯诺依曼机）
• SIMD：单指令流多数据流计算机
• MISD：多指令流单数据流计算机
• MIMD：多指令流多数据流计算机

（2）几种MIMD

• PVP并行向量处理机：多VP（向量处理器）通过交叉开关和多个SM（共享内存）相连
• SMP对称多处理机：多P/C（商品微处理器）通过交叉开关/总线和多个SM（共享内存）相连
• MPP大规模并行处理机：处理节点有商品微处理器+LM（分布式本地内存），节点间通过高带宽低延迟定制网络互联，异步MIMD，多个进程有自己的地址空间，通过消息传递机制通信
• COW工作站机群：节点是完整操作系统的工作站，且有磁盘
• DSM分布共享存储处理机：高速缓存目录DIR确保缓存一致性，将物理分布式LM组成逻辑共享SM从而提供统一地址的编程空间

注：对称指所有处理器都能同等地访问I/O很同样的运行程序（如OS和I/O服务程序），而非对称主从式是仅有主处理器运行OS和控制访问I/O并监控从处理器执行

2.并行计算机访存模型

UMA（Uniform Memory Access）均匀存储访问：物理存储器被所有处理器均匀共享，所有处理器对所有SM访存时间相同，每台处理器可带有高速私有缓存，外围设备共享。
- NUMA非均匀存储访问：共享的SM是由物理分布式的LM逻辑构成，处理器访存时间不一样，访问LM或CSM（群内共享存储器）内存储器比访问GSM（群间共享存储器）快
- COMA（Cache-Only MA）全高速缓存存储访问：NUMA的特例、全高速缓存实现
- CC-NUMA（Coherent-Cache NUMA）高速缓存一致性NUMA：NUMA＋高速缓存一致性协议
- NORMA（No-Remote MA）非远程存储访问：无SM，所有LM私有，通过消息传递通信

3.Cache一致性协议

• 监听总线协议：总线连接通信，写无效和写更新策略
• 基于目录的协议：目录记录共享数据缓存状态，读缺失时查看目录D，写更新时通知目录D

4.其他并行计算概念

衡量并行计算机性能单位：

• PFLOPS：每秒1千万亿 (=10^15) 次的浮点运算
• TFLOPS：每秒1万亿 (=10^12) 次的浮点运算
• GFLOPS：每秒10亿 (=10^9) 次的浮点运算

TOP500前500名超级计算机排名指标(GFLOPS)：
Rmax：Maximal LINPACK(Linear system package) performance achieved
Rpeak：Theoretical peak performance

三、并行计算性能评测

1、机器级的性能测评

主要包括CPU和存储器的某些基本性能指标，并行通信开销以及机器的成本、价格和性能价格比等。
（1）CPU和存储器的某些基本性能指标

• 工作负载：指计算操作的数目，通常可用执行时间、所执行的指令数目和所完成的浮点运算数3个物理量来度量。

• 并行执行时间：$T_n=T_{comput}+T_{paro}+T_{comm}$。其中$T_{comput}$ 为计算时间，$T_{paro}$ 为并行开销时间，包括进程管理（生成、切换和结束等）、组操作（进程组的生成与消亡等）时间，和进程查询（询问进程的标志、等级和组大小等）时间；$T_{comm}$ 为相互通信时间，包括同步（如路障、临界区、事件等）、通信（如点到点通信、整体通信、读/写共享变量等）时间和聚操作（如归约和前缀运算等）时间。

• 存储器的层次结构：各层存储器的容量、延迟和带宽。

• 存储器带宽的估算

（2）通信开销$T_{comm}$

（3）机器的成本、价格和性能价格比

2、算法级的性能测评

包括加速、效率和可扩放性等。

（1）并行系统的加速（比）
指对于一个给定的应用，并行算法（或并行程序）的执行速度相对于串行算法（或串行程序）的执行速度加快了多少倍。

加速比性能定律：
Amdahl加速定律（固定计算负载）
Gustafson定律（适用于可扩放问题）
Sun&Ni定律（受限于存储器）

约定：

• p是处理器数
• 问题规模W=程序中串行分量Ws+可并行部分Wh
• f为串行部分比例，f=Ws/W
• S为加速比

（2）扩放性
最朴素的含义是在确定的应用背景下，计算机系统（或算法或编程等）性能随处理器数的增加而按比例提高的能力。被广泛用来描述并行算法能否有效利用可扩充的处理器的能力。

3、程序级的性能测评

基本测试程序、数学库测试程序和并行测试程序等

四、并行计算理论基础：并行算法设计

1、并行算法

并行算法是一些可同时执行的诸进行的集合，这些进程互相作用和协调动作从而实现给定问题的求解。
从不同角度分类成数值计算和非数值计算的并行算法；同步、异步和分布的并行算法；共享存储和分布存储的并行算法；确定的和随机的并行算法等。

2、并行算法的表达

（1）par-do

n个节点并行完成for循环（每个节点不同，和i相关）：

for i = 1 to n par-do
    ...
endfor

（2）for all

所有节点都执行相同语句：

for all Pi, where 0 <= i <= k do
    ...
endfor

3、并行算法的复杂度

通常分析以下指标：

• 运行时间t(n)：求解问题的时间，包括计算时间和通信时间
• 处理器数目p(n)：求解给定问题所用的处理器数目
• 成本c(n)：c(n) = t(n)*p(n)
• 成本最优——并行算法的成本在数量级上等于最坏情况下串行求解次问题所需要的执行步数
• 工作量W(n)：并行算法完成的总的操作数量
• 工作量最优：功耗低、环保

并行算法的WT表示——Brent定理：

令W(n)是并行算法A在运行时间T(n)内所执行的运算量，则A使用p台处理器可在t(n) = O(W(n)/p + T(n))时间内执行完毕

4、并行计算模型

1.PRAM模型（SIMD-SM）

PRAM（Parallel Random Access Machine）并行随机存取机器，是一种抽象并行计算模型，它假设：

• 存在容量无限大的SM
• 有限或无限个功能相同的处理器，且均有简单算术运算和逻辑判断功能
• 任何时刻各处理器可通过SM交换数据

根据并发访问机制，又分为：

• 不允许同时读和同时写的PRAM-EREW
• 允许同时读但不允许同时写的PRAM-CREW
• 允许同时读和同时写的PRAM-CRCW

PRAM-CRCW又分为：

• 只允许同时写相同的数CPRAM-CRCW
• 只允许优先处理器先写PPRAM-CRCW
• 允许任意处理器自由写APRAM-CRCW

PRAM优点：

• 适合并行算法表达、分析和比较
• 使用简单，屏蔽了通信、存储管理、进程同步等并行细节
• 易于修改算法设计以适应不同并行机

PRAM缺点：

• PRAM是同步模型，同步锁费时
• 不适用于MIMD和DM
• 假设任何处理器可在单位时间内访问任何处理单元而不考虑竞争和带宽，不现实

2.异步APRAM模型（MIMD-SM）

异步APRAM模型假设：

• 每个处理器有LM、局部时钟、局部程序
• 处理器通信经过SM
• 无全局时钟，各处理器异步执行各自指令
• 处理器之间的指令相互依赖关系必须显式加入同步障
• 一条指令可以在非确定但有界时间内完成

指令类型有：

• 全局读：从SM读到LM
• 局部操作：LM操作存入LM
• 全局写：LM写入SM
• 同步：各处理器需等到其他处理器到达后才能继续执行

APRAM比PRAM更加接近实际并行机

3.BSP模型（MIMD-DM）

BSP（Bulk Synchronous Parallel）大同步并行机（APRAM算作轻量）是一个分布式存储的MIMD模型，它的计算由若干全局同步分开的、周期为L的超级步组成，各超级步中处理器做LM操作并通过选路器接收和发送消息；然后做一次全局检查，以确定该超级步是否已经完成（块内异步并行，块间显式同步）

参数：处理器数p、选路器吞吐率g、全局同步间隔L、一个超级步中一个处理器至多发送或接收h条消息

4.LogP模型：MIMD-DM，点到点通讯

LogP模型是分布式存储、点到点通信的MIMD模型

LogP采取隐式同步，而不显式同步障

5、并行算法基本设计策略

1 、串改并
发掘和利用现有串行算法中的并行性，直接将串行算法改造为并行算法

最常用的设计思路但并不普适，好的串行算法一般无法并行化（数值串行算法可以）

快速排序的串改并
SISD上串行执行，最坏情况下O(n^2)，理想情况下O(nlogn)

思路：将O(n)的划分（Partition）并行化是关键，算法：

• 构造一棵二叉排序树，主元是根
• 小于等于主元素处于左子树，大于主元素处于右子树
• 左右子树均是二叉排序树

在CRCW模型上，用伪代码描述如下：

//A[1...n]排序用n个处理器，处理器i中存有A[i]
//f[i]中存有其元素是主元素的处理器号
//LC[1...n]和RC[1...n]分别记录给定主元素的左儿子和右儿子

for each processor i par-do
    root = i    //所有处理器竞争，只有一个写入root
    f[i] = root //所有处理器都把root写入自己的f[i]
    LC[i] = RC[i] = n + 1   //初始值
endfor

repeat for each processor i != root do
    if (A[i] < A[f[i]]) || (A[i] == A[f[i]] && i < f[i])
        LC[f[i]] = i    //并发写，所有满足条件的i只有一个写入LC[f[i]]作为左子树的根，也就是下一次循环的主元素
        if i == LC[f[i]] then
            exit        //若当前处理器写入则什么也不做
        else
            f[i] = LC[f[i]] //若当前处理器没有写入，那么它只能当LC[f[i]]的字节点了
        endif
    else
        RC[f[i]] = i    //并发写，所有满足条件的i只有一个写入RC[f[i]]作为右子树的根，也就是下一次循环的主元素
        if i == RC[f[i]] then
            exit        //若当前处理器写入则什么也不做
        else
            f[i] = RC[f[i]] //若当前处理器没有写入，那么它只能当RC[f[i]]的字节点了
        endif
    endif
endrepeat

每次迭代构造一层排序二叉树花费O(1)时间，在基本平衡的情况下树高O(logn)，则算法复杂度为O(logn)

2 、全新设计
从问题本身描述出发,不考虑相应的串行算法,设计 一个全新的并行算法

有向环K着色算法的并行化设计
问题：有向环顶点着色，一共K种颜色，相邻顶点不允许同色

串行算法（3着色）：交替2种颜色，若顶点是奇数则需要用到第3种颜色（难以并行化）

SIMD-EREW上的并行K着色算法：

//初始随机着色为c[i]，每个顶点着色不同
//输出着色方案为nc[i]
for i = 1 to n par-do
    k = c[i]和c[i的后继]的最低的不同二进制位
    nc[i] = 2 * k + c[i]的二进制第k位
endfor

O(1)时间完成，需要算法正确性证明

3、借用法
找出求解问题和某个已解决问题之间的联系，改造或利用已知算法应用到求解问题上

利用矩阵乘法求所有点对间最短路径
d[k][i][j]表示从vi到vj**至多**经过k-1个中间顶点时的最短路径，d[k][i][j] ＝ min{d[k/2][i][l]+d[k/2][l][j]}

那么可以用矩阵乘法的改进（乘变加，求和变min）做logn次矩阵乘法即可

思路是这样，具体伪代码略，需要O(n^3)个节点，时间复杂度为O(log^2(n))

6、并行算法常用设计技术

1、划分设计技术

使用划分法把问题求解分成两步：

• 把给定问题划分成p个几乎等尺寸的子问题
• 用p台处理器并行求解子问题

（1）均匀划分(PSRS排序)

长度为n的待处理序列均匀划分给p个处理器，每个处理器处理n/p个元素

MIMD机器上的PSRS排序算法：

(1)均匀划分:将n个元素A[1..n]均匀划分成p段，分配给p个处理器
(2)局部排序:pi调用串行排序算法对A[(i-1)n/p+1..in/p]排序
(3)选取样本:pi从其有序子序列A[(i-1)n/p+1..in/p]中选取p个样本元素
(4)样本排序:用一台处理器对p2个样本元素进行串行排序
(5)选择主元:用一台处理器从排好序的样本序列中选取p-1个主元,并播送给其他pi 
(6)主元划分:pi按主元将有序段A[(i-1)n/p+1..in/p]划分成p段
(7)全局交换:各处理器将其有序段按段号交换到对应的处理器中（一定保证均匀划分？？）
(8)归并排序:各处理器对接收到的元素进行归并排序

（2）方根划分(Valiant归并排序)

长度为n的待处理序列，取i*sqrt(n)为划分元，将元素划分成若干段交给处理器处理

SIMD-CREW机器上的Valiant归并排序：

(1)方根划分：把A和B（长n有序段）的第i*sqrt(n)元素作为划分元，把A和B分成若干段
(2)段间比较：A中的所有划分元和B中的所有划分元比较，确定A中划分元应插入B中的哪一段
(3)段内比较：A中的划分元和B中相应段比较并确定插入位置，这些插入位置又将B重新划分成了若干段
(4)段组合并：插入A划分元后，又得到若干有序段组需要归并，递归直到有一组（A）的长度为0

使用n个处理器可以在O(loglogn)内完成

（3）对数划分(并行归并排序)

取i*logn作为划分元划分

定义位序rank(x,X)为X中小于等于x的元素个数，则有PRAM-CREW上的对数划分：

//非降有序的A={a[1],...,a[n]}和B={b[1],...,b[m]}归并
//假设log(m)和k=m/log(m)均为整数
j[0]=0
j[k]=n

for i = 1 to k - 1 par-do
    j[i] = rank(b[ilogm], A)
endfor

for i = 1 to k - 1 par-do
    Bi = {b[ilogm+1], ..., b[(i+1)logm]}
    Ai = {a[j[i]+1], ..., a[j[i+1]]}
endfor

//将原问题转化为子序列组Bi和Ai的归并，那么同样可以递归调用完成整个序列的归并
//对数划分保证Bi和Ai中的元素均大于Bi-1和Ai-1中的元素

（4）功能划分( (m,n)-选择)

功能划分是根据特定问题而把序列分成p个等长组，每组符合问题特性的一种划分方法

(m,n)-选择问题是将长为n的序列中选取前m个较小的元素，利用功能划分来实现并行化的(m,n)-选择问题求解：

(1)功能划分：将A划分成g=n/m组，每组含m个元素
(2)局部排序：使用Batcher排序网络将各组并行排序
(3)两两比较：将所排序的各组两两比较，从而形成MIN序列
(4)排序-比较：对各个MIN序列，重复执行第(2)和第(3)步，直至选出m个最小者

2、分治设计技术

分治将复杂问题划分成较小规模特性相同的子问题，且子问题类型和原问题类型相同，通常用递归完成分治算法

3、平衡树设计技术

以树的叶结点为输入，中间结点为处理结点，由叶向根或由根向叶逐层进行并行处理

4、倍增设计技术

递归调用时，所要处理数据之间的距离逐步加倍, 经过k步后即可完成距离为2^k的所有数据的计算

5、流水线技术

将算法路程分成p个前后衔接的任务片段，一个任务片段完成之后其后继任务片段可以立即开始，那么久可以引入流水线的思想来处理多条数据

五、并行计算软件支撑：并行编程

1、并行程序设计基础

1.1 并行语言构造方法

• 库例程：MPI、Pthreads
• 扩展串行语言：Fortran90
• 加编译注释构造：OpenMP

1.2 并行性问题

可利用SPMD来伪造MPMD

需要运行MPMD：parbegin S1 S2 S3 parend
可以改造成SPMD：

for i = 1 to 3 par-do
    if i == 1 then S1
    else if i == 2 then S2
    else if i == 3 then S3
endfor

那么并行扩展至需要支持SPMD即可

1.3 交互/通信

（1）交互类型

• 通信：进程间传数（共享变量、消息传递、参数传递［父进程传给子进程］）
• 同步：进程间相互等待或继续执行的操作（原子同步、控制同步［同步障、临界区］、数据同步［锁、condition、监控程序和事件］）
• 聚合：将分进程所计算的的结果整合起来（规约、扫描）

（2）交互方式

• 同步交互：所有参与者全部到达后继续执行
• 异步交互：任意进程到达后不必等待其他进程即可继续执行

（3）交互模式

按编译时是否能确定交互模式可分为静态的、动态的

按多少发送者和接收者：

• 一对一：点到点通信
• 一对多：广播、散播
• 多对一：收集、规约
• 多对多：全交换、扫描、置换、移位

1.4 并行编程风范

• 相并行：BSP模式，程序由一组超级步组成，步内各自并行计算，步间通信同步
• 主从并行：主进程串行执行并且协调任务，子进程计算任务，需要划分设计并结合相并行
• 分治并行：父进程把负载分割并指派给子进程，难以平衡负载
• 流水线并行：进程划分成流水线，依次依赖，数据开始流动
• 工作池并行：进程从工作池中取任务执行

1.5 并行程序设计模型

（1）隐式并行

用串行语言编程，编译器货操作系统自动转化成并行代码

特点：语义简单、可以执行好、易调试易验证、but效率低

（2）数据并行

SIMD模型，包括数据选路和局部计算，特点：但现场、松散同步、常用聚合操作

数据并行计算π：

long i,j,t,N=100000;
double local[N], temp[N], pi, w;
w = 1.0/N;

forall (i = 0; i < N; i++) {
    local[i] = (i + 0.5) * w;
    temp[i] = 4.0 / (1.0 + local[i] * local[i]);
}

pi = reduce(temp, +);

（3）消息传递

MPP、COW自然模型，MPI广泛应用：多线程异步并行、地址空间分开、常用SPMD形式编码

MPI消息传递计算π：

#define N 100000
main(){
    double local=0.0,pi,w,temp=0.0;
    long i,taskid,numtask;

    w=1.0/N;

    MPI_Init(&argc,&argv);
    MPI_Comm_rank(MPI_COMM_WORLD,&taskid);
    MPI_Comm_Size(MPI_COMM_WORLD,&numtask);

    for (i = taskid; i < N; i=i + numtask){
        temp = (i+0.5)*w;
        local = 4.0/(1.0+temp*temp)+local; 
    }

    MPI_Reduce(&local,&pi,1,MPI_Double,MPI_MAX,0, MPI_COMM_WORLD);

    if (taskid == 0) printf(“pi is %f \n”,pi*w); 

    MPI_Finalize() ;
}

（4）共享变量

PVP、SMP、DSM自然模型，多线程异步，显示同步而隐式通信

OpenMP使用共享变量计算π：

#define N 100000
main(){
    double local,pi=0.0,w;
    long i;

    w=1.0/N;

    #pragma parallel
    #pragma shared(pi, w)
    #pragma local(i, local) 
    {
        #pragma parallel for (i = 0; i < N; i++)
        {
            local = (i + 0.5) * w;
            local = 4.0 / (1.0 + local * local);
        }
        #pragma critical
        {
            pi = pi + local
        }
    }
}

六、共享存储系统并行编程——OpenMP编程

OpenMP是FORK-JOIN模型，主线程串行执行，直到编译制导并行域出现

七、分布存储系统并行编程——MPI编程

MPI是一种消息传递接口的标准，适用于分布式存储系统的编程模型

与OpenMP不同的是，MPI是多进程的并行模式，运行时需要在外部指定开启进程数，并且是用SPMD的编程风格去模拟MPMD的编程风格（用进程号区别），不会FORK-JOIN而是通过消息传递同步

八、GPU体系结构及编程——CUDA编程