url: https://hpc-tutorials.llnl.gov/mpi/exercise_3/

Exercise 3

使用 MPI 实现并行的 “data decomposition on an array”。

为了方便理解，先看串行版本：

#include <stdio.h>
#include <stdlib.h>
#define	ARRAYSIZE	20000000
float	data[ARRAYSIZE]; 	/* the initial array */

int main(int argc, char *argv[])
{
	int     i; 			/* loop variable */

	printf("Starting serial array example...\n");
	printf("Using array of %d floats. Requires %ld bytes\n",ARRAYSIZE,sizeof(data));

	/* Initialize the array */
	printf("Initializing array...\n");
	for(i=0; i<ARRAYSIZE; i++) 
  	data[i] =  i * 1.0;

	/* Do a simple value assignment to each of the array elements */
	printf("Performing computation on array elements...\n");
	for(i=1; i < ARRAYSIZE; i++)
   	data[i] = data[i] + i * 1.0;

	/* Print a few sample results */
	printf("Sample results\n");
	printf("   data[1]=%e\n",  data[1]);
	printf("   data[100]=%e\n",  data[100]);
	printf("   data[1000]=%e\n",  data[1000]);
	printf("   data[10000]=%e\n",  data[10000]);
	printf("   data[100000]=%e\n",  data[100000]);
	printf("   data[1000000]=%e\n",  data[1000000]);
	printf("\nAll Done!\n");
}

编译运行：

gcc test.c
./a.out

预期输出：

Starting serial array example...
Using array of 20000000 floats. Requires 80000000 bytes
Initializing array...
Performing computation on array elements...
Sample results
   data[1]=2.000000e+00
   data[100]=2.000000e+02
   data[1000]=2.000000e+03
   data[10000]=2.000000e+04
   data[100000]=2.000000e+05
   data[1000000]=2.000000e+06

All Done!

接下来是实现的并行版本：

#include "mpi.h"
#include <stdio.h>
#include <stdlib.h>
#define  ARRAYSIZE	20000000  // 定义数组大小：2000万个元素
#define  MASTER		0         // 定义主进程编号为0

double  data[ARRAYSIZE];      // 全局数组，所有进程共享（在MPI中每个进程有独立副本）

