kernel_3D.cu

#ifndef _KERNEL_3D_H_
#define _KERNEL_3D_H_

#include <stdio.h>
#include "shmem_operations_3D.cu"
#include "prng.cu"
#include "reduce.cu"
#include "config.h"

__device__ float corona(data *gmem, data *sdata, data node, float n,
                        float gamma, float lambda, unsigned int x,
                        unsigned int y, unsigned int z, int size) {

  float corona = 0.0f;

  for (int i = 0; i < N; i++) {

    corona += (sdata[index3D(x + 1, y, z)].vector[i] +
               sdata[index3D(x - 1, y, z)].vector[i] +
               sdata[index3D(x, y + 1, z)].vector[i] +
               sdata[index3D(x, y - 1, z)].vector[i] +
               sdata[index3D(x, y, z + 1)].vector[i] +
               sdata[index3D(x, y, z - 1)].vector[i]);
  }

  return corona;
}

__device__ float newcorona(data *gmem, data *sdata, data node, float n,
                           float gamma, float lambda, unsigned int x,
                           unsigned int y, unsigned int z, float corona,
                           int size) {

  float newcorona = 0.0f;

  for (int i = 0; i < N; i++) {

    newcorona =
        -lambda * (2.0f * (sdata[index3D(x - 1, y - 1, z)].vector[i] +
                           sdata[index3D(x - 1, y + 1, z)].vector[i] +
                           sdata[index3D(x + 1, y + 1, z)].vector[i] +
                           sdata[index3D(x + 1, y - 1, z)].vector[i] +
                           sdata[index3D(x, y - 1, z - 1)].vector[i] +
                           sdata[index3D(x, y + 1, z - 1)].vector[i] +
                           sdata[index3D(x, y - 1, z + 1)].vector[i] +
                           sdata[index3D(x, y + 1, z + 1)].vector[i] +
                           sdata[index3D(x - 1, y, z - 1)].vector[i] +
                           sdata[index3D(x + 1, y, z - 1)].vector[i] +
                           sdata[index3D(x - 1, y, z + 1)].vector[i] +
                           sdata[index3D(x + 1, y, z + 1)].vector[i]) -
                   12.0f * corona + (sdata[index3D(x + 2, y, z)].vector[i] +
                                     sdata[index3D(x - 2, y, z)].vector[i] +
                                     sdata[index3D(x, y + 2, z)].vector[i] +
                                     sdata[index3D(x, y - 2, z)].vector[i] +
                                     sdata[index3D(x, y, z + 2)].vector[i] +
                                     sdata[index3D(x, y, z - 2)].vector[i]));
  }
  return newcorona;
}

__device__ float delta_E(data *gmem, data *sdata, data node, float n,
                         float gamma, float lambda, unsigned int x,
                         unsigned int y, unsigned int z, float corona,
                         float newcorona, float _equation_element1,
                         float _equation_element2, int size) {

  float c = 0.0f;
  float d = 0.0f;
  float sum = 0.0f;

  for (int i = 0; i < N; i++) {

    c +=
        (sdata[index3D(x, y, z)].vector[i] * sdata[index3D(x, y, z)].vector[i]);

    d += (node.vector[i] * node.vector[i]);

    sum = -node.vector[i] + sdata[index3D(x, y, z)].vector[i];
  }

  float e = (d - c) * (_equation_element1 + _equation_element2 * (d + c)) +
            (sum) * (corona + newcorona);

  return e;
}

