1 // Copyright (c) 2017-2024, 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-tools.h> 11 #include <cuda.h> 12 #include <cuda_runtime.h> 13 #include <stdbool.h> 14 #include <stddef.h> 15 #include <string.h> 16 17 #include "../cuda/ceed-cuda-common.h" 18 #include "../cuda/ceed-cuda-compile.h" 19 #include "ceed-cuda-ref.h" 20 21 //------------------------------------------------------------------------------ 22 // Compile restriction kernels 23 //------------------------------------------------------------------------------ 24 static inline int CeedElemRestrictionSetupCompile_Cuda(CeedElemRestriction rstr) { 25 Ceed ceed; 26 bool is_deterministic; 27 char *restriction_kernel_source; 28 const char *restriction_kernel_path; 29 CeedInt num_elem, num_comp, elem_size, comp_stride; 30 CeedRestrictionType rstr_type; 31 CeedElemRestriction_Cuda *impl; 32 33 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 34 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 35 CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type)); 36 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 37 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 38 CeedCallBackend(CeedElemRestrictionGetCompStride(rstr, &comp_stride)); 39 if (rstr_type == CEED_RESTRICTION_POINTS) { 40 CeedCallBackend(CeedElemRestrictionGetMaxPointsInElement(rstr, &elem_size)); 41 } else { 42 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 43 } 44 is_deterministic = impl->d_l_vec_indices != NULL; 45 46 // Compile CUDA kernels 47 switch (rstr_type) { 48 case CEED_RESTRICTION_STRIDED: { 49 bool has_backend_strides; 50 CeedInt strides[3] = {1, num_elem * elem_size, elem_size}; 51 52 CeedCallBackend(CeedElemRestrictionHasBackendStrides(rstr, &has_backend_strides)); 53 if (!has_backend_strides) { 54 CeedCallBackend(CeedElemRestrictionGetStrides(rstr, strides)); 55 } 56 57 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-restriction-strided.h", &restriction_kernel_path)); 58 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source -----\n"); 59 CeedCallBackend(CeedLoadSourceToBuffer(ceed, restriction_kernel_path, &restriction_kernel_source)); 60 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source Complete! -----\n"); 61 CeedCallBackend(CeedCompile_Cuda(ceed, restriction_kernel_source, &impl->module, 6, "RSTR_ELEM_SIZE", elem_size, "RSTR_NUM_ELEM", num_elem, 62 "RSTR_NUM_COMP", num_comp, "RSTR_STRIDE_NODES", strides[0], "RSTR_STRIDE_COMP", strides[1], "RSTR_STRIDE_ELEM", 63 strides[2])); 64 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "StridedNoTranspose", &impl->ApplyNoTranspose)); 65 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "StridedTranspose", &impl->ApplyTranspose)); 66 } break; 67 case CEED_RESTRICTION_POINTS: { 68 const char *offset_kernel_path; 69 char **file_paths = NULL; 70 CeedInt num_file_paths = 0; 71 72 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-restriction-at-points.h", &restriction_kernel_path)); 73 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source -----\n"); 74 CeedCallBackend(CeedLoadSourceAndInitializeBuffer(ceed, restriction_kernel_path, &num_file_paths, &file_paths, &restriction_kernel_source)); 75 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-restriction-offset.h", &offset_kernel_path)); 76 CeedCallBackend(CeedLoadSourceToInitializedBuffer(ceed, offset_kernel_path, &num_file_paths, &file_paths, &restriction_kernel_source)); 77 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source Complete! -----\n"); 78 CeedCallBackend(CeedCompile_Cuda(ceed, restriction_kernel_source, &impl->module, 6, "RSTR_ELEM_SIZE", elem_size, "RSTR_NUM_ELEM", num_elem, 79 "RSTR_NUM_COMP", num_comp, "RSTR_NUM_NODES", impl->num_nodes, "RSTR_COMP_STRIDE", comp_stride, 80 "USE_DETERMINISTIC", is_deterministic ? 1 : 0)); 81 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OffsetNoTranspose", &impl->ApplyNoTranspose)); 82 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "AtPointsTranspose", &impl->ApplyTranspose)); 83 } break; 84 case CEED_RESTRICTION_STANDARD: { 85 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-restriction-offset.h", &restriction_kernel_path)); 86 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source -----\n"); 87 CeedCallBackend(CeedLoadSourceToBuffer(ceed, restriction_kernel_path, &restriction_kernel_source)); 88 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source Complete! -----\n"); 89 CeedCallBackend(CeedCompile_Cuda(ceed, restriction_kernel_source, &impl->module, 6, "RSTR_ELEM_SIZE", elem_size, "RSTR_NUM_ELEM", num_elem, 90 "RSTR_NUM_COMP", num_comp, "RSTR_NUM_NODES", impl->num_nodes, "RSTR_COMP_STRIDE", comp_stride, 91 "USE_DETERMINISTIC", is_deterministic ? 1 : 0)); 92 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OffsetNoTranspose", &impl->ApplyNoTranspose)); 93 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OffsetTranspose", &impl->ApplyTranspose)); 94 } break; 95 case CEED_RESTRICTION_ORIENTED: { 96 const char *offset_kernel_path; 97 char **file_paths = NULL; 98 CeedInt num_file_paths = 0; 99 100 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-restriction-oriented.h", &restriction_kernel_path)); 101 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source -----\n"); 102 CeedCallBackend(CeedLoadSourceAndInitializeBuffer(ceed, restriction_kernel_path, &num_file_paths, &file_paths, &restriction_kernel_source)); 103 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-restriction-offset.h", &offset_kernel_path)); 104 CeedCallBackend(CeedLoadSourceToInitializedBuffer(ceed, offset_kernel_path, &num_file_paths, &file_paths, &restriction_kernel_source)); 105 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source Complete! -----\n"); 106 CeedCallBackend(CeedCompile_Cuda(ceed, restriction_kernel_source, &impl->module, 6, "RSTR_ELEM_SIZE", elem_size, "RSTR_NUM_ELEM", num_elem, 107 "RSTR_NUM_COMP", num_comp, "RSTR_NUM_NODES", impl->num_nodes, "RSTR_COMP_STRIDE", comp_stride, 108 "USE_DETERMINISTIC", is_deterministic ? 1 : 0)); 109 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OrientedNoTranspose", &impl->ApplyNoTranspose)); 110 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OffsetNoTranspose", &impl->ApplyUnsignedNoTranspose)); 111 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OrientedTranspose", &impl->ApplyTranspose)); 112 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OffsetTranspose", &impl->ApplyUnsignedTranspose)); 113 // Cleanup 114 CeedCallBackend(CeedFree(&offset_kernel_path)); 115 for (CeedInt i = 0; i < num_file_paths; i++) CeedCallBackend(CeedFree(&file_paths[i])); 116 CeedCallBackend(CeedFree(&file_paths)); 117 } break; 118 case CEED_RESTRICTION_CURL_ORIENTED: { 119 const char *offset_kernel_path; 120 char **file_paths = NULL; 121 CeedInt num_file_paths = 0; 122 123 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-restriction-curl-oriented.h", &restriction_kernel_path)); 124 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source -----\n"); 125 CeedCallBackend(CeedLoadSourceAndInitializeBuffer(ceed, restriction_kernel_path, &num_file_paths, &file_paths, &restriction_kernel_source)); 126 CeedCallBackend(CeedGetJitAbsolutePath(ceed, "ceed/jit-source/cuda/cuda-ref-restriction-offset.h", &offset_kernel_path)); 127 CeedCallBackend(CeedLoadSourceToInitializedBuffer(ceed, offset_kernel_path, &num_file_paths, &file_paths, &restriction_kernel_source)); 128 CeedDebug256(ceed, CEED_DEBUG_COLOR_SUCCESS, "----- Loading Restriction Kernel Source Complete! -----\n"); 129 CeedCallBackend(CeedCompile_Cuda(ceed, restriction_kernel_source, &impl->module, 6, "RSTR_ELEM_SIZE", elem_size, "RSTR_NUM_ELEM", num_elem, 130 "RSTR_NUM_COMP", num_comp, "RSTR_NUM_NODES", impl->num_nodes, "RSTR_COMP_STRIDE", comp_stride, 131 "USE_DETERMINISTIC", is_deterministic ? 1 : 0)); 132 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "CurlOrientedNoTranspose", &impl->ApplyNoTranspose)); 133 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "CurlOrientedUnsignedNoTranspose", &impl->ApplyUnsignedNoTranspose)); 134 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OffsetNoTranspose", &impl->ApplyUnorientedNoTranspose)); 135 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "CurlOrientedTranspose", &impl->ApplyTranspose)); 136 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "CurlOrientedUnsignedTranspose", &impl->ApplyUnsignedTranspose)); 137 CeedCallBackend(CeedGetKernel_Cuda(ceed, impl->module, "OffsetTranspose", &impl->ApplyUnorientedTranspose)); 138 // Cleanup 139 CeedCallBackend(CeedFree(&offset_kernel_path)); 140 for (CeedInt i = 0; i < num_file_paths; i++) CeedCallBackend(CeedFree(&file_paths[i])); 141 CeedCallBackend(CeedFree(&file_paths)); 142 } break; 143 } 144 CeedCallBackend(CeedFree(&restriction_kernel_path)); 145 CeedCallBackend(CeedFree(&restriction_kernel_source)); 146 return CEED_ERROR_SUCCESS; 147 } 148 149 //------------------------------------------------------------------------------ 150 // Core apply restriction code 151 //------------------------------------------------------------------------------ 152 static inline int CeedElemRestrictionApply_Cuda_Core(CeedElemRestriction rstr, CeedTransposeMode t_mode, bool use_signs, bool use_orients, 153 CeedVector u, CeedVector v, CeedRequest *request) { 154 Ceed ceed; 155 CeedRestrictionType rstr_type; 156 const CeedScalar *d_u; 157 CeedScalar *d_v; 158 CeedElemRestriction_Cuda *impl; 159 160 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 161 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 162 CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type)); 163 164 // Assemble kernel if needed 165 if (!impl->module) { 166 CeedCallBackend(CeedElemRestrictionSetupCompile_Cuda(rstr)); 167 } 168 169 // Get vectors 170 CeedCallBackend(CeedVectorGetArrayRead(u, CEED_MEM_DEVICE, &d_u)); 171 if (t_mode == CEED_TRANSPOSE) { 172 // Sum into for transpose mode, e-vec to l-vec 173 CeedCallBackend(CeedVectorGetArray(v, CEED_MEM_DEVICE, &d_v)); 174 } else { 175 // Overwrite for notranspose mode, l-vec to e-vec 176 CeedCallBackend(CeedVectorGetArrayWrite(v, CEED_MEM_DEVICE, &d_v)); 177 } 178 179 // Restrict 180 if (t_mode == CEED_NOTRANSPOSE) { 181 // L-vector -> E-vector 182 CeedInt elem_size; 183 184 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 185 const CeedInt block_size = elem_size < 1024 ? (elem_size > 32 ? elem_size : 32) : 1024; 186 const CeedInt grid = CeedDivUpInt(impl->num_nodes, block_size); 187 188 switch (rstr_type) { 189 case CEED_RESTRICTION_STRIDED: { 190 void *args[] = {&d_u, &d_v}; 191 192 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyNoTranspose, grid, block_size, args)); 193 } break; 194 case CEED_RESTRICTION_POINTS: 195 case CEED_RESTRICTION_STANDARD: { 196 void *args[] = {&impl->d_offsets, &d_u, &d_v}; 197 198 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyNoTranspose, grid, block_size, args)); 199 } break; 200 case CEED_RESTRICTION_ORIENTED: { 201 if (use_signs) { 202 void *args[] = {&impl->d_offsets, &impl->d_orients, &d_u, &d_v}; 203 204 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyNoTranspose, grid, block_size, args)); 205 } else { 206 void *args[] = {&impl->d_offsets, &d_u, &d_v}; 207 208 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyUnsignedNoTranspose, grid, block_size, args)); 209 } 210 } break; 211 case CEED_RESTRICTION_CURL_ORIENTED: { 212 if (use_signs && use_orients) { 213 void *args[] = {&impl->d_offsets, &impl->d_curl_orients, &d_u, &d_v}; 214 215 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyNoTranspose, grid, block_size, args)); 216 } else if (use_orients) { 217 void *args[] = {&impl->d_offsets, &impl->d_curl_orients, &d_u, &d_v}; 218 219 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyUnsignedNoTranspose, grid, block_size, args)); 220 } else { 221 void *args[] = {&impl->d_offsets, &d_u, &d_v}; 222 223 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyUnorientedNoTranspose, grid, block_size, args)); 224 } 225 } break; 226 } 227 } else { 228 // E-vector -> L-vector 229 const bool is_deterministic = impl->d_l_vec_indices != NULL; 230 const CeedInt block_size = 32; 231 const CeedInt grid = CeedDivUpInt(impl->num_nodes, block_size); 232 233 switch (rstr_type) { 234 case CEED_RESTRICTION_STRIDED: { 235 void *args[] = {&d_u, &d_v}; 236 237 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args)); 238 } break; 239 case CEED_RESTRICTION_POINTS: { 240 if (!is_deterministic) { 241 void *args[] = {&impl->d_offsets, &impl->d_points_per_elem, &d_u, &d_v}; 242 243 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args)); 244 } else { 245 void *args[] = {&impl->d_l_vec_indices, &impl->d_t_indices, &impl->d_points_per_elem, &impl->d_t_offsets, &d_u, &d_v}; 246 247 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args)); 248 } 249 } break; 250 case CEED_RESTRICTION_STANDARD: { 251 if (!is_deterministic) { 252 void *args[] = {&impl->d_offsets, &d_u, &d_v}; 253 254 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args)); 255 } else { 256 void *args[] = {&impl->d_l_vec_indices, &impl->d_t_indices, &impl->d_t_offsets, &d_u, &d_v}; 257 258 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args)); 259 } 260 } break; 261 case CEED_RESTRICTION_ORIENTED: { 262 if (use_signs) { 263 if (!is_deterministic) { 264 void *args[] = {&impl->d_offsets, &impl->d_orients, &d_u, &d_v}; 265 266 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args)); 267 } else { 268 void *args[] = {&impl->d_l_vec_indices, &impl->d_t_indices, &impl->d_t_offsets, &impl->d_orients, &d_u, &d_v}; 269 270 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args)); 271 } 272 } else { 273 if (!is_deterministic) { 274 void *args[] = {&impl->d_offsets, &d_u, &d_v}; 275 276 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyUnsignedTranspose, grid, block_size, args)); 277 } else { 278 void *args[] = {&impl->d_l_vec_indices, &impl->d_t_indices, &impl->d_t_offsets, &d_u, &d_v}; 279 280 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyUnsignedTranspose, grid, block_size, args)); 281 } 282 } 283 } break; 284 case CEED_RESTRICTION_CURL_ORIENTED: { 285 if (use_signs && use_orients) { 286 if (!is_deterministic) { 287 void *args[] = {&impl->d_offsets, &impl->d_curl_orients, &d_u, &d_v}; 288 289 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args)); 290 } else { 291 void *args[] = {&impl->d_l_vec_indices, &impl->d_t_indices, &impl->d_t_offsets, &impl->d_curl_orients, &d_u, &d_v}; 292 293 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyTranspose, grid, block_size, args)); 294 } 295 } else if (use_orients) { 296 if (!is_deterministic) { 297 void *args[] = {&impl->d_offsets, &impl->d_curl_orients, &d_u, &d_v}; 298 299 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyUnsignedTranspose, grid, block_size, args)); 300 } else { 301 void *args[] = {&impl->d_l_vec_indices, &impl->d_t_indices, &impl->d_t_offsets, &impl->d_curl_orients, &d_u, &d_v}; 302 303 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyUnsignedTranspose, grid, block_size, args)); 304 } 305 } else { 306 if (!is_deterministic) { 307 void *args[] = {&impl->d_offsets, &d_u, &d_v}; 308 309 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyUnorientedTranspose, grid, block_size, args)); 310 } else { 311 void *args[] = {&impl->d_l_vec_indices, &impl->d_t_indices, &impl->d_t_offsets, &d_u, &d_v}; 312 313 CeedCallBackend(CeedRunKernel_Cuda(ceed, impl->ApplyUnorientedTranspose, grid, block_size, args)); 314 } 315 } 316 } break; 317 } 318 } 319 320 if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED) *request = NULL; 321 322 // Restore arrays 323 CeedCallBackend(CeedVectorRestoreArrayRead(u, &d_u)); 324 CeedCallBackend(CeedVectorRestoreArray(v, &d_v)); 325 return CEED_ERROR_SUCCESS; 326 } 327 328 //------------------------------------------------------------------------------ 329 // Apply restriction 330 //------------------------------------------------------------------------------ 331 static int CeedElemRestrictionApply_Cuda(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector v, CeedRequest *request) { 332 return CeedElemRestrictionApply_Cuda_Core(rstr, t_mode, true, true, u, v, request); 333 } 334 335 //------------------------------------------------------------------------------ 336 // Apply unsigned restriction 337 //------------------------------------------------------------------------------ 338 static int CeedElemRestrictionApplyUnsigned_Cuda(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 339 CeedRequest *request) { 340 return CeedElemRestrictionApply_Cuda_Core(rstr, t_mode, false, true, u, v, request); 341 } 342 343 //------------------------------------------------------------------------------ 344 // Apply unoriented restriction 345 //------------------------------------------------------------------------------ 346 static int CeedElemRestrictionApplyUnoriented_Cuda(CeedElemRestriction rstr, CeedTransposeMode t_mode, CeedVector u, CeedVector v, 347 CeedRequest *request) { 348 return CeedElemRestrictionApply_Cuda_Core(rstr, t_mode, false, false, u, v, request); 349 } 350 351 //------------------------------------------------------------------------------ 352 // Get offsets 353 //------------------------------------------------------------------------------ 354 static int CeedElemRestrictionGetOffsets_Cuda(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt **offsets) { 355 CeedElemRestriction_Cuda *impl; 356 CeedRestrictionType rstr_type; 357 358 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 359 CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type)); 360 switch (mem_type) { 361 case CEED_MEM_HOST: 362 *offsets = rstr_type == CEED_RESTRICTION_POINTS ? impl->h_offsets_at_points : impl->h_offsets; 363 break; 364 case CEED_MEM_DEVICE: 365 *offsets = rstr_type == CEED_RESTRICTION_POINTS ? impl->d_offsets_at_points : impl->d_offsets; 366 break; 367 } 368 return CEED_ERROR_SUCCESS; 369 } 370 371 //------------------------------------------------------------------------------ 372 // Get orientations 373 //------------------------------------------------------------------------------ 374 static int CeedElemRestrictionGetOrientations_Cuda(CeedElemRestriction rstr, CeedMemType mem_type, const bool **orients) { 375 CeedElemRestriction_Cuda *impl; 376 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 377 378 switch (mem_type) { 379 case CEED_MEM_HOST: 380 *orients = impl->h_orients; 381 break; 382 case CEED_MEM_DEVICE: 383 *orients = impl->d_orients; 384 break; 385 } 386 return CEED_ERROR_SUCCESS; 387 } 388 389 //------------------------------------------------------------------------------ 390 // Get curl-conforming orientations 391 //------------------------------------------------------------------------------ 392 static int CeedElemRestrictionGetCurlOrientations_Cuda(CeedElemRestriction rstr, CeedMemType mem_type, const CeedInt8 **curl_orients) { 393 CeedElemRestriction_Cuda *impl; 394 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 395 396 switch (mem_type) { 397 case CEED_MEM_HOST: 398 *curl_orients = impl->h_curl_orients; 399 break; 400 case CEED_MEM_DEVICE: 401 *curl_orients = impl->d_curl_orients; 402 break; 403 } 404 return CEED_ERROR_SUCCESS; 405 } 406 407 //------------------------------------------------------------------------------ 408 // Get offset for padded AtPoints E-layout 409 //------------------------------------------------------------------------------ 410 static int CeedElemRestrictionGetAtPointsElementOffset_Cuda(CeedElemRestriction rstr, CeedInt elem, CeedSize *elem_offset) { 411 CeedInt layout[3]; 412 413 CeedCallBackend(CeedElemRestrictionGetELayout(rstr, layout)); 414 *elem_offset = 0 * layout[0] + 0 * layout[1] + elem * layout[2]; 415 return CEED_ERROR_SUCCESS; 416 } 417 418 //------------------------------------------------------------------------------ 419 // Destroy restriction 420 //------------------------------------------------------------------------------ 421 static int CeedElemRestrictionDestroy_Cuda(CeedElemRestriction rstr) { 422 Ceed ceed; 423 CeedElemRestriction_Cuda *impl; 424 425 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 426 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 427 if (impl->module) { 428 CeedCallCuda(ceed, cuModuleUnload(impl->module)); 429 } 430 CeedCallBackend(CeedFree(&impl->h_offsets_owned)); 431 CeedCallCuda(ceed, cudaFree((CeedInt *)impl->d_offsets_owned)); 432 CeedCallCuda(ceed, cudaFree((CeedInt *)impl->d_t_offsets)); 433 CeedCallCuda(ceed, cudaFree((CeedInt *)impl->d_t_indices)); 434 CeedCallCuda(ceed, cudaFree((CeedInt *)impl->d_l_vec_indices)); 435 CeedCallBackend(CeedFree(&impl->h_orients_owned)); 436 CeedCallCuda(ceed, cudaFree((bool *)impl->d_orients_owned)); 437 CeedCallBackend(CeedFree(&impl->h_curl_orients_owned)); 438 CeedCallCuda(ceed, cudaFree((CeedInt8 *)impl->d_curl_orients_owned)); 439 CeedCallBackend(CeedFree(&impl->h_offsets_at_points_owned)); 440 CeedCallCuda(ceed, cudaFree((CeedInt *)impl->d_offsets_at_points_owned)); 441 CeedCallBackend(CeedFree(&impl->h_points_per_elem_owned)); 442 CeedCallCuda(ceed, cudaFree((CeedInt *)impl->d_points_per_elem_owned)); 443 CeedCallBackend(CeedFree(&impl)); 444 return CEED_ERROR_SUCCESS; 445 } 446 447 //------------------------------------------------------------------------------ 448 // Create transpose offsets and indices 449 //------------------------------------------------------------------------------ 450 static int CeedElemRestrictionOffset_Cuda(const CeedElemRestriction rstr, const CeedInt elem_size, const CeedInt *indices) { 451 Ceed ceed; 452 bool *is_node; 453 CeedSize l_size; 454 CeedInt num_elem, num_comp, num_nodes = 0; 455 CeedInt *ind_to_offset, *l_vec_indices, *t_offsets, *t_indices; 456 CeedRestrictionType rstr_type; 457 CeedElemRestriction_Cuda *impl; 458 459 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 460 CeedCallBackend(CeedElemRestrictionGetData(rstr, &impl)); 461 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 462 CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type)); 463 CeedCallBackend(CeedElemRestrictionGetLVectorSize(rstr, &l_size)); 464 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 465 const CeedInt size_indices = num_elem * elem_size; 466 467 // Count num_nodes 468 CeedCallBackend(CeedCalloc(l_size, &is_node)); 469 470 for (CeedInt i = 0; i < size_indices; i++) is_node[indices[i]] = 1; 471 for (CeedInt i = 0; i < l_size; i++) num_nodes += is_node[i]; 472 impl->num_nodes = num_nodes; 473 474 // L-vector offsets array 475 CeedCallBackend(CeedCalloc(l_size, &ind_to_offset)); 476 CeedCallBackend(CeedCalloc(num_nodes, &l_vec_indices)); 477 for (CeedInt i = 0, j = 0; i < l_size; i++) { 478 if (is_node[i]) { 479 l_vec_indices[j] = i; 480 ind_to_offset[i] = j++; 481 } 482 } 483 CeedCallBackend(CeedFree(&is_node)); 484 485 // Compute transpose offsets and indices 486 const CeedInt size_offsets = num_nodes + 1; 487 488 CeedCallBackend(CeedCalloc(size_offsets, &t_offsets)); 489 CeedCallBackend(CeedMalloc(size_indices, &t_indices)); 490 // Count node multiplicity 491 for (CeedInt e = 0; e < num_elem; ++e) { 492 for (CeedInt i = 0; i < elem_size; ++i) ++t_offsets[ind_to_offset[indices[elem_size * e + i]] + 1]; 493 } 494 // Convert to running sum 495 for (CeedInt i = 1; i < size_offsets; ++i) t_offsets[i] += t_offsets[i - 1]; 496 // List all E-vec indices associated with L-vec node 497 for (CeedInt e = 0; e < num_elem; ++e) { 498 for (CeedInt i = 0; i < elem_size; ++i) { 499 const CeedInt lid = elem_size * e + i; 500 const CeedInt gid = indices[lid]; 501 502 t_indices[t_offsets[ind_to_offset[gid]]++] = lid; 503 } 504 } 505 // Reset running sum 506 for (int i = size_offsets - 1; i > 0; --i) t_offsets[i] = t_offsets[i - 1]; 507 t_offsets[0] = 0; 508 509 // Copy data to device 510 // -- L-vector indices 511 CeedCallCuda(ceed, cudaMalloc((void **)&impl->d_l_vec_indices, num_nodes * sizeof(CeedInt))); 512 CeedCallCuda(ceed, cudaMemcpy((CeedInt *)impl->d_l_vec_indices, l_vec_indices, num_nodes * sizeof(CeedInt), cudaMemcpyHostToDevice)); 513 // -- Transpose offsets 514 CeedCallCuda(ceed, cudaMalloc((void **)&impl->d_t_offsets, size_offsets * sizeof(CeedInt))); 515 CeedCallCuda(ceed, cudaMemcpy((CeedInt *)impl->d_t_offsets, t_offsets, size_offsets * sizeof(CeedInt), cudaMemcpyHostToDevice)); 516 // -- Transpose indices 517 CeedCallCuda(ceed, cudaMalloc((void **)&impl->d_t_indices, size_indices * sizeof(CeedInt))); 518 CeedCallCuda(ceed, cudaMemcpy((CeedInt *)impl->d_t_indices, t_indices, size_indices * sizeof(CeedInt), cudaMemcpyHostToDevice)); 519 520 // Cleanup 521 CeedCallBackend(CeedFree(&ind_to_offset)); 522 CeedCallBackend(CeedFree(&l_vec_indices)); 523 CeedCallBackend(CeedFree(&t_offsets)); 524 CeedCallBackend(CeedFree(&t_indices)); 525 return CEED_ERROR_SUCCESS; 526 } 527 528 //------------------------------------------------------------------------------ 529 // Create restriction 530 //------------------------------------------------------------------------------ 531 int CeedElemRestrictionCreate_Cuda(CeedMemType mem_type, CeedCopyMode copy_mode, const CeedInt *offsets, const bool *orients, 532 const CeedInt8 *curl_orients, CeedElemRestriction rstr) { 533 Ceed ceed, ceed_parent; 534 bool is_deterministic; 535 CeedInt num_elem, num_comp, elem_size; 536 CeedRestrictionType rstr_type; 537 CeedElemRestriction_Cuda *impl; 538 539 CeedCallBackend(CeedElemRestrictionGetCeed(rstr, &ceed)); 540 CeedCallBackend(CeedGetParent(ceed, &ceed_parent)); 541 CeedCallBackend(CeedIsDeterministic(ceed_parent, &is_deterministic)); 542 CeedCallBackend(CeedElemRestrictionGetNumElements(rstr, &num_elem)); 543 CeedCallBackend(CeedElemRestrictionGetNumComponents(rstr, &num_comp)); 544 CeedCallBackend(CeedElemRestrictionGetElementSize(rstr, &elem_size)); 545 CeedCallBackend(CeedElemRestrictionGetType(rstr, &rstr_type)); 546 // Use max number of points as elem size for AtPoints restrictions 547 if (rstr_type == CEED_RESTRICTION_POINTS) { 548 CeedInt max_points = 0; 549 550 for (CeedInt i = 0; i < num_elem; i++) { 551 max_points = CeedIntMax(max_points, offsets[i + 1] - offsets[i]); 552 } 553 elem_size = max_points; 554 } 555 const CeedInt size = num_elem * elem_size; 556 557 CeedCallBackend(CeedCalloc(1, &impl)); 558 impl->num_nodes = size; 559 CeedCallBackend(CeedElemRestrictionSetData(rstr, impl)); 560 561 // Set layouts 562 { 563 bool has_backend_strides; 564 CeedInt layout[3] = {1, size, elem_size}; 565 566 CeedCallBackend(CeedElemRestrictionSetELayout(rstr, layout)); 567 if (rstr_type == CEED_RESTRICTION_STRIDED) { 568 CeedCallBackend(CeedElemRestrictionHasBackendStrides(rstr, &has_backend_strides)); 569 if (has_backend_strides) { 570 CeedCallBackend(CeedElemRestrictionSetLLayout(rstr, layout)); 571 } 572 } 573 } 574 575 // Pad AtPoints indices 576 if (rstr_type == CEED_RESTRICTION_POINTS) { 577 CeedSize offsets_len = elem_size * num_elem, at_points_size = num_elem + 1; 578 CeedInt max_points = elem_size, *offsets_padded, *points_per_elem; 579 580 CeedCheck(mem_type == CEED_MEM_HOST, ceed, CEED_ERROR_BACKEND, "only MemType Host supported when creating AtPoints restriction"); 581 CeedCallBackend(CeedMalloc(offsets_len, &offsets_padded)); 582 CeedCallBackend(CeedMalloc(num_elem, &points_per_elem)); 583 for (CeedInt i = 0; i < num_elem; i++) { 584 CeedInt num_points = offsets[i + 1] - offsets[i]; 585 586 points_per_elem[i] = num_points; 587 at_points_size += num_points; 588 // -- Copy all points in element 589 for (CeedInt j = 0; j < num_points; j++) { 590 offsets_padded[i * max_points + j] = offsets[offsets[i] + j] * num_comp; 591 } 592 // -- Replicate out last point in element 593 for (CeedInt j = num_points; j < max_points; j++) { 594 offsets_padded[i * max_points + j] = offsets[offsets[i] + num_points - 1] * num_comp; 595 } 596 } 597 CeedCallBackend(CeedSetHostCeedIntArray(offsets, copy_mode, at_points_size, &impl->h_offsets_at_points_owned, &impl->h_offsets_at_points_borrowed, 598 &impl->h_offsets_at_points)); 599 CeedCallCuda(ceed, cudaMalloc((void **)&impl->d_offsets_at_points_owned, at_points_size * sizeof(CeedInt))); 600 CeedCallCuda(ceed, cudaMemcpy((CeedInt **)impl->d_offsets_at_points_owned, impl->h_offsets_at_points, at_points_size * sizeof(CeedInt), 601 cudaMemcpyHostToDevice)); 602 impl->d_offsets_at_points = (CeedInt *)impl->d_offsets_at_points_owned; 603 604 // -- Use padded offsets for the rest of the setup 605 offsets = (const CeedInt *)offsets_padded; 606 copy_mode = CEED_OWN_POINTER; 607 CeedCallBackend(CeedElemRestrictionSetAtPointsEVectorSize(rstr, at_points_size * num_comp)); 608 609 // -- Points per element 610 CeedCallBackend(CeedSetHostCeedIntArray(points_per_elem, CEED_OWN_POINTER, num_elem, &impl->h_points_per_elem_owned, 611 &impl->h_points_per_elem_borrowed, &impl->h_points_per_elem)); 612 CeedCallCuda(ceed, cudaMalloc((void **)&impl->d_points_per_elem_owned, num_elem * sizeof(CeedInt))); 613 CeedCallCuda(ceed, 614 cudaMemcpy((CeedInt **)impl->d_points_per_elem_owned, impl->h_points_per_elem, num_elem * sizeof(CeedInt), cudaMemcpyHostToDevice)); 615 impl->d_points_per_elem = (CeedInt *)impl->d_points_per_elem_owned; 616 } 617 618 // Set up device offset/orientation arrays 619 if (rstr_type != CEED_RESTRICTION_STRIDED) { 620 switch (mem_type) { 621 case CEED_MEM_HOST: { 622 CeedCallBackend(CeedSetHostCeedIntArray(offsets, copy_mode, size, &impl->h_offsets_owned, &impl->h_offsets_borrowed, &impl->h_offsets)); 623 CeedCallCuda(ceed, cudaMalloc((void **)&impl->d_offsets_owned, size * sizeof(CeedInt))); 624 CeedCallCuda(ceed, cudaMemcpy((CeedInt *)impl->d_offsets_owned, impl->h_offsets, size * sizeof(CeedInt), cudaMemcpyHostToDevice)); 625 impl->d_offsets = (CeedInt *)impl->d_offsets_owned; 626 if (is_deterministic) CeedCallBackend(CeedElemRestrictionOffset_Cuda(rstr, elem_size, offsets)); 627 } break; 628 case CEED_MEM_DEVICE: { 629 CeedCallBackend(CeedSetDeviceCeedIntArray_Cuda(ceed, offsets, copy_mode, size, &impl->d_offsets_owned, &impl->d_offsets_borrowed, 630 (const CeedInt **)&impl->d_offsets)); 631 CeedCallBackend(CeedMalloc(size, &impl->h_offsets_owned)); 632 CeedCallCuda(ceed, cudaMemcpy((CeedInt *)impl->h_offsets_owned, impl->d_offsets, size * sizeof(CeedInt), cudaMemcpyDeviceToHost)); 633 impl->h_offsets = impl->h_offsets_owned; 634 if (is_deterministic) CeedCallBackend(CeedElemRestrictionOffset_Cuda(rstr, elem_size, offsets)); 635 } break; 636 } 637 638 // Orientation data 639 if (rstr_type == CEED_RESTRICTION_ORIENTED) { 640 switch (mem_type) { 641 case CEED_MEM_HOST: { 642 CeedCallBackend(CeedSetHostBoolArray(orients, copy_mode, size, &impl->h_orients_owned, &impl->h_orients_borrowed, &impl->h_orients)); 643 CeedCallCuda(ceed, cudaMalloc((void **)&impl->d_orients_owned, size * sizeof(bool))); 644 CeedCallCuda(ceed, cudaMemcpy((bool *)impl->d_orients_owned, impl->h_orients, size * sizeof(bool), cudaMemcpyHostToDevice)); 645 impl->d_orients = impl->d_orients_owned; 646 } break; 647 case CEED_MEM_DEVICE: { 648 CeedCallBackend(CeedSetDeviceBoolArray_Cuda(ceed, orients, copy_mode, size, &impl->d_orients_owned, &impl->d_orients_borrowed, 649 (const bool **)&impl->d_orients)); 650 CeedCallBackend(CeedMalloc(size, &impl->h_orients_owned)); 651 CeedCallCuda(ceed, cudaMemcpy((bool *)impl->h_orients_owned, impl->d_orients, size * sizeof(bool), cudaMemcpyDeviceToHost)); 652 impl->h_orients = impl->h_orients_owned; 653 } break; 654 } 655 } else if (rstr_type == CEED_RESTRICTION_CURL_ORIENTED) { 656 switch (mem_type) { 657 case CEED_MEM_HOST: { 658 CeedCallBackend(CeedSetHostCeedInt8Array(curl_orients, copy_mode, 3 * size, &impl->h_curl_orients_owned, &impl->h_curl_orients_borrowed, 659 &impl->h_curl_orients)); 660 CeedCallCuda(ceed, cudaMalloc((void **)&impl->d_curl_orients_owned, 3 * size * sizeof(CeedInt8))); 661 CeedCallCuda(ceed, 662 cudaMemcpy((CeedInt8 *)impl->d_curl_orients_owned, impl->h_curl_orients, 3 * size * sizeof(CeedInt8), cudaMemcpyHostToDevice)); 663 impl->d_curl_orients = impl->d_curl_orients_owned; 664 } break; 665 case CEED_MEM_DEVICE: { 666 CeedCallBackend(CeedSetDeviceCeedInt8Array_Cuda(ceed, curl_orients, copy_mode, 3 * size, &impl->d_curl_orients_owned, 667 &impl->d_curl_orients_borrowed, (const CeedInt8 **)&impl->d_curl_orients)); 668 CeedCallBackend(CeedMalloc(3 * size, &impl->h_curl_orients_owned)); 669 CeedCallCuda(ceed, 670 cudaMemcpy((CeedInt8 *)impl->h_curl_orients_owned, impl->d_curl_orients, 3 * size * sizeof(CeedInt8), cudaMemcpyDeviceToHost)); 671 impl->h_curl_orients = impl->h_curl_orients_owned; 672 } break; 673 } 674 } 675 } 676 677 // Register backend functions 678 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "Apply", CeedElemRestrictionApply_Cuda)); 679 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyUnsigned", CeedElemRestrictionApplyUnsigned_Cuda)); 680 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "ApplyUnoriented", CeedElemRestrictionApplyUnoriented_Cuda)); 681 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetOffsets", CeedElemRestrictionGetOffsets_Cuda)); 682 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetOrientations", CeedElemRestrictionGetOrientations_Cuda)); 683 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetCurlOrientations", CeedElemRestrictionGetCurlOrientations_Cuda)); 684 if (rstr_type == CEED_RESTRICTION_POINTS) { 685 CeedCallBackend( 686 CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "GetAtPointsElementOffset", CeedElemRestrictionGetAtPointsElementOffset_Cuda)); 687 } 688 CeedCallBackend(CeedSetBackendFunction(ceed, "ElemRestriction", rstr, "Destroy", CeedElemRestrictionDestroy_Cuda)); 689 return CEED_ERROR_SUCCESS; 690 } 691 692 //------------------------------------------------------------------------------ 693