#include "mpi.h"      // knjižnica MPI
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>



#define EPSILON  0.001   /* natancnost */
#define N 20 //velikost kvadrata


int main(int argc, char *argv[])
{
    clock_t startclock, endclock;
    double elapsed;
    int my_rank;    // rank (oznaka) procesa
    int num_of_processes;  // število procesov
    int source;      // rank pošiljatelja
    int destination;    // rank sprejemnika
    int tag = 50;    // zaznamek sporoèila
    char   ch;
    int i, x, y;
    double diff;
    int skupnifinish;
    double temp;
    int dodatni;
    int finish=1;
    int *sizeproc;
    int *velikost;
    double mean;
    double **T;
    double **Tnew;
    int size;
    MPI_Status status;    // status sprejema
    int stevilo_trakov=0;
    int wherestart=0;
    int whereend=0;
    MPI_Datatype  DelnaMatrika;

    int n;
    double h;
    double x1, y1, rezultat, v, koncnaVsota=0;
    double start, stop;

    num_of_processes=1;
    MPI_Init(&argc, &argv);  // inicializacija MPI okolja
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);  // poizvedba po ranku procesa
    MPI_Comm_size(MPI_COMM_WORLD, &num_of_processes);  // poizvedba po številu procesov
    startclock=clock();
    finish=1;
    stevilo_trakov=0;
    wherestart=0;
    whereend=0;
    sizeproc = (int *)malloc(num_of_processes * sizeof(int));
    velikost = (int *)malloc(num_of_processes * sizeof(int));
    T = (double **)malloc(sizeof(double *) * N);
    Tnew = (double **)malloc(sizeof(double *) * N);
    for(i=0; i<N; i++)
    {
        T[i] = (double *)malloc(sizeof(double) * N);
        Tnew[i] = (double *)malloc(sizeof(double) * N);
    }

    /* inicializacija robnih toèk */
    for(i=0; i < N; i++)
    {
        Tnew[N-1][i] = 0.0;
        Tnew[i][0] = Tnew[i][N-1] = Tnew[0][i] = 100.0;
        T[N-1][i] = 0.0;
        T[i][0] = T[i][N-1] = T[0][i] = 100.0;
    }
    /* inicializacija toèk v sredini */
    mean = 0.0;
    for(i=0; i < N; i++)
    {
        mean += T[i][0] + T[i][N-1] + T[0][i] + T[N-1][i];
    }
    mean /= (4.0 * N);
    for(y=1; y < N-1; y++)
    {
        for(x=1; x < N-1; x++)
        {
            T[y][x] = mean;
            Tnew[y][x] = mean;
        }
    }



    // ustvarimo nov tip mpi datatype, array velikosti N*stevilo_trakov
    MPI_Type_contiguous(N, MPI_DOUBLE, &DelnaMatrika);
    MPI_Type_commit(&DelnaMatrika);
    temp=N/num_of_processes;
    stevilo_trakov=floor(temp);
    size=stevilo_trakov;
    dodatni=N%num_of_processes;
    wherestart = my_rank*size;
    if(dodatni>0)
    {
        if(my_rank<dodatni)size++;
        if(dodatni>my_rank) wherestart=wherestart+my_rank;
        else wherestart=wherestart+dodatni;
    }
    whereend = wherestart+size;
    if(my_rank==0)wherestart++;
    if(my_rank+1==num_of_processes)whereend--;
    fflush(stdout);

    MPI_Gather( &wherestart, 1, MPI_INT, sizeproc, 1, MPI_INT, 0, MPI_COMM_WORLD);
    for(i=0;i<num_of_processes;i++)
    {
        if(i+1!=num_of_processes)
        {
            velikost[i]=sizeproc[i+1]-sizeproc[i];
        }
        else
        {
            velikost[i]=N-sizeproc[i]-1;
        }
        while(1)
        {
            diff = 0;
            if(finish>0)
            {
                for(x=wherestart; x <whereend; x++)
                {
                    for(y=1; y < N-1; y++)
                    {
                        Tnew[x][y] = (T[x][y-1]+T[x][y+1]+T[x-1][y]+T[x+1][y])/4.0;
                        if( fabs(Tnew[x][y]-T[x][y]) > diff )
                        {
                            diff = fabs(Tnew[x][y]-T[x][y]);
                        }
                    }
                }
            }
            if(my_rank>0)
            {
                MPI_Send(Tnew[wherestart], 1, DelnaMatrika, my_rank-1,0, MPI_COMM_WORLD);
            }
            if(my_rank+1<num_of_processes)
            {
                MPI_Send(Tnew[whereend-1], 1, DelnaMatrika, my_rank+1,0, MPI_COMM_WORLD);
            }
            if(my_rank>0)
            {
                MPI_Recv(Tnew[wherestart-1], 1, DelnaMatrika, my_rank-1, 0, MPI_COMM_WORLD, &status);
            }
            if(my_rank+1<num_of_processes)
            {
                MPI_Recv(Tnew[whereend], 1, DelnaMatrika, my_rank+1, 0, MPI_COMM_WORLD, &status);
            }
            if(diff<EPSILON)finish=0;
            MPI_Allreduce (&finish, &skupnifinish, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD );
            for(x=wherestart-1; x < whereend+1; x++)
            {
                for(y=1; y < N-1; y++)
                {
                    T[x][y] = Tnew[x][y];
                }
            }
            if( skupnifinish == 0) break;
        }


        if(my_rank>0)
        {
            for(i=wherestart; i < whereend; i++)
            {
                MPI_Send(T[i], 1, DelnaMatrika, 0,0, MPI_COMM_WORLD);
            }
        }
        else
        {
            for(x=1; x < num_of_processes; x++)
            {
                for(i=sizeproc[x]; i <sizeproc[x]+velikost[x]; i++)
                {
                    MPI_Recv(T[i], 1, DelnaMatrika, x, 0, MPI_COMM_WORLD, &status);
                }
            }
        }

        if(my_rank==0)
        {
            for(y=0; y < N; y++)
            {
                for(x=0; x < N; x++)
                {
                    printf(" %3.0f", T[y][x]);
                }
                printf("\n");
            }
        }
        printf("\n");
        /* sproscanje pomnilnika */
        for(i=0; i<N; i++)
        {
            free(T[i]);
            free(Tnew[i]);
        }
        free(T);
        free(Tnew);
        if(my_rank==0)
        {
            endclock=clock();
            elapsed = ((double) (endclock - startclock) * 1000) / CLOCKS_PER_SEC;
            printf("Za kvadrat velikosti : %d je bilo potrebno %d milisekund\n", N, (int)elapsed);
        }
    }
    MPI_Finalize();
    return 0;
}