//2次元ラプラス方程式を差分法で解く
//OpenMP版
#include <stdio.h>
#include <string.h>
#include <sys/time.h>
#include <math.h>
//#include <omp.h>

#define SIZE	(1024*3)
#define EPS		1.0e-1

double u[SIZE+2][SIZE+2], uu[SIZE+2][SIZE+2];

double second()
{
	struct timeval tv;
	gettimeofday(&tv, NULL);
	return tv.tv_sec + tv.tv_usec / 1.0e6;
}

void init(double u[SIZE+2][SIZE+2])
{
	int i, j;
	
	memset(u, 0, sizeof(double) * (SIZE + 2) * (SIZE + 2));
	
	for(i=1; i<=SIZE; i++)
		for(j=1; j<=SIZE; j++)
			u[i][j] = sin((i - 1) / (SIZE * M_PI)) + cos((j - 1) / (SIZE * M_PI));
}

int main()
{
	double start, time, err;
	int	i, j, iter=0;
	
	init(u);	init(uu);
	
	start = second();
#pragma omp parallel private(i, j)
	do {
#pragma omp master
		{ iter++; }
#pragma omp for
		for(i=1; i<=SIZE; i++)
			for(j=1; j<=SIZE; j++)
				uu[i][j] = u[i][j];
#pragma omp for nowait
		for(i=1; i<=SIZE; i++)
			for(j=1; j<=SIZE; j++)
				u[i][j] = (uu[i-1][j] + uu[i+1][j] + uu[i][j-1] + uu[i][j+1]) / 4.0;
#pragma omp single
		{ err = 0.0; }
#pragma omp for reduction(+:err)
		for(i=1; i<=SIZE; i++)
			for(j=1; j<=SIZE; j++)
				err += (uu[i][j] - u[i][j]) * (uu[i][j] - u[i][j]);
	} while(err > EPS);
	time = second() - start;
	
	printf("time=%f[sec], iteration=%d, perf=%f[MFLOPS]\n",
			time, iter, (8.0 * SIZE * SIZE * iter) / time / 1.0e6);
	
	return 0;
}
//EOF