/* mpi_max.c SJ */
/* Example program for maximum finding */
/* Compile: mpicc -lm -o mpi_max [-DVERBOSE] mpi_max.c */
/* Usage mpirun -np 00 mpi_max -N input_size -K #of_iter */

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <math.h>
#include <limits.h>
#include <mpi.h>


int reduceMax(int *bl, int bs,  MPI_Comm comm);
int dcMax(int *input, int n, int pid, int p, MPI_Comm comm);


int main(int argc, char **argv)
{
	int PID,  /* process-ID */
		P;    /* number of processes */
	/* MPI_Status tmpstat; */

	int N = 10000;
	int k, K = 10;
	int blockSize;
	int *block = NULL, *input = NULL;

	int result, a;

	double start_ustime, used_ustime;

	MPI_Init(&argc, &argv);
	MPI_Comm_rank(MPI_COMM_WORLD, &PID);
	MPI_Comm_size(MPI_COMM_WORLD, &P);

	/* read command line parameters */

	for (a = 1; a < argc; a++) {
		if (! strcmp(argv[a], "-N") && argv[a+1] != NULL &&
			atoi(argv[a+1]) > 0)
			N = atoi(argv[a+1]);

		if (! strcmp(argv[a], "-K") && argv[a+1] != NULL &&
			atoi(argv[a+1]) > 0)
			K = atoi(argv[a+1]);
	}

	#ifdef VERBOSE
	if (PID == 0) printf("N = %d, P = %d, K = %d.\n", N, P, K);
	#endif

	/* local part of input */
	blockSize = (int)ceil((double)N/P);
	N = P*blockSize; /* make input size divisible */
	block = (int*)malloc(blockSize * sizeof(int));

	/* generate input */
	if (PID == 0) {
		int i, r, m = -1;
		input = (int*)malloc(N * sizeof(int));

		for (i = 0; i < N; i++) {
			r = rand();
			if (r > m)
				m = r;
			input[i] = r;
		}
		printf("Input: N=%d, P=%d, max=%d.\n", N, P, m);
	}

	start_ustime = MPI_Wtime();

	/* distribute input to all processes */
	for (k = 0; k < K; k++)
		MPI_Scatter(input, blockSize, MPI_INT,
					block, blockSize, MPI_INT,
					0, MPI_COMM_WORLD);
	
	used_ustime = MPI_Wtime() - start_ustime;
	if (PID == 0) printf("Scatter in %f us.\n",
						 1000000.0 * used_ustime / K);

	MPI_Barrier(MPI_COMM_WORLD);

	/* test maximum */
	start_ustime = MPI_Wtime();

	for (k = 0; k < K; k++)
		result = reduceMax(block, blockSize, MPI_COMM_WORLD);
	
	used_ustime = MPI_Wtime() - start_ustime;
	if (PID == 0) printf("reduceMax=%d in %.1f us.\n",
						 result, used_ustime * 1000000.0 / K);

	MPI_Barrier(MPI_COMM_WORLD);

	/* test another maximum */
	start_ustime = MPI_Wtime();

	for (k = 0; k < K; k++)
		result = dcMax(input, N, PID, P, MPI_COMM_WORLD);

	used_ustime = MPI_Wtime() - start_ustime;
	if (PID == 0) printf("dcMax=%d in %.1f us.\n",
						 result, used_ustime * 1000000.0 / K);



	free(block);
	if (PID == 0) free(input);
	MPI_Finalize();

	return 0;

}	/* main() */

/* maximum locally, then use MPI_Allreduce() */
int reduceMax(int *bl, int bs,  MPI_Comm comm) {

	int result, i, lmax = INT_MIN;
	for (i = 0; i < bs; i++)
		if (bl[i] > lmax)
			lmax = bl[i];

	MPI_Allreduce(&lmax, &result, 1, MPI_INT, MPI_MAX, comm);

	return result;
}


/* maximum in divide&conquer, divides processes to half,
   new 0 send half of data to new 1 */
/* Not really efficient algorithm, just to demonstrate recursive
   MPI program */

/* input only in pid 0, result returned only by pid0, others return 0 */
int dcMax(int *input, int n, int pid, int p, MPI_Comm comm) {

	int half;
	int *rcvbuf = NULL;

	MPI_Status tmpstat;
	MPI_Comm newcomm;
	int newpid, newp, result;

	if (p == 1 && pid == 0) {
		/* one process only */
		int i, lmax = INT_MIN;

		for (i = 0; i < n; i++)
			if (input[i] > lmax)
				lmax = input[i];

		#ifdef VERBOSE
		int gpid;
		MPI_Comm_rank(MPI_COMM_WORLD, &gpid);
		printf("Process %d has %d inputs, max = %d.\n", gpid, n, lmax);
		#endif

		return lmax;
	} else {
		if (n == 0 || p == 1)
			/* no more input (for me) */
			return INT_MIN;
	}
	
	/* n > 1, p > 1 */

	half = n/2;

	/* pid0 sends half of input to pid1 */
	if (pid == 0)
		MPI_Send(input + half, n - half, MPI_INT, 1, 0, comm);
	else if (pid == 1) {
		rcvbuf = (int*)malloc(sizeof(int)*(n - half));
		MPI_Recv(rcvbuf,  n - half, MPI_INT, 0, 0, comm, &tmpstat);
	}

	/* all processes split to 2 groups, pid0 and pid1 will lead
	   the new groups */

	MPI_Comm_split(comm, pid%2, 0, &newcomm);

	/* get new p, pid */
	MPI_Comm_rank(newcomm, &newpid);
	MPI_Comm_size(newcomm, &newp);

	/* make recursive call according to the group */
	if (pid%2 == 0)
		/* group of the old pid0 */
		result = dcMax(input, n/2, newpid, newp, newcomm);
	else
		/* group of the old pid1 */
		result = dcMax(rcvbuf, n - half, newpid, newp, newcomm);

	/* combine the results, old pid1 sends result to pid0 */
	if (pid == 1) {
		MPI_Send(&result, 1, MPI_INT, 0, 0, comm);
		result = 0;
		free(rcvbuf);
	} else if (pid == 0) {
		int result1;
		MPI_Recv(&result1, 1, MPI_INT, 1, 0, comm, &tmpstat);
		if (result1 > result)
			result = result1;
	} else
		result = 0;	/* others won't need result */

	return result;

} /* dcMax() */