/**
@file "heat_openmp_2.c"
@brief Paralelno resenje Heat Equation pomocu OpenMP. Verzija 2.

       Resenje je znatno jednostavnije nego u slucaju sa MPI, iz prostog razloga
       sto se matrica toplote nalazi u deljenoj memoriji, pa se resenje svodi na 
       trivijalnu podelu posla izmedju niti, tako da svaka nit racuna nove vrednosti
       za odredjeni broj redova matrice.
       
       Za razliku od verzije 1 (fajl heat_openmp_1.c), ova verzija je efikasnija,
	   jer se tim niti-radnika kreira samo jednom, a ne iznova u svakoj iteraciji
	   petlje sto je drasticno kvarilo perfrmanse.	

*/
/******************************************************************************/

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

/******************************************************************************/

/**
@def MIN_MATRIX_DIM
@brief Miminalna dimenzija matrice.
*/
#define MIN_MATRIX_DIM 6

/**
@def MAX_MATRIX_DIM
@brief Maksimalna dimenzija matrice.
*/
#define MAX_MATRIX_DIM 1000

/**
@def MIN_NUM_OF_STEPS
@brief Minimalni broj koraka za izracunavanje konacne toplote
*/
#define MIN_NUM_OF_STEPS 1

/**
@def MAX_NUM_OF_STEPS
@brief Maksimalni broj koraka za izracunavanje konacne toplote
*/
#define MAX_NUM_OF_STEPS 1000000

/**
@def MAX_NUM_OF_THREADS
@brief Maksimalni broj niti
*/
#define MAX_NUM_OF_THREADS 10

/**
@def ErrorDimension
@brief Greska dimenzije matrice
*/
#define ErrorDimension 1

/**
@def ErrorNumOfSteps
@brief Greska u zadatom broju koraka
*/
#define ErrorNumOfSteps 2

/**
@def ErrorFileOpen
@brief Greska prilikom otvaranja fajla za ucitavanje podataka
*/
#define ErrorFileOpen 3

/**
@def ErrorNumerOfInputParameters
@brief Greska prilikom otvaranja fajla za ucitavanje podataka
*/
#define ErrorNumerOfInputParameters 4

/**
@def Cx
@brief Konstanta koja se koristi pri izracunavanju narednost stanja
Pogledati formulu za izracunavanje narednog stanja
*/
#define Cx 0.1

/**
@def Cy
@brief Konstanta koja se koristi pri izracunavanju narednog stanja
Pogledati formulu za izracunavanje narednog stanja
*/
#define Cy 0.1

/*
double test[6][6] = {0, 0, 0, 0, 0, 0,
                     0, 2, 2, 2, 2, 0,
                     0, 2, 8, 8, 2, 0,
                     0, 2, 8, 8, 2, 0,
		     0, 2, 2, 2, 2, 0,
		     0, 0, 0, 0, 0, 0};
*/

/******************************************************************************/

/**
@brief Matrica toplote iz koje se citaju vrednosti u parnoj iteraciji, a upisuju u neparnoj.
*/
static double matrix0[MAX_MATRIX_DIM][MAX_MATRIX_DIM];

/**
@brief Matrica toplote iz koje se citaju vrednosti u neparnoj iteraciji, a upisuju u parnoj.
*/
static double matrix1[MAX_MATRIX_DIM][MAX_MATRIX_DIM];

/**
@brief Pokazivac na matricu iz koje treba citati u tekucoj iteraciji.
*/
double (*readMatrix)[MAX_MATRIX_DIM][MAX_MATRIX_DIM];

/**
@brief Pokazivac na matricu u koju treba upisivati u tekucoj iteraciji.
*/
double (*writeMatrix)[MAX_MATRIX_DIM][MAX_MATRIX_DIM];

/**
@brief Pokazivac na matricu u koju sadrzi rezultat.
*/
double (*resultMatrix)[MAX_MATRIX_DIM][MAX_MATRIX_DIM];


/**
@brief Dimenzija matrice 
*/
unsigned int matrixDimension;

