/*
 * Sequential version.
 * No paralleliization used.
 */

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

double **compute_jacobi(int n, double init_value, double threshold, borders b);

int main(int argc, char* argv[]) {
  int numprocs;
  int config_loaded;
  configuration config;
  borders b;
  double startwtime = 0.0, endwtime;
  double **x;

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

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

  config_loaded = load_config(&config);
  if (config_loaded != 0) {
    return 1;
  }

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

  startwtime = MPI_Wtime();
  x = compute_jacobi(config.n, config.init_value, config.threshold, b);

  endwtime = MPI_Wtime();
  printf("Wall clock time: %fs\n", endwtime - startwtime);

  destroy_matrix(x, config.n + 2);
  MPI_Finalize();

  return 0;
}

double **compute_jacobi(int n, double init_value, double threshold, borders b) {
  double **x;
  double max_diff, new_x;
  int i, j;

  /* Initialize boundary regions */
  x = create_matrix(n + 2, n + 2);
  for (i = 1; i <= n; i++) {
    x[0][i] = b.north;
    x[n + 1][i] = b.south;
    x[i][0] = b.west;
    x[i][n + 1] = b.east;
  }
  /* Initialize the rest of the matrix */
  for (i = 1; i <= n; i++) {
    for (j = 1; j <= n; j++) {
      x[i][j] = init_value;
    }
  }
  /* Iterative refinement of x until values converge */
  do {
    max_diff = 0;
    for (i = 1; i <= n; 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;
      }
    }
  } while (max_diff > threshold);
  return x;
}