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