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, CeedInt start, CeedInt stop, CeedTransposeMode tmode, 21 CeedTransposeMode lmode, CeedVector u, CeedVector v, CeedRequest *request) { 22 int ierr; 23 CeedElemRestriction_Ref *impl; 24 ierr = CeedElemRestrictionGetData(r, (void *)&impl); CeedChk(ierr);; 25 const CeedScalar *uu; 26 CeedScalar *vv; 27 CeedInt nelem, elemsize, nnodes, ncomp, voffset; 28 ierr = CeedElemRestrictionGetNumElements(r, &nelem); CeedChk(ierr); 29 ierr = CeedElemRestrictionGetElementSize(r, &elemsize); CeedChk(ierr); 30 ierr = CeedElemRestrictionGetNumNodes(r, &nnodes); CeedChk(ierr); 31 ierr = CeedElemRestrictionGetNumComponents(r, &ncomp); 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_1(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, start, stop, tmode, lmode, u, 95 v, request); 96 } 97 98 static int CeedElemRestrictionApply_Ref_Core_8(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, start, stop, tmode, lmode, u, 102 v, request); 103 } 104 105 static int CeedElemRestrictionApply_Ref(CeedElemRestriction r, 106 CeedTransposeMode tmode, 107 CeedTransposeMode lmode, CeedVector u, 108 CeedVector v, CeedRequest *request) { 109 int ierr; 110 CeedInt numblk, blksize; 111 ierr = CeedElemRestrictionGetNumBlocks(r, &numblk); CeedChk(ierr); 112 ierr = CeedElemRestrictionGetBlockSize(r, &blksize); CeedChk(ierr); 113 114 if (blksize == 1) 115 return CeedElemRestrictionApply_Ref_Core_1(r, 0, numblk, tmode, lmode, u, 116 v, request); 117 else if (blksize == 8) 118 return CeedElemRestrictionApply_Ref_Core_8(r, 0, numblk, tmode, lmode, u, 119 v, request); 120 else 121 // LCOV_EXCL_START 122 return CeedElemRestrictionApply_Ref_Core(r, blksize, 0, numblk, tmode, 123 lmode, u, v, request); 124 // LCOV_EXCL_STOP 125 } 126 127 static int CeedElemRestrictionApplyBlock_Ref(CeedElemRestriction r, 128 CeedInt block, CeedTransposeMode tmode, CeedTransposeMode lmode, 129 CeedVector u, CeedVector v, CeedRequest *request) { 130 int ierr; 131 CeedInt blksize; 132 ierr = CeedElemRestrictionGetBlockSize(r, &blksize); CeedChk(ierr); 133 134 if (blksize == 1) 135 return CeedElemRestrictionApply_Ref_Core_1(r, block, block+1, tmode, lmode, 136 u, v, request); 137 else if (blksize == 8) 138 return CeedElemRestrictionApply_Ref_Core_8(r, block, block+1, tmode, lmode, 139 u, v, request); 140 else 141 // LCOV_EXCL_START 142 return CeedElemRestrictionApply_Ref_Core(r, blksize, block, block+1, tmode, 143 lmode, u, v, request); 144 // LCOV_EXCL_STOP 145 } 146 147 static int CeedElemRestrictionDestroy_Ref(CeedElemRestriction r) { 148 int ierr; 149 CeedElemRestriction_Ref *impl; 150 ierr = CeedElemRestrictionGetData(r, (void *)&impl); CeedChk(ierr); 151 152 ierr = CeedFree(&impl->indices_allocated); CeedChk(ierr); 153 ierr = CeedFree(&impl); CeedChk(ierr); 154 return 0; 155 } 156 157 int CeedElemRestrictionCreate_Ref(CeedMemType mtype, CeedCopyMode cmode, 158 const CeedInt *indices, CeedElemRestriction r) { 159 int ierr; 160 CeedElemRestriction_Ref *impl; 161 CeedInt elemsize, nelem; 162 ierr = CeedElemRestrictionGetNumElements(r, &nelem); CeedChk(ierr); 163 ierr = CeedElemRestrictionGetElementSize(r, &elemsize); CeedChk(ierr); 164 Ceed ceed; 165 ierr = CeedElemRestrictionGetCeed(r, &ceed); CeedChk(ierr); 166 167 if (mtype != CEED_MEM_HOST) 168 // LCOV_EXCL_START 169 return CeedError(ceed, 1, "Only MemType = HOST supported"); 170 // LCOV_EXCL_STOP 171 ierr = CeedCalloc(1,&impl); CeedChk(ierr); 172 switch (cmode) { 173 case CEED_COPY_VALUES: 174 ierr = CeedMalloc(nelem*elemsize, &impl->indices_allocated); 175 CeedChk(ierr); 176 memcpy(impl->indices_allocated, indices, 177 nelem * elemsize * sizeof(indices[0])); 178 impl->indices = impl->indices_allocated; 179 break; 180 case CEED_OWN_POINTER: 181 impl->indices_allocated = (CeedInt *)indices; 182 impl->indices = impl->indices_allocated; 183 break; 184 case CEED_USE_POINTER: 185 impl->indices = indices; 186 } 187 188 ierr = CeedElemRestrictionSetData(r, (void *)&impl); CeedChk(ierr); 189 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", 190 CeedElemRestrictionApply_Ref); CeedChk(ierr); 191 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyBlock", 192 CeedElemRestrictionApplyBlock_Ref); 193 CeedChk(ierr); 194 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", 195 CeedElemRestrictionDestroy_Ref); CeedChk(ierr); 196 return 0; 197 } 198