1 // Copyright (c) 2017-2022, 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 <stddef.h> 12 13 #include "../cuda/ceed-cuda-common.h" 14 #include "../cuda/ceed-cuda-compile.h" 15 #include "ceed-cuda-gen-operator-build.h" 16 #include "ceed-cuda-gen.h" 17 18 //------------------------------------------------------------------------------ 19 // Destroy operator 20 //------------------------------------------------------------------------------ 21 static int CeedOperatorDestroy_Cuda_gen(CeedOperator op) { 22 CeedOperator_Cuda_gen *impl; 23 CeedCallBackend(CeedOperatorGetData(op, &impl)); 24 CeedCallBackend(CeedFree(&impl)); 25 return CEED_ERROR_SUCCESS; 26 } 27 28 static int Waste(int threads_per_sm, int warp_size, int threads_per_elem, int elems_per_block) { 29 int useful_threads_per_block = threads_per_elem * elems_per_block; 30 // round up to nearest multiple of warp_size 31 int block_size = CeedDivUpInt(useful_threads_per_block, warp_size) * warp_size; 32 int blocks_per_sm = threads_per_sm / block_size; 33 return threads_per_sm - useful_threads_per_block * blocks_per_sm; 34 } 35 36 // Choose the least wasteful block size constrained by blocks_per_sm of max_threads_per_block. 37 // 38 // The x and y part of block[] contains per-element sizes (specified on input) while the z part is number of elements. 39 // 40 // Problem setting: we'd like to make occupancy high with relatively few inactive threads. CUDA (cuOccupancyMaxPotentialBlockSize) can tell us how 41 // many threads can run. 42 // 43 // Note that full occupancy sometimes can't be achieved by one thread block. 44 // For example, an SM might support 1536 threads in total, but only 1024 within a single thread block. 45 // So cuOccupancyMaxPotentialBlockSize may suggest a block size of 768 so that two blocks can run, versus one block of 1024 will prevent a second 46 // block from running. The cuda-gen kernels are pretty heavy with lots of instruction-level parallelism (ILP) so we'll generally be okay with 47 // relatively low occupancy and smaller thread blocks, but we solve a reasonably general problem here. Empirically, we find that blocks bigger than 48 // about 256 have higher latency and worse load balancing when the number of elements is modest. 49 // 50 // 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). 51 // They also have a lot of __syncthreads(), which is another point against excessively large thread blocks. 52 // Suppose I have elements with 7x7x7 quadrature points. 53 // This will loop over the last dimension, so we have 7*7=49 threads per element. 54 // Suppose we have two elements = 2*49=98 useful threads. 55 // CUDA schedules in units of full warps (32 threads), so 128 CUDA hardware threads are effectively committed to that block. 56 // Now suppose cuOccupancyMaxPotentialBlockSize returned 352. 57 // 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 58 // hardware threads). 59 // 60 // If instead, we had packed 3 elements, we'd have 3*49=147 useful threads occupying 160 slots, and could schedule two blocks. 61 // Alternatively, we could pack a single block of 7 elements (2*49=343 useful threads) into the 354 slots. 62 // The latter has the least "waste", but __syncthreads() over-synchronizes and it might not pay off relative to smaller blocks. 63 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], 64 int *grid) { 65 const int threads_per_sm = blocks_per_sm * max_threads_per_block; 66 const int threads_per_elem = block[0] * block[1]; 67 int elems_per_block = 1; 68 int waste = Waste(threads_per_sm, warp_size, threads_per_elem, 1); 69 for (int i = 2; i <= CeedIntMin(max_threads_per_block / threads_per_elem, num_elem); i++) { 70 int i_waste = Waste(threads_per_sm, warp_size, threads_per_elem, i); 71 // We want to minimize waste, but smaller kernels have lower latency and less __syncthreads() overhead so when a larger block size has the same 72 // waste as a smaller one, go ahead and prefer the smaller block. 73 if (i_waste < waste || (i_waste == waste && threads_per_elem * i <= 128)) { 74 elems_per_block = i; 75 waste = i_waste; 76 } 77 } 78 // 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 79 // check before setting the z-dimension size of the block. 80 block[2] = CeedIntMin(elems_per_block, max_threads_z); 81 *grid = CeedDivUpInt(num_elem, elems_per_block); 82 return CEED_ERROR_SUCCESS; 83 } 84 85 // callback for cuOccupancyMaxPotentialBlockSize, providing the amount of dynamic shared memory required for a thread block of size threads. 86 static size_t dynamicSMemSize(int threads) { return threads * sizeof(CeedScalar); } 87 88 //------------------------------------------------------------------------------ 89 // Apply and add to output 90 //------------------------------------------------------------------------------ 91 static int CeedOperatorApplyAdd_Cuda_gen(CeedOperator op, CeedVector input_vec, CeedVector output_vec, CeedRequest *request) { 92 Ceed ceed; 93 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 94 Ceed_Cuda *cuda_data; 95 CeedCallBackend(CeedGetData(ceed, &cuda_data)); 96 CeedOperator_Cuda_gen *data; 97 CeedCallBackend(CeedOperatorGetData(op, &data)); 98 CeedQFunction qf; 99 CeedQFunction_Cuda_gen *qf_data; 100 CeedCallBackend(CeedOperatorGetQFunction(op, &qf)); 101 CeedCallBackend(CeedQFunctionGetData(qf, &qf_data)); 102 CeedInt num_elem, num_input_fields, num_output_fields; 103 CeedCallBackend(CeedOperatorGetNumElements(op, &num_elem)); 104 CeedOperatorField *op_input_fields, *op_output_fields; 105 CeedCallBackend(CeedOperatorGetFields(op, &num_input_fields, &op_input_fields, &num_output_fields, &op_output_fields)); 106 CeedQFunctionField *qf_input_fields, *qf_output_fields; 107 CeedCallBackend(CeedQFunctionGetFields(qf, NULL, &qf_input_fields, NULL, &qf_output_fields)); 108 CeedEvalMode eval_mode; 109 CeedVector vec, output_vecs[CEED_FIELD_MAX] = {NULL}; 110 111 // Creation of the operator 112 CeedCallBackend(CeedOperatorBuildKernel_Cuda_gen(op)); 113 114 // Input vectors 115 for (CeedInt i = 0; i < num_input_fields; i++) { 116 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 117 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 118 data->fields.inputs[i] = NULL; 119 } else { 120 // Get input vector 121 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 122 if (vec == CEED_VECTOR_ACTIVE) vec = input_vec; 123 CeedCallBackend(CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->fields.inputs[i])); 124 } 125 } 126 127 // Output vectors 128 for (CeedInt i = 0; i < num_output_fields; i++) { 129 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 130 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 131 data->fields.outputs[i] = NULL; 132 } else { 133 // Get output vector 134 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 135 if (vec == CEED_VECTOR_ACTIVE) vec = output_vec; 136 output_vecs[i] = vec; 137 // Check for multiple output modes 138 CeedInt index = -1; 139 for (CeedInt j = 0; j < i; j++) { 140 if (vec == output_vecs[j]) { 141 index = j; 142 break; 143 } 144 } 145 if (index == -1) { 146 CeedCallBackend(CeedVectorGetArray(vec, CEED_MEM_DEVICE, &data->fields.outputs[i])); 147 } else { 148 data->fields.outputs[i] = data->fields.outputs[index]; 149 } 150 } 151 } 152 153 // Get context data 154 CeedCallBackend(CeedQFunctionGetInnerContextData(qf, CEED_MEM_DEVICE, &qf_data->d_c)); 155 156 // Apply operator 157 void *opargs[] = {(void *)&num_elem, &qf_data->d_c, &data->indices, &data->fields, &data->B, &data->G, &data->W}; 158 const CeedInt dim = data->dim; 159 const CeedInt Q_1d = data->Q_1d; 160 const CeedInt P_1d = data->max_P_1d; 161 const CeedInt thread_1d = CeedIntMax(Q_1d, P_1d); 162 int max_threads_per_block, min_grid_size; 163 CeedCallCuda(ceed, cuOccupancyMaxPotentialBlockSize(&min_grid_size, &max_threads_per_block, data->op, dynamicSMemSize, 0, 0x10000)); 164 int block[3] = 165 { 166 thread_1d, 167 dim < 2 ? 1 : thread_1d, 168 -1, 169 }, 170 grid; 171 CeedChkBackend(BlockGridCalculate(num_elem, min_grid_size / cuda_data->device_prop.multiProcessorCount, max_threads_per_block, 172 cuda_data->device_prop.maxThreadsDim[2], cuda_data->device_prop.warpSize, block, &grid)); 173 CeedInt shared_mem = block[0] * block[1] * block[2] * sizeof(CeedScalar); 174 CeedCallBackend(CeedRunKernelDimShared_Cuda(ceed, data->op, grid, block[0], block[1], block[2], shared_mem, opargs)); 175 176 // Restore input arrays 177 for (CeedInt i = 0; i < num_input_fields; i++) { 178 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_input_fields[i], &eval_mode)); 179 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 180 } else { 181 CeedCallBackend(CeedOperatorFieldGetVector(op_input_fields[i], &vec)); 182 if (vec == CEED_VECTOR_ACTIVE) vec = input_vec; 183 CeedCallBackend(CeedVectorRestoreArrayRead(vec, &data->fields.inputs[i])); 184 } 185 } 186 187 // Restore output arrays 188 for (CeedInt i = 0; i < num_output_fields; i++) { 189 CeedCallBackend(CeedQFunctionFieldGetEvalMode(qf_output_fields[i], &eval_mode)); 190 if (eval_mode == CEED_EVAL_WEIGHT) { // Skip 191 } else { 192 CeedCallBackend(CeedOperatorFieldGetVector(op_output_fields[i], &vec)); 193 if (vec == CEED_VECTOR_ACTIVE) vec = output_vec; 194 // Check for multiple output modes 195 CeedInt index = -1; 196 for (CeedInt j = 0; j < i; j++) { 197 if (vec == output_vecs[j]) { 198 index = j; 199 break; 200 } 201 } 202 if (index == -1) { 203 CeedCallBackend(CeedVectorRestoreArray(vec, &data->fields.outputs[i])); 204 } 205 } 206 } 207 208 // Restore context data 209 CeedCallBackend(CeedQFunctionRestoreInnerContextData(qf, &qf_data->d_c)); 210 211 return CEED_ERROR_SUCCESS; 212 } 213 214 //------------------------------------------------------------------------------ 215 // Create operator 216 //------------------------------------------------------------------------------ 217 int CeedOperatorCreate_Cuda_gen(CeedOperator op) { 218 Ceed ceed; 219 CeedCallBackend(CeedOperatorGetCeed(op, &ceed)); 220 CeedOperator_Cuda_gen *impl; 221 222 CeedCallBackend(CeedCalloc(1, &impl)); 223 CeedCallBackend(CeedOperatorSetData(op, impl)); 224 225 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", CeedOperatorApplyAdd_Cuda_gen)); 226 CeedCallBackend(CeedSetBackendFunction(ceed, "Operator", op, "Destroy", CeedOperatorDestroy_Cuda_gen)); 227 return CEED_ERROR_SUCCESS; 228 } 229 230 //------------------------------------------------------------------------------ 231