/*
 * OpenMP version.
 */

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

float *compute_jacobi(int n, float init_value, float threshold, borders b, int *iterations) {
  float *x;
  float *new_x;
  float *tmp_x;
  float max_diff, new_value;
  int i, j;
  int nb = n + 2; // n plus the border

  /* Initialize boundary regions */
  x = create_sa_matrix(nb, nb);
  new_x = create_sa_matrix(nb, nb);
  for (i = 0; i < nb; i++) {
    x[IDX(nb, 0, i)] = b.north;
    x[IDX(nb, n + 1, i)] = b.south;
    x[IDX(nb, i, 0)] = b.west;
    x[IDX(nb, i, n + 1)] = b.east;
    new_x[IDX(nb, 0, i)] = b.north;
    new_x[IDX(nb, n + 1, i)] = b.south;
    new_x[IDX(nb, i, 0)] = b.west;
    new_x[IDX(nb, i, n + 1)] = b.east;
  }
  /* Initialize the rest of the matrix */
  for (i = 1; i <= n; i++) {
    for (j = 1; j <= n; j++) {
      x[IDX(nb, i, j)] = init_value;
    }
  }
  /* Iterative refinement of x until values converge */
  *iterations = 0;
  do {
    max_diff = 0;
    #pragma omp parallel for \
      reduction (max:max_diff) \
      default(none) private(new_value, j) shared(x, new_x, n, nb)
    for (i = 1; i <= n; i++) {
      for (j = 1; j <= n; j++) {
        new_value = 0.25 * (x[IDX(nb, i - 1, j)] + x[IDX(nb, i, j + 1)] + x[IDX(nb, i + 1, j)] + x[IDX(nb, i, j - 1)]);
        max_diff = (float) fmaxf(max_diff, fabs(new_value - x[IDX(nb, i, j)]));
        new_x[IDX(nb, i, j)] = new_value;
      }
    }
    tmp_x = new_x;
    new_x = x;
    x = tmp_x;
    (*iterations)++;
  } while (max_diff > threshold);
  return x;
}