xref: /libCEED/include/ceed/jit-source/cuda/cuda-ref-restriction-at-points.h (revision 4aa5fbb4c55e771431ba815a9dddd623e4e3362a)
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 /// @file
9 /// Internal header for CUDA offset element restriction kernels
10 
11 #include <ceed.h>
12 
13 //------------------------------------------------------------------------------
14 // E-vector -> L-vector, standard (with offsets)
15 //------------------------------------------------------------------------------
16 #if !USE_DETERMINISTIC
17 extern "C" __global__ void AtPointsTranspose(const CeedInt *__restrict__ indices, const CeedInt *__restrict__ points_per_elem,
18                                              const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v) {
19   for (CeedInt node = blockIdx.x * blockDim.x + threadIdx.x; node < RSTR_NUM_ELEM * RSTR_ELEM_SIZE; node += blockDim.x * gridDim.x) {
20     const CeedInt ind      = indices[node];
21     const CeedInt loc_node = node % RSTR_ELEM_SIZE;
22     const CeedInt elem     = node / RSTR_ELEM_SIZE;
23 
24     if (loc_node >= points_per_elem[elem]) continue;
25     for (CeedInt comp = 0; comp < RSTR_NUM_COMP; comp++) {
26       atomicAdd(v + ind + comp * RSTR_COMP_STRIDE, u[loc_node + comp * RSTR_ELEM_SIZE * RSTR_NUM_ELEM + elem * RSTR_ELEM_SIZE]);
27     }
28   }
29 }
30 #else
31 extern "C" __global__ void AtPointsTranspose(const CeedInt *__restrict__ l_vec_indices, const CeedInt *__restrict__ t_indices,
32                                              const CeedInt *__restrict__ points_per_elem, const CeedInt *__restrict__ t_offsets,
33                                              const CeedScalar *__restrict__ u, CeedScalar *__restrict__ v) {
34   CeedScalar value[RSTR_NUM_COMP];
35 
36   for (CeedInt i = blockIdx.x * blockDim.x + threadIdx.x; i < RSTR_NUM_NODES; i += blockDim.x * gridDim.x) {
37     const CeedInt ind     = l_vec_indices[i];
38     const CeedInt range_1 = t_offsets[i];
39     const CeedInt range_N = t_offsets[i + 1];
40 
41     for (CeedInt comp = 0; comp < RSTR_NUM_COMP; comp++) value[comp] = 0.0;
42 
43     for (CeedInt j = range_1; j < range_N; j++) {
44       const CeedInt t_ind    = t_indices[j];
45       const CeedInt loc_node = t_ind % RSTR_ELEM_SIZE;
46       const CeedInt elem     = t_ind / RSTR_ELEM_SIZE;
47 
48       if (loc_node >= points_per_elem[elem]) continue;
49       for (CeedInt comp = 0; comp < RSTR_NUM_COMP; comp++) {
50         value[comp] += u[loc_node + comp * RSTR_ELEM_SIZE * RSTR_NUM_ELEM + elem * RSTR_ELEM_SIZE];
51       }
52     }
53 
54     for (CeedInt comp = 0; comp < RSTR_NUM_COMP; comp++) v[ind + comp * RSTR_COMP_STRIDE] += value[comp];
55   }
56 }
57 #endif
58