__device__ void shmem_metropolis_3D(data *gmem, data *sdata,
                                    unsigned int location, unsigned int *seed,
                                    unsigned int size, float n, float gamma,
                                    float lambda, float eq_el1, float eq_el2,
                                    unsigned int x, unsigned int y,
                                    unsigned int z) {

  register data node;

  register float _equation_element1 = 3.0f + 0.5f * n + 21.0f * lambda;
  register float _equation_element2 = 0.0416666f * gamma;

  float __corona = corona(gmem, sdata, node, n, gamma, lambda, x, y, z, size);
  float __newcorona =
      newcorona(gmem, sdata, node, n, gamma, lambda, x, y, z, __corona, size);

  for (int i = 0; i < 4; i++) {

    for (int i = 0; i < N; i++) {

      node.vector[i] = sdata[index3D(x, y, z)].vector[i] +
                       0.4f * HybridTaus(seed[0], seed[1], seed[2], seed[3]) -
                       0.2f;
    }

    register float dE =
        delta_E(gmem, sdata, node, n, gamma, lambda, x, y, z, __corona,
                __newcorona, _equation_element1, _equation_element2, size);

    if (dE < 0) {

      sdata[index3D(x, y, z)] = node;
    }

    if (HybridTaus(seed[0], seed[1], seed[2], seed[3]) < __expf(-dE)) {

      sdata[index3D(x, y, z)] = node;
    }
  }
}

__device__ void Metropolis(data *gmem, data *sdata, unsigned int location,
                           unsigned int *register_seed, unsigned int size,
                           float n, float gamma, float lambda, float eq_el1,
                           float eq_el2, unsigned int x, unsigned int y,
                           unsigned int z, unsigned int y1, unsigned int y2) {

  for (int i = 0; i < 50; ++i) {

#define metropolis_3D(x, y, z)                                                 \
  shmem_metropolis_3D(gmem, sdata, location, register_seed, size, n, gamma,    \
                      lambda, eq_el1, eq_el2, x, y, z);                        \
  __syncthreads();

    metropolis_3D(2 * threadIdx.x, y1, z);
    metropolis_3D(2 * threadIdx.x + 1, y1, z);
    metropolis_3D(2 * threadIdx.x, y2, z);
    metropolis_3D(2 * threadIdx.x + 1, y2, z);
    metropolis_3D(2 * threadIdx.x, y1 + 1, z);
    metropolis_3D(2 * threadIdx.x + 1, y1 + 1, z);
    metropolis_3D(2 * threadIdx.x, y2 + 1, z);
    metropolis_3D(2 * threadIdx.x + 1, y2 + 1, z);
    metropolis_3D(2 * threadIdx.x, y1, z + 1);
    metropolis_3D(2 * threadIdx.x + 1, y1, z + 1);
    metropolis_3D(2 * threadIdx.x, y2, z + 1);
    metropolis_3D(2 * threadIdx.x + 1, y2, z + 1);
    metropolis_3D(2 * threadIdx.x, y1 + 1, z + 1);
    metropolis_3D(2 * threadIdx.x + 1, y1 + 1, z + 1);
    metropolis_3D(2 * threadIdx.x, y2 + 1, z + 1);
    metropolis_3D(2 * threadIdx.x + 1, y2 + 1, z + 1);
  }
}

