并行程序设计及ScaLAPACK函数库应用简介 报告人：罗正平 导师：肖炳甲研究员 Report for Friday Symposium (1/29/2010) 2017/9/11

并行计算--高性能计算 并行计算（Parallel Computing） 是指在并行计算机上，将一个应用分解成多个子任务，分配给不同处理器，各个处理器之间相互协同，并行地执行子任务，从而达到加快求解的速度或提高求解应用问题规模的目的。 高性能计算（High Performance Computing，HPC) 超级计算（Super Computing） 2017/9/11

为什么需要并行计算机 问题: 科学和工程问题的数值模拟与仿真 计算密集 数据密集 网络密集 三种混合 要求:在合理的时限内完成计算任务 秒级 制造业 分钟级 短时天气预报(当天) 小时级 中期天气预报(3~10日) 尽可能快 长期天气预报(气候) 可计算湍流模拟 2017/9/11

并行计算的基本条件 并行计算机 应用问题必须具有并行度 并行编程 2017/9/11

并行计算机 由多个计算单元组成，运算速度快、存储容量大、可靠性高的计算机系统。 也称巨型机或超级计算机。 Roadrunner 1.026PetaFlops 本页图片是IBM的最新军用超级计算机“Roadrunner”，每秒计算能力达到了1.026PetaFlops，是目前最强的 IBM BlueGene/L的两倍还多。Roadrunner一共拥有116640颗计算核心，由三种不同的处理器组成，包括12960颗改进版的IBM Cell，以及少量的AMD Opteron。Roadrunner超级计算机的布线总长达到了57公里，功率为3.9兆瓦，占地约6000平方英尺，总重超过500000磅。由IBM和洛斯阿拉莫斯(Los Alamos)国家实验室技术人员共同开发和组装。其计算能力形象的描述是：如果全球60亿人每人都拿上一台手持计算器，然后一天24小时不停地计算，如此不停地继续工作64年，其总工作量仅相当于Roadrunner一天的工作量。（美国能源部下属国家核安全管理局(NNSA)局长托马斯·阿戈蒂诺(Thomas Agostino) ） 2017/9/11

并行计算机体系结构 对称多处理机SMP（Symmetric Multiprocessor） 分布式共享存储并行计算机DSM（Distributed Shared Memory） 大规模并行计算机MPP（Massively Parallel Processing） 2017/9/11

机群系统（cluster） 每个结点都是一个完整的计算机；每个结点包含数个微处理器，结点 内部采取共享存储 各个结点之间由高性能互联网络链接 网络接口和I/O总线松耦合链接 每个结点都有完整的操作系统 2017/9/11

Cluster based on SMP -- MPP IBM刀片机群（HS21） Cluster based on SMP -- MPP 2017/9/11

并行程序开发策略 应用问题的并行度，也就是应用问题可以分解为多个可并行执行的子任务。 现有的串行源代码 有目的稍许修改源代码 自动并行化 并行应用程序 (a)自动并行化 开发并行库 重新链接 (b)并行库 作重大修改 编译器支持并行化 (c)重新编写并行代码 2017/9/11

并行程序设计 共享变量模型（Shared Variable） 消息传递模型（Message Passing） 2017/9/11

共享变量（Shared Variable） OpenMP OpenMP是基于线程的并行编程模型，使用Fork-Join并行执行模型。所有的OpenMP程序开始于一个单独的主线程，主线程串行执行直到遇到并行域才开始并行执行。 所有的OpenMP并行化，都是通过使用嵌入到C/C++或Fortran源代码中的编译制导语句来达到的。 Fork Join 主线程 并行域 #include “omp.h” int main(int argc, char* argv[]) { int nthreads, tid; int nprocs; char buf[32]; /* Fork a team of threads */ # pragma omp parallel private(nthreds,tid) /* Obtain and print thread id */ tid=omp_get_thread_num(); printf(“Hello World from OMP thread %d\n”,tid); /* Only master thread does this */ if(tip==0){ nthreads=omp_get_num_threads(); printf(“Number of threads %d\n”,nthreads); } } return 0; 2017/9/11

