#include <omp.h>
#include <stdio.h>
#include <stdlib.h>

#define MAXT 20
#define CT 0.335

#define END 20
#define N  5000
float mt_1[N][N]; // matriz com o estado da chapa em t-1
float mt[N][N];   // idem em t

int main (int argc, char *argv[]) 
{
	int nthreads, i, j, step;
	FILE *f;
	float *pt_1;
	float *pt;
 
	nthreads = strtol(argv[1], NULL, 10);
	# pragma omp parallel for num_threads(nthreads)
	for(i=0; i<N; i++) {
		for(j=0; j<N; j++)
			if ( (i== 0) || (j== 0) || ( i==N-1) || (j==N-1) ) 
				mt_1[i][j] = 200;
			else 
				mt_1[i][j] = 20;

	}
	pt_1 = &mt_1[0][0]; 
	pt = &mt[0][0];
	# pragma omp parallel num_threads(nthreads) default(none) \
				private(i,j,step) shared(pt_1,pt,mt_1,mt)
	for(step=0; step < END; step++) {
		#pragma omp for
		for(i=1; i<N-1; i++) {
			for(j=1; j<N-1; j++)
				pt[i*N+j]= CT*(pt_1[(i-1)*N+j] + pt_1[(i+1)*N+j] + pt_1[i*N+j+1] + pt_1[i*N+j-1])/4;
		}
		# pragma omp master
		{
			float *tmp;
  			tmp = pt_1;
  			pt_1 = pt;
  			pt = tmp;
  		}
  	}
  	
	f= fopen("simulacao.bin", "w");
	fwrite(pt_1, sizeof(float), N*N, f);
	fclose(f);

	return 0;
 }