/* parallelpascalstriangle.c */

/* © Copyright 2003 Dauger Research, Inc. 
    
Our lawyers made us say this: 
DISCLAIMER: We provide the following on an "AS IS" basis. Use it at your own risk.  */

#include <stdio.h> /* standard I/O routines */
#include <stdlib.h> /* standard library routines, such as memory */
#include <math.h> /* math libraries */

#include "mpi.h" /* MPI library */

#define myType		unsigned long

/* Routines performing elemental work. */
void addright(myType *inputArray, long size); 
void addright(myType *inputArray, long size) 
	{
	if (inputArray) {
		long i=size-1; 
		while (i--) 
			inputArray[i+1]+=inputArray[i]; 
		}
	}

void addleft(myType *inputArray, long size); 
void addleft(myType *inputArray, long size) 
	{
	if (inputArray) {
		long i; 
		for(i=0; i<size-1; i++) 
			inputArray[i]+=inputArray[i+1]; 
		}
	}

void iterationloop(myType *theArray, long arraySize, long lastIteration, short verbose, int idproc, int nproc); 
void iterationloop(myType *theArray, long arraySize, long lastIteration, short verbose, int idproc, int nproc)
	{
	
	if (theArray) if (arraySize>0) {
		long index; 
		long leftIDProc=idproc-1, rightIDProc=idproc+1; 
		MPI_Request request; MPI_Status status; 
		
		/* get to know my (processor's) neighbors */
		if (leftIDProc<0) 
			leftIDProc=nproc-1; 
		if (rightIDProc>=nproc) 
			rightIDProc=0; 
		
		/* loop over indicies */
		for(index=0; index<lastIteration; index++) {
			/* call the left or right routine for each index, reusing array */
			if (1&index) {
				addleft(theArray, arraySize); 
				
				MPI_Irecv(&theArray[arraySize-1], 1, MPI_LONG, rightIDProc, 
					index, MPI_COMM_WORLD, &request); 
				MPI_Send(&theArray[0], 1, MPI_LONG, leftIDProc, index, MPI_COMM_WORLD); 
				MPI_Wait(&request, &status); 
				}
			else {
				addright(theArray, arraySize); 
				
				MPI_Irecv(&theArray[0], 1, MPI_LONG, leftIDProc, 
					index, MPI_COMM_WORLD, &request); 
				MPI_Send(&theArray[arraySize-1], 1, MPI_LONG, rightIDProc, index, MPI_COMM_WORLD); 
				MPI_Wait(&request, &status); 
				}
			
			if (verbose) {long i; 
				for(i=0; i<arraySize; i++) {
					if (1&verbose) 
						printf("%d\t", theArray[i]); 
					else 
						printf("%c", theArray[i]&1?'*':' '); 
					}
				printf("\n"); 
				}
			}
		
		}
	
	}
	
void computeloop(myType **theArray, long binomialExponent, int idproc, int nproc); 
void computeloop(myType **theArray, long binomialExponent, int idproc, int nproc)
	{
	myType *myArray; 
	long arraySize=((binomialExponent+nproc)/nproc)+1; 
	
	/* allocate an array to hold all the results */
	myArray=malloc(sizeof(myType)*arraySize); 
	
	if (myArray) {
		long i; 
		
		/* clear the array */
		for(i=arraySize; i--; )
			myArray[i]=0; 
		
		/* seed the array */
		if (idproc==((binomialExponent>>1)/(arraySize-1))) 
			myArray[(binomialExponent>>1)-idproc*(arraySize-1)]=1; 
		
		/* propagate */
		iterationloop(myArray, arraySize, binomialExponent, 2, idproc, nproc); 
		
		/* proc 0's has a bigger array so we can gather all the data at the end */
		{myType *bigArray=idproc?NULL:malloc(sizeof(myType)*(binomialExponent+nproc)); 
		
		/* gathers the data from the other arrays ... */
		MPI_Gather(myArray, arraySize-1, MPI_LONG, 
			bigArray, arraySize-1, MPI_LONG, /* ... to proc 0 */ 0, MPI_COMM_WORLD); 
		
		if (!idproc) {
			free(myArray); 
			myArray=bigArray; 
			}
		
		}
		
		}
	
	/* return the array */
	*theArray=myArray; 
	}
	
int ppinit(int argc, char *argv[], int *idproc, int *nproc); 
void ppexit(void); 

/* BinomialExponent specifies how far to propagate */

#define BinomialExponent			30L

int main(int argc, char *argv[])
	{
	/* main copies of the sum and the array */
	myType *theArray=NULL; 
	int idproc, nproc, ierr; 
	
	/* initialize parallel processing */
	ierr = ppinit(argc, argv, &idproc, &nproc); 
	
	if (ierr) return ierr; /* stop right there if there's a problem */
	
	printf("I'm processor #%d in a %d-processor cluster.\n", idproc, nproc); 
	
	printf("Beginning computation...\n"); 
	
	computeloop(&theArray, BinomialExponent, idproc, nproc); 
	
	if (theArray) {/* error checking */
		
		if (!idproc) {/* only processor 0 */
			FILE *fp; 
		
			/* save the array into a data file */
			fp=fopen("output", "w"); 
			if (fp) {
				printf("writing array...\n"); 
				
				fwrite(theArray, sizeof(myType), BinomialExponent+1, fp); 
				
				fclose(fp); 
				}
			}
		
		if (1) 	{long i; 
			for(i=0; i<BinomialExponent+1; i++) 
				printf("%d\t", theArray[i]); 
			printf("\n"); 
			}
		
		/* clean up after myself */
		free(theArray); 
		
		}
	else 
		printf("memory allocation failure\n"); 
	
	/* only proc 0 pauses for user exit */
	if (!idproc) {
		printf("press return to continue\n"); 
		getc(stdin); 
		}
	
	ppexit(); 
	
	return 0; 
	}

#ifdef __MWERKS__
/* only for Metrowerks CodeWarrior */
#include <SIOUX.h>
#endif 

int ppinit(int argc, char *argv[], int *idproc, int *nproc)
	{
/* this subroutine initializes parallel processing
   idproc = processor id
   nproc = number of real or virtual processors obtained */
   int ierr;
   *nproc=0; 

/* initialize the MPI execution environment */

   ierr = MPI_Init(&argc, &argv);
   if (!ierr) {

		/* determine the rank of the calling process in the communicator */
		ierr = MPI_Comm_rank(MPI_COMM_WORLD, idproc);

		/* determine the size of the group associated with a communicator */
		ierr = MPI_Comm_size(MPI_COMM_WORLD, nproc);

#ifdef __MWERKS__
/* only for Metrowerks CodeWarrior */
		SIOUXSettings.asktosaveonclose=0; 
		SIOUXSettings.autocloseonquit=1;
#endif 
		
		}
	
	return ierr; 
	}

void ppexit(void)
	{
/* this subroutine terminates parallel processing */
	int ierr;
	
/* waits until everybody gets here */
	ierr = MPI_Barrier(MPI_COMM_WORLD); 
	
/* terminate MPI execution environment */
	ierr = MPI_Finalize();
	
	}