高性能计算-雅可比算法-MPI重复非阻塞优化(7)

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

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

//并行-重复非阻塞-笛卡尔64个进程块
/*
优化：通信接口分步优化 MPI_Start
要求：
1、仍然采取二维行列同时分块的方式，对数据进行区域分解。
2、使用笛卡尔虚拟拓扑的相关接口，进行二维进程网格阵列构建、邻居进程编号获取等操作。
3、使用64个进程，每个进程初始化一个二维子块，并负责该子块的Jacobi迭代计算。
4、使用重复非阻塞通信接口。
思路:使用虚拟进程和笛卡尔，每个进程都对上下左右块通信
*/
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 top=0,bottom=0,left=0,right=0;	//标记每个block实际数据的边界
	//int ltBID=0,rtBID=1,lbBID=2,rbBID=3;//标记四个角落位置的进程
	MPI_Status arr_status[2+S*2]={0};		//
	MPI_Request* arr_requestS = calloc(2+S*2,sizeof(MPI_Request));	//发送请求
	MPI_Request* arr_requestR = calloc(2+S*2,sizeof(MPI_Request));	//接收请求
	int dims[2]={N/S,N/S};	//8,8
	int periods[2] = {0,0};	//每个维度中不定期
	MPI_Comm cartcomm;
	int pid;				//当前进程ID
	int size=0;
	int coords[2] ={0};		//当前进程笛卡尔坐标
	int nbrs[4] ={0};		//上下左右进程的id
	enum DIR{UP,DOWN,LEFT,RIGHT};

	MPI_Init(&argc,&argv);
	MPI_Comm_size(MPI_COMM_WORLD,&size);
	MPI_Comm_rank(MPI_COMM_WORLD,&pid);
	if(size == B)
	{
		MPI_Cart_create(MPI_COMM_WORLD,2,dims,periods,0,&cartcomm);
		MPI_Cart_coords(cartcomm,pid,2,coords);

		//初始化
		for(int i=0; i<BS; i++)
		{
			for(int j=0; j<BS; j++)
			{
				rows[i][j] = 0.0;
				rows2[i][j] = 0.0;
			}
		}
		//确定上下左右进程的id
		MPI_Cart_shift(cartcomm,0,1,&nbrs[UP],&nbrs[DOWN]);
		MPI_Cart_shift(cartcomm,1,1,&nbrs[LEFT],&nbrs[RIGHT]);
		//有效数据边界初始化
		if(pid%RB == 0)//left
		{
			for(int i=1;i<BS-1;i++)
				rows[i][1] = 8.0;
		}
		if(pid%RB==RB-1)//right
		{
			for(int i=1;i<BS-1;i++)
				rows[i][BS-2] = 8.0;
		}
		if(pid>=0 && pid<RB)//top
		{
			for(int i=1;i<BS-1;i++)
				rows[1][i] = 8.0;
		}
		if(pid<B && pid>=B-RB)//bottom
		{
			for(int i=1;i<BS-1;i++)
				rows[BS-2][i] = 8.0;
		}


		//建立通信连接
		for(int i=0;i<4;i++)
		{
			if(nbrs[i] == -1)
				nbrs[i] = MPI_PROC_NULL;
		}
		//上下左右
		MPI_Recv_init(&rows[0][1],S,MPI_FLOAT,nbrs[UP],0,cartcomm,&arr_requestR[0]);
		MPI_Recv_init(&rows[BS-1][1],S,MPI_FLOAT,nbrs[DOWN],0,cartcomm,&arr_requestR[1]);
		for(int i=1,k=2;i<BS-1;i++,k+=2)
		{
			MPI_Recv_init(&rows[i][0],1,MPI_FLOAT,nbrs[LEFT],0,cartcomm,&arr_requestR[k]);
			MPI_Recv_init(&rows[i][BS-1],1,MPI_FLOAT,nbrs[RIGHT],0,cartcomm,&arr_requestR[k+1]);
		}
		
		MPI_Send_init(&rows[1][1],S,MPI_FLOAT,nbrs[UP],0,cartcomm,&arr_requestS[0]);
		MPI_Send_init(&rows[BS-2][1],S,MPI_FLOAT,nbrs[DOWN],0,cartcomm,&arr_requestS[1]);
		for(int i=1,k=2;i<BS-1;i++,k+=2)
		{
			MPI_Send_init(&rows[i][1],1,MPI_FLOAT,nbrs[LEFT],0,cartcomm,&arr_requestS[k]);
			MPI_Send_init(&rows[i][BS-2],1,MPI_FLOAT,nbrs[RIGHT],0,cartcomm,&arr_requestS[k+1]);
		}

		//块内需要计算数据的边界索引
		int rbegin,rend;	//块内起始 终止列号
		int cbegin,cend;	//块内列起始 终止列号
		rbegin = (pid>=0 && pid<RB)?2:1;
		rend = (pid<B && pid>=B-RB)?BS-3:BS-2;
		cbegin =(pid%RB == 0)?2:1;
		cend = (pid%RB == RB-1)?BS-3:BS-2;

		//迭代
		for(int step=0; step<T; step++)
		{
			//每个进程都完成收发数据才能计算
			RequestStart(2+S*2,arr_requestR);
			RequestStart(2+S*2,arr_requestS);
			MPI_Waitall(2+S*2,arr_requestR,arr_status);
			MPI_Waitall(2+S*2,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=1,m=0;i<BS-1;i++,m++)
		{
			for(int j=1,n=0;j<BS-1;j++,n++)
				temprows[m][n] = rows[i][j];
		}
		MPI_Barrier(MPI_COMM_WORLD);
		MPI_Gather(temprows,S*S,MPI_FLOAT,temprows1,S*S,MPI_FLOAT,0,MPI_COMM_WORLD);

		//对数据重新整理
		//遍历temprows1
		int index=0;
		for(int rb=0;rb<RB;rb++)//块行索引
		{
			for(int cb=0;cb<RB;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(2+S*2,arr_requestR);
		RequestFree(2+S*2,arr_requestS);
	}
	else if(pid==0)
	{
		printf("parameter:should -n %d\n",B);
	}
	free(arr_requestS);
	free(arr_requestR);
	MPI_Finalize();
	return 0;
}