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