高性能计算-雅可比算法MPI通信优化(5)

1. 雅可比算法原理：如下图对方阵非边界元素求上下左右元素的均值，全部计算元素的数值计算完成后更新矩阵，进行下一次迭代。
   
   

2. 测试目标：用MPI实现对8*8方阵雅可比算法迭代并行计算，用重复非阻塞的通信方式

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

#define N 8		//方阵行列数
#define B 4		//并行进程数
#define S N/(B/2)	//分块方阵的大小
#define BS S+1	//块包含交换数据的方阵大小
#define T 2		//迭代次数

//并行-重复非阻塞-4*4分块并行计算，共4个块，每个计算块包含其他块计算数据的块大小为 5*5
/*
优化：通信接口分步优化 MPI_Start
*/
void printRows(int pid,float rows[BS][BS])
{
	printf("result in %d\n",pid);
	for(int i=0;i<BS;i++)
	{
		for(int j=0;j<BS;j++)
			printf("%.3f\t",rows[i][j]);
		printf("\n");
	}
}

void RequestStart(int count,MPI_Request arr_request[])
{
	for(int i=0;i<count;i++)
		MPI_Start(&arr_request[i]);
}

void RequestFree(int count,MPI_Request arr_request[])
{
	for(int i=0;i<count;i++)
		MPI_Request_free(&arr_request[i]);
}

