第3章 MPI 分布式内存编程 摘录

在消息传递程序中,运行在一个核-内存的程序通常称为进程。两个进程通过调用函数来进行通信:一个进程调用发送函数,另一个调用接收函数。

我们把使用消息传递的实现称为消息传递接口(Message - Passing Interface, MPI)。MPI定义了一个可以被程序调用的函数库。

涉及多于两个进程的‘全局’通信函数,这些函数称为集合通信。

编写消息传递程序的基本问题:数据的分割和分布式内存系统的I/O。

不让每个进程都进行I/O,而是指派其中一个进程来负责输出,而其他进程向它发送要打印的消息。

3.1 预备知识

3.1.1 编译和执行

linux系统下,编译指令是:

mpicc -g -Wall -o mpi_hello mpi_hello.cpp

执行命令:

mpiexec -n numbe_of_process ./mpi_hello

Hello_world 的MPI版本;

void mpiHello()
{
	char greeting[MAX_STRING];
	int comm_sz;
	int my_rank;

	MPI_Init(NULL, NULL);
	MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);
	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

	MPI_Status status;

	if (my_rank != 0)
	{
		sprintf(greeting, "Greetings from process %d of %d", my_rank, comm_sz);
		MPI_Send(greeting, strlen(greeting) + 1, MPI_CHAR, 0, 0, MPI_COMM_WORLD);
	}
	else
	{
		printf("Greetings from process %d of %d!\n", my_rank, comm_sz);

		for (int q = 1; q < comm_sz; q++)
		{
			
			MPI_Recv(greeting, MAX_STRING, MPI_CHAR, q, 0, MPI_COMM_WORLD, &status);
			printf("%s\n", greeting);
			
			//检测接收到数据量
			//int count = 0;
			//MPI_Get_count(&status, MPI_CHAR, &count);
			//cout << count << endl;
		}
        // 接受任何发送给0进程的消息
		//MPI_Recv(greeting, MAX_STRING, MPI_CHAR, MPI_ANY_SOURCE, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
	}

	MPI_Finalize();
}

3.1.2 MPI程序

mpi.h头文件包括了MPI函数的原型、宏定义、类型定义等,它还包括了编译MPI程序所需要的全部定义与声明。

所有MPI定义的标识符都由字符串MPI_开始。下划线后第一个字母大写,代表函数名和MPI定义的类型。MPI定义的宏和常量的所有字母都大写。

3.1.3 MPI_Init和MPI_Finalize

调用MPI_Init是为了告诉MPI系统进行所有必要的初始化设置。例如,系统需要为消息缓冲区分配存储空间,为进程指定进程号。从经验上来看,在程序调用MPI_Init之前,不应该调用其他MPI函数。

int MPI_Init(
        int*    argc_p    /*in/out*/
        char*** argv_p    /*in/out*/
);
// 当程序不使用这些参数时,可以将其设置为NULL

MPI_Init返回一个int型错误码。

int MPI_Finalize(void)

调用MPI_Finalize是为了告知系统MPI已经使用完毕,为MPI分配的任何资源可以释放了。一般而言,在调用MPI_Finalize之后,就不应该在调用任何MPI函数了。

一个典型的MPI程序框架是;

...
#include<mpi.h>
...

