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

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

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

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

int main(int argc, char* argv[]) {
  int config_loaded;
  configuration config;
  borders b;

  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;

  compute_jacobi(config.n, config.init_value, config.threshold, b);

  return 0;
}

void compute_jacobi(int n, double init_value, double threshold, borders b) {
  double x[n + 2][n + 2];
  double max_diff, new_x;
  int i, j;

  /* Initialize boundary regions */
  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);
}