#include "mg.h"
#include "mg_dim.h"

#define PSINV_DIM       PSINV_DIM_GPU
#define RESID_DIM       RESID_DIM_GPU
#define RPRJ3_DIM       RPRJ3_DIM_GPU
#define INTERP_1_DIM    INTERP_1_DIM_GPU
#define NORM2U3_DIM     NORM2U3_DIM_GPU
#define COMM3_1_DIM     COMM3_1_DIM_GPU
#define COMM3_2_DIM     COMM3_2_DIM_GPU
#define COMM3_3_DIM     COMM3_3_DIM_GPU
#define ZERO3_DIM       ZERO3_DIM_GPU

__kernel
void kernel_norm2u3(__global double *r,
                    const int n1, const int n2, const int n3,
                    __global double *res_sum,
                    __global double *res_max,
                    __local double *scratch_sum,
                    __local double *scratch_max)
{
  int i3 = get_global_id(1) + 1;
  int i2 = get_group_id(0) + 1;
  int i1 = get_local_id(0) + 1;

  double s = 0.0;
  double my_rnmu = 0.0;
  double a;
  while (i1 < n1-1) {
    double r321 = r[i3*n2*n1+i2*n1+i1];
    //s = s + pow(r321, 2.0);
    s = s + r321 * r321;
    a = fabs(r321);
    my_rnmu = (a > my_rnmu) ? a : my_rnmu;

    i1 += get_local_size(0);
  }
  int lid = get_local_id(0);
  scratch_sum[lid] = s;
  scratch_max[lid] = my_rnmu;

  barrier(CLK_LOCAL_MEM_FENCE);

  // reduction in a work-group
  for (int i = get_local_size(0) / 2; i > 0; i >>= 1) {
    if (lid < i) {
      scratch_sum[lid] += scratch_sum[lid + i];
      scratch_max[lid] = (scratch_max[lid] > scratch_max[lid + i])
                       ? scratch_max[lid] : scratch_max[lid + i];
    }
    barrier(CLK_LOCAL_MEM_FENCE);
  }
  if (lid == 0) {
    int idx = get_group_id(1) * get_num_groups(0) + get_group_id(0);
    res_sum[idx] = scratch_sum[0];
    res_max[idx] = scratch_max[0];
  }
}