int main(int argc, char* argv[]){
...
/* NO MPI calls before


}

3.1.4 通信子、MPI_Comm_size和MPI_Comm_rank

在MPI中,通信子(communicator)指的是一组可以互相发送消息的进程集合。MPI_Init的其中一个目的,是在用户启动程序时,定义由用户启动的所有进程组成的通信子。这个通信子称为MPI_COMM_WORLD。

int MPI_Comm_size(MPI_Comm comm, int* comm_sz_p);

int MPI_Comm_rank(MPI_Comm comm, int* my_rank_p);

第一个参数是通信子,它所属的类型是MPI为通信子定义的特殊类型:MPI_Comm。MPI_Comm_size函数在它的第二个参数里返回通信子的进程数;MPI_Comm_rank函数在它的第二个函数里返回正在调用进程在通信子中的进程号。在MPI_COMM_WORLD中经常用参数comm_sz表示进程的数量,用参数my_rank来表示进程号。

3.1.5 SPMD程序(Single Program, Multiple Data)

编写单个程序,让不同进程产生不同动作,简单地让进程按照它们的进程号来匹配程序分支。

3.1.6 通信

int MPI_Send(
    void* msg_buf_p, 
    int msg_size, 
    MPI_Datatype msg_type, 
    int dest, 
    int tag, 
    MPI_Comm communicator);

第一个参数msg_buf_p是一个指向包含消息内容的内存块的指针。

第二个msg_size和第三个参数msg_type指定了要发送的数据量和数据类型。MPI的数据类型MPI_INT是指C++语言中的int类型,其他类型类似。

发送消息的大小必须小于或等于缓冲区的大小。

第四个参数dest指定要接收信息的进程的进程号。

第五个参数tag是个非负int类型,用于区分看上去完全一样的信息。比如用于计算或IO输出。

第六个参数communicator是通信子。所有涉及通信的MPI函数都有一个通信子参数,通信子最重要的目的之一是指定通信范围。

通信子指的是一组可用互相发送消息的进程的集合。反过来,一个通信子中的进程所发送的消息不能被另一个通信子中的进程所接收。因此,通信子这一特性可用于复杂程序。

3.1.8 MPI_Recv

int MPI_Recv(
    void* msg_buf_p, 
    int buf_size, 
    MPI_Datatype nuf_type, 
    int source, 
    int tag, 
    MPI_Comm communicator, 
    MPI_Status*, status_p);

前三个参数指定了用于接收信息的内存:msg_buf_p指向内存块,buf_size指定了内存块中要存储对象的数量,buf_type说明了对象的类型。

参数tag要与发送消息的参数tag相匹配,参数communicator必须域发送进程所用的通信子相匹配。简短介绍statue_p,在大部分情况下,调用函数不使用该参数,使用MPI_STATUS_IGNORE就行。

3.1.9 消息匹配

假定q号进程调用了MPI_Send函数:

MPI_Send(send_buf_p, send_buf_size, send_buf_type, dest, tag, send_comm);

并且假定r号进程调用了MPI_Recv函数;

MPI_Recv(recv_buf_p, recv_buf_size, recv_buf_type, source, RECV_tag, recv_comm, MPI_STATUS_IGNORE);

若q号进程发送的消息要被r号进程接收,则必须满足;

recvc_comm= send_comm;

recv_tag = send_tag;

dest = r;

src = r;

并且send_buf_p / recv_buf_p, send_buf_size / recv_buf_size, send_type / recv_type必须指定兼容的缓冲区。一般需要满足;

recv_type = send_type;

recv_buf_size >= send_buf_size;

那么由q号进程发送的消息就可用被r号进程成功地接收。

一个进程可接收多个进程发送的消息,接收进程并不知道其他进程发送消息的顺序。如果分配给每个进程的工作所要耗费的时间是无法预测的,那么0号进程就无法知道其他进程完成工作的顺序。如果只是简单的按照进程号顺序接收结果(for 循环)。但是,如果comm_s - 1 号进程先完成工作,但是他必须等待其他进程完成。

为避免该问题,MPI提供一个特殊的的常量MPI_ANY_SOURCE,可用传递给MPI_Recv,0号进程可以按照进程完成工作的顺序来接收结果了。

类似地,一个进程也有可能接收多条来自另一个进程的有着不同标签的消息,并且接收进程并不知道消息发送的顺序。MPI提供常量MPI_ANY_TAG传递给MPI_Recv的参数tag。

当使用上述MPI_ANY_?的通配符时,有几点需要强调
        1. 只有接收者可以使用通配符。发送者必须指定一个进程与一个非负整数标签。因此,MPI使用所谓的“推(push)”通信机制,而不是“拉(pull)”通信机制。
        2. 通信子参数没有通配符。发送者与接收者都必须指定通信子。

3.1.10 status_p参数

接收者可以在不知道以下信息的情况下接收信息:
        1. 消息中的数据量;
        2. 消息的发送者;
        3. 消息的标签;

那么接收者是如何找出这些值的? MPI_Recv的最后一个参数类型为MPI_Status*类型。MPI_Status是一个至少有3个成员的结构,MPI_Source,MPI_Tag,MPI_Error。

MPI_Status status, 通过将&status传递给MPI_Recv函数并调用它后,可以通过检查以下两个成员来确定发送者和标签:status.MPI_Source与status.MPI_Tag。

MPI_Get_count(MPI_Status* status, MPI_Datatype recv_type, int* count);

使用MPI_Get_count()找回接收数据量,会返回count参数接收到的元素数量。count值并不能简单地作为MPI_Status变量地成员直接访问,因为它取决于接收数据的类型。如果不是很必要获得count值,没有必要为得到该值而浪费一次计算。

3.1.11 MPI_Send和MPI_Recv的语义

发送进程组装消息,例如,它为实际要发送的数据添加”信封“信息。”信封“信息包括目标进程的进程号、发送进程的进程号、标签、通信子,以及数据类型和数量等信息。一旦消息组装完毕,有两种可能性:发送进程可以缓冲消息,也可以阻塞(block)。如果它缓冲消息,则MPI系统将会把消息(数据和信封)放置在它自己的内部存储器里,并返回MPI_Send的调用。

另一方面,如果系统发送阻塞,那么它将一直等待,直到可以开始发送消息,并不立即返回对MPI_Send的调用。因此,如果使用MPI_Send,当函数返回时,实际上并不知道消息是否以及发送出去。

MPI_Send的精确行为是由MPI实现所决定的。典型的实现方法有一个默认的消息”截止“大小(cutoff message size)。如果一条消息的大小小于“截止”大小,它将被缓冲;如果大于截止大小,那么MPI_Send函数将被阻塞。

MPI_Recv函数总是阻塞的,直到接收一条匹配的消息。因此,当MPI_Recv函数调用返回时,就知道有一条消息以及存储在接收缓冲区了(除非产生了错误)。

MPI要求消息是不可超越的(nonovertaking)。即如果q号进程发送了两条消息给r号进程,那么q号进程发送的第一条消息必须在第二条消息之前。但是,如果消息是来自于不同进程时,消息的到达顺序是没有限制的。

3.1.12 潜在的陷阱

MPI_Recv的语义会导致MPI编程中的一个潜在陷阱。如果一个进程试图接收消息,但没有相匹配的消息,那么该进程将会永远阻塞在哪里,即进程悬挂。因此,在设计程序时,需要保证每条接收都有一条相匹配的发送。注意进程号与标签匹配。

简单地说,如果调用MPI_Send/MPI_Recv发生了阻塞,并且没有相匹配地接收/发送,那么发送/接收进程就悬挂起来。

3.2 用MPI来实现梯形积分法

3.2.1 梯形积分

h = (b - a)/ n;

3.2.2 并行梯形积分法

可以用四个基本步骤来设计一个并行程序:
        1. 将问题的解决方案划分多个任务;
        2. 在任务间识别需要的通信通道;
        3. 在任务聚合成复合任务;
        4. 在核上分配复合任务;

对于梯形积分法,可以识别出两种任务:一种是获取单个矩形区域的面积,另一个是计算这些面积之和。然后可以通过通信通道将这些每个第一种任务与第二种任务相连接。

如何将聚合任务并将其分配到核? 梯形的数目越多,计算精度越高。因此小梯形的数目将超过核的数量,需要将梯形区域面积的计算聚合成组。一个很自然的方法是将区间[a, b]分成comm_sz 个梯形核所有comm_sz个子空间上应用梯形积分法。最后,利用0号进程将这些梯形面积的估计值都加起来,伪代码如下:

Get a, b, n;
h = (b - a) / n;
local_n = n / comm_sz;
local_b = local_a + my_rank*local_n*h;
local_intergral = Trap(local_a, local_b, local_n, h);

if (my_rank != 0)
    Send local_integral to process 0;
else /* my_rank */
    total_integral = local_integral;
    for (proc = 1; proc < comm_sz; pro++) {
        Receive local_integral from proc;
        total_integral += local_integral;
    }
}
if (my_rank == 0)
    print result;

