diff --git a/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp b/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp index 99da342f66..7ce356bd8a 100644 --- a/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp +++ b/backends/cuda-gen/ceed-cuda-gen-operator-build.cpp @@ -29,11 +29,14 @@ static int CeedOperatorBuildKernelData_Cuda_gen(Ceed ceed, CeedInt num_input_fie CeedQFunctionField *qf_input_fields, CeedInt num_output_fields, CeedOperatorField *op_output_fields, CeedQFunctionField *qf_output_fields, CeedInt *max_P_1d, CeedInt *Q_1d, CeedInt *dim, bool *is_tensor, bool *use_3d_slices) { + bool valid_tensor_status = false; + // Find dim, P_1d, Q_1d *max_P_1d = 0; *Q_1d = 0; *dim = 0; *is_tensor = true; + for (CeedInt i = 0; i < num_input_fields; i++) { CeedBasis basis; @@ -44,14 +47,17 @@ 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)); + CeedCheck(!valid_tensor_status || *is_tensor == is_field_tensor, ceed, CEED_ERROR_BACKEND, "Fields must be compatible"); + valid_tensor_status = true; + *is_tensor = *is_tensor && is_field_tensor; + 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 +73,17 @@ 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)); + CeedCheck(!valid_tensor_status || *is_tensor == is_field_tensor, ceed, CEED_ERROR_BACKEND, "Fields must be compatible"); + valid_tensor_status = true; + *is_tensor = *is_tensor && is_field_tensor; + 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 +133,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 +156,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 +194,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 +219,13 @@ 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"; 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 +233,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 +259,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 +280,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 +311,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 +327,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 +357,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 +373,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 +395,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 +413,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 +429,45 @@ static int CeedOperatorBuildKernelBasis_Cuda_gen(std::ostringstream &code, CeedO break; case CEED_EVAL_INTERP: if (is_at_points) { + std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d"; + 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"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", r_c" << var_suffix << ");\n"; } else { + 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 << "*" << 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"; } break; case CEED_EVAL_GRAD: if (is_at_points) { + std::string function_name = (dim == 1 ? "Interp" : "InterpTensor") + std::to_string(dim) + "d"; + 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"; + 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*" << Q_name << "];\n"; + code << " " << function_name << "(data, r_e" << var_suffix << ", s_B" + << var_suffix << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n"; } else { + std::string function_name = "GradNonTensor"; + 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"; + code << " " << function_name << "(data, r_e" << var_suffix + << ", s_G" << var_suffix << ", r_q" << var_suffix << ");\n"; } break; case CEED_EVAL_WEIGHT: { @@ -446,11 +475,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"; 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 +498,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 +556,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; @@ -675,7 +723,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; @@ -777,7 +825,7 @@ static int CeedOperatorBuildKernelQFunction_Cuda_gen(std::ostringstream &code, C if (dim != 3 || is_at_points || use_3d_slices) { code << "1"; } else { - code << "Q_1d"; + code << Q_name; } code << ", inputs, outputs);\n"; @@ -950,11 +998,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"; @@ -1014,7 +1068,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"; @@ -1032,13 +1086,13 @@ extern "C" int CeedOperatorBuildKernel_Cuda_gen(CeedOperator op) { // 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 @@ -1125,17 +1179,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 +1197,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..cffe9c3893 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,53 @@ 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; + + for (CeedInt i = 0; i < num_input_fields; i++) { + CeedBasis basis; + + CeedCallBackend(CeedOperatorFieldGetBasis(op_input_fields[i], &basis)); + if (basis != CEED_BASIS_NONE) { + const char *resource; + char *resource_root; + Ceed basis_ceed; + + 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) { + const char *resource; + char *resource_root; + Ceed basis_ceed; + + 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) { 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 to large non-tensor bases"); CeedCallBackend(CeedOperatorGetFallback(op, &op_fallback)); CeedCallBackend(CeedOperatorApplyAdd(op_fallback, input_vec, output_vec, request)); return CEED_ERROR_SUCCESS; @@ -132,7 +170,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,8 +269,9 @@ 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));