消息传递（Data Parallel） MPI 互连网络 …… Include “mpif.h” 其他参数说明 Call MPI_Init() Call MPI_Comm_rank() Call MPI_Comm_size() 建立进程拓扑结构和新的通信器 应用程序体： 1.计算控制程序体； 2.进程间通信。 Call MPI_Finalize() end P/M 互连网络 MPI …… 基于MPI的并行程序，是有一个或多个彼此通过调用库函数进行消息收发通信的进程所组成。 绝大部分MPI实现中，一组固定的进程在程序初始化时生成，可以执行SPMD或MPMD模式。 2017/9/11

MPI 数据交换 进入MPI系统 2017/9/11 退出MPI系统 MPI函数 进程0 进程1 进程2 进程3 init() program example1 include “mpif.h” !! MPI系统头文件 integer status(MPI_STATUS_SIZE),my_rank,p,source,dest,tag,ierr,data c-------进入MPI系统 call MPI_Init(ierr) call MPI_Comm_rank(MPI_COMM_WORLD,my_rank,ierr) call MPI_Comm_size(MPI_COMM_WORLD,p,ierr) c-------数据交换 data=0 tag= 5 source= my_rank-1 if(source.eq.-1) source=p-1 dest =my_rank+1 if(dest.eq.p) dest=0 if(my_rank.eq.0) then call MPI_Send(data,1,MPI_INTEGER,dest,tag,MPI_COMM_WORLD,ierr) call MPI_Recv(data,1,MPI_INTEGER,source,tag,MPI_COMM_WORLD,status,ierr) else data=data+1 endif c-------广播数据 call MPI_Bcast(data,1,MPI_INTEGER,0,MPI_COMM_WORLD,ierr) c------打印输出 if(data.eq.p-1) then print *,”Successful, data=”,data print *,”Failure, data=”,data call MPI_Finalize(ierr) end init() 进入MPI系统 comm_rank() myrank=0 myrank=1 myrank=2 myrank=3 comm_size() p=4 p=4 p=4 p=4 data=0 data=0 data=0 data=0 source=3 source=0 source=1 source=2 dest=1 dest=2 dest=3 dest=0 数据交换 send() recv() recv() recv() recv() data=data+1 Example1：进程0发送一个整数data给进程1；进程1将该数加1后传递给进程2；进程2加1后将该数传递给进程3；以此类推，最后进程p-1将该数传递给进程0，由进程0负责广播给所有进程，并由进程0负责打印输出。 此处以P=4为例实现。 send() data=data+1 send() data=data+1 send() broadcast() send() recv() recv() recv() output “data” 2017/9/11 退出MPI系统 finalize()

ScaLAPACK 并行库开发高性能计算程序基本思想： 并行函数库提供经过优化的通用并行算法代码，用户只需根据需要调用相应函数就可编写并行程序。 ScaLAPACK (Scalable LAPACK)是LAPACK在分布式计算环境中的扩展，主要运行在基于分布式存储和消息传递机制的MIMD计算机以及支持PVM或MPI的集群上。 2017/9/11

ScaLAPACK层次结构 BLAS: Basic Linear Algebra Subprograms PBLAS: Parallel Basic Linear Algebra Subprograms BLACS: Basic Linear Algebra Communication Subprograms LAPACK: Linear Algebra PACKage ScaLAPACK: Scalable Linear AlgebraPACKage 2017/9/11

BLAS简介 BLAS是一组高质量的基本向量矩阵运算子程序。其从结构上分成3部分: Level 1 BLAS: 向量和向量，向量和标量之间的运算； Level 2 BLAS: 向量和矩阵间的运算； Level 3 BLAS：矩阵和矩阵之间的运算。 BLAS支持四种浮点格式运算：单精度实数(REAL)、双精度实数(DOUBLE)、单精度复数(COMPLEX)和双精度复数(COMPLEX*16)。对应子程序的首字母分别为S、D、C和Z。 1. BLAS的主要贡献是将高性能代数计算程序的开发同针对特定机器的性能优化独立开来：代数算法程序的开发者只需运用适当的分块技术将计算程序变成矩阵向量的基本运算并调用相应的BLAS子程序，而不必考虑与计算机体系结构相关的性能优化问题（非常繁杂）。 2. BLAS中除了正文提到的三个部分外，还包括一个辅助子程序XERBLA，用于打印错误信息。 3. 在优化BLAS库中，层次越高的子程序性能改善越大；在使用BLAS库的一个基本原则是：尽可能地使用Level 3 BLAS中的子程序，其次才是Level 2 BLAS。 2017/9/11