在MPI编程中,局部变量只在使用它们的进程内有效。全局变量在所有的进程中有效。

3.3 I/O处理

3.3.1 输出

虽然MPI标准没有指定那些进程可以访问那些I/O设备,但是几乎所有的MPI实现都允许MPI_COMM_WORLD里的所有进程都能访问标准输出stdout和标准错误输出stderr,所以,大部分的MPI实现都允许所有进程执行printf和fprintf。

但是,大部分的MPI实现并不提供对这些I/O设备访问的自动调动。如果多个进程试图写标准输出stdout,那么这些进程的输出顺序是无法预测的,甚至会发生一个进程的输出被另一个进程的输出打断。

第一版: 固定a,b,n

double trapeArea(double left_endpt, double right_endpt, int trap_count, double base_len)
{
	double estimate, x;
	
	estimate = (func(left_endpt) + func(right_endpt)) / 2;
	
	for (int i = 0; i < trap_count; i++)
	{
		x = left_endpt + i * base_len;
		estimate += func(x);
	}
	
	estimate *= base_len;

	return estimate;
}

void mpiTrapeAreaIntegra()
{
	int my_rank, comm_sz, n = 1024, local_n;
	double a = 0, b = 3, h, local_a, local_b;
	double local_int, total_int;

	MPI_Init(NULL, NULL);
	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
	MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);


	h = (b - a) / n;
	local_n = n / comm_sz;

	local_a = a + local_n * h * my_rank;
	local_b = local_a + local_n * h;
	local_int = trapeArea(local_a, local_b, local_n, h);


	if (my_rank != 0)
	{

		MPI_Send(&local_int, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	}
	else
	{
		total_int = local_int;	
		for (int i = 1; i < comm_sz; i++)
		{
			MPI_Recv(&local_int, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
			total_int += local_int;
		}
		//MPI_Recv(&)
	}

	if (my_rank == 0)
	{
		printf("the integral of area is %f\n", total_int);
	}

	MPI_Finalize();

}

3.3.2 输入a,b,n

与输出不同,大部分的MPI实现只允许MPI_COMM_WORLD中的0号进程访问标准输入stdin。因为如果多个进程都能访问输入,那么哪个进程应该得到输入数据的哪个部分呢?

第二版:输入a,b,n

void getInput(int my_rank, int comm_sz, double& a, double& b, int& n)
{
	int src = 0;

	if (my_rank == 0)
	{
		cout << "Enter a, b, and n\n"; //如果使用printf函数,貌似先运行下面三行的cin,然后再是printf
		cin >> a;
		cin >> b;
		cin >> n;
		// 由于getInput函数中a,b, n是控制台输入的,所以只对0号进程可见,于是需要将a,b,n输入值发送给其他进程。否则其他进程中的a,b,n为声明变量时的值。
		for (int dest = 1; dest < comm_sz; dest++)
		{
			MPI_Send(&a, 1, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD);
			MPI_Send(&b, 1, MPI_DOUBLE, dest, 0, MPI_COMM_WORLD);
			MPI_Send(&n, 1, MPI_INT, dest, 0, MPI_COMM_WORLD);
		}
	}
	else
	{
		MPI_Recv(&a, 1, MPI_DOUBLE, src, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		MPI_Recv(&b, 1, MPI_DOUBLE, src, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		MPI_Recv(&n, 1, MPI_INT, src, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
	}
}

void mpiTrapeAreaIntegraGetInput()
{
	int my_rank, comm_sz, n = 0, local_n;
	double a = 0, b = 0, h, local_a, local_b;
	double local_int, total_int;

	MPI_Init(NULL, NULL);
	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
	MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

	getInput(my_rank, comm_sz, a, b, n);
	h = (b - a) / n;
	local_n = n / comm_sz;

	local_a = a + local_n * h * my_rank;
	local_b = local_a + local_n * h;
	local_int = trapeArea(local_a, local_b, local_n, h);

	if (my_rank != 0)
	{
		MPI_Send(&local_int, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	}
	else
	{
		total_int = local_int;	
		for (int i = 1; i < comm_sz; i++)
		{
			MPI_Recv(&local_int, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
			total_int += local_int;
		}
		//MPI_Recv(&)
	}

	if (my_rank == 0)
	{
		printf("the integral of area is %f\n", total_int);
	}

	MPI_Finalize();

}

3.4 集合通信

3.4.1 树形结构通信

 原有方案需要0号进程接收comm_sz - 1 次,并做7次加法。而树形结构通信0号进程只需要3次接收和3次加法,并且其他进程所作的接收与加法都不超过2次。

第三版:广播输入a, b,n(MPI_Bcast)

void getInputBroadcast(int my_rank, int comm_sz, double& a, double& b, int& n)
{
	int src = 0;

	MPI_Datatype input_mpi_t;
	mpiBuildMyDataType(&a, &b, &n, &input_mpi_t);

	if (my_rank == 0)
	{
		cout << "Enter a, b, and n\n"; //如果使用printf函数,貌似先运行下面三行的cin,然后再是printf
		cin >> a;
		cin >> b;
		cin >> n;
	}


	// 由于getInput函数中a,b, n是控制台输入的,所以只对0号进程可见,于是需要将a,b,n输入值发送给其他进程。否则其他进程中的a,b,n为声明变量时的值。
	MPI_Bcast(&a, 1, MPI_DOUBLE, src, MPI_COMM_WORLD);
	MPI_Bcast(&b, 1, MPI_DOUBLE, src, MPI_COMM_WORLD);
	MPI_Bcast(&n, 1, MPI_INT, src, MPI_COMM_WORLD);

	//MPI_Bcast(&a, 1, input_mpi_t, src, MPI_COMM_WORLD);
	MPI_Type_free(&input_mpi_t);
}

void mpiTrapeAreaIntegraGetInput()
{
	int my_rank, comm_sz, n = 0, local_n;
	double a = 0, b = 0, h, local_a, local_b;
	double local_int, total_int;

	MPI_Init(NULL, NULL);
	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
	MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

	getInput(my_rank, comm_sz, a, b, n);
	h = (b - a) / n;
	local_n = n / comm_sz;

	local_a = a + local_n * h * my_rank;
	local_b = local_a + local_n * h;
	local_int = trapeArea(local_a, local_b, local_n, h);

	if (my_rank != 0)
	{
		MPI_Send(&local_int, 1, MPI_DOUBLE, 0, 0, MPI_COMM_WORLD);
	}
	else
	{
		total_int = local_int;	
		for (int i = 1; i < comm_sz; i++)
		{
			MPI_Recv(&local_int, 1, MPI_DOUBLE, i, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
			total_int += local_int;
		}
		//MPI_Recv(&)
	}

	if (my_rank == 0)
	{
		printf("the integral of area is %f\n", total_int);
	}

	MPI_Finalize();

}

3.4.2 MPI_Reduce

在MPI里,涉及通信子中所有进程的通信函数被称为集合通信(collective communication)。

int MPI_Reduce(
    void* input_data_p, 
    void* output_data_p, 
    int count, 
    MPI_Data_type datatype, 
    MPI_Op operator, 
    MPI_Comm comm);

如果count参数大于1,那么MPI_Reduce函数可以应用到数组上。

operator类型是MPI_Op,是一个预定义的MPI类型。其预定义值如下表所示:

运算符值含义运算符值含义
MPI_MAX求最大值MPI_LOR逻辑或
MPI_MIN求最小值MPI_BOR按位或
MPI_SUM求累加和MPI_LXOR逻辑异或
MPI_PROD求累乘积MPI_BXOR按位异或
MPI_LAND逻辑与MPI_MAXLOC求最大值位置
MPI_BAND按位与MPI_MINLOC求最小值位置

第四版:归约梯形和(MPI_Reduce)

void mpiTrapeAreaIntegraMpiReduce()
{
	int my_rank, comm_sz, n = 0, local_n;
	double a = 0, b = 0, h, local_a, local_b;
	double local_int, total_int;

	MPI_Init(NULL, NULL);
	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
	MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

	getInput(my_rank, comm_sz, a, b, n);
	h = (b - a) / n;
	local_n = n / comm_sz;

	local_a = a + local_n * h * my_rank;
	local_b = local_a + local_n * h;
	local_int = trapeArea(local_a, local_b, local_n, h);

	MPI_Reduce(&local_int, &total_int, 1, MPI_DOUBLE, MPI_SUM, 0, MPI_COMM_WORLD);

	if (my_rank == 0)
	{
		printf("the integral of area is %f\n", total_int);
	}

	MPI_Finalize();

}

3.4.3 集合通信与点对点通信

集合通信与点对点通信在多个方面是不同的:
        1. 在通信子中的所有进程都必须调用相同的集合通信函数。例如,试图将一个进程中的MPI_Reduce与另一个进程的MPI_Recv调用相匹配的程序会出错,此时程序会被悬挂或者崩溃。
        2. 每个进程传递给MPI集合通信函数的参数必须是“相容的”。例如,一个进程将0作为dest_process的值传递给参数,而另一个传递的是1,那么对MPI_Reduce调用所产生的结果就是错误的,程序有可能被悬挂或者崩溃。
        3.参数output_data_p只用在dest_process上。然而,所有进程仍需要传递一个与output_data_p相对应的实际参数,尽管它的值只是NULL。
        4.点对点通信是通过标签和通信子来匹配的。集合通信不使用标签,只通过通信子和调用顺序来匹配。内存单元的名字与MPI_Reduce的调用顺序无关,函数调用顺序决定了匹配方式。

最后一个忠告:我们也许会冒风险用同一个缓冲区同时作为输入和输出的MPI_Reduce。

NPI_Reduce(&x, &x, 1, MPI_DOUBLE, MPI_SUN, 0, comm);

程序会挂起或奔溃的原因在于涉及输出参数的别名。两个参数如果指向同一块内存,它们之间存在别名问题。MPI禁止输入或输出参数作为其他参数的别名。

3.4.4 MPI_Allreduce

不难想象这样一种情况,即所有的进程都想到全局求和的结果,以便可以完成一个更大规模的计算。可以通过颠倒整颗树来发布全局和,MPI_Allreduce可以实现该功能。还有另一种方法是蝶形通信。

int MPI_Allreduce(
    void* input_data_p, 
    void* output_data_p, 
    int count, 
    MPI_Datatype datatype, 
    MPI_Op operator, 
    MPI_Comm comm);

与MPI_Reduce参数区别是没有dest_process。

树型结构的全局求和

 

蝶形结构的全局求和

3.4.5 广播

在一个集合通信中,如果属于一个进程的数据被发送到通信子中的所有进程,这样的集合通信就叫做广播(broadcast)。

int MPI_Bcast(
    void* data_p, 
    int count, 
    MPI_Datatype datatype, 
    int source_proc, 
    MPI_Comm comm);

当集合通信函数中的某个参数被标记位输入/输出时,意味着可能某些进程中它是输入参数,某些进程是一个输出参数。

3.4.6 数据分发MPI_Scatter

向量求和案例:计算工作由向量的各个分量求和组成,可能只是指定各个任务求和对应分量。各个任务之间没有通信,向量的并行加法问题就归结为聚合任务以及将它们分配到核上。

将数据分发的三种方式:块划分、循环划分以及块-循环划分

对12个分量的向量划分如下:

进程分量
块划分循环划分块-循环划分
0012303690167
14567147102389
289101125811451011

 可以试着0号进程读入整个向量,但只将分量发送给需要分量的其他进程。MPI_Scatter实现。

int MPI_Scatter(
    void* send_buf_p, 
    int send_count, 
    MPI_Datatype, datatype, 
    void* recv_buf_p, 
    int recv_count, 
    MPI_Datatype datatype, 
    int src_proc, 
    MPI_Comm comm);
void mpiVectorRead(int local_n, int n, double* vector_local, int my_rank, MPI_Comm comm)
{
	double* vector_global = NULL;

	if (my_rank == 0)
	{
		vector_global = new double[n];

		for (int i = 0; i < n; ++i)
		{
			vector_global[i] = static_cast<double>(i);
		}

		MPI_Scatter(vector_global, local_n, MPI_DOUBLE, vector_local, local_n, MPI_DOUBLE, 0, comm);
		delete[] vector_global;
	}
	else
	{
		MPI_Scatter(vector_global, local_n, MPI_DOUBLE, vector_local, local_n, MPI_DOUBLE, 0, comm);
	}
}

void mpiVectorSum(int local_n, double* local_x, double* local_y, double* local_z)
{
	for (int local_i = 0; local_i < local_n; ++local_i)
	{
		local_z[local_i] = local_x[local_i] + local_y[local_i];
	}
}

void mpiVectorPrint(int local_n, int n, double* local_c, int my_rank, MPI_Comm comm)
{
	double* c = NULL;

	if (my_rank == 0)
	{
		c = new double[n];
		MPI_Gather(local_c, local_n, MPI_DOUBLE, c, local_n, MPI_DOUBLE, 0, comm);
		
		cout << "打印输出\n";
		for (int i = 0; i < n; ++i)
		{
			cout << c[i] << endl;
		}
	} 

	else
	{
		MPI_Gather(local_c, local_n, MPI_DOUBLE, c, local_n, MPI_DOUBLE, 0, comm);
	}

	delete[] c;
}

void mpiParallelVectorSum()
{
	int comm_sz, my_rank;
	double* local_a = NULL, * local_b = NULL, * local_c = NULL;
	int n = 0;

	MPI_Init(NULL, NULL);
	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
	MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

	if (my_rank == 0)
	{
		cout << "请输入向量长度\n";
		cin >> n;
	}

	MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);

	int local_n = n / comm_sz;
	local_a = new double[local_n];
	local_b = new double[local_n];
	local_c = new double[local_n];

	mpiVectorRead(local_n, n, local_a, my_rank, MPI_COMM_WORLD);
	mpiVectorRead(local_n, n, local_b, my_rank, MPI_COMM_WORLD);

	mpiVectorSum(local_n, local_a, local_b, local_c);

	mpiVectorPrint(local_n, n, local_c, my_rank, MPI_COMM_WORLD);

	delete[] local_a;
	delete[] local_b;
	delete[] local_c;

	MPI_Finalize();

}

令人惊讶的是send_count 也应该是local_n,因为send_count参数表示的是发送到每个进程的数据量,而不是send_buf_p所指的内存的数据量。需要注意,MPI_Scatter将第一个数据块发送给0号进程,第二块数据发送给1号进程,以此类推。所以MPI_Scatter适用于块划分法。

3.4.8 聚集MPI_Gather

若想看见向量加法的结果,我们需要将向量的所有分量都收集到0号进程。

int MPI_Gather(
    void* send_buf_p, 
    int send_count, 
    MPI_Datatype send_type, 
    void* recv_buf_p, 
    int recv_count, 
    MPI_Datatype recv_type, 
    int dest_proc, 
    MPI_Comm comm);

注意,recv_count是0号进程接收每个进程单独发送到的数据量,而不是所有接收到的数据量的总和。MPI_Gather也是只能应用在块划分中。

3.4.9 全局通信

int MPI_Allgather(
    void* send_buf_p, 
    int send_count, 
    MPI_Datatype send_type, 
    void* recv_buf_p, 
    int recv_count, 
    MPI_Datatype datatype,
    MPI_Comm comm);

MPI_Allgather函数将每个进程的send_buf_p内容串联起来,存储在每个进程的recv_buf_p参数。通常情况下,recv_count与send_count的值相同。

void mpiMatrixRead(int local_n, int n, int m, double* mat_local, int my_rank, MPI_Comm comm)
{
	double* mat_global = NULL;

	if (my_rank == 0)
	{
		mat_global = new double[n * m];

		for (int i = 0; i < n; ++i)
		{
			for (int j = 0; j < m; ++j)
			{
				int k = j + i * m;
				mat_global[k] = 1.;
			}
		}

		MPI_Scatter(mat_global, local_n * m, MPI_DOUBLE, mat_local, local_n * m, MPI_DOUBLE, 0, comm);
		delete[] mat_global;
	}
	else
	{
		MPI_Scatter(mat_global, local_n * m, MPI_DOUBLE, mat_local, local_n * m, MPI_DOUBLE, 0, comm);
	}
}

void mpiMatDMVect(int local_n, int m, double* x, double* mat_local, double* y_local, int my_rank, MPI_Comm comm)
{
	for (int i = 0; i < local_n; ++i)
	{
		y_local[i] = 0;
		for (int j = 0; j < m; ++j)
		{
			y_local[i] += mat_local[i * m + j] * x[j];
		}
	}
}

void mpiParallelMatDotMulVect()
{
	int comm_sz, my_rank;
	double* mat_local = NULL, * x = NULL, * y_local = NULL;
	int n = 0, m = 0, local_n;

	MPI_Init(NULL, NULL);
	MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);
	MPI_Comm_size(MPI_COMM_WORLD, &comm_sz);

	if (my_rank == 0)
	{
		cout << "请输入矩阵n,m\n";
		cin >> n;
		cin >> m;
	}

	MPI_Bcast(&m, 1, MPI_INT, 0, MPI_COMM_WORLD);
	MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);

	// 矩阵与向量分发
	local_n = n / comm_sz;
	mat_local = new double[local_n * m];
	x = new double[m];
	y_local = new double[local_n];

	if (my_rank == 0)
	{
		for (int j = 0; j < m; ++j)
		{
			x[j] = static_cast<double>(j);
		}
	}

	MPI_Bcast(x, m, MPI_DOUBLE, 0, MPI_COMM_WORLD);

	mpiMatrixRead(local_n, n, m, mat_local, my_rank, MPI_COMM_WORLD);

	// 矩阵与向量相乘
	mpiMatDMVect(local_n, m, x, mat_local, y_local, my_rank, MPI_COMM_WORLD);

	// 向量收集
	mpiVectorPrint(local_n, n, y_local, my_rank, MPI_COMM_WORLD);

	delete[] mat_local;
	delete[] x;
	delete[] y_local;

	MPI_Finalize();

}

3.5 MPI的派生数据类型

在几乎所有的分布式内存系统中,通信比本地计算开销大很多。

如果减少发送的消息数量,就能提供程序的性能。

MPI提供了3个方法来整合多条数据消息:不同通信函数的count参数,派生数据类型,以及MPI_Pack/Pack函数。

在MPI中,通过同时存储数据项的类型以及它们在内存中的相对位置,派生数据类型可以用于表示内存中数据项的任意集合。其主要思想是:如果发送数据的函数知道数据项的类型以及在内存中数据项集合的相对位置,就可以在数据项被发送出去之前在内存中将数据项聚集起来。

一个派生数据类型是由一系列的MPI的基本数据类型和每个数据类型的偏移所组成的。

派生数据函数:

int MPI_Type_creat_struct(
    int count, 
    int array_of_blocklengths[], 
    MPI_Aint array_of_displacements[], 
    MPI_Datatype array_of_types[], 
    MPI_Datatype new_type_p);

参数count指的是数据类型中元素的个数,在此例中,为3.

参数array_of_blocklengths允许单独的数据项可能是一个数组或者子数组。在此例为{1, 1, 1}.

参数array_ofdisplacements指定了距离消息起始位置的偏移量。在此例为{0, 16, 24}。未找到位置,可以用MPI_Get_address函数。

int MPI_Get_address(void* location_p, MPI_Aint* address_p);
MPI_Aint a_addr, b_addr, c_addr;

MPI_Get_address(&a, &a_addr);
array_of_displacements[0] = 0;
MPI_Get_address(&b, &b_addr);
array_of_displacements[1] = b_addr - a_addr;
MPI_Get_address(&c, &c_addr);
array_of_displacements[2] = c_addr - b_addr;

MPI_Get_address返回的是location_p所指向的内存单元的地址。

完成以上工作后,可以通过以下代码建立新的数据类型;

MPI_Datatype input_mpi_t;

...
MPI_Type_struct(3, array_of_blocklengths, array_of_displacements, array_of_types, &input_mpi_t);

在使用通信函数中的input_mpi_t之前,必须先用一个MPI_Type_commit函数调用去指定它:

int MPI_Type_commit(MPI_Datatype* new_mpi_t_p);

在构造新数据类型过程中,MPI实现可能要在内部分配额外的存储空间。因此,当我们使用新的数据时,可以用一个函数调用区释放额外的存储空间。

int MPI_Type_free(MPI_Datatype* old_mpi_t_p);

假设在梯形积分法的例子中,a,b,n在内存的位置如下所示:

变量地址
a24
b40
n48

派生的数据类型 {(MPI_DOUBLE, 0),(MPI_DOUBLE, 16),(MPI_INT, 24)}

MPI_Pack每次复制一块要发送的数据到用户提供的缓冲区,该缓冲区既可以接收数据也可以发送数据。当接收数据后,可以使用MPI_Unpack将接收缓冲区的数据解包。

int MPI_Pack(
    void* in_buf, 
    int in_buf_count, 
    MPI_Datatype datatype, 
    void* pack_buf, 
    int pack_buf_sz, 
    int* position_p, 
    MPI_Comm comm);

可以用下面代码来打包梯形积分法程序的输入数据;

char pack_buf[100];
int position = 0;

// 打包
if (my_rank == 0 )
{
    MPI_Pack(&a, 1, MPI_DOUBLE, pack_buf, 100, &position, comm);
    MPI_Pack(&b, 1, MPI_DOUBLE, pack_buf, 100, &position, comm);
    MPI_Pack(&n, 1, MPI_INT, pack_buf, 100, &position, comm);
}

// 广播
MPI_Bcast(pack_buf, 100, MPI_PACKED, 0, comm);

char unpack_buf[100];

// 解压
MPI_Unpack(pack_buf, 100, &position, unpack_buf, MPI_PACKED, comm);

关键是position参数,当调用MPI_Pack时,该参数应该指向pack_buf中第一个可访问元素的位置;当MPI_Pack返回时,该参数指向pack_buf在数据打包后第一个可以被访问的位置。

int MPI_Unpack(
    void* pack_buf, 
    int pack_buf_sz, 
    int* position_p, 
    void* out_buf, 
    int out_buf_count, 
    MPI_Datatype datatype, 
    MPI_Comm comm);

函数MPI_Unpack从position = 0的位置一块块解压数据。

3.6 MPI程序的性能评估

3.6.1 计时

MPI_Wtime,返回从过去某一个时刻开始所经过的秒数;

double MPI_Wtime(void);

计时代码如下:

double start, finish;
...
start = MPI_Wtime();
/* code to be timed */
...
finish = MPI_Wtime();
printf("%d Proc > Elasped time : %f", my_rank, finish - start);

C语言中clock()函数返回的是CPU时间(包括用户代码、库函数以及系统调用函数所消耗的时间),但它并不包括空闲时间。

并行时间取决于”最慢“进程花费时间。计时时应该让所有进程都处于同一时间起点,MPI_Barrier函数实现该要求。

int MPI_Barrier(MPI_Comm comm);

MPI_Wtime返回墙上时间,包括空闲等待时间。

3.6.2 结果

当增加进程数时,运行时间会在一个阶段内减少。然而,达到某点后运行时间后会开始变得很慢。

并行程序会将串行程序的工作分配到各个进程中,但又会增加额外的开销,在MPI程序中,并行计算的开销一般来自于通信,它还同时受到问题集的规模核进程数的影响。

3.6.4 可扩展性

若程序可以在不增加问题规模的前提下维持稳定效率,那么此程序是强可扩展性的;当问题规模增加时,通过增加进程数来维持程序效率的,称为弱可扩展性。

3.7 并行排序算法

设总共有n个键值,p = comm_sz个进程,给每个进程分配n/p个键值,假设能够除断。开始是乱序的。然而在算法结束后要求:
        1. 每个进程上的键值应该以升序的方式存储;
        2. 若0<=q<r<p,则分配的q的每一个键值都应该小于等于分配给r进程的每一个键值。

3.7.1 简单的串行排序算法

void Bubble_sort(int a[], int n)
{
    int list_length, i, temp;
    for (list_length = ; list_length >= 2; list_length--)
        for (i = 0; i < list_length; i++)
        {
            if (a[i] > a[i + 1])
            {
                temp = a[n];
                a[n] = a[i + 1];
                a[n + 1] = temp;
            }
        }

}

冒泡排序的一个变种是奇偶交换排序。关键在于去耦的比较-交换。

该算法由一系列阶段组成,阶段分2种类型。

在偶数阶段,比较-交换由下数对执行:

        (a[0], a[1]),(a[2], a[3]),(a[4], a[5])...

在奇数阶段,比较-交换由下数对执行:

        (a[1], a[2]),(a[3], a[4]),(a[5], a[6])...

定理:设A是由n个键值的列表,作为奇偶交换排序算法的输入,那么经过n个阶段后,A能够排好。

串行的奇偶交换排序算法实现:

void Odd_even_sort(int a[], int n)
{
    int phase, i, temp;
    for (phase = 0; phase < n; phase++)
    {  
        if (phase % 2 == 0)
        {
            for (i = 0; i < n; i+=2)
            {
               if (a[i] > a[i + 1])
                {
                   temp = a[n];
                   a[n] = a[i + 1];
                   a[n + 1] = temp;
                 }
             }
        }
        else
        {
            for (i = 1; i < n - 1; i+=2)
            {
               if (a[i] > a[i + 1])
                {
                   temp = a[n];
                   a[n] = a[i + 1];
                   a[n + 1] = temp;
                 }
             }
        }
    }
}

3.7.2 并行奇偶交换排序算法

根据foster方法;

        任务:在阶段j结束时确定a[i]的值;

        通信:确定a[i]值的任务需要与其他a[i+1]或a[i-1]的任务进行通信,同时,在阶段j结束时,a[i]的值需要用来在阶段j+1结束时确定a[i]的值。

一般而言不会有p等于n的情况,不可能有几千上百万核。另外,即使有这样的处理器,每一次比较-交换时,发送和接收信息导致的额外时间开销会严重降低程序的运行时间。记住通信开销远比“局部”计算开销大。

所以给每个进程n/p个键值,让每个进程先进行内部排序,然后在进行进程间通信-比较排序。

时间进程
0123
开始15, 11, 9, 163, 14, 8, 74, 6, 12, 105, 2, 13, 1
局部排序后9, 11, 15, 163, 7, 8, 144, 6, 10, 121, 2, 5, 13
阶段0后3, 7, 8, 911, 14, 15, 161, 2, 4, 56, 10, 12, 13
阶段1后3, 7, 8, 91, 2, 4, 511, 14, 15, 166, 10, 12, 13
阶段2后1, 2, 3, 45, 7, 8, 96, 10, 11, 1213, 14, 15, 16
阶段3后1, 2, 3, 45, 6, 7, 89, 10, 11, 1213, 14, 15, 16

 伪代码如下所示:

Sort local keys:
for (phase = 0; phase < comm_sz; phase++){
    partner = Compute_partner(phase, my_rank);
    if (I'm not idle){
    Send my keys to partner;
    Receive keys from partner;
    if (my_rank < partner)
        Keep smaller keys;
    else
        Keep larger keys;
    }
}

MPI_PROC_NULL是MPI库定义的一个常量。在点对点通信中,将它作为源进程或目标进程的进程号,此时,调用函数后会直接返回,不会产生任何通信。

3.7.3 MPI程序的安全性

许多MPI函数都设置了使系统从缓冲到阻塞间切换的阈值,即相对较小的消息就交由MPI_Send缓冲,但对于大型数据就选择阻塞模式。如果每个进程都阻塞在MPI_Send上,则没有进程会去调用MPI_Recv,此时程序就会死锁或挂起,每个进程都在等待一个不会发生的事件发生。

依赖MPI提供的缓冲机制是不安全的,可能在运行其他输入集时导致程序奔溃或挂起,引出一些问题:
        1. 什么样的程序是安全的?
        2. 怎末修改并行奇偶交换排序程序的通信过程,使其安全?

针对第一个问题,MPI标准提供了MPI_Ssend函数。额外的s代表同步(synchronize),MPI_Ssend保证了直到对应的接收开始前,发送端一直阻塞。将MPI_Ssend替换MPI_Send来检查程序是否安全。

int MPI_Ssend(void* msg_buf_p, int count, MPI_Datatype datatype, int dest, int tag, MPI_Comm comm);

针对第二个问题,解决方案是重构通信。造成程序不安全的最大因素在于多个进程先同时发送消息,在同时接收消息。奇偶排序配对进程之间的数据就是一个例子。另一个是“环状”传递,每个进程q向q+1发送消息,进程comm_sz - 1向0号进程发送消息。

MPI_Send(msg, size, MPI_INT, (my_rank + 1) % my_rank, 0, omm);
MPI_Recv(new_msg, size, MPI_INT, (my_rank + comm_sz - 1) % comm_sz, 0, omm);

重构这两个通信函数,使一些进程先接收消息在发送消息,代码如下:

if (my_rank % 2 == 0)
{
    MPI_Send(msg, size, MPI_INT, (my_rank + 1) % my_rank, 0, omm);
    MPI_Recv(new_msg, size, MPI_INT, (my_rank + comm_sz - 1) % comm_sz, 0, omm);
}
else
{
    MPI_Recv(new_msg, size, MPI_INT, (my_rank + comm_sz - 1) % comm_sz, 0, omm);
    MPI_Send(msg, size, MPI_INT, (my_rank + 1) % my_rank, 0, omm);
}

如果comm_sz是偶数,这样改的效果很好。

MPI提供了自己的调度通信的方法,该函数为MPI_Sendrecv;调用一次该函数,它会分别执行一次阻塞式消息发送和一次消息接收。

int MPI_Sendrecv(
    void* send_buf_p, 
    int send_buf_size, 
    MPI_Datatype send_datatype, 
    int dest, 
    int send_tag, 
    void* recv_buf_p, 
    int recv_size,
    MPI_Datatype recv_datatype, 
    int source, 
    int recv_tag, 
    MPI_Comm comm, 
    MPI_Status* status_p);

如果发送和接收使用的是同一缓冲区,MPI库还提供了MPI_Sendrecv_replace函数,

int MPI_Sendrecv_replace(
    void* buf_p,
    int buf_size,
    MPI_Datatype buf_type,
    int dest,
    int send_tag,
    int source,
    int recv_tag,
    MPI_Comm comm,
    MPI_Status* status_p
)

然后就是保留最小n/p个键值了,代码如下:

void Merge_low(int my_keys[], int recv_keys[], int temp_keys[], int local_n)
{
    int m_i = r_i = t_i = 0;
    while (t_i < local_n)
    {
        if (my_keys[m_i] < recv_keys[r_i])
        {
            temp_keys[t_i] = my_keys[m_i];
            t_i++, m_i++;
        }
        else
        {
            temp_keys[t_i] = recv_keys[r_i];
            t_i++, r_i++;
        }
    }
    
    for (m_i = 0; m_i < local_n; m_i++)
        my_keys[m_i] = temp_keys[m_i];
}

3.8 小结

消息传递接口(MPI),是一个可以被C,C++,Fortan程序调用的函数库。许多系统使用mpicc编译MPI程序,并通过mpiexec运行它们。C程序需要包含mpi.h头文件才能使用MPI库定义的函数与宏。

MPI_Init函数建立MPI程序,它应该最先被调用。若程序不使用argc和argv,可以直接传入NULL。用完MPI后,记得调用MPI_Finalize函数结束程序。

在MPI中,一个通信子是一组进程的集合,该集合中的进程之间可以相互发送消息。MPI程序启动后,MPI创建由所有进程组成的通信子,称为MPI_COMM_WORLD。

许多并行程序使用单程序多数据流方法(SPMD),通过根据不同的进程号转移到不同的分支语句,使得运行一个程序就能获得运行多个不同效果的效果。

要想从一个MPI进程发送数据给另一个进程,可以调用MPI_Send函数,而MPI_Recv函数是用来接受消息的。

MPI_Recv是阻塞的,直到消息接受到(或者发生错误)前,该函数不会返回。MPI_Send的行为则由MPI的定义来实现,它可以阻塞或缓冲发送消息。当它阻塞时,在相应的接收启动前它不会返回;如果缓冲消息。MPI会将其复制到它私有的存储空间,一旦复制完成,MPI_Send就会返回。

编写MPI程序时,非常重要的一点是区分全局变量和局部变量。局部变量只作用在定义它的那个进程,而全局变量作用全部进程。

MPI程序允许所有的进程将结果打印到stdout和stderr中。但是一般只允许进程0来读取stdin输入的数据。

集合通信涉及一个通信子中的所有进程。

MPI_Send和MPI_Recv常常被称为点对点通信。

两个常见的集合通信函数MPI_Reduce和MPI_Allreduce。MPI_Reduce存储全局操作的结果到指定进程,MPI_Allreduce将结果存储到通信子中的全部进程中。

有些类似于MPI_Reduce的MPI函数,会将同样的参数传递到输入和输出缓冲区中,这种现象称为参数别名,MPI显式禁止了输出参数和其他参数相混淆。

MPI_Bcast从单个进程向同一个通信子中的所有其他进程发送消息。

MPI_Scatter指定进程在各个进程间分配一个数组的元素。

MPI_Gather将收集各个进程间的元素到指定进程。

MPI_Allgarher收集所有的元素发送给所有的进程。

MPI_Barrier用来同步进程,在同一个通信子的所有进程调用该函数前,所有调用MPI_Barrier的进程都不能返回。

在编写MPI程序时,如何将数据分割并分配到各个进程是关键点。一般有块划分,循环划分,块-循环划分三种方式。

与只涉及CPU与主存操作相比,发送消息的开销是昂贵的。而且,用尽可能少的次数发送尽可能多的消息可以节省开销,所以将多个消息合并为一个消息发送是降低开销的好方法。

MPI有降低消息发送开销三种方法:通信函数中的count参数,派生数据类型和MPI_Pack/Unpack函数。

计算并行程序运行时间时,我们一般关注经过的总时间或叫“墙上时间”,即运行一段代码所需要的时间,它包括用户级代码、库函数、用户代码调用系统函数的运行时间和空闲时间,一般建议多次运行后报告最短的结果。

统计完时间后,可以用加速比和效率来评估程序的性能。

并行开销是并行程序运行时间的一部分,指的是不能被串行程序执行的额外时间(例如,进程之间通信)。

如果MPI程序的正确行为取决为MPI_Send正在缓冲的输入,则它是不安全的。

若进程双方都阻塞并等待永不发生的事件发生,则进程就死锁了。

MPI程序能通过互相替换每个MPI_Send和MPI_Recv函数来检测是否安全。

MPI_Ssend使用MPI_Send一样的参数,但它一直阻塞直到对应的接收端开启。额外的s表明同步。如果使用MPI_Ssend的MPI程序根据要求的输入和通信规模正确地完成了运行,则该程序是安全的。

一个不安全的MPI程序有多种方法变成安全,程序员可以调度MPI_Send和MPI_Recv使某些进程(偶数序列号的进程)先调用MPI_Send,其他进程(奇数序列号的进程)先调用MPI_Recv。另外,可以使用MPI_Sendrecv或者MPI_Sendrecv_replace,这些函数个发送和接收一条消息,它们保证各自程序不会奔溃或死锁。

  • 1
    点赞
  • 22
    收藏
    觉得还不错? 一键收藏
  • 0
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值