diff --git a/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp b/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp index 99da342f66..d19eedd491 100644 --- a/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp +++ b/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp @@ -34,6 +34,7 @@ static int CeedOperatorBuildKernelData_Cuda_gen(Ceed ceed, CeedInt num_input_fie *Q_1d = 0; *dim = 0; *is_tensor = true; + for (CeedInt i = 0; i < num_input_fields; i++) { CeedBasis basis; @@ -44,14 +45,15 @@ static int CeedOperatorBuildKernelData_Cuda_gen(Ceed ceed, CeedInt num_input_fie // Collect dim, P_1d, and Q_1d CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); - CeedCheck(is_field_tensor, ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); *is_tensor = *is_tensor && is_field_tensor; - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); + if (is_field_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &field_P_1d)); *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); CeedCheck(*dim == 0 || field_dim == *dim, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *dim = field_dim; - CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); + if (is_field_tensor) CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); + else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &field_Q_1d)); CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *Q_1d = field_Q_1d; } @@ -67,14 +69,15 @@ static int CeedOperatorBuildKernelData_Cuda_gen(Ceed ceed, CeedInt num_input_fie // Collect dim, P_1d, and Q_1d CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); - CeedCheck(is_field_tensor, ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); *is_tensor = *is_tensor && is_field_tensor; - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); + if (is_field_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &field_P_1d)); *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); CeedCheck(*dim == 0 || field_dim == *dim, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *dim = field_dim; - CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); + if (is_field_tensor) CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); + else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &field_Q_1d)); CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *Q_1d = field_Q_1d; } @@ -124,10 +127,10 @@ static int CeedOperatorBuildKernelData_Cuda_gen(Ceed ceed, CeedInt num_input_fie // Setup fields //------------------------------------------------------------------------------ static int CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, CeedInt i, CeedOperatorField op_field, - CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool is_at_points, + CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool is_tensor, bool is_at_points, bool use_3d_slices) { std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); - std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; + std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix, Q_name = is_tensor ? "Q_1d" : "Q"; std::string option_name = (is_input ? "inputs" : "outputs"); CeedEvalMode eval_mode = CEED_EVAL_NONE; CeedInt elem_size = 0, num_comp = 0, P_1d = 0; @@ -147,7 +150,8 @@ static int CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream &code, C CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); if (basis != CEED_BASIS_NONE) { CeedCallBackend(CeedBasisGetData(basis, &basis_data)); - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d)); } CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); @@ -184,7 +188,7 @@ static int CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream &code, C if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; else data->B.outputs[i] = basis_data->d_interp_1d; } - code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; + code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_name << "*" << Q_name << "];\n"; code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; break; case CEED_EVAL_GRAD: @@ -209,13 +213,15 @@ static int CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream &code, C if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; else data->B.outputs[i] = basis_data->d_interp_1d; } - code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; - code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; + if (is_tensor) { + code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_name << "*" << Q_name << "];\n"; + code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; + } if (is_at_points) break; // No G mat for AtPoints if (use_3d_slices) { if (is_input) data->G.inputs[i] = basis_data->d_collo_grad_1d; else data->G.outputs[i] = basis_data->d_collo_grad_1d; - code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * Q_1d << "];\n"; + code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_name << "*" << Q_name << "];\n"; code << " LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; } else { bool has_collo_grad = basis_data->d_collo_grad_1d; @@ -223,11 +229,12 @@ static int CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream &code, C if (is_input) data->G.inputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; else data->G.outputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; if (has_collo_grad) { - code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * Q_1d << "];\n"; + code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_name << "*" << Q_name << "];\n"; code << " LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; } else { - code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * P_1d << "];\n"; - code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; + code << " __shared__ CeedScalar s_G" << var_suffix << "[" << P_name << "*" << Q_name << (is_tensor ? "" : "*dim") << "];\n"; + code << " LoadMatrix<" << P_name << ", " << Q_name << (is_tensor ? "" : "*dim") << ">(data, G." << option_name << "[" << i << "], s_G" + << var_suffix << ");\n"; } } break; @@ -248,9 +255,9 @@ static int CeedOperatorBuildKernelFieldData_Cuda_gen(std::ostringstream &code, C //------------------------------------------------------------------------------ static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, CeedInt i, CeedInt dim, CeedInt field_input_buffer[], CeedOperatorField op_field, CeedQFunctionField qf_field, - CeedInt Q_1d, bool is_input, bool is_at_points, bool use_3d_slices) { + CeedInt Q_1d, bool is_input, bool is_tensor, bool is_at_points, bool use_3d_slices) { std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); - std::string P_name = "P_1d" + var_suffix; + std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix; CeedEvalMode eval_mode = CEED_EVAL_NONE; CeedInt elem_size = 0, num_comp = 0, P_1d = 0; CeedSize l_size; @@ -269,7 +276,8 @@ static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, } CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); if (basis != CEED_BASIS_NONE) { - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d)); } CeedCallBackend(CeedBasisDestroy(&basis)); CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); @@ -299,8 +307,8 @@ static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); code << " // CompStride: " << comp_stride << "\n"; data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; - code << " ReadLVecStandard" << dim << "d(data, l_size" - << var_suffix << ", elem, indices.inputs[" << i << "], d" << var_suffix << ", r_e" << var_suffix << ");\n"; + code << " ReadLVecStandard" << (is_tensor ? dim : 1) << "d(data, l_size" << var_suffix << ", elem, indices.inputs[" << i << "], d" << var_suffix << ", r_e" << var_suffix << ");\n"; break; } case CEED_RESTRICTION_STRIDED: { @@ -315,8 +323,8 @@ static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); } code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; - code << " ReadLVecStrided" << dim << "d(data, elem, d" << var_suffix << ", r_e" << var_suffix << ");\n"; + code << " ReadLVecStrided" << (is_tensor ? dim : 1) << "d(data, elem, d" << var_suffix << ", r_e" << var_suffix << ");\n"; break; } case CEED_RESTRICTION_POINTS: { @@ -345,8 +353,8 @@ static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); code << " // CompStride: " << comp_stride << "\n"; data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets; - code << " WriteLVecStandard" << dim << "d(data, l_size" - << var_suffix << ", elem, indices.outputs[" << i << "], r_e" << var_suffix << ", d" << var_suffix << ");\n"; + code << " WriteLVecStandard" << (is_tensor ? dim : 1) << "d(data, l_size" << var_suffix << ", elem, indices.outputs[" << i << "], r_e" << var_suffix << ", d" << var_suffix << ");\n"; break; } case CEED_RESTRICTION_STRIDED: { @@ -361,8 +369,8 @@ static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); } code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; - code << " WriteLVecStrided" << dim << "d(data, elem, r_e" << var_suffix << ", d" << var_suffix << ");\n"; + code << " WriteLVecStrided" << (is_tensor ? dim : 1) << "d(data, elem, r_e" << var_suffix << ", d" << var_suffix << ");\n"; break; } case CEED_RESTRICTION_POINTS: @@ -383,10 +391,10 @@ static int CeedOperatorBuildKernelRestriction_Cuda_gen(std::ostringstream &code, // Basis //------------------------------------------------------------------------------ static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedOperator_Cuda_gen *data, CeedInt i, CeedInt dim, - CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, + CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool is_tensor, bool is_at_points, bool use_3d_slices) { std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); - std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; + std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix, Q_name = is_tensor ? "Q_1d" : "Q"; CeedEvalMode eval_mode = CEED_EVAL_NONE; CeedInt elem_size = 0, num_comp = 0, P_1d = 0; CeedElemRestriction elem_rstr; @@ -401,7 +409,8 @@ static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedO CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); if (basis != CEED_BASIS_NONE) { - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d)); } CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); @@ -416,29 +425,45 @@ static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedO break; case CEED_EVAL_INTERP: if (is_at_points) { - code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; - code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d"; + + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", r_c" << var_suffix << ");\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; - code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; + std::string function_name = is_tensor ? ((dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d") : "InterpNonTensor"; + + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", r_q" << var_suffix << ");\n"; } break; case CEED_EVAL_GRAD: if (is_at_points) { - code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; - code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d"; + + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", r_c" << var_suffix << ");\n"; } else if (use_3d_slices) { + std::string function_name = (dim > 1 ? "InterpTensor" : "Interp") + std::to_string(dim) + "d"; + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; - code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", r_q" << var_suffix << ");\n"; + } else if (is_tensor) { + bool is_collocated = dim == 3 && Q_1d >= P_1d; + std::string function_name = (dim == 1 ? "Grad" : (is_collocated ? "GradTensorCollocated" : "GradTensor")) + std::to_string(dim) + "d"; + + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << (dim >= 3 ? Q_name : "1") << "];\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; - code << " Grad" << (dim > 1 ? "Tensor" : "") << (dim == 3 && Q_1d >= P_1d ? "Collocated" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix << ", r_q" - << var_suffix << ");\n"; + std::string function_name = "GradNonTensor"; + + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; + code << " " << function_name << "(data, r_e" << var_suffix + << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n"; } break; case CEED_EVAL_WEIGHT: { @@ -446,11 +471,12 @@ static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedO code << " // Nothing to do AtPoints\n"; } else { CeedBasis_Cuda_shared *basis_data; + std::string function_name = is_tensor ? ((dim == 1 ? "Weight" : "WeightTensor") + std::to_string(dim) + "d") : "WeightNonTensor"; - code << " CeedScalar r_q" << var_suffix << "[" << Q_name << "];\n"; + code << " CeedScalar r_q" << var_suffix << "[" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; CeedCallBackend(CeedBasisGetData(basis, &basis_data)); data->W = basis_data->d_q_weight_1d; - code << " Weight" << (dim > 1 ? "Tensor" : "") << dim << "d<" << Q_name << ">(data, W, r_q" << var_suffix << ");\n"; + code << " " << function_name << "<" << Q_name << ">(data, W, r_q" << var_suffix << ");\n"; } break; } @@ -468,25 +494,42 @@ static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedO case CEED_EVAL_INTERP: code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; if (is_at_points) { - code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_c" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d"; + + code << " " << function_name << "(data, r_c" << var_suffix << ", s_B" + << var_suffix << ", r_e" << var_suffix << ");\n"; } else { - code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + std::string function_name = + is_tensor ? ((dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d") : "InterpTransposeNonTensor"; + + code << " " << function_name << "(data, r_q" << var_suffix << ", s_B" + << var_suffix << ", r_e" << var_suffix << ");\n"; } break; case CEED_EVAL_GRAD: code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; if (is_at_points) { - code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_c" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d"; + + code << " " << function_name << "(data, r_c" << var_suffix << ", s_B" + << var_suffix << ", r_e" << var_suffix << ");\n"; } else if (use_3d_slices) { - code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d"; + + code << " " << function_name << "(data, r_q" << var_suffix << ", s_B" + << var_suffix << ", r_e" << var_suffix << ");\n"; + } else if (is_tensor) { + bool is_collocated = dim == 3 && Q_1d >= P_1d; + std::string function_name = + (dim == 1 ? "GradTranspose" : (is_collocated ? "GradTransposeTensorCollocated" : "GradTransposeTensor")) + std::to_string(dim) + "d"; + + code << " " << function_name << "(data, r_q" << var_suffix << ", s_B" + << var_suffix << ", s_G" << var_suffix << ", r_e" << var_suffix << ");\n"; } else { - code << " GradTranspose" << (dim > 1 ? "Tensor" : "") << (dim == 3 && Q_1d >= P_1d ? "Collocated" : "") << dim << "d(data, r_q" << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix - << ", r_e" << var_suffix << ");\n"; + std::string function_name = "GradTransposeNonTensor"; + + code << " " << function_name << "(data, r_q" << var_suffix + << ", s_G" << var_suffix << ", r_e" << var_suffix << ");\n"; } break; // LCOV_EXCL_START @@ -509,8 +552,9 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C CeedInt num_input_fields, CeedOperatorField *op_input_fields, CeedQFunctionField *qf_input_fields, CeedInt num_output_fields, CeedOperatorField *op_output_fields, CeedQFunctionField *qf_output_fields, - std::string qfunction_name, CeedInt Q_1d, bool is_at_points, bool use_3d_slices) { - std::string Q_name = "Q_1d"; + std::string qfunction_name, CeedInt Q_1d, bool is_tensor, bool is_at_points, + bool use_3d_slices) { + std::string Q_name = is_tensor ? "Q_1d" : "Q"; CeedEvalMode eval_mode = CEED_EVAL_NONE; CeedElemRestriction elem_rstr; @@ -526,25 +570,25 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C if (is_at_points) { code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "];\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; } break; case CEED_EVAL_INTERP: if (is_at_points) { // Accumulator for point data - code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; - code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "; i++) {\n"; + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n"; + code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "; i++) {\n"; code << " r_c" << var_suffix << "[i] = 0.0;\n"; code << " }\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; } break; case CEED_EVAL_GRAD: if (is_at_points) { // Accumulator for point data - code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "*dim];\n"; - code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "; i++) {\n"; + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "*dim];\n"; + code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "; i++) {\n"; code << " r_c" << var_suffix << "[i] = 0.0;\n"; code << " }\n"; } else if (use_3d_slices) { @@ -554,7 +598,7 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C code << " r_q" << var_suffix << "[i] = 0.0;\n"; code << " }\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; } break; case CEED_EVAL_WEIGHT: @@ -675,7 +719,7 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); } code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; - code << " ReadEVecSliceStrided3d(data, elem, q, d" << var_suffix << ", r_s" << var_suffix << ");\n"; } else { CeedSize l_size = 0; @@ -774,10 +818,10 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C // Apply QFunction code << "\n // -- Apply QFunction\n"; code << " " << qfunction_name << "(ctx, "; - if (dim != 3 || is_at_points || use_3d_slices) { + if (dim != 3 || is_at_points || use_3d_slices || !is_tensor) { code << "1"; } else { - code << "Q_1d"; + code << Q_name; } code << ", inputs, outputs);\n"; @@ -950,11 +994,17 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { } // Load basis source files - // TODO: Add non-tensor, AtPoints - code << "// Tensor basis source\n"; - code << "#include \n\n"; - code << "// AtPoints basis source\n"; - code << "#include \n\n"; + if (is_tensor) { + code << "// Tensor basis source\n"; + code << "#include \n\n"; + } else { + code << "// Non-tensor basis source\n"; + code << "#include \n\n"; + } + if (is_at_points) { + code << "// AtPoints basis source\n"; + code << "#include \n\n"; + } code << "// CodeGen operator source\n"; code << "#include \n\n"; @@ -966,7 +1016,7 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { // Define CEED_Q_VLA code << "\n#undef CEED_Q_VLA\n"; - if (dim != 3 || is_at_points || use_3d_slices) { + if (dim != 3 || is_at_points || use_3d_slices || !is_tensor) { code << "#define CEED_Q_VLA 1\n\n"; } else { code << "#define CEED_Q_VLA " << Q_1d << "\n\n"; @@ -990,7 +1040,7 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { code << "// d_[in,out]_i: CeedVector device array\n"; code << "// r_[in,out]_e_i: Element vector register\n"; code << "// r_[in,out]_q_i: Quadrature space vector register\n"; - code << "// r_[in,out]_c_i: AtPoints Chebyshev coefficents register\n"; + code << "// r_[in,out]_c_i: AtPoints Chebyshev coefficients register\n"; code << "// r_[in,out]_s_i: Quadrature space slice vector register\n"; code << "// \n"; code << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n"; @@ -1014,7 +1064,7 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { } code << " const CeedInt dim = " << dim << ";\n"; - code << " const CeedInt Q_1d = " << Q_1d << ";\n"; + code << " const CeedInt " << (is_tensor ? "Q_1d" : "Q") << " = " << Q_1d << ";\n"; if (is_at_points) { code << " const CeedInt max_num_points = " << max_num_points << ";\n"; code << " const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n"; @@ -1027,18 +1077,18 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { code << " data.t_id_y = threadIdx.y;\n"; code << " data.t_id_z = threadIdx.z;\n"; code << " data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n"; - code << " data.slice = slice + data.t_id_z*T_1D" << (dim > 1 ? "*T_1D" : "") << ";\n"; + code << " data.slice = slice + data.t_id_z*T_1D" << ((!is_tensor || dim == 1) ? "" : "*T_1D") << ";\n"; // Initialize constants, and matrices B and G code << "\n // Input field constants and basis data\n"; for (CeedInt i = 0; i < num_input_fields; i++) { - CeedCallBackend( - CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, i, op_input_fields[i], qf_input_fields[i], Q_1d, true, is_at_points, use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, i, op_input_fields[i], qf_input_fields[i], Q_1d, true, is_tensor, + is_at_points, use_3d_slices)); } code << "\n // Output field constants and basis data\n"; for (CeedInt i = 0; i < num_output_fields; i++) { - CeedCallBackend( - CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, i, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_at_points, use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelFieldData_Cuda_gen(code, data, i, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_tensor, + is_at_points, use_3d_slices)); } // Loop over all elements @@ -1056,7 +1106,7 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); - max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * elem_size); + max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_tensor && (dim >= 3) ? elem_size : 1)); CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); } for (CeedInt i = 0; i < num_output_fields; i++) { @@ -1066,7 +1116,7 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); - max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * elem_size); + max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_tensor && (dim >= 3) ? elem_size : 1)); CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); } code << " // Scratch restriction buffer space\n"; @@ -1125,17 +1175,17 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { // ---- Restriction CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, f, dim, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f], - Q_1d, true, is_at_points, use_3d_slices)); + Q_1d, true, is_tensor, is_at_points, use_3d_slices)); // ---- Basis action - CeedCallBackend( - CeedOperatorBuildKernelBasis_Cuda_gen(code, data, f, dim, op_input_fields[f], qf_input_fields[f], Q_1d, true, is_at_points, use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, f, dim, op_input_fields[f], qf_input_fields[f], Q_1d, true, is_tensor, + is_at_points, use_3d_slices)); } // -- Q function CeedCallBackend(CeedOperatorBuildKernelQFunction_Cuda_gen(code, data, dim, max_num_points, num_input_fields, op_input_fields, qf_input_fields, - num_output_fields, op_output_fields, qf_output_fields, qfunction_name, Q_1d, is_at_points, - use_3d_slices)); + num_output_fields, op_output_fields, qf_output_fields, qfunction_name, Q_1d, is_tensor, + is_at_points, use_3d_slices)); // -- Output basis and restriction code << "\n // -- Output field basis action and restrictions\n"; @@ -1143,12 +1193,12 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { code << " // ---- Output field " << i << "\n"; // ---- Basis action - CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, i, dim, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_at_points, - use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelBasis_Cuda_gen(code, data, i, dim, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_tensor, + is_at_points, use_3d_slices)); // ---- Restriction CeedCallBackend(CeedOperatorBuildKernelRestriction_Cuda_gen(code, data, i, dim, NULL, op_output_fields[i], qf_output_fields[i], Q_1d, false, - is_at_points, use_3d_slices)); + is_tensor, is_at_points, use_3d_slices)); } // Close loop and function diff --git a/backends/cuda-gen/ceed-cuda-gen-operator.c b/backends/cuda-gen/ceed-cuda-gen-operator.c index 570e8735c1..43d388e293 100644 --- a/backends/cuda-gen/ceed-cuda-gen-operator.c +++ b/backends/cuda-gen/ceed-cuda-gen-operator.c @@ -11,6 +11,7 @@ #include #include #include +#include #include "../cuda/ceed-cuda-common.h" #include "../cuda/ceed-cuda-compile.h" @@ -98,7 +99,7 @@ static size_t dynamicSMemSize(int threads) { return threads * sizeof(CeedScalar) // Apply and add to output //------------------------------------------------------------------------------ static int CeedOperatorApplyAdd_Cuda_gen(CeedOperator op, CeedVector input_vec, CeedVector output_vec, CeedRequest *request) { - bool is_at_points; + bool is_at_points, is_tensor; Ceed ceed; Ceed_Cuda *cuda_data; CeedInt num_elem, num_input_fields, num_output_fields; @@ -110,16 +111,62 @@ static int CeedOperatorApplyAdd_Cuda_gen(CeedOperator op, CeedVector input_vec, CeedOperatorField *op_input_fields, *op_output_fields; CeedOperator_Cuda_gen *data; - // Check for tensor-product bases + // Check for shared bases + CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); { - bool has_tensor_bases; + bool has_shared_bases = true, is_all_tensor = true, is_all_nontensor = true; + + for (CeedInt i = 0; i < num_input_fields; i++) { + CeedBasis basis; + + CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); + if (basis != CEED_BASIS_NONE) { + bool is_tensor = true; + const char *resource; + char *resource_root; + Ceed basis_ceed; + + CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); + is_all_tensor &= is_tensor; + is_all_nontensor &= !is_tensor; + CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed)); + CeedCallBackend(CeedGetResource(basis_ceed, &resource)); + CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root)); + has_shared_bases &= !strcmp(resource_root, "/gpu/cuda/shared"); + CeedCallBackend(CeedFree(&resource_root)); + CeedCallBackend(CeedDestroy(&basis_ceed)); + } + CeedCallBackend(CeedBasisDestroy(&basis)); + } - CeedCallBackend(CeedOperatorHasTensorBases(op, &has_tensor_bases)); - // -- Fallback to ref if not all bases are tensor-product - if (!has_tensor_bases) { + for (CeedInt i = 0; i < num_output_fields; i++) { + CeedBasis basis; + + CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); + if (basis != CEED_BASIS_NONE) { + bool is_tensor = true; + const char *resource; + char *resource_root; + Ceed basis_ceed; + + CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); + is_all_tensor &= is_tensor; + is_all_nontensor &= !is_tensor; + + CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed)); + CeedCallBackend(CeedGetResource(basis_ceed, &resource)); + CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root)); + has_shared_bases &= !strcmp(resource_root, "/gpu/cuda/shared"); + CeedCallBackend(CeedFree(&resource_root)); + CeedCallBackend(CeedDestroy(&basis_ceed)); + } + CeedCallBackend(CeedBasisDestroy(&basis)); + } + // -- Fallback to ref if not all bases are shared + if (!has_shared_bases || (!is_all_tensor && !is_all_nontensor)) { CeedOperator op_fallback; - CeedDebug256(CeedOperatorReturnCeed(op), CEED_DEBUG_COLOR_SUCCESS, "Falling back to /gpu/cuda/ref CeedOperator due to non-tensor bases"); + CeedDebug256(CeedOperatorReturnCeed(op), CEED_DEBUG_COLOR_SUCCESS, "Falling back to /gpu/cuda/ref CeedOperator due unsupported bases"); CeedCallBackend(CeedOperatorGetFallback(op, &op_fallback)); CeedCallBackend(CeedOperatorApplyAdd(op_fallback, input_vec, output_vec, request)); return CEED_ERROR_SUCCESS; @@ -132,7 +179,6 @@ static int CeedOperatorApplyAdd_Cuda_gen(CeedOperator op, CeedVector input_vec, CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); CeedCallBackend(CeedQFunctionGetData(qf, &qf_data)); CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); - CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); // Creation of the operator @@ -232,11 +278,19 @@ static int CeedOperatorApplyAdd_Cuda_gen(CeedOperator op, CeedVector input_vec, const CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); int max_threads_per_block, min_grid_size, grid; + CeedCallBackend(CeedOperatorHasTensorBases(op, &is_tensor)); CeedCallCuda(ceed, cuOccupancyMaxPotentialBlockSize(&min_grid_size, &max_threads_per_block, data->op, dynamicSMemSize, 0, 0x10000)); - int block[3] = {thread_1d, dim < 2 ? 1 : thread_1d, -1}; + int block[3] = {thread_1d, ((!is_tensor || dim == 1) ? 1 : thread_1d), -1}; - CeedCallBackend(BlockGridCalculate(num_elem, min_grid_size / cuda_data->device_prop.multiProcessorCount, max_threads_per_block, - cuda_data->device_prop.maxThreadsDim[2], cuda_data->device_prop.warpSize, block, &grid)); + if (is_tensor) { + CeedCallBackend(BlockGridCalculate(num_elem, min_grid_size / cuda_data->device_prop.multiProcessorCount, max_threads_per_block, + cuda_data->device_prop.maxThreadsDim[2], cuda_data->device_prop.warpSize, block, &grid)); + } else { + CeedInt elems_per_block = CeedIntMin(cuda_data->device_prop.maxThreadsDim[2], CeedIntMax(512 / thread_1d, 1)); + + grid = num_elem / elems_per_block + (num_elem % elems_per_block > 0); + block[2] = elems_per_block; + } CeedInt shared_mem = block[0] * block[1] * block[2] * sizeof(CeedScalar); CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->op, grid, block[0], block[1], block[2], shared_mem, opargs)); diff --git a/backends/cuda-ref/ceed-cuda-ref-operator.c b/backends/cuda-ref/ceed-cuda-ref-operator.c index fe873c74b5..319d10145f 100644 --- a/backends/cuda-ref/ceed-cuda-ref-operator.c +++ b/backends/cuda-ref/ceed-cuda-ref-operator.c @@ -132,7 +132,7 @@ static int CeedOperatorSetupFields_Cuda(CeedQFunction qf, CeedOperator op, bool // Input CEED_VECTOR_ACTIVE // Output CEED_VECTOR_ACTIVE without CEED_EVAL_NONE // Input CEED_VECTOR_NONE with CEED_EVAL_WEIGHT - // Input passive vectorr with CEED_EVAL_NONE and strided restriction with CEED_STRIDES_BACKEND + // Input passive vector with CEED_EVAL_NONE and strided restriction with CEED_STRIDES_BACKEND CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &l_vec)); is_active = l_vec == CEED_VECTOR_ACTIVE; CeedCallBackend(CeedVectorDestroy(&l_vec)); diff --git a/backends/hip-gen/ceed-hip-gen-operator-build.cpp b/backends/hip-gen/ceed-hip-gen-operator-build.cpp index bd56a77350..1df23c3d7e 100644 --- a/backends/hip-gen/ceed-hip-gen-operator-build.cpp +++ b/backends/hip-gen/ceed-hip-gen-operator-build.cpp @@ -61,6 +61,7 @@ static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fiel *Q_1d = 0; *dim = 0; *is_tensor = true; + for (CeedInt i = 0; i < num_input_fields; i++) { CeedBasis basis; @@ -71,14 +72,15 @@ static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fiel // Collect dim, P_1d, and Q_1d CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); - CeedCheck(is_field_tensor, ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); *is_tensor = *is_tensor && is_field_tensor; - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); + if (is_field_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &field_P_1d)); *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); CeedCheck(*dim == 0 || field_dim == *dim, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *dim = field_dim; - CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); + if (is_field_tensor) CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); + else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &field_Q_1d)); CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *Q_1d = field_Q_1d; } @@ -94,14 +96,15 @@ static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fiel // Collect dim, P_1d, and Q_1d CeedCallBackend(CeedBasisIsTensor(basis, &is_field_tensor)); - CeedCheck(is_field_tensor, ceed, CEED_ERROR_BACKEND, "Backend does not implement operators with non-tensor basis"); *is_tensor = *is_tensor && is_field_tensor; - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); + if (is_field_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &field_P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &field_P_1d)); *max_P_1d = CeedIntMax(*max_P_1d, field_P_1d); CeedCallBackend(CeedBasisGetDimension(basis, &field_dim)); CeedCheck(*dim == 0 || field_dim == *dim, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *dim = field_dim; - CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); + if (is_field_tensor) CeedCallBackend(CeedBasisGetNumQuadraturePoints1D(basis, &field_Q_1d)); + else CeedCallBackend(CeedBasisGetNumQuadraturePoints(basis, &field_Q_1d)); CeedCheck(*Q_1d == 0 || field_Q_1d == *Q_1d, ceed, CEED_ERROR_BACKEND, "Quadrature spaces must be compatible"); *Q_1d = field_Q_1d; } @@ -151,9 +154,10 @@ static int CeedOperatorBuildKernelData_Hip_gen(Ceed ceed, CeedInt num_input_fiel // Setup fields //------------------------------------------------------------------------------ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt i, CeedOperatorField op_field, - CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool is_at_points, bool use_3d_slices) { + CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool is_tensor, bool is_at_points, + bool use_3d_slices) { std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); - std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; + std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix, Q_name = is_tensor ? "Q_1d" : "Q"; std::string option_name = (is_input ? "inputs" : "outputs"); CeedEvalMode eval_mode = CEED_EVAL_NONE; CeedInt elem_size = 0, num_comp = 0, P_1d = 0; @@ -173,7 +177,8 @@ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, Ce CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); if (basis != CEED_BASIS_NONE) { CeedCallBackend(CeedBasisGetData(basis, &basis_data)); - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d)); } CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); @@ -210,7 +215,7 @@ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, Ce if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; else data->B.outputs[i] = basis_data->d_interp_1d; } - code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; + code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_name << "*" << Q_name << "];\n"; code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; break; case CEED_EVAL_GRAD: @@ -235,13 +240,15 @@ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, Ce if (is_input) data->B.inputs[i] = basis_data->d_interp_1d; else data->B.outputs[i] = basis_data->d_interp_1d; } - code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_1d * Q_1d << "];\n"; - code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; + if (is_tensor) { + code << " __shared__ CeedScalar s_B" << var_suffix << "[" << P_name << "*" << Q_name << "];\n"; + code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, B." << option_name << "[" << i << "], s_B" << var_suffix << ");\n"; + } if (is_at_points) break; // No G mat for AtPoints if (use_3d_slices) { if (is_input) data->G.inputs[i] = basis_data->d_collo_grad_1d; else data->G.outputs[i] = basis_data->d_collo_grad_1d; - code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * Q_1d << "];\n"; + code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_name << "*" << Q_name << "];\n"; code << " LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; } else { bool has_collo_grad = basis_data->d_collo_grad_1d; @@ -249,11 +256,12 @@ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, Ce if (is_input) data->G.inputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; else data->G.outputs[i] = has_collo_grad ? basis_data->d_collo_grad_1d : basis_data->d_grad_1d; if (has_collo_grad) { - code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * Q_1d << "];\n"; + code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_name << "*" << Q_name << "];\n"; code << " LoadMatrix<" << Q_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; } else { - code << " __shared__ CeedScalar s_G" << var_suffix << "[" << Q_1d * P_1d << "];\n"; - code << " LoadMatrix<" << P_name << ", " << Q_name << ">(data, G." << option_name << "[" << i << "], s_G" << var_suffix << ");\n"; + code << " __shared__ CeedScalar s_G" << var_suffix << "[" << P_name << "*" << Q_name << (is_tensor ? "" : "*dim") << "];\n"; + code << " LoadMatrix<" << P_name << ", " << Q_name << (is_tensor ? "" : "*dim") << ">(data, G." << option_name << "[" << i << "], s_G" + << var_suffix << ");\n"; } } break; @@ -274,9 +282,9 @@ static int CeedOperatorBuildKernelFieldData_Hip_gen(std::ostringstream &code, Ce //------------------------------------------------------------------------------ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt i, CeedInt dim, CeedInt field_input_buffer[], CeedOperatorField op_field, CeedQFunctionField qf_field, - CeedInt Q_1d, bool is_input, bool is_at_points, bool use_3d_slices) { + CeedInt Q_1d, bool is_input, bool is_tensor, bool is_at_points, bool use_3d_slices) { std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); - std::string P_name = "P_1d" + var_suffix; + std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix; CeedEvalMode eval_mode = CEED_EVAL_NONE; CeedInt elem_size = 0, num_comp = 0, P_1d = 0; CeedSize l_size; @@ -295,7 +303,8 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, } CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); if (basis != CEED_BASIS_NONE) { - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d)); } CeedCallBackend(CeedBasisDestroy(&basis)); CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); @@ -325,8 +334,8 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); code << " // CompStride: " << comp_stride << "\n"; data->indices.inputs[i] = (CeedInt *)rstr_data->d_offsets; - code << " ReadLVecStandard" << dim << "d(data, l_size" - << var_suffix << ", elem, indices.inputs[" << i << "], d" << var_suffix << ", r_e" << var_suffix << ");\n"; + code << " ReadLVecStandard" << (is_tensor ? dim : 1) << "d(data, l_size" << var_suffix << ", elem, indices.inputs[" << i << "], d" << var_suffix << ", r_e" << var_suffix << ");\n"; break; } case CEED_RESTRICTION_STRIDED: { @@ -341,8 +350,8 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); } code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; - code << " ReadLVecStrided" << dim << "d(data, elem, d" << var_suffix << ", r_e" << var_suffix << ");\n"; + code << " ReadLVecStrided" << (is_tensor ? dim : 1) << "d(data, elem, d" << var_suffix << ", r_e" << var_suffix << ");\n"; break; } case CEED_RESTRICTION_POINTS: { @@ -371,8 +380,8 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, CeedCallBackend(CeedElemRestrictionGetCompStride(elem_rstr, &comp_stride)); code << " // CompStride: " << comp_stride << "\n"; data->indices.outputs[i] = (CeedInt *)rstr_data->d_offsets; - code << " WriteLVecStandard" << dim << "d(data, l_size" - << var_suffix << ", elem, indices.outputs[" << i << "], r_e" << var_suffix << ", d" << var_suffix << ");\n"; + code << " WriteLVecStandard" << (is_tensor ? dim : 1) << "d(data, l_size" << var_suffix << ", elem, indices.outputs[" << i << "], r_e" << var_suffix << ", d" << var_suffix << ");\n"; break; } case CEED_RESTRICTION_STRIDED: { @@ -387,8 +396,8 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); } code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; - code << " WriteLVecStrided" << dim << "d(data, elem, r_e" << var_suffix << ", d" << var_suffix << ");\n"; + code << " WriteLVecStrided" << (is_tensor ? dim : 1) << "d(data, elem, r_e" << var_suffix << ", d" << var_suffix << ");\n"; break; } case CEED_RESTRICTION_POINTS: @@ -409,10 +418,10 @@ static int CeedOperatorBuildKernelRestriction_Hip_gen(std::ostringstream &code, // Basis //------------------------------------------------------------------------------ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt i, CeedInt dim, - CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, + CeedOperatorField op_field, CeedQFunctionField qf_field, CeedInt Q_1d, bool is_input, bool is_tensor, bool is_at_points, bool use_3d_slices) { std::string var_suffix = (is_input ? "_in_" : "_out_") + std::to_string(i); - std::string P_name = "P_1d" + var_suffix, Q_name = "Q_1d"; + std::string P_name = (is_tensor ? "P_1d" : "P") + var_suffix, Q_name = is_tensor ? "Q_1d" : "Q"; CeedEvalMode eval_mode = CEED_EVAL_NONE; CeedInt elem_size = 0, num_comp = 0, P_1d = 0; CeedElemRestriction elem_rstr; @@ -427,7 +436,8 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); CeedCallBackend(CeedOperatorFieldGetBasis(op_field, &basis)); if (basis != CEED_BASIS_NONE) { - CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + if (is_tensor) CeedCallBackend(CeedBasisGetNumNodes1D(basis, &P_1d)); + else CeedCallBackend(CeedBasisGetNumNodes(basis, &P_1d)); } CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_field, &eval_mode)); @@ -442,29 +452,45 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp break; case CEED_EVAL_INTERP: if (is_at_points) { - code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; - code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d"; + + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", r_c" << var_suffix << ");\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; - code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; + std::string function_name = is_tensor ? ((dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d") : "InterpNonTensor"; + + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", r_q" << var_suffix << ");\n"; } break; case CEED_EVAL_GRAD: if (is_at_points) { - code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; - code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_c" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d"; + + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", r_c" << var_suffix << ");\n"; } else if (use_3d_slices) { + std::string function_name = (dim > 1 ? "InterpTensor" : "Interp") + std::to_string(dim) + "d"; + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; - code << " Interp" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", r_q" << var_suffix << ");\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", r_q" << var_suffix << ");\n"; + } else if (is_tensor) { + bool is_collocated = dim == 3 && Q_1d >= P_1d; + std::string function_name = (dim == 1 ? "Grad" : (is_collocated ? "GradTensorCollocated" : "GradTensor")) + std::to_string(dim) + "d"; + + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << (dim >= 3 ? Q_name : "1") << "];\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; - code << " Grad" << (dim > 1 ? "Tensor" : "") << (dim == 3 && Q_1d >= P_1d ? "Collocated" : "") << dim << "d(data, r_e" << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix << ", r_q" - << var_suffix << ");\n"; + std::string function_name = "GradNonTensor"; + + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim];\n"; + code << " " << function_name << "(data, r_e" << var_suffix + << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n"; } break; case CEED_EVAL_WEIGHT: { @@ -472,11 +498,12 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp code << " // Nothing to do AtPoints\n"; } else { CeedBasis_Hip_shared *basis_data; + std::string function_name = is_tensor ? ((dim == 1 ? "Weight" : "WeightTensor") + std::to_string(dim) + "d") : "WeightNonTensor"; - code << " CeedScalar r_q" << var_suffix << "[" << Q_name << "];\n"; + code << " CeedScalar r_q" << var_suffix << "[" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; CeedCallBackend(CeedBasisGetData(basis, &basis_data)); data->W = basis_data->d_q_weight_1d; - code << " Weight" << (dim > 1 ? "Tensor" : "") << dim << "d<" << Q_name << ">(data, W, r_q" << var_suffix << ");\n"; + code << " " << function_name << "<" << Q_name << ">(data, W, r_q" << var_suffix << ");\n"; } break; } @@ -494,25 +521,42 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp case CEED_EVAL_INTERP: code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; if (is_at_points) { - code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_c" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d"; + + code << " " << function_name << "(data, r_c" << var_suffix << ", s_B" + << var_suffix << ", r_e" << var_suffix << ");\n"; } else { - code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + std::string function_name = + is_tensor ? ((dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d") : "InterpTransposeNonTensor"; + + code << " " << function_name << "(data, r_q" << var_suffix << ", s_B" + << var_suffix << ", r_e" << var_suffix << ");\n"; } break; case CEED_EVAL_GRAD: code << " CeedScalar *r_e" << var_suffix << " = r_e_scratch;\n"; if (is_at_points) { - code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_c" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d"; + + code << " " << function_name << "(data, r_c" << var_suffix << ", s_B" + << var_suffix << ", r_e" << var_suffix << ");\n"; } else if (use_3d_slices) { - code << " InterpTranspose" << (dim > 1 ? "Tensor" : "") << dim << "d(data, r_q" << var_suffix << ", s_B" << var_suffix << ", r_e" << var_suffix << ");\n"; + std::string function_name = (dim == 1 ? "InterpTranspose" : "InterpTransposeTensor") + std::to_string(dim) + "d"; + + code << " " << function_name << "(data, r_q" << var_suffix << ", s_B" + << var_suffix << ", r_e" << var_suffix << ");\n"; + } else if (is_tensor) { + bool is_collocated = dim == 3 && Q_1d >= P_1d; + std::string function_name = + (dim == 1 ? "GradTranspose" : (is_collocated ? "GradTransposeTensorCollocated" : "GradTransposeTensor")) + std::to_string(dim) + "d"; + + code << " " << function_name << "(data, r_q" << var_suffix << ", s_B" + << var_suffix << ", s_G" << var_suffix << ", r_e" << var_suffix << ");\n"; } else { - code << " GradTranspose" << (dim > 1 ? "Tensor" : "") << (dim == 3 && Q_1d >= P_1d ? "Collocated" : "") << dim << "d(data, r_q" << var_suffix << ", s_B" << var_suffix << ", s_G" << var_suffix - << ", r_e" << var_suffix << ");\n"; + std::string function_name = "GradTransposeNonTensor"; + + code << " " << function_name << "(data, r_q" << var_suffix + << ", s_G" << var_suffix << ", r_e" << var_suffix << ");\n"; } break; // LCOV_EXCL_START @@ -534,9 +578,9 @@ static int CeedOperatorBuildKernelBasis_Hip_gen(std::ostringstream &code, CeedOp static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, CeedOperator_Hip_gen *data, CeedInt dim, CeedInt max_num_points, CeedInt num_input_fields, CeedOperatorField *op_input_fields, CeedQFunctionField *qf_input_fields, CeedInt num_output_fields, CeedOperatorField *op_output_fields, - CeedQFunctionField *qf_output_fields, std::string qfunction_name, CeedInt Q_1d, bool is_at_points, - bool use_3d_slices) { - std::string Q_name = "Q_1d"; + CeedQFunctionField *qf_output_fields, std::string qfunction_name, CeedInt Q_1d, bool is_tensor, + bool is_at_points, bool use_3d_slices) { + std::string Q_name = is_tensor ? "Q_1d" : "Q"; CeedEvalMode eval_mode = CEED_EVAL_NONE; CeedElemRestriction elem_rstr; @@ -552,25 +596,25 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce if (is_at_points) { code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "];\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; } break; case CEED_EVAL_INTERP: if (is_at_points) { // Accumulator for point data - code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "];\n"; - code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "; i++) {\n"; + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "];\n"; + code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "; i++) {\n"; code << " r_c" << var_suffix << "[i] = 0.0;\n"; code << " }\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << Q_name << "];\n"; + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; } break; case CEED_EVAL_GRAD: if (is_at_points) { // Accumulator for point data - code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "*dim];\n"; - code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim > 2 ? Q_name : "1") << "; i++) {\n"; + code << " CeedScalar r_c" << var_suffix << "[num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "*dim];\n"; + code << " for (CeedInt i = 0; i < num_comp" << var_suffix << "*" << (dim >= 3 ? Q_name : "1") << "; i++) {\n"; code << " r_c" << var_suffix << "[i] = 0.0;\n"; code << " }\n"; } else if (use_3d_slices) { @@ -580,7 +624,7 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce code << " r_q" << var_suffix << "[i] = 0.0;\n"; code << " }\n"; } else { - code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << Q_name << "];\n"; + code << " CeedScalar r_q" << var_suffix << "[num_comp" << var_suffix << "*dim*" << (is_tensor && (dim >= 3) ? Q_name : "1") << "];\n"; } break; case CEED_EVAL_WEIGHT: @@ -701,7 +745,7 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce CeedCallBackend(CeedElemRestrictionGetStrides(elem_rstr, strides)); } code << " // Strides: {" << strides[0] << ", " << strides[1] << ", " << strides[2] << "}\n"; - code << " ReadEVecSliceStrided3d(data, elem, q, d" << var_suffix << ", r_s" << var_suffix << ");\n"; } else { CeedSize l_size = 0; @@ -717,6 +761,7 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce code << " ReadEVecSliceStandard3d(data, l_size" << var_suffix << ", elem, q, indices.inputs[" << i << "], d" << var_suffix << ", r_s" << var_suffix << ");\n"; } + CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); break; case CEED_EVAL_INTERP: code << " CeedScalar r_s" << var_suffix << "[num_comp" << var_suffix << "];\n"; @@ -799,10 +844,10 @@ static int CeedOperatorBuildKernelQFunction_Hip_gen(std::ostringstream &code, Ce // Apply QFunction code << "\n // -- Apply QFunction\n"; code << " " << qfunction_name << "(ctx, "; - if (dim != 3 || is_at_points || use_3d_slices) { + if (dim != 3 || is_at_points || use_3d_slices || !is_tensor) { code << "1"; } else { - code << "Q_1d"; + code << Q_name; } code << ", inputs, outputs);\n"; @@ -962,11 +1007,17 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { } // Load basis source files - // TODO: Add non-tensor, AtPoints - code << "// Tensor basis source\n"; - code << "#include \n\n"; - code << "// AtPoints basis source\n"; - code << "#include \n\n"; + if (is_tensor) { + code << "// Tensor basis source\n"; + code << "#include \n\n"; + } else { + code << "// Non-tensor basis source\n"; + code << "#include \n\n"; + } + if (is_at_points) { + code << "// AtPoints basis source\n"; + code << "#include \n\n"; + } code << "// CodeGen operator source\n"; code << "#include \n\n"; @@ -978,7 +1029,7 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { // Define CEED_Q_VLA code << "\n#undef CEED_Q_VLA\n"; - if (dim != 3 || is_at_points || use_3d_slices) { + if (dim != 3 || is_at_points || use_3d_slices || !is_tensor) { code << "#define CEED_Q_VLA 1\n\n"; } else { code << "#define CEED_Q_VLA " << Q_1d << "\n\n"; @@ -1002,7 +1053,7 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { code << "// d_[in,out]_i: CeedVector device array\n"; code << "// r_[in,out]_e_i: Element vector register\n"; code << "// r_[in,out]_q_i: Quadrature space vector register\n"; - code << "// r_[in,out]_c_i: AtPoints Chebyshev coefficents register\n"; + code << "// r_[in,out]_c_i: AtPoints Chebyshev coefficients register\n"; code << "// r_[in,out]_s_i: Quadrature space slice vector register\n"; code << "// \n"; code << "// s_B_[in,out]_i: Interpolation matrix, shared memory\n"; @@ -1026,7 +1077,7 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { } code << " const CeedInt dim = " << dim << ";\n"; - code << " const CeedInt Q_1d = " << Q_1d << ";\n"; + code << " const CeedInt " << (is_tensor ? "Q_1d" : "Q") << " = " << Q_1d << ";\n"; if (is_at_points) { code << " const CeedInt max_num_points = " << max_num_points << ";\n"; code << " const CeedInt coords_comp_stride = " << coords_comp_stride << ";\n"; @@ -1039,18 +1090,18 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { code << " data.t_id_y = threadIdx.y;\n"; code << " data.t_id_z = threadIdx.z;\n"; code << " data.t_id = threadIdx.x + threadIdx.y*blockDim.x + threadIdx.z*blockDim.y*blockDim.x;\n"; - code << " data.slice = slice + data.t_id_z*T_1D" << (dim > 1 ? "*T_1D" : "") << ";\n"; + code << " data.slice = slice + data.t_id_z*T_1D" << ((!is_tensor || dim == 1) ? "" : "*T_1D") << ";\n"; // Initialize constants, and matrices B and G code << "\n // Input field constants and basis data\n"; for (CeedInt i = 0; i < num_input_fields; i++) { - CeedCallBackend( - CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_input_fields[i], qf_input_fields[i], Q_1d, true, is_at_points, use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_input_fields[i], qf_input_fields[i], Q_1d, true, is_tensor, + is_at_points, use_3d_slices)); } code << "\n // Output field constants and basis data\n"; for (CeedInt i = 0; i < num_output_fields; i++) { - CeedCallBackend( - CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_at_points, use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelFieldData_Hip_gen(code, data, i, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_tensor, + is_at_points, use_3d_slices)); } // Loop over all elements @@ -1068,7 +1119,7 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_input_fields[i], &elem_rstr)); CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); - max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * elem_size); + max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_tensor && (dim >= 3) ? elem_size : 1)); CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); } for (CeedInt i = 0; i < num_output_fields; i++) { @@ -1078,7 +1129,7 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { CeedCallBackend(CeedOperatorFieldGetElemRestriction(op_output_fields[i], &elem_rstr)); CeedCallBackend(CeedElemRestrictionGetNumComponents(elem_rstr, &num_comp)); CeedCallBackend(CeedElemRestrictionGetElementSize(elem_rstr, &elem_size)); - max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * elem_size); + max_rstr_buffer_size = CeedIntMax(max_rstr_buffer_size, num_comp * (is_tensor && (dim >= 3) ? elem_size : 1)); CeedCallBackend(CeedElemRestrictionDestroy(&elem_rstr)); } code << " // Scratch restriction buffer space\n"; @@ -1137,17 +1188,17 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { // ---- Restriction CeedCallBackend(CeedOperatorBuildKernelRestriction_Hip_gen(code, data, f, dim, field_rstr_in_buffer, op_input_fields[f], qf_input_fields[f], Q_1d, - true, is_at_points, use_3d_slices)); + true, is_tensor, is_at_points, use_3d_slices)); // ---- Basis action - CeedCallBackend( - CeedOperatorBuildKernelBasis_Hip_gen(code, data, f, dim, op_input_fields[f], qf_input_fields[f], Q_1d, true, is_at_points, use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelBasis_Hip_gen(code, data, f, dim, op_input_fields[f], qf_input_fields[f], Q_1d, true, is_tensor, + is_at_points, use_3d_slices)); } // -- Q function CeedCallBackend(CeedOperatorBuildKernelQFunction_Hip_gen(code, data, dim, max_num_points, num_input_fields, op_input_fields, qf_input_fields, - num_output_fields, op_output_fields, qf_output_fields, qfunction_name, Q_1d, is_at_points, - use_3d_slices)); + num_output_fields, op_output_fields, qf_output_fields, qfunction_name, Q_1d, is_tensor, + is_at_points, use_3d_slices)); // -- Output basis and restriction code << "\n // -- Output field basis action and restrictions\n"; @@ -1155,12 +1206,12 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { code << " // ---- Output field " << i << "\n"; // ---- Basis action - CeedCallBackend( - CeedOperatorBuildKernelBasis_Hip_gen(code, data, i, dim, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_at_points, use_3d_slices)); + CeedCallBackend(CeedOperatorBuildKernelBasis_Hip_gen(code, data, i, dim, op_output_fields[i], qf_output_fields[i], Q_1d, false, is_tensor, + is_at_points, use_3d_slices)); // ---- Restriction CeedCallBackend(CeedOperatorBuildKernelRestriction_Hip_gen(code, data, i, dim, NULL, op_output_fields[i], qf_output_fields[i], Q_1d, false, - is_at_points, use_3d_slices)); + is_tensor, is_at_points, use_3d_slices)); } // Close loop and function @@ -1173,7 +1224,7 @@ extern "C" int CeedOperatorBuildKernel_Hip_gen(CeedOperator op) { // Compile CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); - CeedCallBackend(BlockGridCalculate_Hip_gen(dim, num_elem, data->max_P_1d, Q_1d, block_sizes)); + CeedCallBackend(BlockGridCalculate_Hip_gen(is_tensor ? dim : 1, num_elem, data->max_P_1d, Q_1d, block_sizes)); CeedCallBackend(CeedCompile_Hip(ceed, code.str().c_str(), &data->module, 2, "T_1D", block_sizes[0], "BLOCK_SIZE", block_sizes[0] * block_sizes[1] * block_sizes[2])); CeedCallBackend(CeedGetKernel_Hip(ceed, data->module, operator_name.c_str(), &data->op)); diff --git a/backends/hip-gen/ceed-hip-gen-operator.c b/backends/hip-gen/ceed-hip-gen-operator.c index d3455ac137..a2a6ccd1f1 100644 --- a/backends/hip-gen/ceed-hip-gen-operator.c +++ b/backends/hip-gen/ceed-hip-gen-operator.c @@ -35,7 +35,7 @@ static int CeedOperatorDestroy_Hip_gen(CeedOperator op) { // Apply and add to output //------------------------------------------------------------------------------ static int CeedOperatorApplyAdd_Hip_gen(CeedOperator op, CeedVector input_vec, CeedVector output_vec, CeedRequest *request) { - bool is_at_points; + bool is_at_points, is_tensor; Ceed ceed; CeedInt num_elem, num_input_fields, num_output_fields; CeedEvalMode eval_mode; @@ -46,16 +46,62 @@ static int CeedOperatorApplyAdd_Hip_gen(CeedOperator op, CeedVector input_vec, C CeedOperatorField *op_input_fields, *op_output_fields; CeedOperator_Hip_gen *data; - // Check for tensor-product bases + // Check for shared bases + CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); { - bool has_tensor_bases; + bool has_shared_bases = true, is_all_tensor = true, is_all_nontensor = true; + + for (CeedInt i = 0; i < num_input_fields; i++) { + CeedBasis basis; + + CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); + if (basis != CEED_BASIS_NONE) { + bool is_tensor = true; + const char *resource; + char *resource_root; + Ceed basis_ceed; + + CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); + is_all_tensor &= is_tensor; + is_all_nontensor &= !is_tensor; + CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed)); + CeedCallBackend(CeedGetResource(basis_ceed, &resource)); + CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root)); + has_shared_bases &= !strcmp(resource_root, "/gpu/hip/shared"); + CeedCallBackend(CeedFree(&resource_root)); + CeedCallBackend(CeedDestroy(&basis_ceed)); + } + CeedCallBackend(CeedBasisDestroy(&basis)); + } - CeedCallBackend(CeedOperatorHasTensorBases(op, &has_tensor_bases)); - // -- Fallback to ref if not all bases are tensor-product - if (!has_tensor_bases) { + for (CeedInt i = 0; i < num_output_fields; i++) { + CeedBasis basis; + + CeedCallBackend(CeedOperatorFieldGetBasis(op_output_fields[i], &basis)); + if (basis != CEED_BASIS_NONE) { + bool is_tensor = true; + const char *resource; + char *resource_root; + Ceed basis_ceed; + + CeedCallBackend(CeedBasisIsTensor(basis, &is_tensor)); + is_all_tensor &= is_tensor; + is_all_nontensor &= !is_tensor; + + CeedCallBackend(CeedBasisGetCeed(basis, &basis_ceed)); + CeedCallBackend(CeedGetResource(basis_ceed, &resource)); + CeedCallBackend(CeedGetResourceRoot(basis_ceed, resource, ":", &resource_root)); + has_shared_bases &= !strcmp(resource_root, "/gpu/hip/shared"); + CeedCallBackend(CeedFree(&resource_root)); + CeedCallBackend(CeedDestroy(&basis_ceed)); + } + CeedCallBackend(CeedBasisDestroy(&basis)); + } + // -- Fallback to ref if not all bases are shared + if (!has_shared_bases || (!is_all_tensor && !is_all_nontensor)) { CeedOperator op_fallback; - CeedDebug256(CeedOperatorReturnCeed(op), CEED_DEBUG_COLOR_SUCCESS, "Falling back to /gpu/hip/ref CeedOperator due to non-tensor bases"); + CeedDebug256(CeedOperatorReturnCeed(op), CEED_DEBUG_COLOR_SUCCESS, "Falling back to /gpu/hip/ref CeedOperator due to unsupported bases"); CeedCallBackend(CeedOperatorGetFallback(op, &op_fallback)); CeedCallBackend(CeedOperatorApplyAdd(op_fallback, input_vec, output_vec, request)); return CEED_ERROR_SUCCESS; @@ -67,7 +113,6 @@ static int CeedOperatorApplyAdd_Hip_gen(CeedOperator op, CeedVector input_vec, C CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); CeedCallBackend(CeedQFunctionGetData(qf, &qf_data)); CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); - CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); // Creation of the operator @@ -158,10 +203,19 @@ static int CeedOperatorApplyAdd_Hip_gen(CeedOperator op, CeedVector input_vec, C const CeedInt Q_1d = data->Q_1d; const CeedInt P_1d = data->max_P_1d; const CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); - CeedInt block_sizes[3]; - CeedCallBackend(BlockGridCalculate_Hip_gen(dim, num_elem, P_1d, Q_1d, block_sizes)); - if (dim == 1) { + CeedCallBackend(CeedOperatorHasTensorBases(op, &is_tensor)); + CeedInt block_sizes[3] = {thread_1d, ((!is_tensor || dim == 1) ? 1 : thread_1d), -1}; + + if (is_tensor) { + CeedCallBackend(BlockGridCalculate_Hip_gen(is_tensor ? dim : 1, num_elem, P_1d, Q_1d, block_sizes)); + } else { + CeedInt elems_per_block = 64 * thread_1d > 256 ? 256 / thread_1d : 64; + + elems_per_block = elems_per_block > 0 ? elems_per_block : 1; + block_sizes[2] = elems_per_block; + } + if (dim == 1 || !is_tensor) { CeedInt grid = num_elem / block_sizes[2] + ((num_elem / block_sizes[2] * block_sizes[2] < num_elem) ? 1 : 0); CeedInt sharedMem = block_sizes[2] * thread_1d * sizeof(CeedScalar); diff --git a/backends/hip-ref/ceed-hip-ref-operator.c b/backends/hip-ref/ceed-hip-ref-operator.c index 6f1119084b..67a3533ee2 100644 --- a/backends/hip-ref/ceed-hip-ref-operator.c +++ b/backends/hip-ref/ceed-hip-ref-operator.c @@ -131,7 +131,7 @@ static int CeedOperatorSetupFields_Hip(CeedQFunction qf, CeedOperator op, bool i // Input CEED_VECTOR_ACTIVE // Output CEED_VECTOR_ACTIVE without CEED_EVAL_NONE // Input CEED_VECTOR_NONE with CEED_EVAL_WEIGHT - // Input passive vectorr with CEED_EVAL_NONE and strided restriction with CEED_STRIDES_BACKEND + // Input passive vector with CEED_EVAL_NONE and strided restriction with CEED_STRIDES_BACKEND CeedCallBackend(CeedOperatorFieldGetVector(op_fields[i], &l_vec)); is_active = l_vec == CEED_VECTOR_ACTIVE; CeedCallBackend(CeedVectorDestroy(&l_vec)); diff --git a/doc/sphinx/source/releasenotes.md b/doc/sphinx/source/releasenotes.md index 816187fd1f..a9f765a2fe 100644 --- a/doc/sphinx/source/releasenotes.md +++ b/doc/sphinx/source/releasenotes.md @@ -24,6 +24,7 @@ On this page we provide a summary of the main API changes, new features and exam - Enable `#pragma once` for all JiT source; remove duplicate includes in JiT source string before compilation. - Allow user to set additional compiler options for CUDA and HIP JiT. Specifically, directories set with `CeedAddJitSourceRoot(ceed, "foo/bar")` will be used to set `-Ifoo/bar` and defines set with `CeedAddJitDefine(ceed, "foo=bar")` will be used to set `-Dfoo=bar`. +- Added non-tensor basis support to code generation backends `/gpu/cuda/gen` and `/gpu/hip/gen`. ### Examples