MPI的基础知识
MPI(Message Passing Interface 11)英文直译就是基于消息传递的用户界面。
什么是MPI
- MPI是一个库，而不是一门编程语言。MPI库中包含了可以被Fortran/C调用的一系列子例程和函数，MPI库函数/子例程遵守所有对库函数/子例程的调用规则，和一般的函数/过程没有什么区别 。MPI程序的源代码被C或Fortran的编译器编译之后，再与与MPI库链接，就生成了可执行的MPI程序。
- MPI是一种标准或规范的代表，而不特指某一个对它的具体实现（如openMPI，MPICH等）。迄今为止所有的并行计算机制造商都提供对MPI的支持 。
- MPI是一种消息传递的编程模型。最终目的是完成进程间的通信。
MPI的目的
- 较高的通信性能
- 较好的程序可移植性
- 强大的功能
MPI的语言绑定
目前标准的MPI库都提供Fortran和C的API。
- Fortran — 主要用于科学与工程计算
- C — 主要用于系统和应用程序的开发
目前主要的 MPI 实现
以下两个执行MPI标准的库，均为免费开源。
MPI的历史
- 1992年 — MPI 1.0版本问世
- 1997年 — MPI 2.0版本出现
- 2012年 — MPI 3.0版本出现
- 2018年 — MPI 4.0版本出现，目前最新版本
关于MPI的一些名称和概念
- 进程（ process）进程与程序相联， 程序一旦在操作系统中运行即成为进程。进程拥有独立的执行环境（内存、寄存器、程序计数器等），是操作系统中独立存在的可执行的基本程序单位， 串行应用程序编译形成的可执行代码，分为“指令”和“数据”两个部分，并在程序执行时独立地申请和占有内存空间，且所有计算均局限于该内存空间。多个进程可同时存在于单机内同一操作系统，操作系统负责调度分时共享处理机资源（ CPU、内存、存储、外设等），进程间相互独立（内存空间不相交） 。进程间可以相互交换信息， 例如数据交换、同步等待，消息是这些交换信息的基本单位，消息传递是指这些信息在进程间的相互交换，是实现进程间通信的唯一方式 。
- 进程组（ process group） 指MPI 程序的全部进程集合的一个有序子集，且进程组中每个进程被赋于一个在该组中唯一的序号(rank)， 用于在该组中标识该进程，序号的取值范围是[0,进程数- 1] 。
- 通信器（ communicator）通信器也叫通信域，可以理解为一类进程的集合即一个进程组，且在该进程组，进程间可以相互通信。任何MPI通信函数均必须在某个通信器内发生。 MPI系统提供省缺的通信器MPI_COMM_WORLD，所有启动的MPI进程通过调用函数MPI_Init()包含在该通信器内；各进程通过函数 MPI_Comm_size()获取通信器包含的的MPI进程个数。MPI系统在一个MPI程序运行时会自动创建两个通信器：一个称为MPI_COMM_WORLD，它包含 MPI 程序中所有进程，另一个称为 MPI_COMM_SELF，它指单个进程自己所构成的通信器。
- 进程序号（ rank） 用来在一个进程组或通信器中标识一个进程，MPI程序中的进程由进程组或通信器序号唯一确定， 序号相对于进程组或通信器而言（假设np个处理器，标号0…np-1）。同一个进程在不同的进程组或通信器中可以有不同的序号，进程的序号是在进程组或通信器被创建时赋予的。MPI 系统提供了一个特殊的进程序号MPI_PROC_NULL，它代表空进程(不存在的进程)， 与MPI_PROC_NULL 间的通信实际上没有任何作用 。
- 消息（ message） 分为数据（ data）和包装（ envelope） 两个部分。 包装由接收进程序号/发送进程序号、消息标号和通信器三部分组成，数据包含用户将要传递的内容 。
- 通信（ communication） 指在进程之间进行消息的收发、同步等操作。
- MPI对象 MPI系统内部定义的数据结构，包括数据类型（如MPI_INT）、通信器（ MPI_Comm）、通信请求（ MPI_Request）等，它们对用户不透明。在Fortran语言中，所有MPI对象均必须说明为整型变量INTEGER 。
- MPI联接器（ handles） 联接MPI对象的具体变量，用户可以通过它访问和参与相应MPI对象的具体操作。例如， MPI系统内部提供的通信器MPI_COMM_WORLD。在Fortran语言中，所有MPI联接器均必须说明为整型变量INTEGER。
一个最简单的MPI程序 Hello world!
运行一个MPI程序的过程可以概括为：首先，将MPI程序源代码编译，与MPI库函数/子例程链接之后生成可执行的MPI程序；接下来，向操作系统提交MPI可执行程序，操作系统再将可执行程序复制N份交给N个处理器；然后，每个处理器各自运行自己的MPI程序，运行过程会调用MPI库函数子例程，得到该节点上的进程my_id，根据进程号程序执行不同的任务，并适时地进行不同进程间的消息传递，以完成并行的任务。
Hello world MPI程序的Fortran语言实现
PROGRAM hello_world
! or the older form: INCLUDE ’mpif.h’ !引入头文件
USE MPI ! 调用MPI库
IMPLICIT NONE
INTEGER :: my_id, proc_num, ierr, name_len
CHARACTER(LEN=MPI_MAX_PROCESSOR_NAME) :: processor_name
!MPI_MAX_PROCESSOR_NAME是MPI预定义的宏，即MPI所允许的机器名字的最大长度
CALL MPI_INIT(ierr) !MPI初始化
CALL MPI_COMM_RANK(MPI_COMM_WORLD,my_id,ierr) !得到当前正在运行的进程的标识号
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,proc_num,ierr) !得到所有参加运算的进程的个数
CALL MPI_GET_PROCESSOR_NAME(processor_name, name_len, ierr) !得到运行本进程的机器的名称
WRITE(*,10) my_id,proc_num,processor_name
10 FORMAT('Hello World! Process ',I2,' of ',I2,' on ', 20A)
CALL MPI_FINALIZE(ierr) !终止MPI
END PROGRAM
MPI程序的编译和运行
Fortran 语言的MPI程序由 mpif90 来编译，用法与一般的的 Fortran 编译器（如 Gfortran ）类似，向系统提交MPI程序可由 mpirun 命令完成，-np 参数表示该程序的进程数。
mpif90 -c hello_world.f90
mpif90 -o hello_world hello_world.o
mpirun -np 5 hello_world
运行结果如下图
MPI的一些惯例
- MPI 的所有常量、变量与函数 /过程均以MPI_ 开头。
- MPI 的 C 语言接口为函数， Fortran 接口为子程序，且对应接口的名称相同。
- 在C 程序中，所有常数的定义除下划线外一律由大写字母组成，在函数和数据类型定义中, 接MPI_ 之后的第一个字母大写，其余全部为小写字母，即MPI_Xxxx_xxx形式。
- 对于FORTRAN 程序， MPI 函数全部以过程方式调用，一般全用大写字母表示，即MPI_XXXX_XXX形式（ Fortran不区分大小写）。
- 除 MPI_WTIME 和 MPI_WTICK 外，所有C函数调用之后都将返回一个错误信息码，而 MPI 的所有 FORTRAN 子程序中都有一个哑元参数（ ierr）代表调用错误码，错误代码返回值为0表示函数/过程调用成功，返回值为非0则表示调用失败。
- MPI 是按进程组(Process Group) 方式工作的，所有MPI 程序在开始时均被认为是在通信器 MPI_COMM_WORLD 所拥有的进程组中工作，之后用户可以根据自己的需要，建立其它的进程组。
- 所有MPI 的通信一定要在通信器(communicator) 中进行。
- Fortran的数组下标是以1开始，而C语言的数组是以0开始。
- 由于C语言的函数调用机制是值传递，所以MPI的所有C函数中的输出参数用的都是指针 。
MPI基本的用法
MPI 的六个基本接口
MPI有上百个可调用的函数（子例程），具体调用语法可查阅 Open MPI 的说明文档，实际编写MPI程序中最常使用的MPI函数（子例程）的个数是十分有限的 ，并且理论上讲MPI所有的通信功能可以用它的6个基本的函数（子例程）来实现。以下以Fortran语言的MPI实现为例，说明MPI的6个基本的接口：
- MPI初始化
CALL MPI_INIT(ierr)
! OUT INTEGER :: ierr -返回值为0表示成功，返回值为非0表示出错
- MPI结束
CALL MPI_FINALIZE(ierr)
! OUT INTEGER :: ierr
- 得到当前进程在给定通信域（也叫通信子）中的编号
CALL MPI_COMM_RANK(comm,my_id,ierr)
! IN INTEGER :: comm -给定的通信域，可以为MPI预定义的通信域MPI_COMM_WORLD（包含所有进程），也可以为其他自定义通信域
! OUT INTEGER :: my_id -当前节点的进程号
! OUT INTEGER :: ierr -给定的通信域
- 得到给定通信域内的进程数
call MPI_COMM_SIZE(comm,proc_num,ierr)
! IN INTEGER :: comm -给定的通信域
! OUT INTEGER :: proc_num -给定通信域内的进程数
! OUT INTEGER :: ierr
- 消息发送
CALL MPI_SEND(buf,count,datatype,dest,tag,comm,ierr)
! IN <TYPE> :: buf -发送数据的起始地址
! IN INTEGER :: count -发送数据的个数
! IN <handle> :: data_type -发送数据的类型，可以是MPI预定义的数据类型，也可以是用户的自定义的新的数据结构（非连续）
! IN INTEGER :: dest -目的进程的编号
! IN INTEGER :: tag -此次发送的标签
! IN INTEGER :: comm -给定的通信域
! OUT INTEGER :: ierr
说明：MPI_SEND将count个datatype数据类型的数据发送到目的进程， 目的进程在通信域中的标识号是dest，本次发送的消息标志是tag，使用这一标志就可以把本次发送的消息和本进程向同一目的进程发送的其它消息区别开来 ， 其中datatype数据类型可以是MPI的预定义类型 也可以是用户自定义的类型 （后文介绍）。
- 消息接收
CALL MPI_RECV(buf,count,datatype,source,tag,comm,status,ierr)
! OUT <TYPE> :: buf -接受数据的起始地址
! IN INTEGER :: count -接受数据的个数
! IN <handle> :: data_type -接收数据的类型，可以是MPI预定义的数据类型，也可以是用户的自定义的新的数据结构（非连续）
! IN INTEGER :: source -所接受数据的进程的编号，可以为MPI_ANY_SOURCE（MPI预定义常量），用以匹配任意消息来源
! IN INTEGER :: tag -此次接收的标签，可以为MPI_ANY_TAG（MPI预定义常量）用以匹配任意消息标签
! IN INTEGER :: comm -给定的通信域
! OUT INTEGER :: status(MPI_STATUS_SIZE) 返回状态和调试信息
! OUT INTEGER :: ierr
说明：MPI_RECV从指定的进程source接收消息，发送消息的数据类型datatype和消息标识tag要与接收进程指定的datatype和tag相一致，接收消息的起始地址为buf，接收消息的长度必须不小于发送消息的长度，这是因为如果接收到的数据过大，接收缓冲区会发生溢出错误，如果一个长度短于count的消息到达，那么只有相应于这个消息的那些地址被修改，count可以是零，这种情况下消息的数据部分是空的，其中datatype数据类型可以是MPI的预定义类型 也可以是用户自定义的类型 （后文介绍），status是包含MPI_STATUS_SIZE个整数的数组，status(MPI_SOURCE)，status(MPI_TAG)和status(MPI_ERROR)分别表示发送数据的进程标识，发送数据使用tag标识和该接收操作返回的错误代码 。
MPI 预定义数据类型
MPI系统中数据的发送和接收都是基于数据类型进行的。数据类型可以是MPI系统预定义的，称为原始数据类型；也可以是用户在原始数据类型的基础上自己定义的新的数据类型。MPI 的基本数据类型定义与相应的Fortran 的数据类型对照关系如下:
| MPI Datatype | Fortran Datatype |
|---|---|
| MPI_INTEGER | INTEGER |
| MPI_REAL | REAL |
| MPI_DOUBLE_PRECISION | DOUBLE PRECISION |
| MPI_COMPLEX | COMPLEX |
| MPI_LOGICAL | LOGICAL |
| MPI_CHARACTER | CHARACTER(1) |
| MPI_BYTE（由一个字节组成 ） | 无类型与之对应 |
| MPI_PACKED | 无类型与之对应 |
在MPI中，消息传递要做到数据类型的匹配。数据类型匹配有两方面的含义：一方面，宿主语言（Fortran或C）的类型和通信操作所指定的类型相匹配，比如在 Fortran 中声明为 INTEGER 类型的变量在发送和接收时要使用 MPI_INTEGER 与之相对应，声明为REAL类型的变量在发送和接收时要使用 MPI_REAL 与之相对应 前面部分给出的 Fortran 与 MPI 类型的对应关系就是为了满足类型匹配的要求；另一方面，发送方和接收方的类型相匹配，即发送方用 MPI_INTEGER 则接收方也必须使用 MPI_INTEGER，发送方用 MPI_REAL 则接收方也必须用MPI_REAL。但是，对于 MPI 提供的 MPI_BYTE 和 MPI_PACKED 它们可以和任何以字节为单位的存储相匹配，包含这些字节的类型是任意的，MPI_BYTE用于不加修改地传送内存中的二进制值， MPI_PACK用于数据的打包和解包（MPI_UNPACK ）。
MPI消息的组成
MPI 消息包括信封和数据两个部分。信封指出了发送或接收消息的对象及相关信息， 而数据是本消息将要传递的内容。信封和数据又分别包括三个部分，可以用一个三元组来表示： \(信封 <源/目,标识,通信域>\\ 数据 <起始地址,数据个数,数据类型>\) 以MPI_SEND和MPI_RECV为例，下图分别给出了它们的信封和数据部分 ：
死锁现象
编写MPI程序如果通信调用的顺序使用的不当很容易造成死锁。若所有进程都首先进行接收数据操作，那么所有进程都会等待从其他进程接收数据，而没有进程进行数据发送操作，因此程序会停滞，从而发生死锁。若所有的进程都首先进行发送数据操作，由于发送数据系统提供缓存区，若系统缓存区不足，则发送数据将会失败，显然把所有进程的发送数据操作都写在接收操作的前面也有可能是不安全的。因此MPI程序的发送必须合理的安排顺序，否则将会导致死锁，捆绑式的发送和接收操作 MPI_SENDRECV 也可以解决这个问题。
其他常用的MPI函数/子例程
计时函数
在MPI程序中经常会用到时间函数，比如用来统计程序运行的时间，或根据时间的不同选取不同的随机数种子，或根据时间的不同对程序的执行进行控制等。 MPI提供了时间函数可供调用 。
MPI_WTIME() !DOUBLE PRECISION
MPI_WTICK() !DOUBLE PRECISION
MPI_WTIME（）返回一个用浮点数表示的秒数， 它表示从过去某一时刻到调用时刻所经历的时间 ，MPI_WTICK返回MPI_WTIME的精度单位是秒 。
获取主机名和MPI版本号
在实际使用MPI编写并行程序的过程中，经常要将一些中间结果或最终的结果输出到程序自己创建的文件中，对于在不同机器上的进程 常希望输出的文件名包含该机器名，或者是需要根据不同的机器执行不同的操作，这样仅仅靠进程标识rank是不够的 MPI为此提供了一个专门的调用，使各个进程在运行时可以动态得到该进程所运行机器的名字 。另外MPI也提供了可以获得MPI版本号的调用。
CALL MPI_GET_PROCESSOR_NAME(processor_name,name_len,ierr)
! OUT CHARACTER(LEN=MPI_MAX_PROCESSOR_NAME) :: processor_name -主机名
! OUT INTEGER :: name_len -名长度
! OUT INTEGER :: ierr
CALL MPI_GET_VERSION(version, sub_version,ierr)
! OUT INTEGER :: version -主版本号
! OUT INTEGER :: sub_version -子版本号
! OUT INTEGER :: ierr
检测MPI是否初始化及异常终止MPI
MPI 程序中唯一一个可以用在 MPI_INIT 之前的MPI调用就是 MPI_INITALIZED，它的功能就是判断 MPI_INIT 是否已经执行。在编写 MPI 程序的过程中若发现已出现无法恢复的严重错误因而只好退出MPI程序的执行，MPI 提供了 MPI_ABORT 调用,并且在退出时可以返回给调用环境一个错误码.
CALL MPI_INITIALIZED(flag, ierr)
! OUT LOGICAL :: flag - MPI 是否初始化
! OUT INTEGER :: ierr
CALL MPI_ABORT(comm, errorcode, ierr)
! IN INTEGER :: comm -要退出的通信域
! IN INTEGER :: errorcode -返回到所嵌环境的错误码
! OUT INTEGER :: ierr
同步MPI过程
MPI_BARRIER用于一个通信域中所有进程的同步，调用MPI_BARRIER时进程将处于等待状态，直到通信子中所有进程 都调用了该函数后才继续执行。
MPI_BARRIER(comm, ierr)
! OUT INTEGER :: comm -通信域
! OUT INTEGER :: ierr
捆绑发送和接收
MPI提供了捆绑发送和接收操作，可以在一条MPI语句中同时实现向其它进程的数据发送和从其它进程接收数据操作。捆绑发送和接收操作把发送一个消息到一个目的地和从另一个进程接收一个消息合并到一个调用中，源和目的可以是相同的。捆绑发送接收操作虽然在语义上等同于一个发送操作和一个接收操作的结合，但是它可以有效地避免由于单独进行发送或接收操作时，由于次序的错误而造成的死锁。这是因为该操作由通信系统来实现，系统会优化通信次序从而有效地避免不合理的通信次序，最大限度避免死锁的产生。捆绑发送接收操作是不对称的，即一个由捆绑发送接收调用发出的消息可以被一个普通接收操作接收，一个捆绑发送接收调用可以接收一个普通发送操作发送的消息。
CALL MPI_SENDRECV(sendbuf, sendcount, sendtype, dest, sendtag,
recvbuf, recvcount, rectype, source, recvtag, comm,
status, ierr)
! IN <type> :: sendbuf -发送数据的起始地址
! IN INTEGER :: sendcount, sendtype, dest, sendtag -发送数据的个数，类型，目标进程号，消息标识
! IN INTEGER :: recvcount, recvtype, source, recvtag -接收数据的个数，类型，目标进程号，消息标识
! IN INTEGER :: comm -通信域
! OUT INTEGER :: STATUS(MPI_STATUS_SIZE) -返回状态
! OUT INTEGER :: ierr
! OUT <type> :: recvbuf -接收数据的起始地址
一个与 MPI_SENDRECV 类似的操作是 MPI_SENDRECV_REPLACE，它与MPI_SENDRECV的不同就在于其发送数据的地址与接收数据的地址相同。这一调用的执行结果是，将指定地址中原来的数据被传递给指定的目的进程，然后该地址上的数据又被从指定进程接收到的相应类型的数据所取代。
CALL MPI_SENDRECV_REPLACE(buf, count, type, dest, sendtag, source, recvtag, comm, status, ierr)
! INOUT <type> :: buf -发送以及接收数据的起始地址
! IN INTEGER :: count, type -发送以及接收数据的个数，类型
! IN INTEGER :: dest, sendtag -发送数据的目标进程号，消息标识
! IN INTEGER :: source, recvtag -接收数据目标进程号，消息标识
! IN INTEGER :: comm -通信域
! OUT INTEGER :: STATUS(MPI_STATUS_SIZE) -返回状态
! OUT INTEGER :: ierr
! OUT <type> :: recvbuf -接收数据的起始地址
简单的MPI并行程序
计算n!
PROGRAM sum_n
USE MPI
! or the older form: INCLUDE ’mpif.h’
IMPLICIT NONE
INTEGER :: my_id, proc_num, ierr
INTEGER, DIMENSION(MPI_STATUS_SIZE) :: satus
INTEGER :: sum_global, sum_local
INTEGER :: n_start, n_global, i
REAL(KIND=8) :: time_start,time_end,time_total
CALL MPI_INIT(ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD,my_id,ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,proc_num,ierr)
time_start = MPI_WTIME()
n_global = 10000
sum_local = 0
n_start = my_id+1
DO i=n_start,n_global,proc_num
sum_local = sum_local+i
END DO
WRITE(*,*) " "
WRITE(*,*) "==========================================="
WRITE(*,10) my_id,proc_num
10 FORMAT('Process ',I2,' of ',I2)
WRITE(*,*) "The local caculation finished."
WRITE(*,*) "n_start=", n_start, "n_end=", i-proc_num
WRITE(*,*) "sum_local=", sum_local
WRITE(*,*) "==========================================="
CALL MPI_REDUCE(sum_local,sum_global,1,MPI_INTEGER,MPI_SUM,0,MPI_COMM_WORLD,ierr)
time_end = MPI_WTIME()
time_total = time_end-time_start
IF( my_id == 0 )THEN
WRITE(*,*) " "
WRITE(*,*) "==========================================="
WRITE(*,*) "The global caculation finished."
WRITE(*,*) "The process is ", my_id
WRITE(*,*) "sum_local=", sum_global
WRITE(*,*) "Total time is", time_total
WRITE(*,*) "==========================================="
END IF
CALL MPI_FINALIZE(ierr)
END PROGRAM
接收相邻进程的数据
PROGRAM send_recv_id
USE MPI
! or the older form: INCLUDE ’mpif.h’
IMPLICIT NONE
INTEGER :: my_id, proc_num, tag, ierr
INTEGER, DIMENSION(MPI_STATUS_SIZE) :: status
REAL(KIND=8) :: time_start,time_end,time_total
INTEGER :: dest, source
INTEGER :: my_data, left_data, right_data
CALL MPI_INIT(ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD,my_id,ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,proc_num,ierr)
time_start = MPI_WTIME()
my_data = my_id+1000
! send data to left process
tag = 10
IF (my_id-1<0)THEN
dest = my_id-1+proc_num
ELSE
dest = my_id-1
END IF
IF (my_id+1>=proc_num)THEN
source = my_id+1-proc_num
ELSE
source = my_id+1
END IF
! CALL MPI_SEND(my_data,1,MPI_INTEGER,dest,tag,MPI_COMM_WORLD,ierr)
! CALL MPI_RECV(right_data,1,MPI_INTEGER,source,tag,MPI_COMM_WORLD,status,ierr)
! CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
CALL MPI_SENDRECV(my_data, 1, MPI_INTEGER, dest, tag,&
&right_data, 1, MPI_INTEGER, source, tag, MPI_COMM_WORLD,&
&status, ierr)
! send data to right process
tag = 20
IF (my_id-1<0)THEN
source = my_id-1+proc_num
ELSE
source = my_id-1
END IF
IF (my_id+1>=proc_num)THEN
dest = my_id+1-proc_num
ELSE
dest = my_id+1
END IF
! CALL MPI_SEND(my_data,1,MPI_INTEGER,dest,tag,MPI_COMM_WORLD,ierr)
! CALL MPI_RECV(left_data,1,MPI_INTEGER,source,tag,MPI_COMM_WORLD,status,ierr)
! CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
CALL MPI_SENDRECV(my_data, 1, MPI_INTEGER, dest, tag,&
&left_data, 1, MPI_INTEGER, source, tag, MPI_COMM_WORLD,&
&status, ierr)
! Write resualts
WRITE(*,*) " "
WRITE(*,*) "==========================================="
WRITE(*,10) my_id,proc_num
10 FORMAT('Process ',I2,' of ',I2)
WRITE(*,*) "my_data = ", my_data
WRITE(*,*) "left_data=", left_data, "right_data=", right_data
WRITE(*,*) "==========================================="
time_end = MPI_WTIME()
time_total = time_end-time_start
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
IF( my_id == 0 )THEN
WRITE(*,*) " "
WRITE(*,*) "==========================================="
WRITE(*,*) "Process is ", my_id
WRITE(*,*) "The massage passing is finished."
WRITE(*,*) "Total time is ", time_total
WRITE(*,*) "==========================================="
END IF
CALL MPI_FINALIZE(ierr)
END PROGRAM
积分计算PI
program main
use mpi
implicit none
real(kind=8), parameter :: PI25DT = 3.141592653589793238462643d0
real(kind=8) :: my_pi, pi, h, sum, x, f, a
REAL(KIND=8) :: time_start,time_end,time_total
integer :: n, my_id, proc_num, i, ierr
! function to integrate
f(a) = 4.d0 / (1.d0 + a*a)
call MPI_INIT(ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, my_id, ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, proc_num, ierr)
do
if (my_id == 0) then
write(*,*), "Enter the number of intervals: (0 quits) "
read(*,*) n
endif
time_start = MPI_WTIME()
! broadcast n
call MPI_BCAST(n, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
! check for quit signal
if (n <= 0) then
exit
end if
! calculate the interval size
h = 1.0d0/n
sum = 0.0d0
do i = my_id+1, n, proc_num
x = h * (dble(i) - 0.5d0)
sum = sum + f(x)
enddo
my_pi = h * sum
! collect all the partial sums
call MPI_REDUCE(my_pi, pi, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD, ierr)
time_end = MPI_WTIME()
time_total = time_end-time_start
! node 0 prints the C_elementwer.
if (my_id .eq. 0) then
write(*,*) " "
write(*,*) "==========================================="
write(*,*), "pi is ’, pi, ’ Error is", abs(pi - PI25DT)
WRITE(*,*) "Total time is", time_total
write(*,*) "==========================================="
endif
end do
call MPI_FINALIZE(ierr)
end
计算矩阵乘向量
program main
use mpi
integer :: my_id, proc_num, ierr
integer, dimension(MPI_STATUS_SIZE) :: status
integer :: m_num, n_num
real(kind=8), allocatable, dimension(:,:) :: A
real(kind=8), allocatable, dimension(:) :: B, C, buffer
real(kind=8) :: C_element
integer :: manager, worker_num
integer :: i, j, sent_num, sender, row
call MPI_INIT(ierr)
call MPI_COMM_RANK(MPI_COMM_WORLD, my_id, ierr)
call MPI_COMM_SIZE(MPI_COMM_WORLD, proc_num, ierr)
manager = 0
m_num = 100
n_num = 100
worker_num = proc_num-1
if(worker_num>m_num)then
if (my_id == manager) then
write(*,*) "Process manager : ERROR!"
write(*,*) "The worker process number is larger than colums number of vector C!"
write(*,*) "The program is forced to abort!"
CALL MPI_ABORT(MPI_COMM_WORLD,100,ierr)
end if
end if
allocate(B(n_num),buffer(n_num), STAT=ierr)
if (my_id == manager) then
allocate(A(m_num,n_num),C(m_num), STAT=ierr)
do j = 1,n_num
B(j) = 1
do i = 1,m_num
A(i,j) = i
enddo
enddo
! send b to each worker process
call MPI_BCAST(b, n_num, MPI_DOUBLE_PRECISION, manager, MPI_COMM_WORLD, ierr)
! send a row to each worker process; tag with row number
sent_num = 0
do i = 1, worker_num
buffer(:) = a(i,:)
call MPI_SEND(buffer, n_num, MPI_DOUBLE_PRECISION, i, i, MPI_COMM_WORLD, ierr)
sent_num = sent_num+1
enddo
do i = 1,m_num
call MPI_RECV(C_element, 1, MPI_DOUBLE_PRECISION, &
MPI_ANY_SOURCE, MPI_ANY_TAG, &
MPI_COMM_WORLD, status, ierr)
sender = status(MPI_SOURCE)
row = status(MPI_TAG) ! row is tag value
c(row) = C_element
if (sent_num < m_num) then ! send another row
do j = 1,n_num
buffer(j) = a(sent_num+1,j)
enddo
call MPI_SEND(buffer, n_num, MPI_DOUBLE_PRECISION, &
sender, sent_num+1, MPI_COMM_WORLD, ierr)
sent_num = sent_num+1
else ! Tell sender that there is no more work
call MPI_SEND(MPI_BOTTOM, 0, MPI_DOUBLE_PRECISION, &
sender, 0, MPI_COMM_WORLD, ierr)
endif
enddo
else
! workers receive b, then compute dot products until
! done message received
call MPI_BCAST(b, n_num, MPI_DOUBLE_PRECISION, manager, &
MPI_COMM_WORLD, ierr)
do
call MPI_RECV(buffer, n_num, MPI_DOUBLE_PRECISION, &
manager, MPI_ANY_TAG, MPI_COMM_WORLD, &
status, ierr)
if (status(MPI_TAG) == 0) exit
row = status(MPI_TAG)
C_element = 0.0
do i = 1,n_num
C_element = C_element+buffer(i)*b(i)
enddo
call MPI_SEND(C_element, 1, MPI_DOUBLE_PRECISION, manager, &
row, MPI_COMM_WORLD, ierr)
enddo
endif
if (my_id .eq. manager) then
write(*,*) "============================================"
write(*,*) "This is the manager process ", my_id
do i=1,m_num
write(*,"('C[',I3,']=',F10.1)") i, c(i)
end do
write(*,*) "============================================"
deallocate(A, C, STAT=ierr)
end if
deallocate(B, buffer, STAT=ierr)
call MPI_FINALIZE(ierr)
end
矩阵相乘
program main
use MPI
implicit none
integer :: my_id, dest_id, source_id, proc_num, ierr
integer, dimension(MPI_STATUS_SIZE) :: status
real(kind=8) :: x, y, S_local, S_global
real(kind=8), allocatable, dimension(:,:) :: A, B, B_temp, C
real(kind=8) :: time_start,time_end,time_total
integer :: N, NP, i, j, k, step
time_start = MPI_WTIME()
N = 2000
CALL MPI_INIT(ierr)
CALL MPI_COMM_RANK(MPI_COMM_WORLD,my_id,ierr)
CALL MPI_COMM_SIZE(MPI_COMM_WORLD,proc_num,ierr)
NP=N/proc_num
if (my_id == 0) then
if( mod(N,proc_num)/=0 )then
write(*,*) "Process 0: ERROR!"
write(*,*) "The remainder of the division of N by roc_num should be zero."
write(*,*) "The program is forced to abort!"
CALL MPI_ABORT(MPI_COMM_WORLD,100,ierr)
end if
end if
allocate(A(NP,N),B(N,NP),B_temp(N,NP),C(NP,N),STAT=ierr)
! Initialize A,B,C
do j=1,N
do i=1,NP
x = (dble(my_id*NP+i)-1.0d0)/(dble(N)-1.0d0)
y = (dble(j)-1.0d0)/(dble(N)-1.0d0)
A(i,j) = exp(x)*sin(3.0d0*x)
end do
end do
do j=1,Np
do i=1,N
x = (dble(i)-1.0d0)/(dble(N)-1.0d0)
y = (dble(my_id*NP+j)-1.0d0)/(dble(N)-1.0d0)
B(i,j) = ( x+cos(4.0d0*x) )*(1.0d0+y)
end do
end do
C = 0.0d0
!
do j=1,NP
do i=1,NP
do k=1,N
C(i,my_id*NP+j)=C(i,my_id*NP+j)+A(i,k)*B(k,j)
end do
end do
end do
do step=1,proc_num-1
if(my_id+step<=proc_num-1)then
dest_id = my_id+step
else
dest_id = my_id+step-proc_num
end if
if(my_id-step>=0)then
source_id = my_id-step
else
source_id = my_id-step+proc_num
end if
CALL MPI_SENDRECV(B,N*NP,MPI_DOUBLE_PRECISION,dest_id, step,&
&B_temp,N*NP,MPI_DOUBLE_PRECISION,source_id,step,&
&MPI_COMM_WORLD,status,ierr)
do j=1,NP
do i=1,NP
do k=1,N
C(i,source_id*NP+j)=C(i,source_id*NP+j)+A(i,k)*B_temp(k,j)
end do
end do
end do
CALL MPI_BARRIER(MPI_COMM_WORLD,ierr)
end do
S_local = 0.0d0
do j=1,N
do i=1,NP
S_local = S_local+C(i,j)**2
end do
end do
CALL MPI_REDUCE(S_local,S_global,1,MPI_DOUBLE_PRECISION,MPI_SUM,0,MPI_COMM_WORLD,ierr)
S_global = S_global/N**2
time_end = MPI_WTIME()
time_total = time_end-time_start
if (my_id == 0) then
write(*,*) "===================================="
write(*,*) "Process 0"
write(*,*) "S_global=", S_global
WRITE(*,*) "Total time is", time_total
write(*,*) "===================================="
end if
CALL MPI_FINALIZE(ierr)
end program
不同进程数运行所需时间
| 进程数 | 1 | 2 | 4 | 8 | 16 | 25 |
|---|---|---|---|---|---|---|
| 时间（秒） | 70.65 | 26.76 | 16.01 | 8.30 | 4.45 | 3.51 |