Skip to content

Commit

Permalink
hip - add nontensor shared
Browse files Browse the repository at this point in the history
  • Loading branch information
jeremylt committed Jan 7, 2025
1 parent aa4002a commit 083ac09
Show file tree
Hide file tree
Showing 8 changed files with 492 additions and 6 deletions.
179 changes: 179 additions & 0 deletions backends/hip-shared/ceed-hip-shared-basis.c
Original file line number Diff line number Diff line change
Expand Up @@ -479,6 +479,124 @@ static int CeedBasisApplyAddAtPoints_Hip_shared(CeedBasis basis, const CeedInt n
return CEED_ERROR_SUCCESS;
}

//------------------------------------------------------------------------------
// Apply basis
//------------------------------------------------------------------------------
static int CeedBasisApplyNonTensorCore_Hip_shared(CeedBasis basis, bool apply_add, const CeedInt num_elem, CeedTransposeMode t_mode,
CeedEvalMode eval_mode, CeedVector u, CeedVector v) {
Ceed ceed;
Ceed_Hip *ceed_Hip;
CeedInt dim, num_comp;
const CeedScalar *d_u;
CeedScalar *d_v;
CeedBasis_Hip_shared *data;

CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
CeedCallBackend(CeedGetData(ceed, &ceed_Hip));
CeedCallBackend(CeedBasisGetData(basis, &data));
CeedCallBackend(CeedBasisGetDimension(basis, &dim));
CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));

// Get read/write access to u, v
if (u != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u));
else CeedCheck(eval_mode == CEED_EVAL_WEIGHT, ceed, CEED_ERROR_BACKEND, "An input vector is required for this CeedEvalMode");
if (apply_add) CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_DEVICE, &d_v));
else CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v));

// Apply basis operation
switch (eval_mode) {
case CEED_EVAL_INTERP: {
CeedInt P, Q;
CeedInt block_size = data->block_sizes[0];

CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &Q));
CeedInt thread = CeedIntMax(Q, P);
void *interp_args[] = {(void *)&num_elem, &data->d_interp_1d, &d_u, &d_v};

{
CeedInt elems_per_block = 64 * thread > 256 ? 256 / thread : 64;
elems_per_block = elems_per_block > 0 ? elems_per_block : 1;
CeedInt grid = num_elem / elems_per_block + (num_elem % elems_per_block > 0);
CeedInt shared_mem = elems_per_block * thread * sizeof(CeedScalar);

if (t_mode == CEED_TRANSPOSE) {
CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, apply_add ? data->InterpTransposeAdd : data->InterpTranspose, grid, thread, 1,
elems_per_block, shared_mem, interp_args));
} else {
CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->Interp, grid, thread, 1, elems_per_block, shared_mem, interp_args));
}
}
} break;
case CEED_EVAL_GRAD: {
CeedInt P, Q;
CeedInt block_size = data->block_sizes[1];

CeedCallBackend(CeedBasisGetNumNodes(basis, &P));
CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &Q));
CeedInt thread = CeedIntMax(Qd, Pd);

void *grad_args[] = {(void *)&num_elem, &data->d_interp_1d, &data->d_grad_1d, &d_u, &d_v};

{
CeedInt elems_per_block = 64 * thread > 256 ? 256 / thread : 64;
elems_per_block = elems_per_block > 0 ? elems_per_block : 1;
CeedInt grid = num_elem / elems_per_block + (num_elem % elems_per_block > 0);
CeedInt shared_mem = elems_per_block * thread * sizeof(CeedScalar);

if (t_mode == CEED_TRANSPOSE) {
CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, apply_add ? data->GradTransposeAdd : data->GradTranspose, grid, thread, 1, elems_per_block,
shared_mem, grad_args));
} else {
CeedCallBackend(CeedRunKernelDimShared_Hip(ceed, data->Grad, grid, thread, 1, elems_per_block, shared_mem, grad_args));
}
}
} break;
case CEED_EVAL_WEIGHT: {
CeedInt Q;
CeedInt block_size = data->block_sizes[2];

CeedCheck(data->d_q_weight_1d, ceed, CEED_ERROR_BACKEND, "%s not supported; q_weights_1d not set", CeedEvalModes[eval_mode]);
CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &Q));
void *weight_args[] = {(void *)&num_elem, (void *)&data->d_q_weight_1d, &d_v};

