#include <mpi.h>
#include <stdlib.h>
#include <stdio.h>
#include <omp.h>
#define Type double

#if Type == double
#define MPI_TYPE MPI_DOUBLE
#else
#define MPI_TYPE MPI_FLOAT
#endif

inline Type **Allocate2D(int M, int N);
inline void Free(Type **ptr);
inline void Swap(Type ***A, Type ***B);

int main()
{
	MPI_Init(NULL, NULL);
	int rank, numproc;
	FILE *fp;
	FILE *fp_out;
	MPI_Comm_size(MPI_COMM_WORLD, &numproc);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);

	if (rank == 0)
	{
		fp = fopen("./input_m.txt", "w");
		if (!fp)
		{
			printf("Fail to creat\n");
			exit(0);
		}
	}

	const int N = 1000;
	/*distribute whole matrix along row diretion*/
	int row_start = (rank * N + numproc - 1) / numproc;
	int row_end = ((rank + 1) * N + numproc - 1) / numproc;
	int subsize = row_end - row_start;
	int *count = (int *)malloc(numproc * sizeof(int));
	int *disp = (int *)malloc(numproc * sizeof(int));

	MPI_Gather(&subsize, 1, MPI_INT, count, 1, MPI_INT, 0, MPI_COMM_WORLD);

	if (rank == 0)
	{
		for (int i = 0; i < numproc; i++)
			count[i] *= N;
		disp[0] = 0;
		for (int i = 1; i < numproc; i++)
			disp[i] = disp[i - 1] + count[i - 1];
	}

	Type **A_all = Allocate2D(N, N);
	Type **A = Allocate2D(subsize + 2, N);
	Type **A_tmp = Allocate2D(subsize + 2, N);

	/*set to value*/
	for (int i = 0; i < subsize + 2; i++)
	{
		for (int j = 0; j < N; j++)
		{
			/*set the value lying on left and right boundary to be 8*/
			if (j == 0 || j == N - 1)
				A[i][j] = A_tmp[i][j] = 8;
			else
				A[i][j] = A_tmp[i][j] = 0;
		}
	}

	/*set the value lying on top to be 8*/
	if (rank == 0)
		for (int j = 0; j < N; j++)
			A[1][j] = A_tmp[1][j] = 8;

	/*set the value lying on botoom to be 8*/
	if (rank == numproc - 1)
		for (int j = 0; j < N; j++)
			A[subsize][j] = A_tmp[subsize][j] = 8;

	/*gather input*/
	MPI_Gatherv(A[1], subsize * N, MPI_TYPE, A_all[0], count, disp, MPI_TYPE, 0, MPI_COMM_WORLD);

	if (rank == 0)
	{
		fprintf(fp, "Input data:\n");
		for (int i = 0; i < N; i++)
		{
			fprintf(fp, "\t");
			for (int j = 0; j < N; j++)
				fprintf(fp, "%.2lf ", A_all[i][j]);
			fprintf(fp, "\n");
		}
	}
	MPI_Barrier(MPI_COMM_WORLD);
	/*build relationship*/
	int up = rank - 1;
	int down = rank + 1;
	if (rank == 0)
		up = MPI_PROC_NULL;
	if (rank == numproc - 1)
		down = MPI_PROC_NULL;

	/*iteration*/
	row_start = 1;
	row_end = subsize + 1;
	if (rank == 0)
		row_start = 2;
	if (rank == numproc - 1)
		row_end = subsize;

	double t_start = omp_get_wtime();
	for (int ite = 0; ite < 100; ite++)
	{
		/*send data to upper process and recv form lower process*/
		MPI_Sendrecv(A[1], N, MPI_TYPE, up, 0, A[subsize + 1], N, MPI_TYPE, down, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);
		/*send data to lower process and recv form upper process*/
		MPI_Sendrecv(A[subsize], N, MPI_TYPE, down, 0, A[0], N, MPI_TYPE, up, 0, MPI_COMM_WORLD, MPI_STATUS_IGNORE);

		Swap(&A, &A_tmp);

		for (int i = row_start; i < row_end; i++)
		{
			for (int j = 1; j < N - 1; j++)
			{
				/*update A using value of the adjcent elemets of last iteration*/
				A[i][j] = 0.25 * (A_tmp[i - 1][j] + A_tmp[i + 1][j] + A_tmp[i][j - 1] + A_tmp[i][j + 1]);
			}
		}
	}
	double t_end = omp_get_wtime();
	if (rank == 0)
		printf("Runing time is %lf\n", t_end - t_start);

	/*gather the output result*/
	MPI_Gatherv(A[1], subsize * N, MPI_TYPE, A_all[0], count, disp, MPI_TYPE, 0, MPI_COMM_WORLD);

	if (rank == 0)
	{
		fp_out = fopen("./result_m.txt", "w");
		/*print the result into output file*/
		fprintf(fp_out, "Output data:\n");
		for (int i = 0; i < N; i++)
		{
			fprintf(fp_out, "\t");
			for (int j = 0; j < N; j++)
				fprintf(fp_out, "%.2lf ", A_all[i][j]);
			fprintf(fp_out, "\n");
		}
		fclose(fp_out);
	}
	MPI_Barrier(MPI_COMM_WORLD);
	free(count);
	free(disp);
	Free(A_all);
	Free(A);
	Free(A_tmp);
	MPI_Finalize();
}

//创建二维矩阵，分配二维数组内存空间，存储数据
Type **Allocate2D(int M, int N)
{
    Type **tmp = NULL;
    tmp = malloc(M * sizeof(Type *));

    Type *Data = malloc(M * N * sizeof(Type));

    for (int i = 0; i < M; i++)
        tmp[i] = Data + i * N;

    return tmp;
}

//free
void Free(Type **ptr)
{
    if (ptr[0])
        free(ptr[0]);
    if (ptr)
        free(ptr);
}

//交换地址，在每次迭代开始时，将上一次迭代的结果存储在A_tmp指向的数组里面
void Swap(Type ***A, Type ***B)
{
    Type **tmp = (*A);
    (*A) = (*B);
    (*B) = tmp;
}

加上中文注释