#include <stdlib.h>
#include <stdio.h>
#include <mpi.h>
#define AKYRO Lathos(__LINE__)
void Lathos(int line);

main(int argc, char *argv[])
{
    double     b_pivot;
    double     *computed_x;
    FILE       *fd;
    int        flops;
    int        *have_x;
    double     highest;
    int        i;
    int        j;
    int        k;
    double     multiplier;
    double     **my_a;
    double     *my_b;
    int        my_rank;
    double     num;
    int        num_processes;
    int        order;
    int        *perm_vector;
    double     *pivot_elements;
    int        pivot_row;
    double     *pivot_row_elements;
    FILE       *result_file;
    int        *rp_map;
    int        *sent_yet;
    MPI_Status status;
    int        temp;
    double     total;

    flops = 0;

    /* Ekkinish MPI, pairnei arithmo epeksergastwn kai pairnei kai rank */
    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &num_processes);
    MPI_Comm_rank(MPI_COMM_WORLD, &my_rank);

    /* Anoigma tou arxeiou eisodou. Einai ths morfh matrix.txt
	** Ston pinaka sthn arxh dhlwnoume to plhthos twn grammwn ths eksiswshs
	** kai sthn sunexeia dhlwnoume tis times tou pinaka se morfh an eixame x0,x1,x2,x3
	** prwta oles tis times tou x0, meta tou x1 kai meta tou x2 */
    fd = fopen("matrix.txt", "r");
    fscanf(fd, "%d", &order);

    /* To perm_vector mas dinei thn lista periexomenwn, entos tou my_a,
    ** opou mia dwsmenh grammh ston pinaka einai apothikeymenh
    ** Arxika h grammi 0 einai sto index 0, h grammi 1 sto index 1 , ktl... */
    perm_vector = (int *) malloc(sizeof(int) * order);
    if (perm_vector == NULL) AKYRO;

    for (i = 0; i < order; i++) /* Grammes */
        perm_vector[i] = i;

    rp_map = (int *) malloc(sizeof(int) * order);
    if (rp_map == NULL) AKYRO;

    for (i = 0; i < order; i++)
        rp_map[i] = i % num_processes;

    /*  Desmeysh xwrou tou pinaka sintelestwn kai toy dianysmatos lysewn pou ayto to process einai ypeythino gia apothikeysh. */
    my_a = (double **) malloc(sizeof(double *) * order);
    if (my_a == NULL) AKYRO;

    for (i = 0; i < order; i++)
    {
        if (rp_map[i] == my_rank)
        {
            my_a[i] = (double *) malloc(sizeof(double) * order);
            if (my_a[i] == NULL) AKYRO;
        }
        else my_a[i] = NULL;
    }

    my_b = (double *)  malloc(sizeof(double) * order);
    if (my_b == NULL) AKYRO;

    /* Anagnwsh tou pinaka sintelestwn. */
    for (i = 0; i < order; i++) /* Sthles */
    {
        for (j = 0; j < order; j++) /* Grammes */
        {
            fscanf(fd, "%lf", &num);

            if (rp_map[perm_vector[j]] == my_rank)
                my_a[perm_vector[j]][i] = num;
        }
    }

    /* Anagnwsh tou dianysmatos lysewn. */
    for (i = 0; i < order; i++) /* Grammes */
    {
        fscanf(fd, "%lf", &num);

        if (rp_map[perm_vector[i]] == my_rank)
            my_b[perm_vector[i]] = num;
    }

    /* Ektelesh Gaussian Elimination me partial pivoting. */

    /* Desmeysh xwrou gia ta stoixeia pivot. */
    pivot_elements = (double *) malloc(sizeof(double) * order);
    if (pivot_elements == NULL) AKYRO;

    pivot_row_elements = (double *) malloc(sizeof(double) * order);
    if (pivot_row_elements == NULL) AKYRO;

    for (i = 0; i < order - 1; i++) /* Grammes */
    {
        /* Apostolh twn stoixeiwn pivot stous allous processors. */
        for (j = i; j < order; j++) /* Grammes */
        {
            if (rp_map[perm_vector[j]] == my_rank)
            {
                for (k = 0; k < num_processes; k++)
                {
                    if (k != my_rank)
                    {
                        /* to i einai gia sthles edw. */
                        MPI_Send(&my_a[perm_vector[j]][i], 1, MPI_DOUBLE, k,
                                 j, MPI_COMM_WORLD);
                    }
                }
            }
        }

        /* Symplhrwsh tou pinaka me ta stoixeia pivot. */
        for (j = i; j < order; j++) /* Grammes */
        {
            if (rp_map[perm_vector[j]] == my_rank)
            {
                /* to i einai gia sthles edw. */
                pivot_elements[j] = my_a[perm_vector[j]][i];
                if (pivot_elements[j] < 0.0) pivot_elements[j] *= -1.0;
            }
            else
            {
                MPI_Recv(&pivot_elements[j], 1, MPI_DOUBLE,
                         rp_map[perm_vector[j]], j, MPI_COMM_WORLD,
                         &status);
                if (pivot_elements[j] < 0.0) pivot_elements[j] *= -1.0;
            }
        }

        /* Eyresh tis grammhs pivot. */
        highest = 0.0;

        for (j = i; j < order; j++) /* Grammes */
        {
            if (pivot_elements[j] > highest)
            {
                highest = pivot_elements[j];
                pivot_row = j;
            }
        }

        /* Allagh tis grammhs i me thn grammh pivot_row. */
        temp = perm_vector[i];
        perm_vector[i] = perm_vector[pivot_row];
        perm_vector[pivot_row] = temp;

        if (rp_map[perm_vector[i]] == my_rank)
        {
            for (j = 0; j < num_processes; j++)
            {
                if (j != my_rank)
                {
                    MPI_Send(my_a[perm_vector[i]], order, MPI_DOUBLE,
                             j, i, MPI_COMM_WORLD);
                    MPI_Send(&my_b[perm_vector[i]], 1, MPI_DOUBLE,
                             j, i, MPI_COMM_WORLD);
                }
            }

            /* Antigrafh ths grammhs pivot sto pivot_row_elements. */
            for (j = 0; j < order; j++) /* Sthles */
                pivot_row_elements[j] = my_a[perm_vector[i]][j];

            b_pivot = my_b[perm_vector[i]];
        }
        else
        {
            MPI_Recv(pivot_row_elements, order, MPI_DOUBLE,
                     rp_map[perm_vector[i]], i, MPI_COMM_WORLD,
                     &status);
            MPI_Recv(&b_pivot, 1, MPI_DOUBLE, rp_map[perm_vector[i]],
                     i, MPI_COMM_WORLD, &status);
        }

        /*Diadikasia afaireshs pollaplasiou me skopo h sthlh i na ginei 0*/
        for (j = i + 1; j < order; j++) /* Rows */
        {
            if (rp_map[perm_vector[j]] == my_rank)
            {
                multiplier = my_a[perm_vector[j]][i] /
                pivot_row_elements[i];

                for (k = 0; k < order; k++)
                {
                    my_a[perm_vector[j]][k] -= multiplier *
                    pivot_row_elements[k];
                    flops += 2;
                }

                my_b[perm_vector[j]] -= multiplier * b_pivot;
                flops += 2;
            }
        }
    }

    /* Twra exoume ena anw trigwniko systhma. Mporoume na to lysoume me pros ta pisw antikatastash */

    /* Desmeysh xwrou kai arxikopoihsh me false kathe stoixeio tou pinaka pou deixnei poies times x exoume. */
    have_x = (int *) malloc(sizeof(int) * order);
    if (have_x == NULL) AKYRO;

    for (i = 0; i < order; i++)
        have_x[i] = 0;

    /* Desmeysh xwrou gia th domh dedomenwn pou deixnei se poia processes steilame mia ypologismenh timh x. */
    sent_yet = (int *) malloc(sizeof(int) * num_processes);
    if (sent_yet == NULL) AKYRO;

    /* Desmeysh xwrou gia ton pinaka timwn tou x. */
    computed_x = (double *) malloc(sizeof(double) * order);
    if (computed_x == NULL) AKYRO;

    /* Pros ta pisw antikatastash. */
    for (i = order - 1; i >= 0; i--) /* Grammes */
    {
        if (rp_map[perm_vector[i]] == my_rank)
        {
            total = 0.0;

            for (j = i + 1; j < order; j++) /* times X */
            {
                if (!have_x[j])
                {
                    MPI_Recv(&computed_x[j], 1, MPI_DOUBLE,
                             rp_map[perm_vector[j]], j,
                             MPI_COMM_WORLD, &status);
                    have_x[j] = 1;
                }

                total += computed_x[j] * my_a[perm_vector[i]][j];
                flops += 2;
            }
            if(isnan((my_b[perm_vector[i]] - total) / my_a[perm_vector[i]][i]) != 1)
            {
                computed_x[i] = (my_b[perm_vector[i]] - total) / my_a[perm_vector[i]][i];
                flops++;
                have_x[i] = 1;
            }
            /* Apostolh ayths ths timhs sta alla processors. */
            for (j = 0; j < num_processes; j++)
                sent_yet[j] = 0;

            for (j = i - 1; j >= 0; j--) /* Grammes */
            {
                if (rp_map[perm_vector[j]] != my_rank)
                {
                    if (!sent_yet[rp_map[perm_vector[j]]])
                    {
                        MPI_Send(&computed_x[i], 1, MPI_DOUBLE,
                                 rp_map[perm_vector[j]], i,
                                 MPI_COMM_WORLD);
                        sent_yet[rp_map[perm_vector[j]]] = 1;
                    }
                }
            }
        }
    }

    /* Eggrafh twn ypologismenwn timwn tou x sto disko */
    if (rp_map[perm_vector[0]] == my_rank)
    {
        result_file = fopen("result.txt", "w");
        fprintf(result_file, "Dianysma lyshs:\n\n");

        for (i = 0; i < order; i++)
        {
            fprintf(result_file, "%f\n", computed_x[i]);
        }

        fclose(result_file);
    }

    /* Ektypwsh tou arithmou twn floating point leitoyrgiwn pou ektelese ayti i leitoyrgia. */
    printf("H diadikasia %d ektelese %d floating point leitourgies.\n",
           my_rank, flops);

    /* Termatismos MPI. */
    MPI_Finalize();
}

    /* Synarthsh se periptwsh lathous kai eksodou. */
void Lathos (int line)
{
   printf("Lathos sthn grammh %d\n", line);
   exit(1);
}