/**
@brief Broj niti
*/
int numOfThreads=1;

/**
@brief Broj iteracija.
*/
int numOfSteps;

/**
@brief Broj iteracija.
*/
int chunkSize;


/******************************************************************************/

/**
@brief Ucitava elemente u matricu iz fajla ili sa sa konzole.
@param pHeatMatrix Matrica u koju je potrebno ucitati elemente.
@param pCopyHeatMatrix Matrica kopija 
@param pDim Dimenzija kvadratne matrice. Na osnovu ovoga parametra se odredjuje broj elemenata koji ce se citati iz fajla.
@param stream Fajl otvoren za citanje, iz kojeg se citaju podaci.
*/
void scanfMatrix(double pHeatMatrix[][MAX_MATRIX_DIM], double pCopyHeatMatrix[][MAX_MATRIX_DIM], int pDim, FILE *stream) {
    int i, j;
    double temp=0;

    for (i = 0; i < pDim; i++) {
        for (j = 0; j < pDim; j++) {
	    
	    //pHeatMatrix[i][j] = pCopyHeatMatrix[i][j] = test[i][j];
	    
	    if (stream==NULL)
		scanf("%lf", &temp);
	    else
		fscanf(stream,"%lf", &temp);

	    pHeatMatrix[i][j] = pCopyHeatMatrix[i][j] = temp;
	    
        }
    }
}

/**
@brief Ispis navedene matrice.
@param pHeatMatrix Matrica koja se ispisuje.
@param pRows Broj vrsta matrice
@param pCols Broj kolona matrice
@param pFormat Format koji ce biti koriscen za ispis svakog elementa matrice
*/
void printfMatrix(double pHeatMatrix[][MAX_MATRIX_DIM], int pRows, int pCols, const char *pFormat) {
    int i, j;
    
    for (i = 0; i < pRows; i++) {
        for (j = 0; j < pCols; j++) {
            printf(pFormat, pHeatMatrix[i][j]);
        }
        printf("\n");
    }
}

/**
@brief Provera korektnosti dimenzije matrice toplote.
@param pMatrixDimension Dimenzija matrice
@return Vraca vrednost != 0 ukoliko je zadata dimenzija nekorektna
*/
void ifWrongDimensionAbort(int pMatrixDimension) {
    if (pMatrixDimension < MIN_MATRIX_DIM || pMatrixDimension > MAX_MATRIX_DIM  || pMatrixDimension%2){
	printf("Pogresne dimenzije matrice! Prekidam rad programa...\n");
	exit(ErrorDimension);
    }
}

/**
@brief Provera korektnosti broja niti. Ako zadati broj niti premasuje maksimalno dozvoljeni, koristice se broj niti jednak raspolozivom broju procesora.
@param pNumOfThreads Broj koraka koje treba sprovesti.
*/
void ifWrongNumOfThreadsDefault(int pNumOfThreads) {
    if (pNumOfThreads > MAX_NUM_OF_THREADS){
	numOfThreads = omp_get_num_procs();
	printf("Zadat prevelik broj niti! ");
    }
    printf("Broj niti= %d\n\n\n",numOfThreads);
}

/**
@brief Provera korektnosti broja koraka za izracunavanje konacnog stanja.
@param pNumOfSteps Broj koraka koje treba sprovesti.
@return Vraca vrednost != 0 ukoliko je zadati broj koraka nekorektan.
*/
void ifWrongNumOfStepsAbort(int pNumOfSteps) {
    if (pNumOfSteps < MIN_NUM_OF_STEPS || pNumOfSteps > MAX_NUM_OF_STEPS){
	printf("Pogresnan broj koraka! Prekidam rad programa...");
	exit(ErrorNumOfSteps);
    }
}

/**
@brief Ispis nakon zavrsetka obrade
*/
void printFooter() {
    printf("\n//////////////////  END EXECUTION  ///////////////////////\n");
}

/**
@brief Ispis informacija tokom jedne iteracije
@param pIteration Broj iteracije
*/
void printIteration(int pIteration) {
    printf("\n_ _ _ _ _ _ _ #%d iteracija _ _ _ _ _ _ _ _\n", pIteration); 
}

