//2次元ラプラス方程式を差分法で解く
//Pthread版（問題サイズをスレッド数で割り切れることが前提）
#include <stdio.h>
#include <string.h>
#include <sys/time.h>
#include <math.h>
#include <pthread.h>

#define SIZE		(1024*3)
#define EPS			1.0e-1
#define THREAD_NUM	2	//スレッドの数(1以上指定)
#define SPLITS	(int)(SIZE/THREAD_NUM)

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

typedef struct _thread_arg
{
	int thread_no;
	int myStart;
	pthread_mutex_t *mutex;
}
thread_arg_t;

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));
}

void thread_copy(void *arg)
{
	int i, j;
	
	thread_arg_t* targ = (thread_arg_t *)arg;
	
	for(i=targ->myStart+1; i<=targ->myStart+SPLITS; i++)
		for(j=1; j<=SIZE; j++)
			uu[i][j] = u[i][j];
}

void thread_upd_err(void *arg)
{
	int i, j;
	double lerr;
	
	thread_arg_t* targ = (thread_arg_t *)arg;
	
	for(i=targ->myStart+1; i<=targ->myStart+SPLITS; 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;
	
	lerr = 0.0;
	for(i=targ->myStart+1; i<=targ->myStart+SPLITS; i++)
		for(j=1; j<=SIZE; j++)
			lerr += (uu[i][j] - u[i][j]) * (uu[i][j] - u[i][j]);
	
	pthread_mutex_lock(targ->mutex);
	err += lerr;
	pthread_mutex_unlock(targ->mutex);
}

int main()
{
	double start, time;
	int	i, iter=0;
	
	pthread_t handle[THREAD_NUM];
	thread_arg_t targ[THREAD_NUM];
	pthread_mutex_t mutex;
	
	if( SIZE%THREAD_NUM != 0 ){
		printf("問題サイズをスレッド数で割り切れること\n");
		return 1;
	}
	
	for(i=0; i<THREAD_NUM; i++) {
		targ[i].myStart = SPLITS * i;
		targ[i].mutex = &mutex;
	}
	
	pthread_mutex_init(&mutex, NULL);
	
	init(u);	init(uu);
	
	start = second();
	
	do {
		iter++;
		
		for(i=0; i<THREAD_NUM; i++)
			pthread_create(&handle[i], NULL, (void*)thread_copy, (void*)&targ[i]);
		
		for(i=0; i<THREAD_NUM; i++)
			pthread_join(handle[i], NULL);
		
		err = 0.0;
		for(i=0; i<THREAD_NUM; i++)
			pthread_create(&handle[i], NULL, (void*)thread_upd_err, (void*)&targ[i]);
		
		for(i=0; i<THREAD_NUM; i++)
			pthread_join(handle[i], NULL);
	} 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