int main(int argc,char* argv[])
{
	float rows[BS][BS],rows2[BS][BS],temprows[S][S],temprows1[N][N],finalrows[N][N];	
	
	int pid;
	int top=0,bottom=0,left=0,right=0;	//标记每个block实际数据的边界
	int ltBID=0,rtBID=1,lbBID=2,rbBID=3;//标记四个角落位置的进程
	MPI_Status arr_status[BS]={0};		//
	MPI_Request arr_requestS[BS] = {0};	//发送请求 第0个:行数据请求
	MPI_Request arr_requestR[BS] = {0};	//接收请求 第0个:行数据请求

	MPI_Init(&argc,&argv);
	MPI_Comm_rank(MPI_COMM_WORLD,&pid);
	
	//初始化
	for(int i=0; i<BS; i++)
	{
		for(int j=0; j<BS; j++)
		{
			rows[i][j] = 0.0;
			rows2[i][j] = 0.0;
		}
	}
	//有效数据边界初始化
	if(ltBID==pid || rtBID==pid)
	{
		top = 0;
		bottom = S-1;
	}
	else
	{
		top = 1;
		bottom = S;
	}
	if(ltBID==pid || lbBID==pid)
	{
		left = 0;
		right = S-1;
	}
	else
	{
		left = 1;
		right = S;
	}
	//数据初始化
	if(ltBID==pid || rtBID==pid)
	{
		for(int j=left;j<=right;j++)
			rows[top][j] = 8.0;
	}
	else if(lbBID==pid || rbBID==pid)
	{
		for(int j=left;j<=right;j++)
			rows[bottom][j] = 8.0;
	}
	if(ltBID==pid||lbBID==pid)
	{
		for(int i=top;i<=bottom;i++)
			rows[i][left] = 8.0;
	}
	else if(rtBID==pid || rbBID==pid)
	{
		for(int i=top;i<=bottom;i++)
			rows[i][right] = 8.0;
	}
	
	//建立通信连接
	if(ltBID==pid)
	{
		MPI_Recv_init(&rows[S],S,MPI_FLOAT,pid+2,0,MPI_COMM_WORLD,&arr_requestR[0]);	
		for(int i=top,k=1;i<=bottom;i++,k++)
		{
			MPI_Recv_init(&rows[i][S],1,MPI_FLOAT,pid+1,0,MPI_COMM_WORLD,&arr_requestR[k]);	
			MPI_Send_init(&rows[i][S-1],1,MPI_FLOAT,pid+1,0,MPI_COMM_WORLD,&arr_requestS[k]);
		}
		MPI_Send_init(&rows[S-1],S,MPI_FLOAT,pid+2,0,MPI_COMM_WORLD,&arr_requestS[0]);
	}
	if(rtBID==pid)
	{
		MPI_Recv_init(&rows[S][1],S,MPI_FLOAT,pid+2,0,MPI_COMM_WORLD,&arr_requestR[0]);	
		for(int i=top,k=1;i<=bottom;i++,k++)
		{
			MPI_Recv_init(&rows[i],1,MPI_FLOAT,pid-1,0,MPI_COMM_WORLD,&arr_requestR[k]);	
			MPI_Send_init(&rows[i][1],1,MPI_FLOAT,pid-1,0,MPI_COMM_WORLD,&arr_requestS[k]);
		}
		MPI_Send_init(&rows[S-1][1],S,MPI_FLOAT,pid+2,0,MPI_COMM_WORLD,&arr_requestS[0]);
	}
	if(lbBID==pid)
	{
		MPI_Recv_init(&rows[0],S,MPI_FLOAT,pid-2,0,MPI_COMM_WORLD,&arr_requestR[0]);	
		for(int i=top,k=1;i<=bottom;i++,k++)
		{
			MPI_Recv_init(&rows[i][S],1,MPI_FLOAT,pid+1,0,MPI_COMM_WORLD,&arr_requestR[k]);	
			MPI_Send_init(&rows[i][S-1],1,MPI_FLOAT,pid+1,0,MPI_COMM_WORLD,&arr_requestS[k]);
		}
		MPI_Send_init(&rows[1],S,MPI_FLOAT,pid-2,0,MPI_COMM_WORLD,&arr_requestS[0]);
	}
	if(rbBID==pid)
	{
		MPI_Recv_init(&rows[0][1],S,MPI_FLOAT,pid-2,0,MPI_COMM_WORLD,&arr_requestR[0]);	
		for(int i=top,k=1;i<=bottom;i++,k++)
		{
			MPI_Recv_init(&rows[i],1,MPI_FLOAT,pid-1,0,MPI_COMM_WORLD,&arr_requestR[k]);	
			MPI_Send_init(&rows[i][1],1,MPI_FLOAT,pid-1,0,MPI_COMM_WORLD,&arr_requestS[k]);
		}
		MPI_Send_init(&rows[1][1],S,MPI_FLOAT,pid-2,0,MPI_COMM_WORLD,&arr_requestS[0]);
	}
	//块内需要计算数据的边界索引
	int rbegin,rend;	//块内起始 终止列号
	int cbegin,cend;	//块内列起始 终止列号
	rbegin = 1;
	rend = S-1;
	cbegin = 1;
	cend = S-1;
	//迭代
	for(int step=0; step<T; step++)
	{
		//每个进程都完成收发数据才能计算
		RequestStart(BS,arr_requestR);
		RequestStart(BS,arr_requestS);
		MPI_Waitall(BS,arr_requestR,arr_status);
		MPI_Waitall(BS,arr_requestS,arr_status);

		//计算
		for(int i=rbegin;i<=rend;i++)
		{
			for(int j=cbegin;j<=cend;j++)
				rows2[i][j] =0.25*(rows[i-1][j]+rows[i][j-1]+rows[i][j+1]+rows[i+1][j]);
		}
		//更新
		for(int i=rbegin;i<=rend;i++)
		{
			for(int j=cbegin;j<=cend;j++)
				rows[i][j] = rows2[i][j];
		}
	}
	//打印
	sleep(pid);
	printRows(pid,rows);
	
	//Gather data from all processes
	for(int i=top,m=0;i<=bottom;i++,m++)
	{
		for(int j=left,n=0;j<=right;j++,n++)
			temprows[m][n] = rows[i][j];
	}
	MPI_Barrier(MPI_COMM_WORLD);
	MPI_Gather(temprows,16,MPI_FLOAT,temprows1,16,MPI_FLOAT,0,MPI_COMM_WORLD);
	
	//对数据重新整理
	//遍历temprows1
	int index=0;
	for(int rb=0;rb<2;rb++)//块行索引
	{
		for(int cb=0;cb<2;cb++)//块列索引
		{
			for(int r=0;r<S;r++)
			{
				for(int c=0;c<S;c++)
				{
					finalrows[rb*S+r][cb*S+c] = *((float*)&temprows1+index++);
				}
			}
		}
	}
	if(pid==0)
	{
        fprintf(stderr,"\nResult after gathering data:\n");
        for(int i = 0; i < N; i++)
        {
            for(int j = 0; j < N; j++)
                fprintf(stderr,"%.3f\t", finalrows[i][j]);
            fprintf(stderr,"\n");
        }
        fprintf(stderr,"\n");	
	}
	RequestFree(BS,arr_requestR);
	RequestFree(BS,arr_requestS);
	MPI_Finalize();
	return 0;
}