1 // Copyright (c) 2017-2026, 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 MAGMA tensor basis interpolation in 1D
10 #include "magma-common-tensor.h"
11
12 // macros to abstract access of shared memory and reg. file
13 #define sT(i, j) sT[(j) * P + (i)]
14 #define sTmp(i, j, ldw) sTmp[(j) * (ldw) + (i)]
15
16 ////////////////////////////////////////////////////////////////////////////////
17 // interp basis action (2D)
18 template <typename T, int DIM_U, int DIM_V, int NUM_COMP, int P, int Q, int rU_SIZE, int rV_SIZE>
magma_interp_2d_device(const T * sT,T rU[DIM_U][NUM_COMP][rU_SIZE],T rV[DIM_V][NUM_COMP][rV_SIZE],const int tx,T rTmp,T * swork)19 static __device__ __inline__ void magma_interp_2d_device(const T *sT, T rU[DIM_U][NUM_COMP][rU_SIZE], T rV[DIM_V][NUM_COMP][rV_SIZE], const int tx,
20 T rTmp, T *swork) {
21 // Assumptions
22 // 1. 1D threads of size max(P,Q)
23 // 2. input: rU[DIM_U x NUM_COMP x rU_SIZE] in registers (per thread)
24 // 3. output: rV[DIM_V x NUM_COMP x rV_SIZE] in registers (per thread)
25 // 4. Two products per component
26 // 4.1 Batch P of (1xP) matrices times (PxQ) matrix => Batch P of (1xQ) matrices
27 // 4.2 Batch 1 of (QxP) matrix times (PxQ) matrix => (QxQ) matrix
28 // 5. Each thread computes one row of the output of each product
29 // 6. Sync is recommended before and after the call
30
31 for (int comp = 0; comp < NUM_COMP; comp++) {
32 // 1st product -- Batch P of (1xP) matrices [reg] x (PxQ) [shmem] => Batch P of (1xQ) matrices
33 // the batch output P x (1xQ) is written on the fly to shmem
34 if (tx < P) {
35 const int batchid = tx;
36 const int sld = 1;
37 T *sTmp = swork + batchid * (1 * Q);
38 for (int j = 0; j < Q; j++) {
39 rTmp = 0.0;
40 for (int i = 0; i < P; i++) {
41 rTmp += rU[0][comp][i] * sT(i, j);
42 }
43 sTmp(0, j, sld) = rTmp;
44 }
45 } // end of: if (tx < P)
46 __syncthreads();
47
48 // 2nd product -- Batch 1 of a (QxP) matrix [shmem] x (PxQ) [shmem] => (QxQ) matrix [reg]
49 if (tx < Q) {
50 const int batchid = 0;
51 const int sld = Q;
52 T *sTmp = swork + batchid * (Q * P);
53 for (int j = 0; j < Q; j++) {
54 rTmp = 0.0;
55 for (int i = 0; i < P; i++) {
56 rTmp += sTmp(tx, i, sld) * sT(i, j);
57 }
58 rV[0][comp][j] += rTmp;
59 }
60 }
61 __syncthreads();
62 }
63 }
64
65 ////////////////////////////////////////////////////////////////////////////////
__launch_bounds__(MAGMA_BASIS_BOUNDS (BASIS_MAX_P_Q,MAGMA_MAXTHREADS_2D))66 extern "C" __launch_bounds__(MAGMA_BASIS_BOUNDS(BASIS_MAX_P_Q, MAGMA_MAXTHREADS_2D)) __global__
67 void magma_interpn_2d_kernel(const CeedScalar *dT, const CeedScalar *dU, const int estrdU, const int cstrdU, CeedScalar *dV, const int estrdV,
68 const int cstrdV, const int nelem) {
69 MAGMA_DEVICE_SHARED(CeedScalar, shared_data)
70
71 const int tx = threadIdx.x;
72 const int ty = threadIdx.y;
73 const int elem_id = (blockIdx.x * blockDim.y) + ty;
74
75 if (elem_id >= nelem) return;
76
77 CeedScalar rU[1][BASIS_NUM_COMP][BASIS_P] = {0.0}; // for a non-fused operator BASIS_DIM is always 1
78 CeedScalar rV[1][BASIS_NUM_COMP][BASIS_Q] = {0.0}; // for a non-fused operator BASIS_DIM is always 1
79 CeedScalar rTmp = 0.0;
80
81 // shift global memory pointers by elem stride
82 dU += elem_id * estrdU;
83 dV += elem_id * estrdV;
84
85 // assign shared memory pointers
86 CeedScalar *sT = (CeedScalar *)shared_data;
87 CeedScalar *sTmp = sT + BASIS_P * BASIS_Q;
88 sTmp += ty * (BASIS_P * BASIS_MAX_P_Q);
89
90 // read T
91 if (ty == 0) {
92 read_T_notrans_gm2sm<BASIS_P, BASIS_Q>(tx, dT, sT);
93 }
94
95 // read U -- there is a sync at the end of this function
96 read_U_2d<CeedScalar, BASIS_P, 1, BASIS_NUM_COMP, BASIS_P, 0>(dU, cstrdU, rU, sTmp, tx);
97
98 // no sync needed here -- read_U_2d already syncs at the end
99 magma_interp_2d_device<CeedScalar, 1, 1, BASIS_NUM_COMP, BASIS_P, BASIS_Q, BASIS_P, BASIS_Q>(sT, rU, rV, tx, rTmp, sTmp);
100 __syncthreads();
101
102 // write V
103 write_V_2d<CeedScalar, BASIS_Q, 1, BASIS_NUM_COMP, BASIS_Q, 0>(dV, cstrdV, rV, tx);
104 }
105
106 ////////////////////////////////////////////////////////////////////////////////
__launch_bounds__(MAGMA_BASIS_BOUNDS (BASIS_MAX_P_Q,MAGMA_MAXTHREADS_2D))107 extern "C" __launch_bounds__(MAGMA_BASIS_BOUNDS(BASIS_MAX_P_Q, MAGMA_MAXTHREADS_2D)) __global__
108 void magma_interpt_2d_kernel(const CeedScalar *dT, const CeedScalar *dU, const int estrdU, const int cstrdU, CeedScalar *dV, const int estrdV,
109 const int cstrdV, const int nelem) {
110 MAGMA_DEVICE_SHARED(CeedScalar, shared_data)
111
112 const int tx = threadIdx.x;
113 const int ty = threadIdx.y;
114 const int elem_id = (blockIdx.x * blockDim.y) + ty;
115
116 if (elem_id >= nelem) return;
117
118 CeedScalar rU[1][BASIS_NUM_COMP][BASIS_Q] = {0.0}; // for a non-fused operator BASIS_DIM is always 1
119 CeedScalar rV[1][BASIS_NUM_COMP][BASIS_P] = {0.0}; // for a non-fused operator BASIS_DIM is always 1
120 CeedScalar rTmp = 0.0;
121
122 // shift global memory pointers by elem stride
123 dU += elem_id * estrdU;
124 dV += elem_id * estrdV;
125
126 // assign shared memory pointers
127 CeedScalar *sT = (CeedScalar *)shared_data;
128 CeedScalar *sTmp = sT + BASIS_Q * BASIS_P;
129 sTmp += ty * (BASIS_Q * BASIS_MAX_P_Q);
130
131 // read T
132 if (ty == 0) {
133 read_T_trans_gm2sm<BASIS_Q, BASIS_P>(tx, dT, sT);
134 }
135
136 // read U -- there is a sync at the end of this function
137 read_U_2d<CeedScalar, BASIS_Q, 1, BASIS_NUM_COMP, BASIS_Q, 0>(dU, cstrdU, rU, sTmp, tx);
138
139 // no sync needed here -- read_U_2d already syncs at the end
140 magma_interp_2d_device<CeedScalar, 1, 1, BASIS_NUM_COMP, BASIS_Q, BASIS_P, BASIS_Q, BASIS_P>(sT, rU, rV, tx, rTmp, sTmp);
141 __syncthreads();
142
143 // write V
144 write_V_2d<CeedScalar, BASIS_P, 1, BASIS_NUM_COMP, BASIS_P, 0>(dV, cstrdV, rV, tx);
145 }
146
147 ////////////////////////////////////////////////////////////////////////////////
__launch_bounds__(MAGMA_BASIS_BOUNDS (BASIS_MAX_P_Q,MAGMA_MAXTHREADS_2D))148 extern "C" __launch_bounds__(MAGMA_BASIS_BOUNDS(BASIS_MAX_P_Q, MAGMA_MAXTHREADS_2D)) __global__
149 void magma_interpta_2d_kernel(const CeedScalar *dT, const CeedScalar *dU, const int estrdU, const int cstrdU, CeedScalar *dV, const int estrdV,
150 const int cstrdV, const int nelem) {
151 MAGMA_DEVICE_SHARED(CeedScalar, shared_data)
152
153 const int tx = threadIdx.x;
154 const int ty = threadIdx.y;
155 const int elem_id = (blockIdx.x * blockDim.y) + ty;
156
157 if (elem_id >= nelem) return;
158
159 CeedScalar rU[1][BASIS_NUM_COMP][BASIS_Q] = {0.0}; // for a non-fused operator BASIS_DIM is always 1
160 CeedScalar rV[1][BASIS_NUM_COMP][BASIS_P] = {0.0}; // for a non-fused operator BASIS_DIM is always 1
161 CeedScalar rTmp = 0.0;
162
163 // shift global memory pointers by elem stride
164 dU += elem_id * estrdU;
165 dV += elem_id * estrdV;
166
167 // assign shared memory pointers
168 CeedScalar *sT = (CeedScalar *)shared_data;
169 CeedScalar *sTmp = sT + BASIS_Q * BASIS_P;
170 sTmp += ty * (BASIS_Q * BASIS_MAX_P_Q);
171
172 // read T
173 if (ty == 0) {
174 read_T_trans_gm2sm<BASIS_Q, BASIS_P>(tx, dT, sT);
175 }
176
177 // read U -- there is a sync at the end of this function
178 read_U_2d<CeedScalar, BASIS_Q, 1, BASIS_NUM_COMP, BASIS_Q, 0>(dU, cstrdU, rU, sTmp, tx);
179
180 // no sync needed here -- read_U_2d already syncs at the end
181 magma_interp_2d_device<CeedScalar, 1, 1, BASIS_NUM_COMP, BASIS_Q, BASIS_P, BASIS_Q, BASIS_P>(sT, rU, rV, tx, rTmp, sTmp);
182 __syncthreads();
183
184 // sum into V
185 sum_V_2d<CeedScalar, BASIS_P, 1, BASIS_NUM_COMP, BASIS_P, 0>(dV, cstrdV, rV, tx);
186 }
187