1 // Copyright (c) 2017-2018, Lawrence Livermore National Security, LLC. 2 // Produced at the Lawrence Livermore National Laboratory. LLNL-CODE-734707. 3 // All Rights reserved. See files LICENSE and NOTICE for details. 4 // 5 // This file is part of CEED, a collection of benchmarks, miniapps, software 6 // libraries and APIs for efficient high-order finite element and spectral 7 // element discretizations for exascale applications. For more information and 8 // source code availability see http://github.com/ceed. 9 // 10 // The CEED research is supported by the Exascale Computing Project 17-SC-20-SC, 11 // a collaborative effort of two U.S. Department of Energy organizations (Office 12 // of Science and the National Nuclear Security Administration) responsible for 13 // the planning and preparation of a capable exascale ecosystem, including 14 // software, applications, hardware, advanced system engineering and early 15 // testbed platforms, in support of the nation's exascale computing imperative. 16 17 #include <ceed/ceed.h> 18 #include <ceed/backend.h> 19 #include <stddef.h> 20 #include "ceed-cuda-gen.h" 21 #include "ceed-cuda-gen-operator-build.h" 22 #include "../cuda/ceed-cuda.h" 23 24 //------------------------------------------------------------------------------ 25 // Destroy operator 26 //------------------------------------------------------------------------------ 27 static int CeedOperatorDestroy_Cuda_gen(CeedOperator op) { 28 int ierr; 29 CeedOperator_Cuda_gen *impl; 30 ierr = CeedOperatorGetData(op, &impl); CeedChkBackend(ierr); 31 ierr = CeedFree(&impl); CeedChkBackend(ierr); 32 return CEED_ERROR_SUCCESS; 33 } 34 35 static int Waste(int threads_per_sm, int warp_size, int threads_per_elem, 36 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 = ((useful_threads_per_block + warp_size - 1) / warp_size) * 40 warp_size; 41 int blocks_per_sm = threads_per_sm / block_size; 42 return threads_per_sm - useful_threads_per_block * blocks_per_sm; 43 } 44 45 // Choose the least wasteful block size constrained by blocks_per_sm of 46 // max_threads_per_block. 47 // 48 // The x and y part of block[] contains per-element sizes (specified on input) 49 // while the z part is number of elements. 50 // 51 // Problem setting: we'd like to make occupancy high with relatively few 52 // inactive threads. CUDA (cuOccupancyMaxPotentialBlockSize) can tell us how 53 // many threads can run. 54 // 55 // Note that full occupancy sometimes can't be achieved by one thread block. For 56 // example, an SM might support 1536 threads in total, but only 1024 within a 57 // single thread block. So cuOccupancyMaxPotentialBlockSize may suggest a block 58 // size of 768 so that two blocks can run, versus one block of 1024 will prevent 59 // a second block from running. The cuda-gen kernels are pretty heavy with lots 60 // of instruction-level parallelism (ILP) so we'll generally be okay with 61 // relatvely low occupancy and smaller thread blocks, but we solve a reasonably 62 // general problem here. Empirically, we find that blocks bigger than about 256 63 // have higher latency and worse load balancing when the number of elements is 64 // modest. 65 // 66 // cuda-gen can't choose block sizes arbitrarily; they need to be a multiple of 67 // the number of quadrature points (or number of basis functions). They also 68 // have a lot of __syncthreads(), which is another point against excessively 69 // large thread blocks. Suppose I have elements with 7x7x7 quadrature points. 70 // This will loop over the last dimension, so we have 7*7=49 threads per 71 // element. Suppose we have two elements = 2*49=98 useful threads. CUDA 72 // schedules in units of full warps (32 threads), so 128 CUDA hardware threads 73 // are effectively committed to that block. Now suppose 74 // cuOccupancyMaxPotentialBlockSize returned 352. We can schedule 2 blocks of 75 // size 98 (196 useful threads using 256 hardware threads), but not a third 76 // block (which would need a total of 384 hardware threads). 77 // 78 // If instead, we had packed 3 elements, we'd have 3*49=147 useful threads 79 // occupying 160 slots, and could schedule two blocks. Alternatively, we could 80 // pack a single block of 7 elements (2*49=343 useful threads) into the 354 81 // slots. The latter has the least "waste", but __syncthreads() 82 // over-synchronizes and it might not pay off relative to smaller blocks. 83 static int BlockGridCalculate(CeedInt nelem, int blocks_per_sm, 84 int max_threads_per_block, int warp_size, int block[3], int *grid) { 85 const int threads_per_sm = blocks_per_sm * max_threads_per_block; 86 const int threads_per_elem = block[0] * block[1]; 87 int elems_per_block = 1; 88 int waste = Waste(threads_per_sm, warp_size, threads_per_elem, 1); 89 for (int i=2; 90 i <= CeedIntMin(max_threads_per_block / threads_per_elem, nelem); 91 i++) { 92 int i_waste = Waste(threads_per_sm, warp_size, threads_per_elem, i); 93 // We want to minimize waste, but smaller kernels have lower latency and 94 // less __syncthreads() overhead so when a larger block size has the same 95 // waste as a smaller one, go ahead and prefer the smaller block. 96 if (i_waste < waste || (i_waste == waste && threads_per_elem * i <= 128)) { 97 elems_per_block = i; 98 waste = i_waste; 99 } 100 } 101 block[2] = elems_per_block; 102 *grid = (nelem + elems_per_block - 1) / elems_per_block; 103 return CEED_ERROR_SUCCESS; 104 } 105 106 // callback for cuOccupancyMaxPotentialBlockSize, providing the amount of 107 // dynamic shared memory required for a thread block of size threads. 108 static size_t dynamicSMemSize(int threads) { return threads * sizeof(CeedScalar); } 109 110 //------------------------------------------------------------------------------ 111 // Apply and add to output 112 //------------------------------------------------------------------------------ 113 static int CeedOperatorApplyAdd_Cuda_gen(CeedOperator op, CeedVector invec, 114 CeedVector outvec, CeedRequest *request) { 115 int ierr; 116 Ceed ceed; 117 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 118 Ceed_Cuda *cuda_data; 119 ierr = CeedGetData(ceed, &cuda_data); CeedChkBackend(ierr); 120 CeedOperator_Cuda_gen *data; 121 ierr = CeedOperatorGetData(op, &data); CeedChkBackend(ierr); 122 CeedQFunction qf; 123 CeedQFunction_Cuda_gen *qf_data; 124 ierr = CeedOperatorGetQFunction(op, &qf); CeedChkBackend(ierr); 125 ierr = CeedQFunctionGetData(qf, &qf_data); CeedChkBackend(ierr); 126 CeedInt nelem, numinputfields, numoutputfields; 127 ierr = CeedOperatorGetNumElements(op, &nelem); CeedChkBackend(ierr); 128 ierr = CeedQFunctionGetNumArgs(qf, &numinputfields, &numoutputfields); 129 CeedChkBackend(ierr); 130 CeedOperatorField *opinputfields, *opoutputfields; 131 ierr = CeedOperatorGetFields(op, &opinputfields, &opoutputfields); 132 CeedChkBackend(ierr); 133 CeedQFunctionField *qfinputfields, *qfoutputfields; 134 ierr = CeedQFunctionGetFields(qf, &qfinputfields, &qfoutputfields); 135 CeedChkBackend(ierr); 136 CeedEvalMode emode; 137 CeedVector vec, outvecs[16] = {}; 138 139 // Creation of the operator 140 ierr = CeedCudaGenOperatorBuild(op); CeedChkBackend(ierr); 141 142 // Input vectors 143 for (CeedInt i = 0; i < numinputfields; i++) { 144 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 145 CeedChkBackend(ierr); 146 if (emode == CEED_EVAL_WEIGHT) { // Skip 147 data->fields.in[i] = NULL; 148 } else { 149 // Get input vector 150 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 151 if (vec == CEED_VECTOR_ACTIVE) vec = invec; 152 ierr = CeedVectorGetArrayRead(vec, CEED_MEM_DEVICE, &data->fields.in[i]); 153 CeedChkBackend(ierr); 154 } 155 } 156 157 // Output vectors 158 for (CeedInt i = 0; i < numoutputfields; i++) { 159 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 160 CeedChkBackend(ierr); 161 if (emode == CEED_EVAL_WEIGHT) { // Skip 162 data->fields.out[i] = NULL; 163 } else { 164 // Get output vector 165 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); 166 CeedChkBackend(ierr); 167 if (vec == CEED_VECTOR_ACTIVE) vec = outvec; 168 outvecs[i] = vec; 169 // Check for multiple output modes 170 CeedInt index = -1; 171 for (CeedInt j = 0; j < i; j++) { 172 if (vec == outvecs[j]) { 173 index = j; 174 break; 175 } 176 } 177 if (index == -1) { 178 ierr = CeedVectorGetArray(vec, CEED_MEM_DEVICE, &data->fields.out[i]); 179 CeedChkBackend(ierr); 180 } else { 181 data->fields.out[i] = data->fields.out[index]; 182 } 183 } 184 } 185 186 // Get context data 187 CeedQFunctionContext ctx; 188 ierr = CeedQFunctionGetInnerContext(qf, &ctx); CeedChkBackend(ierr); 189 if (ctx) { 190 ierr = CeedQFunctionContextGetData(ctx, CEED_MEM_DEVICE, &qf_data->d_c); 191 CeedChkBackend(ierr); 192 } 193 194 // Apply operator 195 void *opargs[] = {(void *) &nelem, &qf_data->d_c, &data->indices, 196 &data->fields, &data->B, &data->G, &data->W 197 }; 198 const CeedInt dim = data->dim; 199 const CeedInt Q1d = data->Q1d; 200 const CeedInt P1d = data->maxP1d; 201 const CeedInt thread1d = CeedIntMax(Q1d, P1d); 202 int max_threads_per_block, min_grid_size; 203 CeedChk_Cu(ceed, cuOccupancyMaxPotentialBlockSize(&min_grid_size, 204 &max_threads_per_block, data->op, dynamicSMemSize, 0, 0x10000)); 205 int block[3] = {thread1d, dim < 2 ? 1 : thread1d, -1,}, grid; 206 CeedChkBackend(BlockGridCalculate(nelem, 207 min_grid_size/ cuda_data->deviceProp.multiProcessorCount, max_threads_per_block, 208 cuda_data->deviceProp.warpSize, block, &grid)); 209 CeedInt shared_mem = block[0] * block[1] * block[2] * sizeof(CeedScalar); 210 ierr = CeedRunKernelDimSharedCuda(ceed, data->op, grid, block[0], block[1], 211 block[2], shared_mem, opargs); 212 CeedChkBackend(ierr); 213 214 // Restore input arrays 215 for (CeedInt i = 0; i < numinputfields; i++) { 216 ierr = CeedQFunctionFieldGetEvalMode(qfinputfields[i], &emode); 217 CeedChkBackend(ierr); 218 if (emode == CEED_EVAL_WEIGHT) { // Skip 219 } else { 220 ierr = CeedOperatorFieldGetVector(opinputfields[i], &vec); CeedChkBackend(ierr); 221 if (vec == CEED_VECTOR_ACTIVE) vec = invec; 222 ierr = CeedVectorRestoreArrayRead(vec, &data->fields.in[i]); 223 CeedChkBackend(ierr); 224 } 225 } 226 227 // Restore output arrays 228 for (CeedInt i = 0; i < numoutputfields; i++) { 229 ierr = CeedQFunctionFieldGetEvalMode(qfoutputfields[i], &emode); 230 CeedChkBackend(ierr); 231 if (emode == CEED_EVAL_WEIGHT) { // Skip 232 } else { 233 ierr = CeedOperatorFieldGetVector(opoutputfields[i], &vec); 234 CeedChkBackend(ierr); 235 if (vec == CEED_VECTOR_ACTIVE) vec = outvec; 236 // Check for multiple output modes 237 CeedInt index = -1; 238 for (CeedInt j = 0; j < i; j++) { 239 if (vec == outvecs[j]) { 240 index = j; 241 break; 242 } 243 } 244 if (index == -1) { 245 ierr = CeedVectorRestoreArray(vec, &data->fields.out[i]); 246 CeedChkBackend(ierr); 247 } 248 } 249 } 250 251 // Restore context data 252 if (ctx) { 253 ierr = CeedQFunctionContextRestoreData(ctx, &qf_data->d_c); 254 CeedChkBackend(ierr); 255 } 256 return CEED_ERROR_SUCCESS; 257 } 258 259 //------------------------------------------------------------------------------ 260 // Create operator 261 //------------------------------------------------------------------------------ 262 int CeedOperatorCreate_Cuda_gen(CeedOperator op) { 263 int ierr; 264 Ceed ceed; 265 ierr = CeedOperatorGetCeed(op, &ceed); CeedChkBackend(ierr); 266 CeedOperator_Cuda_gen *impl; 267 268 ierr = CeedCalloc(1, &impl); CeedChkBackend(ierr); 269 ierr = CeedOperatorSetData(op, impl); CeedChkBackend(ierr); 270 271 ierr = CeedSetBackendFunction(ceed, "Operator", op, "ApplyAdd", 272 CeedOperatorApplyAdd_Cuda_gen); CeedChkBackend(ierr); 273 ierr = CeedSetBackendFunction(ceed, "Operator", op, "Destroy", 274 CeedOperatorDestroy_Cuda_gen); CeedChkBackend(ierr); 275 return CEED_ERROR_SUCCESS; 276 } 277 //------------------------------------------------------------------------------ 278