1 // Copyright (c) 2017-2025, Lawrence Livermore National Security, LLC and other CEED contributors. 2 // All Rights Reserved. See the top-level LICENSE and NOTICE files for details. 3 // 4 // SPDX-License-Identifier: BSD-2-Clause 5 // 6 // This file is part of CEED: http://github.com/ceed 7 8 #include <ceed.h> 9 #include <ceed/backend.h> 10 #include <ceed/jit-source/cuda/cuda-types.h> 11 #include <cuda.h> 12 #include <cuda_runtime.h> 13 #include <stddef.h> 14 #include <string.h> 15 16 #include "../cuda/ceed-cuda-common.h" 17 #include "../cuda/ceed-cuda-compile.h" 18 #include "ceed-cuda-gen-operator-build.h" 19 #include "ceed-cuda-gen.h" 20 21 //------------------------------------------------------------------------------ 22 // Destroy operator 23 //------------------------------------------------------------------------------ 24 static int CeedOperatorDestroy_Cuda_gen(CeedOperator op) { 25 Ceed ceed; 26 CeedOperator_Cuda_gen *impl; 27 28 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 29 CeedCallBackend(CeedOperatorGetData(op, &impl)); 30 if (impl->points.num_per_elem) CeedCallCuda(ceed, cudaFree((void **)impl->points.num_per_elem)); 31 CeedCallBackend(CeedFree(&impl)); 32 CeedCallBackend(CeedDestroy(&ceed)); 33 return CEED_ERROR_SUCCESS; 34 } 35 36 static int Waste(int threads_per_sm, int warp_size, int threads_per_elem, int elems_per_block) { 37 int useful_threads_per_block = threads_per_elem * elems_per_block; 38 // round up to nearest multiple of warp_size 39 int block_size = CeedDivUpInt(useful_threads_per_block, warp_size) * warp_size; 40 int blocks_per_sm = threads_per_sm / block_size; 41 return threads_per_sm - useful_threads_per_block * blocks_per_sm; 42 } 43 44 // Choose the least wasteful block size constrained by blocks_per_sm of max_threads_per_block. 45 // 46 // The x and y part of block[] contains per-element sizes (specified on input) while the z part is number of elements. 47 // 48 // Problem setting: we'd like to make occupancy high with relatively few inactive threads. CUDA (cuOccupancyMaxPotentialBlockSize) can tell us how 49 // many threads can run. 50 // 51 // Note that full occupancy sometimes can't be achieved by one thread block. 52 // For example, an SM might support 1536 threads in total, but only 1024 within a single thread block. 53 // So cuOccupancyMaxPotentialBlockSize may suggest a block size of 768 so that two blocks can run, versus one block of 1024 will prevent a second 54 // block from running. The cuda-gen kernels are pretty heavy with lots of instruction-level parallelism (ILP) so we'll generally be okay with 55 // relatively low occupancy and smaller thread blocks, but we solve a reasonably general problem here. Empirically, we find that blocks bigger than 56 // about 256 have higher latency and worse load balancing when the number of elements is modest. 57 // 58 // cuda-gen can't choose block sizes arbitrarily; they need to be a multiple of the number of quadrature points (or number of basis functions). 59 // They also have a lot of __syncthreads(), which is another point against excessively large thread blocks. 60 // Suppose I have elements with 7x7x7 quadrature points. 61 // This will loop over the last dimension, so we have 7*7=49 threads per element. 62 // Suppose we have two elements = 2*49=98 useful threads. 63 // CUDA schedules in units of full warps (32 threads), so 128 CUDA hardware threads are effectively committed to that block. 64 // Now suppose cuOccupancyMaxPotentialBlockSize returned 352. 65 // We can schedule 2 blocks of size 98 (196 useful threads using 256 hardware threads), but not a third block (which would need a total of 384 66 // hardware threads). 67 // 68 // If instead, we had packed 3 elements, we'd have 3*49=147 useful threads occupying 160 slots, and could schedule two blocks. 69 // Alternatively, we could pack a single block of 7 elements (2*49=343 useful threads) into the 354 slots. 70 // The latter has the least "waste", but __syncthreads() over-synchronizes and it might not pay off relative to smaller blocks. 71 static int BlockGridCalculate(CeedInt num_elem, int blocks_per_sm, int max_threads_per_block, int max_threads_z, int warp_size, int block[3], 72 int *grid) { 73 const int threads_per_sm = blocks_per_sm * max_threads_per_block; 74 const int threads_per_elem = block[0] * block[1]; 75 int elems_per_block = 1; 76 int waste = Waste(threads_per_sm, warp_size, threads_per_elem, 1); 77 78 for (int i = 2; i <= CeedIntMin(max_threads_per_block / threads_per_elem, num_elem); i++) { 79 int i_waste = Waste(threads_per_sm, warp_size, threads_per_elem, i); 80 81 // We want to minimize waste, but smaller kernels have lower latency and less __syncthreads() overhead so when a larger block size has the same 82 // waste as a smaller one, go ahead and prefer the smaller block. 83 if (i_waste < waste || (i_waste == waste && threads_per_elem * i <= 128)) { 84 elems_per_block = i; 85 waste = i_waste; 86 } 87 } 88 // In low-order elements, threads_per_elem may be sufficiently low to give an elems_per_block greater than allowable for the device, so we must 89 // check before setting the z-dimension size of the block. 90 block[2] = CeedIntMin(elems_per_block, max_threads_z); 91 *grid = CeedDivUpInt(num_elem, elems_per_block); 92 return CEED_ERROR_SUCCESS; 93 } 94 95 // callback for cuOccupancyMaxPotentialBlockSize, providing the amount of dynamic shared memory required for a thread block of size threads. 96 static size_t dynamicSMemSize(int threads) { return threads * sizeof(CeedScalar); } 97 98 //------------------------------------------------------------------------------ 99 // Apply and add to output 100 //------------------------------------------------------------------------------ 101 static int CeedOperatorApplyAddCore_Cuda_gen(CeedOperator op, CUstream stream, const CeedScalar *input_arr, CeedScalar *output_arr, bool *is_run_good, 102 CeedRequest *request) { 103 bool is_at_points, is_tensor; 104 Ceed ceed; 105 Ceed_Cuda *cuda_data; 106 CeedInt num_elem, num_input_fields, num_output_fields; 107 CeedEvalMode eval_mode; 108 CeedQFunctionField *qf_input_fields, *qf_output_fields; 109 CeedQFunction_Cuda_gen *qf_data; 110 CeedQFunction qf; 111 CeedOperatorField *op_input_fields, *op_output_fields; 112 CeedOperator_Cuda_gen *data; 113 114 // Build the operator kernel 115 CeedCallBackend(CeedOperatorBuildKernel_Cuda_gen(op, is_run_good)); 116 if (!(*is_run_good)) return CEED_ERROR_SUCCESS; 117 118 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 119 CeedCallBackend(CeedGetData(ceed, &cuda_data)); 120 CeedCallBackend(CeedOperatorGetData(op, &data)); 121 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 122 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data)); 123 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 124 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 125 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 126 127 // Input vectors 128 for (CeedInt i = 0; i < num_input_fields; i++) { 129 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 130 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 131 data->fields.inputs[i] = NULL; 132 } else { 133 bool is_active; 134 CeedVector vec; 135 136 // Get input vector 137 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 138 is_active = vec == CEED_VECTOR_ACTIVE; 139 if (is_active) data->fields.inputs[i] = input_arr; 140 else CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->fields.inputs[i])); 141 CeedCallBackend(CeedVectorDestroy(&vec)); 142 } 143 } 144 145 // Output vectors 146 for (CeedInt i = 0; i < num_output_fields; i++) { 147 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 148 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 149 data->fields.outputs[i] = NULL; 150 } else { 151 bool is_active; 152 CeedVector vec; 153 154 // Get output vector 155 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 156 is_active = vec == CEED_VECTOR_ACTIVE; 157 if (is_active) data->fields.outputs[i] = output_arr; 158 else CeedCallBackend(CeedVectorGetArray(vec, CEED_MEM_DEVICE, &data->fields.outputs[i])); 159 CeedCallBackend(CeedVectorDestroy(&vec)); 160 } 161 } 162 163 // Point coordinates, if needed 164 CeedCallBackend(CeedOperatorIsAtPoints(op, &is_at_points)); 165 if (is_at_points) { 166 // Coords 167 CeedVector vec; 168 169 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, NULL, &vec)); 170 CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->points.coords)); 171 CeedCallBackend(CeedVectorDestroy(&vec)); 172 173 // Points per elem 174 if (num_elem != data->points.num_elem) { 175 CeedInt *points_per_elem; 176 const CeedInt num_bytes = num_elem * sizeof(CeedInt); 177 CeedElemRestriction rstr_points = NULL; 178 179 data->points.num_elem = num_elem; 180 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, &rstr_points, NULL)); 181 CeedCallBackend(CeedCalloc(num_elem, &points_per_elem)); 182 for (CeedInt e = 0; e < num_elem; e++) { 183 CeedInt num_points_elem; 184 185 CeedCallBackend(CeedElemRestrictionGetNumPointsInElement(rstr_points, e, &num_points_elem)); 186 points_per_elem[e] = num_points_elem; 187 } 188 if (data->points.num_per_elem) CeedCallCuda(ceed, cudaFree((void **)data->points.num_per_elem)); 189 CeedCallCuda(ceed, cudaMalloc((void **)&data->points.num_per_elem, num_bytes)); 190 CeedCallCuda(ceed, cudaMemcpy((void *)data->points.num_per_elem, points_per_elem, num_bytes, cudaMemcpyHostToDevice)); 191 CeedCallBackend(CeedElemRestrictionDestroy(&rstr_points)); 192 CeedCallBackend(CeedFree(&points_per_elem)); 193 } 194 } 195 196 // Get context data 197 CeedCallBackend(CeedQFunctionGetInnerContextData(qf, CEED_MEM_DEVICE, &qf_data->d_c)); 198 199 // Apply operator 200 void *opargs[] = {(void *)&num_elem, &qf_data->d_c, &data->indices, &data->fields, &data->B, &data->G, &data->W, &data->points}; 201 const CeedInt dim = data->dim; 202 const CeedInt Q_1d = data->Q_1d; 203 const CeedInt P_1d = data->max_P_1d; 204 const CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); 205 int max_threads_per_block, min_grid_size, grid; 206 207 CeedCallBackend(CeedOperatorHasTensorBases(op, &is_tensor)); 208 CeedCallCuda(ceed, cuOccupancyMaxPotentialBlockSize(&min_grid_size, &max_threads_per_block, data->op, dynamicSMemSize, 0, 0x10000)); 209 int block[3] = {thread_1d, ((!is_tensor || dim == 1) ? 1 : thread_1d), -1}; 210 211 if (is_tensor) { 212 CeedCallBackend(BlockGridCalculate(num_elem, min_grid_size / cuda_data->device_prop.multiProcessorCount, is_at_points ? 1 : max_threads_per_block, 213 cuda_data->device_prop.maxThreadsDim[2], cuda_data->device_prop.warpSize, block, &grid)); 214 } else { 215 CeedInt elems_per_block = CeedIntMin(cuda_data->device_prop.maxThreadsDim[2], CeedIntMax(512 / thread_1d, 1)); 216 217 grid = num_elem / elems_per_block + (num_elem % elems_per_block > 0); 218 block[2] = elems_per_block; 219 } 220 CeedInt shared_mem = block[0] * block[1] * block[2] * sizeof(CeedScalar); 221 222 CeedCallBackend(CeedTryRunKernelDimShared_Cuda(ceed, data->op, stream, grid, block[0], block[1], block[2], shared_mem, is_run_good, opargs)); 223 224 // Restore input arrays 225 for (CeedInt i = 0; i < num_input_fields; i++) { 226 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 227 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 228 } else { 229 bool is_active; 230 CeedVector vec; 231 232 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 233 is_active = vec == CEED_VECTOR_ACTIVE; 234 if (!is_active) CeedCallBackend(CeedVectorRestoreArrayRead(vec, &data->fields.inputs[i])); 235 CeedCallBackend(CeedVectorDestroy(&vec)); 236 } 237 } 238 239 // Restore output arrays 240 for (CeedInt i = 0; i < num_output_fields; i++) { 241 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 242 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 243 } else { 244 bool is_active; 245 CeedVector vec; 246 247 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 248 is_active = vec == CEED_VECTOR_ACTIVE; 249 if (!is_active) CeedCallBackend(CeedVectorRestoreArray(vec, &data->fields.outputs[i])); 250 CeedCallBackend(CeedVectorDestroy(&vec)); 251 } 252 } 253 254 // Restore point coordinates, if needed 255 if (is_at_points) { 256 CeedVector vec; 257 258 CeedCallBackend(CeedOperatorAtPointsGetPoints(op, NULL, &vec)); 259 CeedCallBackend(CeedVectorRestoreArrayRead(vec, &data->points.coords)); 260 CeedCallBackend(CeedVectorDestroy(&vec)); 261 } 262 263 // Restore context data 264 CeedCallBackend(CeedQFunctionRestoreInnerContextData(qf, &qf_data->d_c)); 265 266 // Cleanup 267 CeedCallBackend(CeedDestroy(&ceed)); 268 CeedCallBackend(CeedQFunctionDestroy(&qf)); 269 if (!(*is_run_good)) data->use_fallback = true; 270 return CEED_ERROR_SUCCESS; 271 } 272 273 static int CeedOperatorApplyAdd_Cuda_gen(CeedOperator op, CeedVector input_vec, CeedVector output_vec, CeedRequest *request) { 274 bool is_run_good = false; 275 const CeedScalar *input_arr = NULL; 276 CeedScalar *output_arr = NULL; 277 278 // Try to run kernel 279 if (input_vec != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(input_vec, CEED_MEM_DEVICE, &input_arr)); 280 if (output_vec != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArray(output_vec, CEED_MEM_DEVICE, &output_arr)); 281 CeedCallBackend(CeedOperatorApplyAddCore_Cuda_gen(op, NULL, input_arr, output_arr, &is_run_good, request)); 282 if (input_vec != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorRestoreArrayRead(input_vec, &input_arr)); 283 if (output_vec != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorRestoreArray(output_vec, &output_arr)); 284 285 // Fallback on unsuccessful run 286 if (!is_run_good) { 287 CeedOperator op_fallback; 288 289 CeedDebug256(CeedOperatorReturnCeed(op), CEED_DEBUG_COLOR_SUCCESS, "Falling back to /gpu/cuda/ref CeedOperator"); 290 CeedCallBackend(CeedOperatorGetFallback(op, &op_fallback)); 291 CeedCallBackend(CeedOperatorApplyAdd(op_fallback, input_vec, output_vec, request)); 292 } 293 return CEED_ERROR_SUCCESS; 294 } 295 296 static int CeedOperatorApplyAddComposite_Cuda_gen(CeedOperator op, CeedVector input_vec, CeedVector output_vec, CeedRequest *request) { 297 bool is_run_good[CEED_COMPOSITE_MAX] = {false}; 298 CeedInt num_suboperators; 299 const CeedScalar *input_arr = NULL; 300 CeedScalar *output_arr = NULL; 301 Ceed ceed; 302 CeedOperator *sub_operators; 303 304 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 305 CeedCall(CeedCompositeOperatorGetNumSub(op, &num_suboperators)); 306 CeedCall(CeedCompositeOperatorGetSubList(op, &sub_operators)); 307 if (input_vec != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArrayRead(input_vec, CEED_MEM_DEVICE, &input_arr)); 308 if (output_vec != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorGetArray(output_vec, CEED_MEM_DEVICE, &output_arr)); 309 for (CeedInt i = 0; i < num_suboperators; i++) { 310 CeedInt num_elem = 0; 311 312 CeedCall(CeedOperatorGetNumElements(sub_operators[i], &num_elem)); 313 if (num_elem > 0) { 314 cudaStream_t stream = NULL; 315 316 CeedCallCuda(ceed, cudaStreamCreate(&stream)); 317 CeedCallBackend(CeedOperatorApplyAddCore_Cuda_gen(sub_operators[i], stream, input_arr, output_arr, &is_run_good[i], request)); 318 CeedCallCuda(ceed, cudaStreamDestroy(stream)); 319 } 320 } 321 if (input_vec != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorRestoreArrayRead(input_vec, &input_arr)); 322 if (output_vec != CEED_VECTOR_NONE) CeedCallBackend(CeedVectorRestoreArray(output_vec, &output_arr)); 323 CeedCallCuda(ceed, cudaDeviceSynchronize()); 324 325 // Fallback on unsuccessful run 326 for (CeedInt i = 0; i < num_suboperators; i++) { 327 if (!is_run_good[i]) { 328 CeedOperator op_fallback; 329 330 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "Falling back to /gpu/cuda/ref CeedOperator"); 331 CeedCallBackend(CeedOperatorGetFallback(sub_operators[i], &op_fallback)); 332 CeedCallBackend(CeedOperatorApplyAdd(op_fallback, input_vec, output_vec, request)); 333 } 334 } 335 CeedCallBackend(CeedDestroy(&ceed)); 336 return CEED_ERROR_SUCCESS; 337 } 338 339 //------------------------------------------------------------------------------ 340 // Create operator 341 //------------------------------------------------------------------------------ 342 int CeedOperatorCreate_Cuda_gen(CeedOperator op) { 343 bool is_composite; 344 Ceed ceed; 345 CeedOperator_Cuda_gen *impl; 346 347 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 348 CeedCallBackend(CeedCalloc(1, &impl)); 349 CeedCallBackend(CeedOperatorSetData(op, impl)); 350 CeedCall(CeedOperatorIsComposite(op, &is_composite)); 351 if (is_composite) { 352 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAddComposite", CeedOperatorApplyAddComposite_Cuda_gen)); 353 } else { 354 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Cuda_gen)); 355 } 356 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Cuda_gen)); 357 CeedCallBackend(CeedDestroy(&ceed)); 358 return CEED_ERROR_SUCCESS; 359 } 360 361 //------------------------------------------------------------------------------ 362