/**
@file "heat_openmp_1.c"
@brief Paralelno resenje Heat Equation pomocu OpenMP. 
       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.

@author Ognjen Scekic

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

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

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

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

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

/**
@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[5][5] = {0, 0, 0, 0, 0,
                     0, 0, 2, 0, 0,
                     0, 2, 3, 2, 0,
                     0, 0, 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 Dimenzija matrice 
*/
unsigned int matrixDimension;

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

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


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

/**
@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++) {
	    if (stream==NULL)
		scanf("%lf", &temp);
	    else
		fscanf(stream,"%lf", &temp);

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

/**
@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
*/
void ifWrongDimensionAbort(int pMatrixDimension) {
    if (pMatrixDimension < MIN_MATRIX_DIM || pMatrixDimension > MAX_MATRIX_DIM ||
	pMatrixDimension%2 != 1){
	printf("Pogresne dimenzije matrice! Prekidam rad programa...\n");
	exit(ErrorDimension);
    }
}

/**
@brief Provera ispravnosti broja niti.
@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.
*/
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;
    double start, stop;
    FILE *input_file;

    /* Pokazivaci na matrice, pomocu kojih cemo pokazivati iz koje matrice treba citati
       a u koju upisivati u tekucoj iteraciji. */
    double (*readMatrix)[MAX_MATRIX_DIM][MAX_MATRIX_DIM];
    double (*writeMatrix)[MAX_MATRIX_DIM][MAX_MATRIX_DIM];

    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);
    
    start = omp_get_wtime();

    for (step=0; step<numOfSteps; step++){
	if (step%2){ /* ako je neparna iteracija */
	    readMatrix=&matrix1;
	    writeMatrix=&matrix0;
	}else{
	    readMatrix=&matrix0;
	    writeMatrix=&matrix1;
	}

	#pragma omp parallel for shared(readMatrix, writeMatrix) schedule(static)
	for (i=1; i<matrixDimension-1; 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 parallel for. Implicitna barijera i flush operacija. */
	
	
	/* Rezultat je u matrici na koju pokazuje writeMatrix */
	//printIteration(step+1);
        //printfMatrix(writeMatrix, matrixDimension, matrixDimension, "%#6.3g     ");

    }
    stop = omp_get_wtime();
    printIteration(numOfSteps);
    printfMatrix(writeMatrix, matrixDimension, matrixDimension, "%#6.3g     ");
    printf("\n\nVreme izvrsavanja obrade: %g ms\n\n",(stop-start)*1000);
    printFooter();
    return 0;
}

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