int main (int argc, char *argv[])
{
	int   numtasks, taskid, rc, dest, offset, i, j, tag1,
      	tag2, source, chunksize, leftover; 
	double mysum, sum;           // mysum: 局部和, sum: 全局和
	double update(int myoffset, int chunk, int myid);  // 函数声明
	MPI_Status status;

	/***** MPI初始化 *****/
	MPI_Init(&argc, &argv);                    // 初始化MPI环境
	MPI_Comm_size(MPI_COMM_WORLD, &numtasks);  // 获取进程总数
	MPI_Comm_rank(MPI_COMM_WORLD,&taskid);     // 获取当前进程ID
	printf ("MPI task %d has started...  ", taskid);
	chunksize = (ARRAYSIZE / numtasks);        // 计算每个进程处理的数据块大小
	leftover = (ARRAYSIZE % numtasks);         // 计算余数，主进程会处理这些额外数据
	tag2 = 1;  // 消息标签1：用于数据发送
	tag1 = 2;  // 消息标签2：用于偏移量发送

	/***** 主进程逻辑 *****/
	if (taskid == MASTER){

  		/* 初始化数组并计算总和用于验证 */
  		sum = 0;
  		for(i=0; i<ARRAYSIZE; i++) {
    		data[i] =  i * 1.0;        // 数组赋值为索引值
    		sum = sum + data[i];       // 计算初始总和
		}
  		printf("Initialized array sum = %e\n",sum);
  		printf("numtasks= %d  chunksize= %d  leftover= %d\n",numtasks,chunksize,leftover);

  		/* 向其他进程分发数据 - 主进程保留第一部分和剩余元素 */
  		offset = chunksize + leftover;  // 主进程处理前 (chunksize+leftover) 个元素
  		for (dest=1; dest<numtasks; dest++) {
    		MPI_Send(&offset, 1, MPI_INT, dest, tag1, MPI_COMM_WORLD);          // 发送数据偏移量
    		MPI_Send(&data[offset], chunksize, MPI_DOUBLE, dest, tag2, MPI_COMM_WORLD); // 发送数据块
    		printf("Sent %d elements to task %d offset= %d\n",chunksize,dest,offset);
    		offset = offset + chunksize;  // 更新下一个数据块的偏移量
		}

  		/* 主进程处理自己的数据部分 */
  		offset = 0;  // 主进程从数组开头处理
  		mysum = update(offset, chunksize+leftover, taskid);  // 处理数据并计算局部和

  		/* 等待接收其他进程的处理结果 */
  		for (i=1; i<numtasks; i++) {
    		source = i;
    		MPI_Recv(&offset, 1, MPI_INT, source, tag1, MPI_COMM_WORLD, &status);     // 接收偏移量
    		MPI_Recv(&data[offset], chunksize, MPI_DOUBLE, source, tag2,              // 接收处理后的数据
      		MPI_COMM_WORLD, &status);
		}

  		/* 使用归约操作获取最终总和并打印样本结果 */  
  		MPI_Reduce(&mysum, &sum, 1, MPI_DOUBLE, MPI_SUM, MASTER, MPI_COMM_WORLD);

  		printf("Sample results: \n");
  		offset = 0;
  		for (i=0; i<numtasks; i++) {
    		for (j=0; j<5; j++) 
				printf("  %e",data[offset+j]);  // 打印每个进程的前5个元素作为样本
    		printf("\n");
    		offset = offset + chunksize;
		}
  		printf("*** Final sum= %e ***\n",sum);

  	}  /* 主进程逻辑结束 */

	/***** 非主进程（工作进程）逻辑 *****/
	if (taskid > MASTER) {

		/* 从主进程接收分配的数据块 */
		source = MASTER;
		MPI_Recv(&offset, 1, MPI_INT, source, tag1, MPI_COMM_WORLD, &status);     // 接收数据偏移量
		MPI_Recv(&data[offset], chunksize, MPI_DOUBLE, source, tag2,              // 接收数据块
			MPI_COMM_WORLD, &status);

		/* 处理分配的数据部分 */
		mysum = update(offset, chunksize, taskid);  // 计算局部和

		/* 将处理结果发送回主进程 */
		dest = MASTER;
		MPI_Send(&offset, 1, MPI_INT, dest, tag1, MPI_COMM_WORLD);           // 发送偏移量
		MPI_Send(&data[offset], chunksize, MPI_DOUBLE, MASTER, tag2, MPI_COMM_WORLD); // 发送处理后的数据

		/* 参与归约操作，将局部和汇总到主进程 */
		MPI_Reduce(&mysum, &sum, 1, MPI_DOUBLE, MPI_SUM, MASTER, MPI_COMM_WORLD);

	} /* 非主进程逻辑结束 */

	MPI_Finalize();  // 终止MPI环境

}   /* main函数结束 */


/* 数据处理函数：对指定数据段进行操作并计算局部和 */
double update(int myoffset, int chunk, int myid) {
	int i; 
	double mysum;
	/* 对分配的数组元素执行加法操作并计算局部和 */
	mysum = 0;
	for(i=myoffset; i < myoffset + chunk; i++) {
		data[i] = data[i] + (i * 1.0);  // 每个元素加上其索引值
		mysum = mysum + data[i];         // 累加计算局部和
	}
	printf("Task %d mysum = %e\n",myid,mysum);
	return(mysum);
}

编译运行：

mpicc test.c -o test
mpirun -np 4 ./test

预期输出：

MPI task 0 has started...  Initialized array sum = 2.000000e+14
numtasks= 4  chunksize= 5000000  leftover= 0
Sent 5000000 elements to task 1 offset= 5000000
MPI task 1 has started...  Task 1 mysum = 7.500000e+13
Sent 5000000 elements to task 2 offset= 10000000
MPI task 2 has started...  Task 2 mysum = 1.250000e+14
Sent 5000000 elements to task 3 offset= 15000000
Task 0 mysum = 2.500000e+13
MPI task 3 has started...  Task 3 mysum = 1.750000e+14
Sample results: 
  0.000000e+00  2.000000e+00  4.000000e+00  6.000000e+00  8.000000e+00
  1.000000e+07  1.000000e+07  1.000000e+07  1.000001e+07  1.000001e+07
  2.000000e+07  2.000000e+07  2.000000e+07  2.000001e+07  2.000001e+07
  3.000000e+07  3.000000e+07  3.000000e+07  3.000001e+07  3.000001e+07
*** Final sum= 4.000000e+14 ***