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-ref.h" 18 19 static inline int CeedElemRestrictionApply_Ref_Core(CeedElemRestriction r, 20 const CeedInt blksize, const CeedInt ncomp, CeedInt start, CeedInt stop, 21 CeedTransposeMode tmode, CeedTransposeMode lmode, CeedVector u, 22 CeedVector v, CeedRequest *request) { 23 int ierr; 24 CeedElemRestriction_Ref *impl; 25 ierr = CeedElemRestrictionGetData(r, (void *)&impl); CeedChk(ierr);; 26 const CeedScalar *uu; 27 CeedScalar *vv; 28 CeedInt nelem, elemsize, nnodes, voffset; 29 ierr = CeedElemRestrictionGetNumElements(r, &nelem); CeedChk(ierr); 30 ierr = CeedElemRestrictionGetElementSize(r, &elemsize); CeedChk(ierr); 31 ierr = CeedElemRestrictionGetNumNodes(r, &nnodes); CeedChk(ierr); 32 voffset = start*blksize*elemsize*ncomp; 33 34 ierr = CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu); CeedChk(ierr); 35 ierr = CeedVectorGetArray(v, CEED_MEM_HOST, &vv); CeedChk(ierr); 36 // Restriction from lvector to evector 37 // Perform: v = r * u 38 if (tmode == CEED_NOTRANSPOSE) { 39 // No indices provided, Identity Restriction 40 if (!impl->indices) { 41 for (CeedInt e = start*blksize; e < stop*blksize; e+=blksize) 42 for (CeedInt j = 0; j < blksize; j++) 43 for (CeedInt k = 0; k < ncomp*elemsize; k++) 44 vv[e*elemsize*ncomp + k*blksize + j - voffset] 45 = uu[CeedIntMin(e+j,nelem-1)*ncomp*elemsize + k]; 46 } else { 47 // Indices provided, standard or blocked restriction 48 // vv has shape [elemsize, ncomp, nelem], row-major 49 // uu has shape [nnodes, ncomp] 50 for (CeedInt e = start*blksize; e < stop*blksize; e+=blksize) 51 for (CeedInt d = 0; d < ncomp; d++) 52 for (CeedInt i = 0; i < elemsize*blksize; i++) 53 vv[i+elemsize*(d*blksize+ncomp*e) - voffset] 54 = uu[lmode == CEED_NOTRANSPOSE 55 ? impl->indices[i+elemsize*e]+nnodes*d 56 : d+ncomp*impl->indices[i+elemsize*e]]; 57 } 58 } else { 59 // Restriction from evector to lvector 60 // Performing v += r^T * u 61 // No indices provided, Identity Restriction 62 if (!impl->indices) { 63 for (CeedInt e = start*blksize; e < stop*blksize; e+=blksize) 64 for (CeedInt j = 0; j < CeedIntMin(blksize, nelem-e); j++) 65 for (CeedInt k = 0; k < ncomp*elemsize; k++) 66 vv[(e+j)*ncomp*elemsize + k] 67 += uu[e*elemsize*ncomp + k*blksize + j - voffset]; 68 } else { 69 // Indices provided, standard or blocked restriction 70 // uu has shape [elemsize, ncomp, nelem] 71 // vv has shape [nnodes, ncomp] 72 for (CeedInt e = start*blksize; e < stop*blksize; e+=blksize) { 73 for (CeedInt d = 0; d < ncomp; d++) 74 for (CeedInt i = 0; i < elemsize*blksize; i+=blksize) 75 // Iteration bound set to discard padding elements 76 for (CeedInt j = i; j < i+CeedIntMin(blksize, nelem-e); j++) 77 vv[lmode == CEED_NOTRANSPOSE 78 ? impl->indices[j+e*elemsize]+nnodes*d 79 : d+ncomp*impl->indices[j+e*elemsize]] 80 += uu[j+elemsize*(d*blksize+ncomp*e) - voffset]; 81 } 82 } 83 } 84 ierr = CeedVectorRestoreArrayRead(u, &uu); CeedChk(ierr); 85 ierr = CeedVectorRestoreArray(v, &vv); CeedChk(ierr); 86 if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED) 87 *request = NULL; 88 return 0; 89 } 90 91 static int CeedElemRestrictionApply_Ref_Core_11(CeedElemRestriction r, 92 CeedInt start, CeedInt stop, CeedTransposeMode tmode, 93 CeedTransposeMode lmode, CeedVector u, CeedVector v, CeedRequest *request) { 94 return CeedElemRestrictionApply_Ref_Core(r, 1, 1, start, stop, tmode, lmode, 95 u, v, request); 96 } 97 98 static int CeedElemRestrictionApply_Ref_Core_18(CeedElemRestriction r, 99 CeedInt start, CeedInt stop, CeedTransposeMode tmode, 100 CeedTransposeMode lmode, CeedVector u, CeedVector v, CeedRequest *request) { 101 return CeedElemRestrictionApply_Ref_Core(r, 8, 1, start, stop, tmode, lmode, 102 u, v, request); 103 104 } 105 106 static int CeedElemRestrictionApply_Ref_Core_31(CeedElemRestriction r, 107 CeedInt start, CeedInt stop, CeedTransposeMode tmode, 108 CeedTransposeMode lmode, CeedVector u, CeedVector v, CeedRequest *request) { 109 return CeedElemRestrictionApply_Ref_Core(r, 1, 3, start, stop, tmode, lmode, 110 u, v, request); 111 } 112 113 static int CeedElemRestrictionApply_Ref_Core_38(CeedElemRestriction r, 114 CeedInt start, CeedInt stop, CeedTransposeMode tmode, 115 CeedTransposeMode lmode, CeedVector u, CeedVector v, CeedRequest *request) { 116 return CeedElemRestrictionApply_Ref_Core(r, 8, 3, start, stop, tmode, lmode, 117 u, v, request); 118 } 119 120 static int CeedElemRestrictionApply_Ref(CeedElemRestriction r, 121 CeedTransposeMode tmode, 122 CeedTransposeMode lmode, CeedVector u, 123 CeedVector v, CeedRequest *request) { 124 int ierr; 125 CeedInt numblk, ncomp, blksize; 126 ierr = CeedElemRestrictionGetNumBlocks(r, &numblk); CeedChk(ierr); 127 ierr = CeedElemRestrictionGetNumComponents(r, &ncomp); CeedChk(ierr); 128 ierr = CeedElemRestrictionGetBlockSize(r, &blksize); CeedChk(ierr); 129 130 CeedInt idx = -1; 131 if (blksize < 10) 132 idx = 10*ncomp + blksize; 133 switch (idx) { 134 case 11: 135 return CeedElemRestrictionApply_Ref_Core_11(r, 0, numblk, tmode, lmode, u, 136 v, request); 137 break; 138 case 18: 139 return CeedElemRestrictionApply_Ref_Core_18(r, 0, numblk, tmode, lmode, u, 140 v, request); 141 break; 142 case 31: 143 return CeedElemRestrictionApply_Ref_Core_31(r, 0, numblk, tmode, lmode, u, 144 v, request); 145 break; 146 case 38: 147 return CeedElemRestrictionApply_Ref_Core_38(r, 0, numblk, tmode, lmode, u, 148 v, request); 149 break; 150 default: 151 // LCOV_EXCL_START 152 return CeedElemRestrictionApply_Ref_Core(r, blksize, ncomp, 0, numblk, 153 tmode, lmode, u, v, request); 154 // LCOV_EXCL_STOP 155 } 156 } 157 158 static int CeedElemRestrictionApplyBlock_Ref(CeedElemRestriction r, 159 CeedInt block, CeedTransposeMode tmode, CeedTransposeMode lmode, 160 CeedVector u, CeedVector v, CeedRequest *request) { 161 int ierr; 162 CeedInt ncomp, blksize; 163 ierr = CeedElemRestrictionGetNumComponents(r, &ncomp); CeedChk(ierr); 164 ierr = CeedElemRestrictionGetBlockSize(r, &blksize); CeedChk(ierr); 165 166 CeedInt idx = -1; 167 if (blksize < 10) 168 idx = 10*ncomp + blksize; 169 switch (idx) { 170 case 11: 171 return CeedElemRestrictionApply_Ref_Core_11(r, block, block+1, tmode, lmode, 172 u, v, request); 173 break; 174 case 18: 175 return CeedElemRestrictionApply_Ref_Core_18(r, block, block+1, tmode, lmode, 176 u, v, request); 177 break; 178 case 31: 179 return CeedElemRestrictionApply_Ref_Core_31(r, block, block+1, tmode, lmode, 180 u, v, request); 181 break; 182 case 38: 183 return CeedElemRestrictionApply_Ref_Core_38(r, block, block+1, tmode, lmode, 184 u, v, request); 185 break; 186 default: 187 // LCOV_EXCL_START 188 return CeedElemRestrictionApply_Ref_Core(r, blksize, ncomp, block, block+1, 189 tmode, lmode, u, v, request); 190 // LCOV_EXCL_STOP 191 } 192 } 193 194 static int CeedElemRestrictionDestroy_Ref(CeedElemRestriction r) { 195 int ierr; 196 CeedElemRestriction_Ref *impl; 197 ierr = CeedElemRestrictionGetData(r, (void *)&impl); CeedChk(ierr); 198 199 ierr = CeedFree(&impl->indices_allocated); CeedChk(ierr); 200 ierr = CeedFree(&impl); CeedChk(ierr); 201 return 0; 202 } 203 204 int CeedElemRestrictionCreate_Ref(CeedMemType mtype, CeedCopyMode cmode, 205 const CeedInt *indices, CeedElemRestriction r) { 206 int ierr; 207 CeedElemRestriction_Ref *impl; 208 CeedInt elemsize, nelem; 209 ierr = CeedElemRestrictionGetNumElements(r, &nelem); CeedChk(ierr); 210 ierr = CeedElemRestrictionGetElementSize(r, &elemsize); CeedChk(ierr); 211 Ceed ceed; 212 ierr = CeedElemRestrictionGetCeed(r, &ceed); CeedChk(ierr); 213 214 if (mtype != CEED_MEM_HOST) 215 // LCOV_EXCL_START 216 return CeedError(ceed, 1, "Only MemType = HOST supported"); 217 // LCOV_EXCL_STOP 218 ierr = CeedCalloc(1,&impl); CeedChk(ierr); 219 switch (cmode) { 220 case CEED_COPY_VALUES: 221 ierr = CeedMalloc(nelem*elemsize, &impl->indices_allocated); 222 CeedChk(ierr); 223 memcpy(impl->indices_allocated, indices, 224 nelem * elemsize * sizeof(indices[0])); 225 impl->indices = impl->indices_allocated; 226 break; 227 case CEED_OWN_POINTER: 228 impl->indices_allocated = (CeedInt *)indices; 229 impl->indices = impl->indices_allocated; 230 break; 231 case CEED_USE_POINTER: 232 impl->indices = indices; 233 } 234 235 ierr = CeedElemRestrictionSetData(r, (void *)&impl); CeedChk(ierr); 236 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", 237 CeedElemRestrictionApply_Ref); CeedChk(ierr); 238 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyBlock", 239 CeedElemRestrictionApplyBlock_Ref); 240 CeedChk(ierr); 241 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", 242 CeedElemRestrictionDestroy_Ref); CeedChk(ierr); 243 return 0; 244 } 245