#include "bt.h"

//---------------------------------------------------------------------
// this function returns the exact solution at point xi, eta, zeta
//---------------------------------------------------------------------
void exact_solution(double xi, double eta, double zeta, double dtemp[5],
                    __constant double *g_ce);

__kernel void error_norm(__global double *g_u,
                         __constant double *g_ce,
                         __global double *g_rms,
                         __local double *l_rms,
                         int gp0,
                         int gp1,
                         int gp2)
{
  int i, j, k, m, lid;
  double xi, eta, zeta, u_exact[5], add;
  __local double *rms_local;

  k = get_global_id(0) + 1;
  lid = get_local_id(0);
  rms_local = &l_rms[lid * 5];

  for (m = 0; m < 5; m++) {
    rms_local[m] = 0.0;
  }

  if (k < gp2) {
    __global double (*u)[JMAXP+1][IMAXP+1][5] =
      (__global double (*)[JMAXP+1][IMAXP+1][5])g_u;

    zeta = (double)k * dnzm1;
    for (j = 0; j < gp1; j++) {
      eta = (double)j * dnym1;
      for (i = 0; i < gp0; i++) {
        xi = (double)i * dnxm1;
        exact_solution(xi, eta, zeta, u_exact, g_ce);

        for (m = 0; m < 5; m++) {
          add = u[k][j][i][m]-u_exact[m];
          rms_local[m] = rms_local[m] + add*add;
        }
      }
    }
  }

  barrier(CLK_LOCAL_MEM_FENCE);

  if (lid == 0) {
    for (i = 1; i < get_local_size(0); i++) {
      __local double *rms_other = &l_rms[i * 5];
      for (m = 0; m < 5; m++) {
        rms_local[m] += rms_other[m];
      }
    }

    __global double *rms = &g_rms[get_group_id(0) * 5];
    for (m = 0; m < 5; m++) {
      rms[m] = rms_local[m];
    }
  }
}
