/*
 * MPI version with the matrix subdivided by "lines".
 */

#include <stdio.h>
#include <math.h>
#include <mpi.h>
#include "../config/config.h"

#define TAG_BORDER 0

typedef struct borders {
  double north;
  double east;
  double south;
  double west;
} borders;

void compute_jacobi(int n, double init_value, double threshold, borders b);
void print_matrix(int rows, int cols, double x[rows][cols]);

int main(int argc, char* argv[]) {
  int rank, numprocs;
  int n;
  double init_value, threshold;
  double north, south, east, west;
  borders b;
  int config_loaded;
  configuration config;

  MPI_Init(&argc, &argv);

  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
  MPI_Comm_size(MPI_COMM_WORLD, &numprocs);

  if (numprocs != 2) {
    MPI_Abort(MPI_COMM_WORLD, 1);
  }

  if (rank == 0) {
    config_loaded = load_config(&config);
    if (config_loaded != 0) {
      MPI_Abort(MPI_COMM_WORLD, 1);
    }
    n = config.n;
    threshold = config.threshold;
    init_value = config.init_value;
    north = config.north;
    south = config.south;
    east = config.east;
    west = config.west;
  }
  MPI_Bcast(&n, 1, MPI_INT, 0, MPI_COMM_WORLD);
  MPI_Bcast(&init_value, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&threshold, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&north, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&south, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&east, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);
  MPI_Bcast(&west, 1, MPI_DOUBLE, 0, MPI_COMM_WORLD);

  b.north = north;
  b.south = south;
  b.east = east;
  b.west = west;

  if (n != 4) {
    MPI_Abort(MPI_COMM_WORLD, 1);
  }
 
  int rows = 2;
  double x[rows + 2][n + 2];
  double max_diff, global_max_diff, new_x;
  int i, j;
  MPI_Status status;

  /* Initialize the matrix */
  for (i = 0; i < rows + 2; i++) {
    for (j = 1; j <= n; j++) {
      x[i][j] = init_value;
    }
  }
  /* Initialize boundary regions */
  for (i = 0; i < rows + 2; i++) {
    x[i][0] = b.west;
    x[i][n + 1] = b.east;
  }
  if (rank == 0) {
    for (i = 1; i <= n + 1; i++) {
      x[0][i] = b.north;
    }
  } else if (rank == 1){
    for (i = 1; i < n + 1; i++) {
      x[rows + 1][i] = b.south;
    }
  }
  /* Iterative refinement of x until values converge */
  do {
    max_diff = 0;
    global_max_diff = 0;
    for (i = 1; i <= rows; i++) {
      for (j = 1; j <= n; j++) {
        new_x = 0.25 * (x[i - 1][j] + x[i][j + 1] + x[i + 1][j] + x[i][j - 1]);
        max_diff = (double) fmax(max_diff, fabs(new_x - x[i][j]));
        x[i][j] = new_x;
      }
    }
    if (rank == 0) {
      MPI_Send(&x[rows][0], n + 2, MPI_DOUBLE, 1, TAG_BORDER, MPI_COMM_WORLD);
      MPI_Recv(&x[rows + 1][0], n + 2, MPI_DOUBLE, 1, TAG_BORDER, MPI_COMM_WORLD, &status);
    } else {
      MPI_Recv(&x[0][0], n + 2, MPI_DOUBLE, 0, TAG_BORDER, MPI_COMM_WORLD, &status);
      MPI_Send(&x[1][0], n + 2, MPI_DOUBLE, 0, TAG_BORDER, MPI_COMM_WORLD);
    }
    MPI_Allreduce(&max_diff, &global_max_diff, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
  } while (global_max_diff > threshold);

  MPI_Finalize();

  return 0;
}

void print_matrix(int rows, int cols, double x[rows][cols]) {
  int i, j;
  printf("---------------------------------\n");
  for (i = 0; i < rows; i++) {
    for (j = 0; j < cols; j++) {
      printf("%f\t", x[i][j]);
    }
    printf("\n");
  }
  printf("---------------------------------\n");
  fflush(stdout);
}