OpenMPI(一) 消息的发布和接收

之前的OpenMP(Open Multi-Processing)系列是在共享内存(shared memory parallelism)的基础上实现的。而在并行运算中，分布式内存(distributed memory parallelism)也很常见，在这种情况下，每个操作单元独立的拥有一块内存，数据依赖关系将通过直接的交流解决，所以将不存在data race的问题。

在这种情况下，每个进程之间的交流变得尤为重要。而MPI（Message Passing Interface）就是一套关于消息交流接口的标准，OpenMPI则是对MPI的免费开源实现之一。

OpenMPI示例代码一：

#include "stdio.h"
#include "stdlib.h"
#include "mpi.h" // 导入mpi.h

int main( int argc, char **argv)
{
    int numprocs, rank;

    MPI_Init(&argc, &argv); // 初始化
    MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);

    printf("Process %d of %d.\n", rank, numprocs);

    MPI_Finalize(); // 结束
    return EXIT_SUCCESS;
}

注意：OpenMPI的使用必须以MPI_Init()开始，以MPI_Finalize()结束。在这两个函数囊括的范围之外使用MPI语句将会报错。

其中，MPI_COMM_WORLD是一个表示了MPI中所有进程的通信域（communicator）。

MPI_Comm_size(MPI_COMM_WORLD , &numprocs);定义总进程数
MPI_Comm_rank(MPI_COMM_WORLD , &rank);定义进程的rank,范围是[0,总进程数-1）

OpenMPI示例代码二：

再看一段简单的矩阵相加代码：

for(i = 0; i < N; i++){
    c[i]  = a[i] + b[i];
} 

如果我们想用OpenMPI并行化的实现这部分代码，由于这是一个分布式内存系统，会遇到一个问题：

对于矩阵a，假设它存储在一个特定进程（如rank = 0的进程）下的内存里。由于这是分布式的内存系统，其他进程无法直接访问到进程0所在的内存，也就无法直接读取到矩阵a的值。所以需要把矩阵a的值告诉其他进程。同理，矩阵b也会面临相同的问题。

在尝试解决这个问题前，先了解一下OpenMP中用于点对点通信（point to point communication）的几个函数：

1. MPI_Send()函数

以下代码表示rank为0的进程将矩阵a的值发送给其他进程。

if(rank == 0 ){
 for(p = 1; p < numProcs; p++){
    MPI_Send(
        &a[p * localSize], // 指向数据的指针
        localSize, // 发送数据的大小
        MPI_FLOAT, // 数据类型
        p, // 接受者进程的rank
        0, // Tag ; 一般设为0
        MPI_COMM_WORLD // Communicator
    );
 }   
}

在这段代码中，localSize = 矩阵a的大小 / 总进程数，这样就可以把矩阵a平均分配给不同的进程进行计算。

2. MPI_Recv()函数

有了发布信息的进程，那就会有接受信息的进程。其他rank > 0的进程都会接受来自rank为0的进程的消息。

MPI_Status status;
...
if(rank > 0)
{
    MPI_Recv(
        local_a, // 指向数据的指针
        localSize, // 发送过来的数据的大小
        MPI_FLOAT, // 数据类型
        0, // 发布者进程的rank
        0, // Tag 一般设置为0
        MPI_COMM_WORLD, // Communicator
        &status // MPI_Status对象
    );
}

可以通过这张图来加深理解：

其他进程在接受到进程0发送的消息以后，会纷纷在自己的内存中进行矩阵加法。并且把结果通过相同步骤返回给进程0。

注意：

1. 从 p=0开始for循环，因为进程0不需要发送代码给自己。
2. 数据类型包括MPI FLOAT, MPI INT, MPI DOUBLE, MPI CHAR等。
3. 大多数MPI调用成功后会返回MPI_SUCCESS。通过status可以判断是否发生错误。我们也可以在&status处填入MPI_STATUS_IGNORE，表示忽略状态。

完整代码

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>


#define N 100


int main( int argc, char **argv )
{
	int i, p;

	// Initialise MPI and find the total number of processes and the rank of 'this' process.
	int rank, numProcs;
	MPI_Init( &argc, &argv );
	MPI_Comm_size( MPI_COMM_WORLD, &numProcs );
	MPI_Comm_rank( MPI_COMM_WORLD, &rank     );

	// Check that the problem size is a multiple of the number of processes.
	if( N%numProcs )
	{
		if( rank==0 ) printf( "Problem size '%d' is not a multiple of the number of processes '%d.'\n", N, numProcs );
		MPI_Finalize();
		return EXIT_FAILURE;
	}

	// Initialise the full vectors a, b and c on rank 0 only.
	float *a=NULL, *b=NULL, *c=NULL;
	if( rank==0 )
	{
		// Allocate memory for the vectors.
		a = (float*) malloc( N*sizeof(float) );
		b = (float*) malloc( N*sizeof(float) );
		c = (float*) malloc( N*sizeof(float) );
		if( !a || !b || !c )
		{
			printf( "Could not allocate memory for the vectors.\n" );
			MPI_Finalize();
			return EXIT_FAILURE;
		}

		// Put values into vectors a and b only.
		for( i=0; i<N; i++ ) a[i] = b[i] = i+1;
	}

	// Initialise the local vectors on all processes except rank 0. Note that rank 0 can use the
	// same arrays a, b and c as for the full vectors, but it helps keep things easier to understand
	// to use new names for the local variables and arrays.
	int localSize = N / numProcs;
	float *local_a=NULL, *local_b=NULL, *local_c=NULL;
	if( rank>0 )
	{
		local_a = (float*) malloc( localSize*sizeof(float) ),
		local_b = (float*) malloc( localSize*sizeof(float) ),
		local_c = (float*) malloc( localSize*sizeof(float) );
	}

	//
	// Step 1. Rank 0 sends segments of the full vectors a and b to their respective processes;
	// those processes must in turn receive the data.
	//
	if( rank==0 )
	{
		for( p=1; p<numProcs; p++ )
		{
			MPI_Send( &a[p*localSize], localSize, MPI_FLOAT, p, 0, MPI_COMM_WORLD );
			MPI_Send( &b[p*localSize], localSize, MPI_FLOAT, p, 0, MPI_COMM_WORLD );
		}
	}
	else
	{
		MPI_Recv( local_a, localSize, MPI_FLOAT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE );
		MPI_Recv( local_b, localSize, MPI_FLOAT, 0, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE );
	}

	//
	// Step 2. Each process performs vector addition on their local segments of the full arrays.
	//
	if( rank==0 )
		for( i=0; i<localSize; i++ ) c[i] = a[i] + b[i];
	else
		for( i=0; i<localSize; i++ ) local_c[i] = local_a[i] + local_b[i];

	//
	// Step 3. Recombine the full c array back on rank 0.
	//
	if( rank==0 )
	{
		for( p=1; p<numProcs; p++ )
			MPI_Recv( &c[p*localSize], localSize, MPI_FLOAT, p, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE );
	}
	else
	{
		MPI_Send( local_c, localSize, MPI_FLOAT, 0, 0, MPI_COMM_WORLD );
	}

	// Check the answer on rank 0. Only display the first few array elements if N is large.
	if( rank==0 )
	{
		for( i=0; i<(N>20?20:N); i++ )
			printf( "%g + %g = %g.\n", a[i], b[i], c[i] );
		if( N>20 )
			printf( "(remaining elements not displayed)\n" );
	}

	// Clear up.
	if( rank==0 )
	{
		free( a );
		free( b );
		free( c );
	}
	else
	{
		free( local_a );
		free( local_b );
		free( local_c );
	}
	MPI_Finalize();
	return EXIT_SUCCESS;
}