{
const CeedInt opt_elems = block_size / Q;
const CeedInt elems_per_block = opt_elems > 0 ? opt_elems : 1;
const CeedInt grid_size = num_elem / elems_per_block + (num_elem % elems_per_block > 0);

CeedCallBackend(CeedRunKernelDim_Hip(ceed, data->Weight, grid_size, Q, elems_per_block, 1, weight_args));
}
} break;
case CEED_EVAL_NONE: /* handled separately below */
break;
// LCOV_EXCL_START
case CEED_EVAL_DIV:
case CEED_EVAL_CURL:
return CeedError(ceed, CEED_ERROR_BACKEND, "%s not supported", CeedEvalModes[eval_mode]);
// LCOV_EXCL_STOP
}

// Restore vectors, cover CEED_EVAL_NONE
CeedCallBackend(CeedVectorRestoreArray(v, &d_v));
if (eval_mode == CEED_EVAL_NONE) CeedCallBackend(CeedVectorSetArray(v, CEED_MEM_DEVICE, CEED_COPY_VALUES, (CeedScalar *)d_u));
if (eval_mode != CEED_EVAL_WEIGHT) CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u));
CeedCallBackend(CeedDestroy(&ceed));
return CEED_ERROR_SUCCESS;
}

int CeedBasisApplyNonTensor_Hip_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u,
CeedVector v) {
CeedCallBackend(CeedBasisApplyNonTensorCore_Hip_shared(basis, false, num_elem, t_mode, eval_mode, u, v));
return CEED_ERROR_SUCCESS;
}

int CeedBasisApplyAddNonTensor_Hip_shared(CeedBasis basis, const CeedInt num_elem, CeedTransposeMode t_mode, CeedEvalMode eval_mode, CeedVector u,
CeedVector v) {
CeedCallBackend(CeedBasisApplyNonTensorCore_Hip_shared(basis, true, num_elem, t_mode, eval_mode, u, v));
return CEED_ERROR_SUCCESS;
}

//------------------------------------------------------------------------------
// Destroy basis
//------------------------------------------------------------------------------
Expand Down Expand Up @@ -572,3 +690,64 @@ int CeedBasisCreateTensorH1_Hip_shared(CeedInt dim, CeedInt P_1d, CeedInt Q_1d,
}

//------------------------------------------------------------------------------
// Create non-tensor basis
//------------------------------------------------------------------------------
int CeedBasisCreateH1_Hip_shared(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis) {
Ceed ceed;
CeedInt num_comp, q_comp_interp, q_comp_grad;
const CeedInt q_bytes = num_qpts * sizeof(CeedScalar);
CeedBasis_Hip_shared *data;

CeedCallBackend(CeedBasisGetCeed(basis, &ceed));
CeedCallBackend(CeedCalloc(1, &data));

// Check max sizes
CeedCheck(dim <= 3, ceed, CEED_ERROR_BACKEND, "Backend does not implement nontensor bases with dim > 3");
CeedCheck(num_nodes * num_qpts * dim < 52 * 52 * 3, ceed, CEED_ERROR_BACKEND, "Backend does not implement nontensor bases with P * Q this large");

// Copy basis data to GPU
CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_INTERP, &q_comp_interp));
CeedCallBackend(CeedBasisGetNumQuadratureComponents(basis, CEED_EVAL_GRAD, &q_comp_grad));
if (q_weight) {
CeedCallHip(ceed, hipMalloc((void **)&data->d_q_weight_1d, q_bytes));
CeedCallHip(ceed, hipMemcpy(data->d_q_weight_1d, q_weight, q_bytes, hipMemcpyHostToDevice));
}
if (interp) {
const CeedInt interp_bytes = q_bytes * num_nodes * q_comp_interp;

CeedCallHip(ceed, hipMalloc((void **)&data->d_interp_1d, interp_bytes));
CeedCallHip(ceed, hipMemcpy(data->d_interp_1d, interp, interp_bytes, hipMemcpyHostToDevice));
}
if (grad) {
const CeedInt grad_bytes = q_bytes * num_nodes * q_comp_grad;

CeedCallHip(ceed, hipMalloc((void **)&data->d_grad_1d, grad_bytes));
CeedCallHip(ceed, hipMemcpy(data->d_grad_1d, grad, grad_bytes, hipMemcpyHostToDevice));
}

// Compile basis kernels
const char basis_kernel_source[] = "// Non-tensor basis source\n#include <ceed/jit-source/hip/hip-shared-basis-nontensor.h>\n";

