#include "lu.h"

__kernel void rhs(__global double *g_u,
                  __global double *g_rsd,
                  __global double *g_frct,
                  __global double *g_qs,
                  __global double *g_rho_i,
                  int nx,
                  int ny,
                  int nz)
{
  int i, j, k, m;
  double tmp;
  __global double (*u)[ISIZ2/2*2+1][ISIZ1/2*2+1][5];
  __global double (*rsd)[ISIZ2/2*2+1][ISIZ1/2*2+1][5];
  __global double (*frct)[ISIZ2/2*2+1][ISIZ1/2*2+1][5];
  __global double (*qs)[ISIZ2/2*2+1][ISIZ1/2*2+1];
  __global double (*rho_i)[ISIZ2/2*2+1][ISIZ1/2*2+1];

  k = get_global_id(2);
  j = get_global_id(1);
  i = get_global_id(0);
  if (k >= nz || j >= ny || i >= nx) return;

  u     = (__global double (*)[ISIZ2/2*2+1][ISIZ1/2*2+1][5])g_u;
  rsd   = (__global double (*)[ISIZ2/2*2+1][ISIZ1/2*2+1][5])g_rsd;
  frct  = (__global double (*)[ISIZ2/2*2+1][ISIZ1/2*2+1][5])g_frct;
  qs    = (__global double (*)[ISIZ2/2*2+1][ISIZ1/2*2+1])g_qs;
  rho_i = (__global double (*)[ISIZ2/2*2+1][ISIZ1/2*2+1])g_rho_i;

  for (m = 0; m < 5; m++) {
    rsd[k][j][i][m] = - frct[k][j][i][m];
  }
  tmp = 1.0 / u[k][j][i][0];
  rho_i[k][j][i] = tmp;
  qs[k][j][i] = 0.50 * (  u[k][j][i][1] * u[k][j][i][1]
                        + u[k][j][i][2] * u[k][j][i][2]
                        + u[k][j][i][3] * u[k][j][i][3] )
                     * tmp;
}

