/* paralleladder.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 */

/* A routine performing elemental work.
This can be replaced with a much larger routine. */
double kernelroutine(double input); 
double kernelroutine(double input) 
	{
	return (input+1)*(input+1); 
	}

/* HighestIndex specifies how highest index to sample */
#define HighestIndex			10000L

void computeloop(double *theSum, double **theArray, int idproc, int nproc); 
void computeloop(double *theSum, double **theArray, int idproc, int nproc)
	{
	/* local copies of the data to be output */
	double *myArray, mySum=0; 
	/* limit of the loop */
	long loopEnd=(HighestIndex+nproc-1)/nproc; /* just this proc's piece */
	/* this processor's index offset */
	long offset=idproc*loopEnd; 
	
	/* allocate an array to hold all the results */
	myArray=malloc(sizeof(double)*loopEnd); 
	
	if (myArray) {
		long index; 
		
		/* loop over indicies */
		for(index=0; index<loopEnd; index++) {
			/* call the kernel routine for each index, and save into the array */
			myArray[index]=kernelroutine(index+offset); 
			/* sum as we go */
			if (index+offset<HighestIndex) /* limit to the desired indicies */
				mySum+=myArray[index]; 
			
			}
		
		/* proc 0 needs to hold the entire array */
		{double *bigArray=idproc?NULL:malloc(sizeof(double)*loopEnd*nproc);
		
		/* gathers the data from the other arrays ... */
		MPI_Gather(myArray, loopEnd, MPI_DOUBLE, 
			bigArray, loopEnd, MPI_DOUBLE, /* ... to proc 0 */ 0, MPI_COMM_WORLD); 
		
		if (!idproc) {
			free(myArray); 
			myArray=bigArray; 
			}
		}
		
		/* performs a parallel sum across processors and saves the result ... */
		MPI_Reduce(&mySum, &mySum, 1, MPI_DOUBLE, MPI_SUM, /* ... at proc 0 */ 0, MPI_COMM_WORLD); 
		
		}
	
	/* return the sum and the array */
	*theSum=mySum; 
	*theArray=myArray; 
	}
	
int ppinit(int argc, char *argv[], int *idproc, int *nproc); 
void ppexit(void); 

int main(int argc, char *argv[])
	{
	/* main copies of the sum and the array */
	double theSum, *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(&theSum, &theArray, 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(double), HighestIndex, fp); 
				
				fclose(fp); 
				}
			
			printf("the sum is %f\n", theSum); 
			
			}
		
		/* 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;

/* terminate MPI execution environment */
	ierr = MPI_Finalize();
	
	}