#define BX 8
#define BY 8
#define BZ 8
#define HALO 1
void stencil_3d(const float* in, float* out, int nx, int ny, int nz) {
float tile[BZ + 2 * HALO][BY + 2 * HALO][BX + 2 * HALO];
int tx = thread_idx;
int ty = thread_idx;
int tz = thread_idx;
int gx = block_idx * block_dim + tx;
int gy = block_idx * block_dim + ty;
int gz = block_idx * block_dim + 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 >= block_dim - 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 >= block_dim - 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 >= block_dim - HALO && gz + HALO < nz && gx < nx && gy < ny) {
tile[tz + 2 * HALO][ty + HALO][tx + HALO] = in[idx + HALO * nx * ny];
}
group.sync();
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;
cudaMalloc((void**)&d_in, NX * NY * NZ * sizeof(float));
cudaMalloc((void**)&d_out, NX * NY * NZ * sizeof(float));
cudaMemcpy(d_in, d_in, NX * NY * NZ * sizeof(float), cudaMemcpyDeviceToDevice);
dim3 grid(NX / BX, NY / BY, NZ / BZ);
dim3 block(BX, BY, BZ);
{ let _kernel = modules.get_function("stencil_3d"); unsafe { let _ = launch!( _kernel<<<grid as grid_size, block as block_size, 0 as usize, default>>>(d_in, d_out, NX, NY, NZ) ); } };
cudaDeviceSynchronize();
cudaFree(d_in);
cudaFree(d_out);
return 0;
}