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