/*
Author: Eric Thibodeau
Date:	19/07/2006
Course: MGL810
Subject: MPI Parallel calculation of a string vibration
MPI reference: MPI The Complete Reference, Volume 1, The MPI Core, second edition.
MPI_SENDRECV() --> p66.
*/
#include <X11/Xlib.h>
#include <X11/Xutil.h>
#include <X11/Xos.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <stdlib.h>
#include <sys/times.h>
#include <time.h>
#include <mpi.h>

#define X_RESN   800		/* largeur de la fenetre par defaut */
#define Y_RESN   200		/* heuteur de la fenetre par defaut */
#define MAX_ITER 100		/* max iterations par defaut */
				
#define REAL_MIN  0
#define IMAG_MIN -1
#define REAL_MAX  X_RESN
#define IMAG_MAX  1
#define WHITE     16777215
#define BLACK     0

#define timeadd(a,b) a.tv_sec+=b.tv_sec;\
					if ((a.tv_usec+b.tv_usec) > 1000000)\
					{\
						a.tv_usec+=b.tv_usec-1000000;\
 						a.tv_sec++;\
					}\
 					else \
						a.tv_usec+=b.tv_usec;

#define time2double(a) ((double) a.tv_sec + (double) a.tv_usec/1000000)
//TIC_START tags timeofday into start
#define TIC_START(start) gettimeofday(&start,NULL);

//TIC_STOP tags time difference between start TAG and _now_ and adds it to t_store
#define TIC_STOP(t_store,start) gettimeofday(&t_stop,NULL);\
						timeval_subtract(&t_diff,&t_stop,&t_start);\
						timeadd(t_store,t_diff);



/* timeval_subtract extracted from http://www.gnu.org/software/libc/manual/html_mono/libc.html#Elapsed%20Time*/

/* Subtract the `struct timeval' values X and Y,
        storing the result in RESULT.
        Return 1 if the difference is negative, otherwise 0.  
		result = x - y*/

int timeval_subtract (result, x, y)
struct timeval *result, *x, *y;
{
	/* Perform the carry for the later subtraction by updating y. */
	if (x->tv_usec < y->tv_usec) {
		int nsec = (y->tv_usec - x->tv_usec) / 1000000 + 1;
		y->tv_usec -= 1000000 * nsec;
		y->tv_sec += nsec;
	}
	if (x->tv_usec - y->tv_usec > 1000000) {
		int nsec = (x->tv_usec - y->tv_usec) / 1000000;
		y->tv_usec += 1000000 * nsec;
		y->tv_sec -= nsec;
	}

       /* Compute the time remaining to wait.
	tv_usec is certainly positive. */
	result->tv_sec = x->tv_sec - y->tv_sec;
	result->tv_usec = x->tv_usec - y->tv_usec;

	/* Return 1 if result is negative. */
	return x->tv_sec < y->tv_sec;
}

/*
Calcul du déplacement d'une section de corde délimité par les positions en X par
les paramètres start et stop.
Variables d'entrée:
C = Constante d'élasticité de la corde
length = La longueur de la corde (nb de pts sur la corde)
Variables modifiés:
Corde
Corde_h1
Corde_h2
*/
void calc_point(double * Corde, double * Corde_h1, double * Corde_h2, unsigned int length, double C)
{
	unsigned int i;
	double Pos_avant2;
	for(i=1; i<length; i++) //on a besoin de i±1 pour le calcul, 
							//donc commencer à 1 et arrêter à length -1
	{
		Pos_avant2= 2.0 * Corde_h1[i];
		Corde[i]=Pos_avant2-Corde_h2[i]+( C * ( Corde_h1[i-1] - Pos_avant2 + Corde_h1[i-1] ));
	}
	return ;
}

