/*                                                                      */
/* CS 4331: Parallel Programming                        Fall, 2008      */
/*                                                                      */
/* This UPC program computes heat flow in a 1-dimensional wire.         */

#include <stdio.h>
#include <math.h>
#include <time.h>
#include "upc.h"
#include "upc_collective.h"

#define  MAXN (1 << 3)		/* maximum vector length	*/

    /* shared objects must be declared at file scope */

    shared [8] double t[8*THREADS];
    shared double epsilon, my_epsilon, max_epsilon;

int main()

{
    int i, j, s;			
    int base;

    double tmp[2];
    double *tme;
    double e, etmp;

    clock_t start_time, end_time;	/* "wallclock" times		*/

    tme = (double *)&t[MAXN*MYTHREAD];
    for (i=0; i<MAXN; ++i)
	tme[i] = 0.0;

    upc_barrier;
  
    if (MYTHREAD == 0)
    {   t[0] = 100.0;
	t[MAXN*THREADS-1] = 100.0;
	epsilon = 0.1;
    }

	max_epsilon = 1.0;

	upc_barrier;

	/* start the clock */
	start_time = clock();

	base = MYTHREAD*MAXN;

	s = 0;

	/* while (max_epsilon > epsilon) */
	for (j=0; j<9; ++j)
	{
		if (MYTHREAD == 0)
		{	tmp[0] = t[0];
			e = 0.0;
		}
		else
		{	tmp[0] = (t[base-1] + t[base] + t[base+1]) / 3.0;
			e = fabs(t[base]-tmp[0]);
		}
		for (i=base+1; i<base+MAXN-1; ++i)
		{
			tmp[1] = (t[i-1] + t[i] + t[i+1]) / 3.0;
			etmp = fabs(t[i]-tmp[1]);
			if (etmp > e) e = etmp;
			t[i-1] = tmp[0];
			tmp[0] = tmp[1];
		}
		if (MYTHREAD < THREADS-1)
		{
			tmp[1] = (t[base+MAXN-2] + t[base+MAXN-1] + t[base+MAXN]) / 3.0;
			etmp = fabs(t[base+MAXN-1]-tmp[1]);
			if (etmp > e) e = etmp;
			t[base+MAXN-1] = tmp[1];
		}
		upc_barrier;
		t[base+MAXN-2] = tmp[0];

		printf("T=%d  e=%lf\n", MYTHREAD, e);
		/* Convert to shared */
		my_epsilon = e;
		upc_all_reduceD(&max_epsilon, &my_epsilon, UPC_MAX, sizeof(double), 1, NULL, 0);
		if (MYTHREAD==0)
			printf("step %d: E=%lf\n", s, max_epsilon);
		++s;
	}

	upc_barrier;

    if (MYTHREAD == 0)
    {   printf("%d steps  e=%lf\n", s, e);
	for (i=0; i<MAXN*THREADS; ++i)
		printf("%lf ", t[i]);
	printf("\n");
    }

	/* stop the clock and print the time */

  	end_time = clock();

	printf("Length = %d\tElapsed time=%lf\t %lf\n",
		MAXN*THREADS, (end_time - start_time) / (double) CLOCKS_PER_SEC, epsilon);
}