#include <hip/hip_runtime.h>
#define BX 8
#define BY 8
#define BZ 8
#define HALO 1
__global__ void stencil_3d(const float* in, float* out, int nx, int ny, int nz) {
__shared__ float tile[BZ + 2 * HALO][BY + 2 * HALO][BX + 2 * HALO];
int tx = threadIdx.x;
int ty = threadIdx.y;
int tz = threadIdx.z;
int gx = blockIdx.x * blockDim.x + tx;
int gy = blockIdx.y * blockDim.y + ty;
int gz = blockIdx.z * blockDim.z + tz;
int idx = gx + gy * nx + gz * nx * ny;
if (gx < nx && gy < ny && gz < nz) {
tile[tz + HALO][ty + HALO][tx + HALO] = in[idx];
}
if (tx < HALO && gx >= HALO && gy < ny && gz < nz) {
tile[tz + HALO][ty + HALO][tx] = in[idx - HALO];
}
if (tx >= blockDim.x - HALO && gx + HALO < nx && gy < ny && gz < nz) {
tile[tz + HALO][ty + HALO][tx + 2 * HALO] = in[idx + HALO];
}
if (ty < HALO && gy >= HALO && gx < nx && gz < nz) {
tile[tz + HALO][ty][tx + HALO] = in[idx - HALO * nx];
}
if (ty >= blockDim.y - HALO && gy + HALO < ny && gx < nx && gz < nz) {
tile[tz + HALO][ty + 2 * HALO][tx + HALO] = in[idx + HALO * nx];
}
if (tz < HALO && gz >= HALO && gx < nx && gy < ny) {
tile[tz][ty + HALO][tx + HALO] = in[idx - HALO * nx * ny];
}
if (tz >= blockDim.z - HALO && gz + HALO < nz && gx < nx && gy < ny) {
tile[tz + 2 * HALO][ty + HALO][tx + HALO] = in[idx + HALO * nx * ny];
}
__syncthreads();
if (gx < nx && gy < ny && gz < nz) {
float c = tile[tz + HALO][ty + HALO][tx + HALO];
float xp = tile[tz + HALO][ty + HALO][tx + HALO + 1];
float xm = tile[tz + HALO][ty + HALO][tx + HALO - 1];
float yp = tile[tz + HALO][ty + HALO + 1][tx + HALO];
float ym = tile[tz + HALO][ty + HALO - 1][tx + HALO];
float zp = tile[tz + HALO + 1][ty + HALO][tx + HALO];
float zm = tile[tz + HALO - 1][ty + HALO][tx + HALO];
out[idx] = 0.125f * (xp + xm + yp + ym + zp + zm) + c;
}
}
int main(void) {
const int NX = 64;
const int NY = 64;
const int NZ = 64;
float* d_in = nullptr;
float* d_out = nullptr;
hipMalloc((void**)&d_in, NX * NY * NZ * sizeof(float));
hipMalloc((void**)&d_out, NX * NY * NZ * sizeof(float));
hipMemcpy(d_in, d_in, NX * NY * NZ * sizeof(float), cudaMemcpyDeviceToDevice);
dim3 grid(NX / BX, NY / BY, NZ / BZ);
dim3 block(BX, BY, BZ);
hipLaunchKernelGGL(stencil_3d, dim3(grid), dim3(block), 0, 0, d_in, d_out, NX, NY, NZ);
hipDeviceSynchronize();
hipFree(d_in);
hipFree(d_out);
return 0;
}