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 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 blksize, nelem, elemsize, nnodes, ncomp, voffset; 28 ierr = CeedElemRestrictionGetBlockSize(r, &blksize); CeedChk(ierr); 29 ierr = CeedElemRestrictionGetNumElements(r, &nelem); CeedChk(ierr); 30 ierr = CeedElemRestrictionGetElementSize(r, &elemsize); CeedChk(ierr); 31 ierr = CeedElemRestrictionGetNumNodes(r, &nnodes); CeedChk(ierr); 32 ierr = CeedElemRestrictionGetNumComponents(r, &ncomp); CeedChk(ierr); 33 voffset = start*blksize*elemsize*ncomp; 34 35 ierr = CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu); CeedChk(ierr); 36 ierr = CeedVectorGetArray(v, CEED_MEM_HOST, &vv); CeedChk(ierr); 37 // Restriction from lvector to evector 38 // Perform: v = r * u 39 if (tmode == CEED_NOTRANSPOSE) { 40 // No indices provided, Identity Restriction 41 if (!impl->indices) { 42 for (CeedInt e = start*blksize; e < stop*blksize; e+=blksize) 43 for (CeedInt j = 0; j < blksize; j++) 44 for (CeedInt k = 0; k < ncomp*elemsize; k++) 45 vv[e*elemsize*ncomp + k*blksize + j - voffset] 46 = uu[CeedIntMin(e+j,nelem-1)*ncomp*elemsize + k]; 47 } else { 48 // Indices provided, standard or blocked restriction 49 // vv has shape [elemsize, ncomp, nelem], row-major 50 // uu has shape [nnodes, ncomp] 51 for (CeedInt e = start*blksize; e < stop*blksize; e+=blksize) 52 for (CeedInt d = 0; d < ncomp; d++) 53 for (CeedInt i = 0; i < elemsize*blksize; i++) 54 vv[i+elemsize*(d*blksize+ncomp*e) - voffset] 55 = uu[lmode == CEED_NOTRANSPOSE 56 ? impl->indices[i+elemsize*e]+nnodes*d 57 : d+ncomp*impl->indices[i+elemsize*e]]; 58 } 59 } else { 60 // Restriction from evector to lvector 61 // Performing v += r^T * u 62 // No indices provided, Identity Restriction 63 if (!impl->indices) { 64 for (CeedInt e = start*blksize; e < stop*blksize; e+=blksize) 65 for (CeedInt j = 0; j < CeedIntMin(blksize, nelem-e); j++) 66 for (CeedInt k = 0; k < ncomp*elemsize; k++) 67 vv[(e+j)*ncomp*elemsize + k] 68 += uu[e*elemsize*ncomp + k*blksize + j - voffset]; 69 } else { 70 // Indices provided, standard or blocked restriction 71 // uu has shape [elemsize, ncomp, nelem] 72 // vv has shape [nnodes, ncomp] 73 for (CeedInt e = start*blksize; e < stop*blksize; e+=blksize) { 74 for (CeedInt d = 0; d < ncomp; d++) 75 for (CeedInt i = 0; i < elemsize*blksize; i+=blksize) 76 // Iteration bound set to discard padding elements 77 for (CeedInt j = i; j < i+CeedIntMin(blksize, nelem-e); j++) 78 vv[lmode == CEED_NOTRANSPOSE 79 ? impl->indices[j+e*elemsize]+nnodes*d 80 : d+ncomp*impl->indices[j+e*elemsize]] 81 += uu[j+elemsize*(d*blksize+ncomp*e) - voffset]; 82 } 83 } 84 } 85 ierr = CeedVectorRestoreArrayRead(u, &uu); CeedChk(ierr); 86 ierr = CeedVectorRestoreArray(v, &vv); CeedChk(ierr); 87 if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED) 88 *request = NULL; 89 return 0; 90 } 91 92 static int CeedElemRestrictionApply_Ref(CeedElemRestriction r, 93 CeedTransposeMode tmode, 94 CeedTransposeMode lmode, CeedVector u, 95 CeedVector v, CeedRequest *request) { 96 int ierr; 97 CeedInt nblk; 98 ierr = CeedElemRestrictionGetNumBlocks(r, &nblk); CeedChk(ierr); 99 return CeedElemRestrictionApply_Ref_Core(r, 0, nblk, tmode, lmode, u, v, 100 request); 101 } 102 103 static int CeedElemRestrictionApplyBlock_Ref(CeedElemRestriction r, 104 CeedInt block, CeedTransposeMode tmode, CeedTransposeMode lmode, 105 CeedVector u, CeedVector v, CeedRequest *request) { 106 return CeedElemRestrictionApply_Ref_Core(r, block, block+1, tmode, lmode, u, 107 v, request); 108 } 109 110 static int CeedElemRestrictionDestroy_Ref(CeedElemRestriction r) { 111 int ierr; 112 CeedElemRestriction_Ref *impl; 113 ierr = CeedElemRestrictionGetData(r, (void *)&impl); CeedChk(ierr); 114 115 ierr = CeedFree(&impl->indices_allocated); CeedChk(ierr); 116 ierr = CeedFree(&impl); CeedChk(ierr); 117 return 0; 118 } 119 120 int CeedElemRestrictionCreate_Ref(CeedMemType mtype, CeedCopyMode cmode, 121 const CeedInt *indices, CeedElemRestriction r) { 122 int ierr; 123 CeedElemRestriction_Ref *impl; 124 CeedInt elemsize, nelem; 125 ierr = CeedElemRestrictionGetNumElements(r, &nelem); CeedChk(ierr); 126 ierr = CeedElemRestrictionGetElementSize(r, &elemsize); CeedChk(ierr); 127 Ceed ceed; 128 ierr = CeedElemRestrictionGetCeed(r, &ceed); CeedChk(ierr); 129 130 if (mtype != CEED_MEM_HOST) 131 return CeedError(ceed, 1, "Only MemType = HOST supported"); 132 ierr = CeedCalloc(1,&impl); CeedChk(ierr); 133 switch (cmode) { 134 case CEED_COPY_VALUES: 135 ierr = CeedMalloc(nelem*elemsize, &impl->indices_allocated); 136 CeedChk(ierr); 137 memcpy(impl->indices_allocated, indices, 138 nelem * elemsize * sizeof(indices[0])); 139 impl->indices = impl->indices_allocated; 140 break; 141 case CEED_OWN_POINTER: 142 impl->indices_allocated = (CeedInt *)indices; 143 impl->indices = impl->indices_allocated; 144 break; 145 case CEED_USE_POINTER: 146 impl->indices = indices; 147 } 148 149 ierr = CeedElemRestrictionSetData(r, (void *)&impl); CeedChk(ierr); 150 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", 151 CeedElemRestrictionApply_Ref); CeedChk(ierr); 152 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "ApplyBlock", 153 CeedElemRestrictionApplyBlock_Ref); 154 CeedChk(ierr); 155 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", 156 CeedElemRestrictionDestroy_Ref); CeedChk(ierr); 157 return 0; 158 } 159