#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_psinv(__global double *r,
                           __global double *u,
                           __global double *c,
                           int n1,
                           int n2,
                           int n3,
                           int offset)
{
  int i3 = get_global_id(1) + 1;
  int i2 = get_group_id(0) + 1;
  int lid = get_local_id(0);

  __local double r1[M], r2[M];

  int i1;
  for (i1 = lid; i1 < n1; i1 += get_local_size(0)) {
    r1[i1] = r[i3*n2*n1+(i2-1)*n2+i1+offset]
           + r[i3*n2*n1+(i2+1)*n1+i1+offset]
           + r[(i3-1)*n2*n1+i2*n1+i1+offset]
           + r[(i3+1)*n2*n1+i2*n1+i1+offset];
    r2[i1] = r[(i3-1)*n2*n1+(i2-1)*n1+i1+offset]
           + r[(i3-1)*n2*n1+(i2+1)*n1+i1+offset]
           + r[(i3+1)*n2*n1+(i2-1)*n1+i1+offset]
           + r[(i3+1)*n2*n1+(i2+1)*n1+i1+offset];
  }

  barrier(CLK_LOCAL_MEM_FENCE);

  for (i1 = lid+1; i1 < n1-1; i1 += get_local_size(0)) {
    u[i3*n2*n1+i2*n1+i1+offset] = u[i3*n2*n1+i2*n1+i1+offset]
      + c[0] * r[i3*n2*n1+i2*n1+i1+offset]
      + c[1] * ( r[i3*n2*n1+i2*n1+i1-1+offset]
               + r[i3*n2*n1+i2*n1+i1+1+offset]
               + r1[i1] )
      + c[2] * ( r2[i1] + r1[i1-1] + r1[i1+1] );
  }
}