LAPACK简介 LAPACK是建立在BLAS库基础上的线性代数函数库。包含了求解科学与工程计算中常见的数值线性代数计算问题，如线性方程组、线性最小二乘问题、特征值问题和奇异值问题等。还可以实现矩阵分解和条件数估计等相关计算。 LAPACK子程序可以分成3类： (1)驱动程序（driver routines）：用于解决一个完整问题，如线性方程组求解； (2)计算程序（computional routines）：也叫做简单LAPACK子程序，用以完成一个特定的计算任务，如矩阵的LU分解； (3)辅助程序（auxiliary routines）：是被驱动程序和计算程序调用的子程序。主要完成对子块的操作和一些常用的底层计算，如计算矩阵范数等。 LAPACK最初目标是在共享存储向量并行计算机上高效地使用EISPACK和LINPACK。而LINPACK和EISPACK忽视微处理器的多层存储结构，以向量操作的形式调用Level 1 BLAS中的子程序完成基本运算，效率低下。LAPACK对矩阵进行分块，通过分块将许多操作转化为矩阵计算，主要是矩阵乘法，可以直接调用Level 3 BLAS子程序，大大提高计算效率。 2017/9/11

BLACS简介 BLACS是为线性代数设计的消息传递函数库。计算模型由一个一维或二维进程网格构成，每个进程存储矩阵和向量的一些片断。 BLACS是建立在PVM或MPI等底层消息传递函数库基础上的，其目标是为通信提供专用于线性代数的可移植层。 1 2 3 Processor ID, row numbering and column numbering all begin with 0 Have chosen to insert the processors into the grid “by-row” Example: processor 6 has grid coordinates (p,q)=(1,2) 1 2 3 4 5 6 7 1 BLACS二维进程网格 2017/9/11

PBLAS简介 PBLAS即并行BLAS，执行消息传递并且其接口与BLAS基本相似。在基于ScaLAPACK函数库的程序中，进程之间的通信出现在PBLAS内部，从而使ScaLAPACK代码与LAPACK代码相当接近，甚至几乎一样，简化程序设计难度。 PBLAS结构上与BLAS一样分成三部分，其子程序名为BLAS子程序名之前加上P，表示并行。 2017/9/11

ScaLAPACK Local：表示本地组件在一个进程中被调用，其参数只存储在一个进程中。 Global：全局组件是同步的并行程序，其参数包括矩阵和向量，分布在多个进程中。 2017/9/11

ScaLAPACK基本步骤 1 1 2 3 1、初始化进程网格； 2、分配数据到进程网格上； 3、调用ScaLAPACK函数； 1 ScaLAPACK基本步骤 1 2 3 1、初始化进程网格； 2、分配数据到进程网格上； 3、调用ScaLAPACK函数； 4、释放进程网格和退出并行应用。 1 1 a11 a12 a21 a22 a13 a14 a23 a24 a15 a25 a31 a32 a41 a42 a33 a34 a43 a44 a35 a45 a51 a52 a53 a54 a55 a11 a12 a15 a21 a22 a25 a51 a52 a55 a13 a14 a23 a24 a53 a54 a31 a32 a35 a41 a42 a45 a33 a34 a43 a44 2-dimension block-cyclic distribution 1 2017/9/11

Example 16＊16 16＊1 串行：SGEMV －－ LAPACK/BLAS SUBROUTINE SGEMV(TRANS,M,N,ALPHA,A,LDA,X,INCX,BETA,Y,INCY) SGEMV: y := alpha*A*x + beta*y, or y := alpha*A'*x + beta*y, 2017/9/11

Example 2 3 16＊16 16＊1 1 PSGEMV 1－8 9－16 1－8 9－16 1 2017/9/11 1 2 3 16＊16 16＊1 1－8 9－16 1－8 9－16 1 PSGEMV—PBLAS子程序 Subroutine PSGEMV(Trans,M,N,A,IA,JA,DESCA,X,IX,JX,DESCX,INCX,BETA,Y,IY,JY,DESCY,INCY) PSGEMV 2017/9/11

谢谢 2017/9/11