/******************************************************************************/

/**
@brief Glavni program. Upotreba heat_omp.exe ulaznifajl num_threads
@param argc Broj argumenata komandne linije
@param argv Argumenti komandne linije
@return Vraca @b 0 ako je sve uspesno izvrseno
*/
int main(int argc, char *argv[]) {
    
    int step, i, j;
    int startRow, endRow, myID;
    double start, stop;
    FILE *input_file;

    printf("Unesite dimenziju matrice: ");
    scanf("%d", &matrixDimension);
    ifWrongDimensionAbort(matrixDimension);
    
    printf("Unesite broj koraka za izracunavanje krajnje toplote: ");
    scanf("%d", &numOfSteps);
    ifWrongNumOfStepsAbort(numOfSteps);

    if (argc==3) {
	numOfThreads = atoi(argv[1]);
	ifWrongNumOfThreadsDefault(numOfThreads);

	input_file = fopen(argv[2],"r");
	if (input_file==NULL) {
	    printf("Error opening file %s\n", argv[1]);   
	    exit(ErrorFileOpen);
	}
	scanfMatrix(matrix0, matrix1, matrixDimension, input_file);
	fclose(input_file);
    }else 
	if(argc!=1) 
	    exit(ErrorNumerOfInputParameters);
	else { /* argc==1 */
	    
	    printf("Unesite toplotu za svaku celiju: \n");
	    scanfMatrix(matrix0, matrix1, matrixDimension, NULL);
	}
    
    //printfMatrix(matrix0, matrixDimension, matrixDimension, "%#6.3g     ");
    //printf("\n\n");

    omp_set_dynamic(0);
    omp_set_num_threads(numOfThreads);

    chunkSize = matrixDimension / numOfThreads;
    
    start = omp_get_wtime();

    #pragma omp parallel firstprivate(readMatrix, writeMatrix) \
			 private(i, j, step, myID, startRow, endRow) \
		         shared(matrix0, matrix1, resultMatrix, chunkSize, numOfSteps, numOfThreads)
    {
	#ifdef _OPENMP
	myID = omp_get_thread_num();
	startRow = (myID==0)?1:(myID*chunkSize); //jer su po ivicama matrice fiksirane vrednosti, kao granicni uslovi
	endRow = (myID==(numOfThreads-1))?(startRow + chunkSize-2):(startRow + chunkSize-1);
	#endif
    #ifndef _OPENMP
	startRow = 1;
	endRow = matrixDimension-2;
	#endif
	
	for (step=0; step<numOfSteps; step++){
	    if (step%2){ /* ako je neparna iteracija */
		readMatrix=&matrix1;
		writeMatrix=&matrix0;
	    }else{
		readMatrix=&matrix0;
		writeMatrix=&matrix1;
	    }
    	
	    for (i=startRow; i<=endRow; i++){
		for(j=1; j<matrixDimension-1; j++){
		
		(*writeMatrix)[i][j] = (*readMatrix)[i][j] + Cx * ((*readMatrix)[i+1][j] + (*readMatrix)[i-1][j] - 2 * (*readMatrix)[i][j]) 
		    + Cy * ((*readMatrix)[i][j+1] + (*readMatrix)[i][j-1] - 2 * (*readMatrix)[i][j]);
		}
	    }//end for petlje po matrici...
	
	    //Bitno! Na ovaj nacin obezbjedjuje se da se potpuno sracuna citava matrica pre prelaska u novu iteraciju...
	    #pragma omp barrier

	}//end for (step=0...
	#pragma omp single
	resultMatrix = writeMatrix;
    }//end pragma omp parallel

    stop = omp_get_wtime();
    printIteration(numOfSteps);

    //printfMatrix((*resultMatrix), matrixDimension, matrixDimension, "%#6.3g     ");
    printf("\n\nVreme izvrsavanja obrade: %g ms\n\n",(stop-start)*1000);
    printFooter();

    return 0;
}

/******************************************************************************/