/*
Calcul de la moyenne du voisinnage pour permettre la requantification de la corde pour permettre
l'affichage d'une corde plus grande que ce qui peut etre affiché dans la fenêtre présente
*/
/*double * requant(const double *Corde,
				 double *Corde_requant,
				 const unsigned int length,
				 const unsigned int newlength,
				 const unsigned int window)
{
	unsigned int i,j;
	for(i=1; i<length; i++) //on a besoin de i±1 pour le calcul, 
							//donc commencer à 1 et arrêter à length -1
	{
		for (j=0;j<window;j++)
		{
			Corde_requant[i]=
		}
		Pos_avant2=2.0*Corde_h1[i];
		Corde[i]=Pos_avant2-Corde_h2[i]+( C * ( Corde_h1[i-1] - Pos_avant2 + Corde_h1[i-1] ));
	}
	return ;
}*/

/* ---- Fonction principale ------------------------------------------------ */

int main (int argc, char *argv[])
{

	/*Variable pour l'afffichage*/
	Display *display;
	Window win;			/* initialisation de la fenetre */
	GC gc;
	int GUI;			/* Oui ou non pour l'affichage*/
	
	/* Variables partagees*/
	
	/*Variables pour MPI*/
	int rank,size;
	
	unsigned int *DispCorde; 	/*Buffer d'affichage*/
	
	unsigned long  min_color, max_color; 	/* valeurs min et max values des couleurs */
	float scale_x, scale_y;					/* facteurs d'ajustement pour mapper la region a la fenetre */
	unsigned int width = 0, height = 0;		/* taille de la fenetre */

	/* Variables de la corde  */
	unsigned long Length, Window, DispLength;
	
	int maxiter = 0;		/* maximum d'iterations */
		
	/* Autres variables */
	
	char dummy;		/* faire un scanf pour inserer une pause avant de quiter */
	struct tms t;	/* Pour le calcul du temps d'exeution */
	struct timeval t_start,t_stop,t_idle,t_start_c,t_stop_c,t_diff,t_comm,t_disp;
	
	MPI_Init(&argc, &argv);
	MPI_Comm_size(MPI_COMM_WORLD, &size);
	MPI_Comm_rank(MPI_COMM_WORLD, &rank);
	
	t_comm.tv_sec = t_comm.tv_usec = 0;
	t_disp.tv_sec = t_disp.tv_usec = 0;
	t_idle.tv_sec = t_idle.tv_usec = 0;
	
	if(!rank)
	{
		/*Print data header immediately so it shows up on top*/
		fprintf(stdout,"Proc#:\tComm:\t\tDisplay:\tIdle:\t\tCPU_time:\tTotal:\n");
	}
		
	if(size < 2)
	{
		fprintf (stderr, "The Manager refuses to work alone. Start MPI with at least 2 processes!\n");
		MPI_Finalize();
		return 1;
	}

	if (rank == 0)
	{
		unsigned long  white=WHITE, black=BLACK;	/* valeurs des pixels blancs et noirs */
		
		unsigned int    x, y,						/* position de la fenetre */
						border_width,				/* largeur des bords en pixels */
						disp_width, disp_height,	/* taille de l'ecran */
						screen;						/* ecran */
	
		char *window_name = "Mandelbrot", *disp_name = NULL;
		unsigned long valuemask = 0;
		XGCValues values;
		XSizeHints size_hints;
	
		XSetWindowAttributes attr[1];
			
	
		/* Les arguments de la commande s'il y en a */
		/* Utilisation:
			string-mpi [maxiter] [width] [height]
			ou
			string-mpi -h pour afficher le message
		*/
	
		if ((argc > 1) && (strcmp(argv[1], "-h") == 0)) {
			fprintf (stderr, "Usage:  mandelbrot [maxiter] [width] [height] [String length] [GUI 0/1/2]\n");
			exit (EXIT_FAILURE);
		}
		if (argc > 1) 
			maxiter = atoi(argv[1]);
		if (argc > 2)
			width = atoi(argv[2]);
		if (argc > 3)
			height = atoi(argv[3]);
		if (argc > 4)
			length = atoi(argv[4]);
		if (argc > 5)
			GUI = atoi(argv[5]);
		
		maxiter = (maxiter > 0) ? maxiter : MAX_ITER;
		width   = (width   > 0) ? width   : X_RESN;
		height  = (height  > 0) ? height  : Y_RESN;
		
		if (GUI)
		{	/* Connection a Xserver */
			if ( (display = XOpenDisplay (disp_name)) == NULL ) {
				fprintf (stderr, "drawon: cannot connect to X server %s\n",
					XDisplayName (disp_name) );
				exit (EXIT_FAILURE);
			}
				
			/* Initialisation de l'affichage graphique */
		
			screen = DefaultScreen (display);
			disp_width = DisplayWidth (display, screen);
			disp_height = DisplayHeight (display, screen);
		
			x = 0;  y = 0;		/* position de la fenetre  */
		
			border_width = 4;
			win = XCreateSimpleWindow (display, RootWindow (display, screen),
						x, y, width, height, border_width, 
						BlackPixel (display, screen), 
						WhitePixel (display, screen));
		
			size_hints.flags = USPosition|USSize;
			size_hints.x = x;
			size_hints.y = y;
			size_hints.width = width;
			size_hints.height = height;
			size_hints.min_width = 300;
			size_hints.min_height = 300;
		
			XSetNormalHints (display, win, &size_hints);
			XStoreName(display, win, window_name);
		
			/* contexte du graphique */	
			gc = XCreateGC (display, win, valuemask, &values);
			
			/* les valeurs des couleurs blanc et noir */
			white = WhitePixel (display, screen);
			black = BlackPixel (display, screen);
			fprintf (stderr,"white: %lu black: %lu\n",white,black);
			
			XSetBackground (display, gc, white);
			XSetForeground (display, gc, black);
			XSetLineAttributes (display, gc, 1, LineSolid, CapRound, JoinRound);
		
			attr[0].backing_store = Always;
			attr[0].backing_planes = 1;
			attr[0].backing_pixel = white;
		
			XChangeWindowAttributes(display, win, 
					CWBackingStore | CWBackingPlanes | CWBackingPixel,
					attr);
			XMapWindow (display, win);
			XSync(display, 0);
		}

		scale_x = (float) (REAL_MAX - REAL_MIN) / (float) width;
		scale_y = (float) (IMAG_MAX - IMAG_MIN) / (float) height;*/
	}
	TIC_START(t_start)
	MPI_Bcast( &width     , 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
	MPI_Bcast( &height    , 1, MPI_UNSIGNED_LONG, 0, MPI_COMM_WORLD);
	MPI_Bcast( &maxiter   , 1, MPI_INT  , 0, MPI_COMM_WORLD);
	TIC_STOP(t_comm,t_start)
	/* Calcul des facteurs d'ajustement	*/
	
	scale_x = (float) (REAL_MAX - REAL_MIN) / (float) width;
	scale_y = (float) (IMAG_MAX - IMAG_MIN) / (float) height;
	
	/* Calculer et afficher les points */
	{
		unsigned int maxmin;
		MPI_Status	status,rstatus;
		MPI_Request sreq;

		/*Start of actual processing*/
		if(!rank)
		{
			int NotDone=1;
			double * Corde;
			
			Corde = (double*) malloc(length*sizeof(double));
			if(Corde == NULL)
			{
				MPI_Finalize();
				printf("Corde Malloc error, exiting\n");
				return 1;
			}
			do
			{
				/*All slaves start y calculating the first column and sending the result, this "start the ball rolling"*/
				/*int MPI_Probe (int source, int tag, MPI_Comm comm, MPI_Status *status);*/
				TIC_START(t_start)
				MPI_Probe(MPI_ANY_SOURCE,MPI_ANY_TAG,MPI_COMM_WORLD,&status); /*Wait for an available process*/
				TIC_STOP(t_idle,t_start)
				
				TIC_START(t_start)
				MPI_Recv(Corde,length,MPI_DOUBLE,status.MPI_SOURCE,status.MPI_TAG,MPI_COMM_WORLD,&rstatus);
				TIC_STOP(t_comm,t_start)

				/*we send out work immediately so that the display processing can be performed in parallel
				ideally, these would be independant threads*/
				if(i<width)
				{
					i++;
					DataSentOut++;
				}
				else if(i==width)
				{
					i=width+1; /*we are explicitally going over the width bounds to express the termination of the calculation proces*/
				}
				
				TIC_START(t_start)
				MPI_Isend(&i,1,MPI_UNSIGNED_LONG,status.MPI_SOURCE,i,MPI_COMM_WORLD,&sreq);
				TIC_STOP(t_comm,t_start)
				
				if(GUI)
				{
					TIC_START(t_start)
					for (k=0;k<=height;k++)
					{
						XSetForeground (display, gc, DispBuff[k]);
						XDrawPoint (display, win, gc, status.MPI_TAG, k);
					}
					TIC_STOP(t_disp,t_start)
				}
				
			}while( NotDone );
		}
		else
		{
			double * Corde, * Corde_h1, * Corde_h2, ptrDbl;
			unsigned int i=0,j=0;
			int right, left;
			
			length=length/size+1;
			if(rank==1)
			{
				left = MPI_PROC_NULL;
				right= rank+1;
			}
			else if (rank==size)
			{
				right= MPI_PROC_NULL;
				left = rank-1;
			}
			else
			{
				length++;
				left = rank-1;
				right= rank+1;
			}
			Corde    = (double*) malloc(sizeof(double)*length);
			Corde_h1 = (double*) malloc(sizeof(double)*length);
			Corde_h2 = (double*) malloc(sizeof(double)*length);
			if( Corde == NULL || Corde_h1 == NULL || Corde_h2 == NULL )
			{
				MPI_Finalize();
				printf("Corde Malloc error, exiting\n");
				return 1;
			}
			else
			{
				//Initialisation des Cordes
				for (i=0;i<=length;i++) Corde[i] =0.0;
				memcpy((void*)Corde_h1,(void*)Corde,length);
				memcpy((void*)Corde_h2,(void*)Corde,length);
			}
			while(i<=maxiter)
			{
				//Excitation point: dirac
				if(rank == 1 && i == 1)
				{
					Corde[0]=1.0;
				}
				//calc_point(double * Corde, double * Corde_h1, double * Corde_h2, unsigned int length, double C)
				calc_point(Corde, Corde_h1, Corde_h2, length, 0.5);
				
				if(GUI || i == maxiter)
				{
					TIC_START(t_start)
					MPI_Isend(Corde,length,MPI_DOUBLE,0,rank,MPI_COMM_WORLD,&sreq);
					TIC_STOP(t_comm,t_start);
				}
				//send/receive ovelapping points to neighbors
				TIC_START(t_start);
				//right send-receive
				MPI_SENDRECV(Corde[length-1],1,MPI_DOUBLE,right,0,
							 Corde[length]  ,1,MPI_DOUBLE,right,0,
							 MPI_COMM_WORLD,&sreq);
				//left send-receive
				MPI_SENDRECV(Corde[length-1],1,MPI_DOUBLE,left,0,
							 Corde[length]  ,1,MPI_DOUBLE,left,0,
							 MPI_COMM_WORLD,&sreq);
				TIC_STOP(t_comm,t_start);
				
				//switch historical strings around:
				ptrDbl  = Corde_h2;
				Corde_h2= Corde_h1;
				Corde_h1= Corde;
				Corde   = ptrDbl;
				i++;
			}
		}
	}

	if (!rank) /*if Manager*/
	{
		if (GUI)
		{
			XFlush (display);
			if (GUI > 1) do{
				// Attendre la reponse de l'usager, sortir apres
					fprintf (stderr, "Taper q pour quitter le programme.\n");
					fscanf (stdin, "%c", &dummy);
			}while (dummy != 'q');
		}
	}
	free(DispBuff);
	
	/* Calculer le temps d'execution et produire le texte de sortie */
	times(&t);

	{
		double CPU_time=(double) (t.tms_utime + t.tms_stime + t.tms_cutime + t.tms_cstime) / CLK_TCK;
		fprintf(stdout,"%i\t%f\t%f\t%f\t%f\t%f\n",rank,
				time2double(t_comm),
				time2double(t_disp),
				time2double(t_idle),
				CPU_time,
				time2double(t_comm)+time2double(t_disp)+time2double(t_idle)+CPU_time);
	}
	
		
	MPI_Finalize();
	return EXIT_SUCCESS;
}