/* SolarSim.c */
#include "SolarSim.h"

typedef struct task {
	int i, j;
} Task, *TaskPtr;

int sock = 0;
long SIM_STEPS, DT, SAVE;

void
do_work(AState state) {
	int tid, nthreads;
	long i, j, t, ntasks, sim_step;
	int n = state->particles;
	vector *r = state->positions;
	vector *v = state->velocities;
	vector *a = state->accelerations;
	vector *ac;
	real *m = state->masses;
	TaskPtr tasks;

	tasks = malloc( n*(n-1)/2 * sizeof(Task) );
	t = 0;
	for( i = 0; i < n-1; i++ )
		for( j = i+1; j < n; j++ ) {
			tasks[t].i = i;
			tasks[t].j = j;
			t++;
		}
	ntasks = t;

	#pragma omp parallel shared(nthreads)
	#pragma omp master
	nthreads = omp_get_num_threads();

	ac = malloc(nthreads * n * sizeof(vector));
	
	#pragma omp parallel shared(r,v,a,ac,m,n,tasks,ntasks,nthreads,DT,sim_step,SIM_STEPS,SAVE) private(i,j,t,tid) default(none)
	sim_step = 0;
	while( sim_step < SIM_STEPS ) {
		if( (SAVE != 0) && (sim_step % SAVE == 0) ) {
			write_state(sock, state);
		}
		
		#pragma omp for schedule(static,n) //shared(ac,n,nthreads) private(i) default(none) nowait
		for( i = 0; i < n*nthreads; i++ ) {
			ac[i].x = 0;
			ac[i].y = 0;
			ac[i].z = 0;
		}

		// compute accelerations. each thread uses ac[tid] vector
		#pragma omp for schedule(static) //shared(tasks, ntasks, r, v, ac, m, n) private(i,j,t,tid) default(none) nowait
		for( t = 0; t < ntasks; t++ ) {
			i = tasks[t].i;
			j = tasks[t].j;
			tid = omp_get_thread_num();
			compute_acc(r+i, r+j, v+i, v+j, ac+tid*n+i, ac+tid*n+j, m[i], m[j]);
		}

		// add up the accelerations computed by each thread
		#pragma omp parallel for schedule(static,1) //shared(a,ac,n,nthreads) private(i,j) default(none) nowait
		for( i = 0; i < n; i++ ) {
			v_init(a+i,0,0,0);
			for( j = 0; j < nthreads; j++ ) {
				a[i].x += ac[j*n+i].x;
				a[i].y += ac[j*n+i].y;
				a[i].z += ac[j*n+i].z;
			}
		}

		// update positions and speeds
		#pragma omp parallel for schedule(static,1) //shared(r,v,a,n,DT) private(i) default(none) nowait
		for( i = 0; i < n; i++ ) {
			update_pos_vel(r+i, v+i, a+i, DT);
		}

		#pragma omp master
		sim_step ++;

		#pragma omp barrier
	}
	if( SAVE != 0 )
		write_state(sock, state);
}

int
main( int argc, char **argv ) {
	TState state;
	double before, after;
	int i, n=8;

	SIM_STEPS = atol(argv[3]);
	DT = atol(argv[4]);
	SAVE = atol(argv[5]);

	if( SAVE != 0 ) {
		sock = open( argv[2], O_WRONLY | O_CREAT | O_TRUNC, 644 );
	}

//	read_initial( &state, argv[1] );
	rand_data( &state, 100 );

	before = omp_get_wtime();
	do_work(&state);
	after = omp_get_wtime();

	if( SAVE != 0 )
		close(sock);

	#pragma omp parallel
	#pragma omp master
	printf("%d %lf\n", omp_get_num_threads(), after-before);

	return EXIT_SUCCESS;
}