#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include "mpi.h"

#define V 4

//Zakładamy, że liczba procesów p jest równa r^2, gdzie r jest liczbą całkowitą, p <= V^2 oraz V % r = 0

float min(float a, float b)
{
    if(a <= b)
        return a;
    return b;
}

int oblicz_pozycje(int i, int j, int w, int r)
{
    return w * w * ((i / w) * r + (j / w)) + (i % w) * w + (j % w);
}

void oblicz_pozycje_odwrotna(int a, int w, int r, int *i, int *j)
{
    *i = (a / (w * w)) / r * w + (a % (w * w)) / w;
    *j = (a / (w * w)) % r * w + (a % w);
}

int main(int argc, char *argv[])
{
    int i, j, k;
    float macierz[V][V];
    float wynik[V][V];
    float **macierz_lokalna;
    float *macierz_tablica, *macierz_lokalna_tablica;
    float *komunikat_wiersz, *komunikat_kolumna;
    const float nieskonczonosc = 1.0/0.0;
    int rank;
    MPI_Comm comm_cart, comm_wiersz, comm_kolumna;
    int dims[2];
    int periods[] = {1, 1};
    int nr_wiersza, nr_kolumny;
    int coords[2];
    int p, r, wielkosc_lokalna;

    MPI_Init(&argc, &argv);
    MPI_Comm_size(MPI_COMM_WORLD, &p);

    r = (int)floor(sqrt(p));
    dims[0] = dims[1] = r;
    wielkosc_lokalna = V / r;
    macierz_lokalna = (float **)calloc(wielkosc_lokalna, sizeof(float *));
    macierz_tablica = (float *)calloc(V * V, sizeof(float));
    macierz_lokalna_tablica = (float *)calloc(wielkosc_lokalna * wielkosc_lokalna, sizeof(float));
    komunikat_wiersz = (float *)calloc(wielkosc_lokalna, sizeof(float));
    komunikat_kolumna = (float *)calloc(wielkosc_lokalna, sizeof(float));

    for(i = 0; i < wielkosc_lokalna; ++i)
    {
        macierz_lokalna[i] = (float *)calloc(wielkosc_lokalna, sizeof(float));
    }

    MPI_Cart_create(MPI_COMM_WORLD, 2, dims, periods, 1, &comm_cart);
    MPI_Comm_rank(comm_cart, &rank);
    MPI_Cart_coords(comm_cart, rank, 2, coords);
    nr_wiersza = coords[0];
    nr_kolumny = coords[1];
    dims[0] = 0;
    dims[1] = 1;
    MPI_Cart_sub(comm_cart, dims, &comm_wiersz);
    dims[0] = 1;
    dims[1] = 0;
    MPI_Cart_sub(comm_cart, dims, &comm_kolumna);
    if(rank == 0)
    {

/*
 * Ustawienie wartości w macierzy
 * dla macierz[i][i] = 0
 * jeśli istnieje krawędź z i do j to macierz[i][j] = waga krawędzi
 * jeśli krawędź z i do j nie istnieje to macierz[i][j] = nieskończoność
 */
        for(i = 0; i < V; ++i)
        {
            for(j = 0; j < V; ++j)
            {
                if(i == j)
                    macierz[i][i] = 0;
                else
                    macierz[i][j] = nieskonczonosc;
            }
        }

macierz[0][2] = -2;
macierz[1][0] = 4;
macierz[1][2] = 3;
macierz[2][3] = 2;
macierz[3][1] = -1;

/*
macierz[0][1] = 1;
macierz[0][3] = 5;
macierz[1][2] = 2;
macierz[1][3] = 1;
macierz[2][0] = 1;
macierz[2][3] = 3;
macierz[3][1] = 3;
macierz[3][2] = 2;
*/

//Rozesłanie lokalnych macierzy do odpowiednich procesów
        for(i = 0; i < V; ++i)
        {
            for(j = 0; j < V; ++j)
            {
                macierz_tablica[oblicz_pozycje(i, j, wielkosc_lokalna, r)] = macierz[i][j];
            }
        }
    }
    MPI_Scatter(macierz_tablica, wielkosc_lokalna * wielkosc_lokalna, MPI_FLOAT, macierz_lokalna_tablica, wielkosc_lokalna * wielkosc_lokalna, MPI_FLOAT, 0, comm_cart);
    for(i = 0; i < wielkosc_lokalna; ++i)
    {
        for(j = 0; j < wielkosc_lokalna; ++j)
        {
            macierz_lokalna[i][j] = macierz_lokalna_tablica[i * wielkosc_lokalna + j];
        }
    }

//Główna część równoległego algorytmu Floyda
    for(k = 0; k < V; ++k)
    {

//Rozesłanie fragmentów kolumn i wierszy
        if(nr_wiersza == k / wielkosc_lokalna)
        {
            for(i = 0; i < wielkosc_lokalna; ++i)
            {
                komunikat_kolumna[i] = macierz_lokalna[k % wielkosc_lokalna][i];
            }
        }
        if(nr_kolumny == k / wielkosc_lokalna)
        {
            for(i = 0; i < wielkosc_lokalna; ++i)
            {
                komunikat_wiersz[i] = macierz_lokalna[i][k % wielkosc_lokalna];
            }
        }
        MPI_Bcast(komunikat_wiersz, wielkosc_lokalna, MPI_FLOAT, k / wielkosc_lokalna, comm_wiersz);
        MPI_Bcast(komunikat_kolumna, wielkosc_lokalna, MPI_FLOAT, k / wielkosc_lokalna, comm_kolumna);

//Obliczenia w algorytmie
        for(i = 0; i < wielkosc_lokalna; ++i)
        {
            for(j = 0; j < wielkosc_lokalna; ++j)
            {
                macierz_lokalna[i][j] = min(macierz_lokalna[i][j], komunikat_wiersz[i] + komunikat_kolumna[j]);
            }
        }
    }

//Łączenie wyników
    for(i = 0; i < wielkosc_lokalna; ++i)
    {
        for(j = 0; j < wielkosc_lokalna; ++j)
        {
            macierz_lokalna_tablica[i * wielkosc_lokalna + j] = macierz_lokalna[i][j];
        }
    }
    MPI_Gather(macierz_lokalna_tablica, wielkosc_lokalna * wielkosc_lokalna, MPI_FLOAT, macierz_tablica, wielkosc_lokalna * wielkosc_lokalna, MPI_FLOAT, 0, comm_cart);
    if(rank == 0)
    {
        for(k = 0; k < V * V; ++k)
        {
            oblicz_pozycje_odwrotna(k, wielkosc_lokalna, r, &i, &j);
            wynik[i][j] = macierz_tablica[k];
        }

//Wyswietlanie wyniku
        for(i = 0; i < V; ++i)
        {
            for(j = 0; j < V; ++j)
                printf("%.3f\t", wynik[i][j]);
            printf("\n");
        }
    }

    MPI_Finalize();
    return 0;
}