__global__ void Phi_3D(data *g_idata, unsigned int *seed, unsigned int size,
                       float n, float gamma, float lambda, float eq_el1,
                       float eq_el2, unsigned int offsetx, unsigned int offsety,
                       unsigned int offsetz) {

  int row = 2 * (blockIdx.x) * SHMEM_CUBE_SIZE + 2 * threadIdx.x +
            SHMEM_CUBE_SIZE * offsetx;

  int col = 2 * (blockIdx.y) * SHMEM_CUBE_SIZE + 4 * threadIdx.y +
            SHMEM_CUBE_SIZE * offsety;

  int dim = 2 * (blockIdx.z) * SHMEM_CUBE_SIZE + 2 * threadIdx.z +
            SHMEM_CUBE_SIZE * offsetz;

  __shared__ data
  sdata[(4 + SHMEM_CUBE_SIZE) * (4 + SHMEM_CUBE_SIZE) * (4 + SHMEM_CUBE_SIZE)];

  unsigned int location = row * size + col + dim * size * size;
  unsigned int seed_location =
      ((blockIdx.x) * blockDim.x + threadIdx.x) * blockDim.y +
      ((blockIdx.y) * blockDim.y + threadIdx.y) +
      (((blockIdx.z) * blockDim.z + threadIdx.z) * blockDim.x * blockDim.y);

  unsigned int register_seed[4];

  register_seed[0] = seed[4 * seed_location];
  register_seed[1] = seed[4 * seed_location + 1];
  register_seed[2] = seed[4 * seed_location + 2];
  register_seed[3] = seed[4 * seed_location + 3];

  int x = (2 * threadIdx.x);
  int y = (4 * threadIdx.y);
  int z = (2 * threadIdx.z);

  globmem_to_shmem_3D(sdata, g_idata, x, y, z, row, col, dim, location, size);
  __syncthreads();

  int y1 = ((threadIdx.z & 1) ? (((x + 2) & 3 ? y : y + 2))
                              : ((x + 2) & 3 ? y + 2 : y));
  int y2 = ((threadIdx.z & 1) ? (((x + 2) & 3 ? y + 2 : y))
                              : (((x + 2) & 3 ? y : y + 2)));

  Metropolis(g_idata, sdata, location, register_seed, size, n, gamma, lambda,
             eq_el1, eq_el2, x, y, z, y1, y2);

  int first_line_index =
      (x + 2) * (SHMEM_CUBE_SIZE + 4) + (y + 2) + (z + 2) * SHM_STRIDE_Z;
  int second_line_index =
      (x + 3) * (SHMEM_CUBE_SIZE + 4) + (y + 2) + (z + 2) * SHM_STRIDE_Z;

  g_idata[location] = sdata[first_line_index];
  g_idata[location + 1] = sdata[first_line_index + 1];
  g_idata[location + 2] = sdata[first_line_index + 2];
  g_idata[location + 3] = sdata[first_line_index + 3];
  g_idata[(row + 1) * size + col + dim * size * size] =
      sdata[second_line_index];
  g_idata[(row + 1) * size + col + dim * size * size + 1] =
      sdata[second_line_index + 1];
  g_idata[(row + 1) * size + col + dim * size * size + 2] =
      sdata[second_line_index + 2];
  g_idata[(row + 1) * size + col + dim * size * size + 3] =
      sdata[second_line_index + 3];

  g_idata[location + size * size] = sdata[first_line_index + SHM_STRIDE_Z];
  g_idata[location + size * size + 1] =
      sdata[first_line_index + 1 + SHM_STRIDE_Z];
  g_idata[location + size * size + 2] =
      sdata[first_line_index + 2 + SHM_STRIDE_Z];
  g_idata[location + size * size + 3] =
      sdata[first_line_index + 3 + SHM_STRIDE_Z];
  g_idata[(row + 1) * size + col + dim * size * size + size * size] =
      sdata[second_line_index + SHM_STRIDE_Z];
  g_idata[(row + 1) * size + col + dim * size * size + size * size + 1] =
      sdata[second_line_index + 1 + SHM_STRIDE_Z];
  g_idata[(row + 1) * size + col + dim * size * size + size * size + 2] =
      sdata[second_line_index + 2 + SHM_STRIDE_Z];
  g_idata[(row + 1) * size + col + dim * size * size + size * size + 3] =
      sdata[second_line_index + 3 + SHM_STRIDE_Z];

  seed[4 * seed_location] = register_seed[0];
  seed[4 * seed_location + 1] = register_seed[1];
  seed[4 * seed_location + 2] = register_seed[2];
  seed[4 * seed_location + 3] = register_seed[3];
}

template <class T>
__global__ void correlation(T *g_idata, float *g_odata, unsigned int n,
                            int size, int R) {
  float *sdata = SharedMemory<float>();

  unsigned int tid = threadIdx.x;
  unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
  unsigned int x = (i / size);
  unsigned int y = i - (x * size);

  unsigned int x_stride = ((x - R) & (size - 1) * size + y);
  unsigned int y_stride = (x * size + (y + R) & (size - 1));

  sdata[tid] = (i < n) ? (g_idata[i].vector[0] * g_idata[x_stride].vector[0] +
                          (g_idata[i].vector[0] * g_idata[y_stride].vector[0]))
                       : 0;

  __syncthreads();

  for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
    if (tid < s) {
      sdata[tid] += sdata[tid + s];
    }
    __syncthreads();
  }

  if (tid == 0)
    g_odata[blockIdx.x] = sdata[0];
}
#endif