CeedCallBackend(CeedBasisGetNumComponents(basis, &num_comp));
CeedCallBackend(CeedCompile_Hip(ceed, basis_kernel_source, &data->module, 5, "BASIS_Q", num_qpts, "BASIS_P", num_nodes, "T_1D",
CeedIntMax(num_qpts, num_nodes), "BASIS_DIM", dim, "BASIS_NUM_COMP", num_comp));
CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "Interp", &data->Interp));
CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "InterpTranspose", &data->InterpTranspose));
CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "InterpTransposeAdd", &data->InterpTransposeAdd));
CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "Grad", &data->Grad));
CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "GradTranspose", &data->GradTranspose));
CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "GradTransposeAdd", &data->GradTransposeAdd));
CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, "Weight", &data->Weight));

CeedCallBackend(CeedBasisSetData(basis, data));

// Register backend functions
CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Apply", CeedBasisApplyNonTensor_Hip_shared));
CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "ApplyAdd", CeedBasisApplyAddNonTensor_Hip_shared));
CeedCallBackend(CeedSetBackendFunction(ceed, "Basis", basis, "Destroy", CeedBasisDestroy_Hip_shared));
CeedCallBackend(CeedDestroy(&ceed));
return CEED_ERROR_SUCCESS;
}

//------------------------------------------------------------------------------
1 change: 1 addition & 0 deletions backends/hip-shared/ceed-hip-shared.c
Original file line number Diff line number Diff line change
Expand Up @@ -36,6 +36,7 @@ static int CeedInit_Hip_shared(const char *resource, Ceed ceed) {
CeedCallBackend(CeedDestroy(&ceed_ref));

CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "BasisCreateTensorH1", CeedBasisCreateTensorH1_Hip_shared));
CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "BasisCreateH1", CeedBasisCreateH1_Hip_shared));
CeedCallBackend(CeedSetBackendFunction(ceed, "Ceed", ceed, "Destroy", CeedDestroy_Hip));
return CEED_ERROR_SUCCESS;
}
Expand Down
3 changes: 3 additions & 0 deletions backends/hip-shared/ceed-hip-shared.h
Original file line number Diff line number Diff line change
Expand Up @@ -38,3 +38,6 @@ typedef struct {

CEED_INTERN int CeedBasisCreateTensorH1_Hip_shared(CeedInt dim, CeedInt P_1d, CeedInt Q_1d, const CeedScalar *interp_1d, const CeedScalar *grad_1d,
const CeedScalar *q_ref_1d, const CeedScalar *q_weight_1d, CeedBasis basis);

CEED_INTERN int CeedBasisCreateH1_Hip_shared(CeedElemTopology topo, CeedInt dim, CeedInt num_nodes, CeedInt num_qpts, const CeedScalar *interp,
const CeedScalar *grad, const CeedScalar *q_ref, const CeedScalar *q_weight, CeedBasis basis);
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,7 @@
// This file is part of CEED: http://github.com/ceed

/// @file
/// Internal header for CUDA shared memory non-tensor product basis templates
/// Internal header for CUDA shared memory non-tensor basis templates
#include <ceed/types.h>

//------------------------------------------------------------------------------
Expand Down
6 changes: 3 additions & 3 deletions include/ceed/jit-source/cuda/cuda-shared-basis-nontensor.h
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@
#include "cuda-shared-basis-read-write-templates.h"

//------------------------------------------------------------------------------
// Interp kernel by dim
// Interp kernels
//------------------------------------------------------------------------------
extern "C" __global__ void Interp(const CeedInt num_elem, const CeedScalar *c_B, const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
extern __shared__ CeedScalar slice[];
Expand Down Expand Up @@ -96,7 +96,7 @@ extern "C" __global__ void InterpTransposeAdd(const CeedInt num_elem, const Ceed
}

//------------------------------------------------------------------------------
// Grad kernel by dim
// Grad kernels
//------------------------------------------------------------------------------
extern "C" __global__ void Grad(const CeedInt num_elem, const CeedScalar *c_G, const CeedScalar *__restrict__ d_U, CeedScalar *__restrict__ d_V) {
extern __shared__ CeedScalar slice[];
Expand Down Expand Up @@ -179,7 +179,7 @@ extern "C" __global__ void GradTransposeAdd(const CeedInt num_elem, const CeedSc
}

//------------------------------------------------------------------------------
// Weight kernels by dim
// Weight kernel
//------------------------------------------------------------------------------
extern "C" __global__ void Weight(const CeedInt num_elem, const CeedScalar *__restrict__ q_weight, CeedScalar *__restrict__ d_W) {
extern __shared__ CeedScalar slice[];
Expand Down
98 changes: 98 additions & 0 deletions include/ceed/jit-source/hip/hip-shared-basis-nontensor-templates.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,98 @@
// Copyright (c) 2017-2024, Lawrence Livermore National Security, LLC and other CEED contributors.
// All Rights Reserved. See the top-level LICENSE and NOTICE files for details.
//
// SPDX-License-Identifier: BSD-2-Clause
//
// This file is part of CEED: http://github.com/ceed

/// @file
/// Internal header for HIP shared memory non-tensor basis templates
#include <ceed/types.h>

//------------------------------------------------------------------------------
// 1D tensor contraction
//------------------------------------------------------------------------------
template <int NUM_COMP, int P_1D, int Q_1D>
inline __device__ void Contract1d(SharedData_Hip &data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) {
data.slice[data.t_id_x] = *U;
__syncthreads();
*V = 0.0;
if (data.t_id_x < Q_1D) {
for (CeedInt i = 0; i < P_1D; i++) {
*V += B[i + data.t_id_x * P_1D] * data.slice[i]; // Contract x direction
}
}
__syncthreads();
}

//------------------------------------------------------------------------------
// 1D transpose tensor contraction
//------------------------------------------------------------------------------
template <int NUM_COMP, int P_1D, int Q_1D>
inline __device__ void ContractTranspose1d(SharedData_Hip &data, const CeedScalar *U, const CeedScalar *B, CeedScalar *V) {
data.slice[data.t_id_x] = *U;
__syncthreads();
if (data.t_id_x < P_1D) {
for (CeedInt i = 0; i < Q_1D; i++) {
*V += B[data.t_id_x + i * P_1D] * data.slice[i]; // Contract x direction
}
}
__syncthreads();
}

//------------------------------------------------------------------------------
// Interpolate to quadrature points
//------------------------------------------------------------------------------
template <int NUM_COMP, int P, int Q>
inline __device__ void Interp1d(SharedData_Hip &data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B, CeedScalar *__restrict__ r_V) {
for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
Contract1d<NUM_COMP, P, Q>(data, &r_U[comp], c_B, &r_V[comp]);
}
}

//------------------------------------------------------------------------------
// Interpolate transpose
//------------------------------------------------------------------------------
template <int NUM_COMP, int P, int Q>
inline __device__ void InterpTranspose1d(SharedData_Hip &data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_B,
CeedScalar *__restrict__ r_V) {
for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
r_V[comp] = 0.0;
ContractTranspose1d<NUM_COMP, P, Q>(data, &r_U[comp], c_B, &r_V[comp]);
}
}

//------------------------------------------------------------------------------
// Derivatives at quadrature points
//------------------------------------------------------------------------------
template <int NUM_COMP, int P, int Q>
inline __device__ void Grad1d(SharedData_Hip &data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_G,
CeedScalar *__restrict__ r_V) {
for (CeedInt dim = 0; dim < DIM; dim++) {
for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
Contract1d<NUM_COMP, P, Q>(data, &r_U[comp], &c_G[dim * P * Q], &r_V[comp + dim * NUM_COMP]);
}
}
}

//------------------------------------------------------------------------------
// Derivatives transpose
//------------------------------------------------------------------------------
template <int NUM_COMP, int P, int Q>
inline __device__ void GradTranspose1d(SharedData_Hip &data, const CeedScalar *__restrict__ r_U, const CeedScalar *c_G,
CeedScalar *__restrict__ r_V) {
for (CeedInt comp = 0; comp < NUM_COMP; comp++) r_V[comp] = 0.0;
for (CeedInt dim = 0; dim < DIM; dim++) {
for (CeedInt comp = 0; comp < NUM_COMP; comp++) {
ContractTranspose1d<NUM_COMP, P, Q>(data, &r_U[comp + dim * NUM_COMP], &c_G[dim * P * Q], &r_V[comp]);
}
}
}

//------------------------------------------------------------------------------
// Quadrature weights
//------------------------------------------------------------------------------
template <int Q>
inline __device__ void Weight1d(SharedData_Hip &data, const CeedScalar *__restrict__ q_weight_1d, CeedScalar *w) {
*w = (data.t_id_x < Q) ? q_weight_1d[data.t_id_x] : 0.0;
}
Loading

0 comments on commit 083ac09

Please sign in to comment.