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 int CeedElemRestrictionApply_Ref(CeedElemRestriction r, 20 CeedTransposeMode tmode, 21 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 nblk, blksize, nelem, elemsize, ndof, ncomp; 29 ierr = CeedElemRestrictionGetNumBlocks(r, &nblk); CeedChk(ierr); 30 ierr = CeedElemRestrictionGetBlockSize(r, &blksize); CeedChk(ierr); 31 ierr = CeedElemRestrictionGetNumElements(r, &nelem); CeedChk(ierr); 32 ierr = CeedElemRestrictionGetElementSize(r, &elemsize); CeedChk(ierr); 33 ierr = CeedElemRestrictionGetNumDoF(r, &ndof); CeedChk(ierr); 34 ierr = CeedElemRestrictionGetNumComponents(r, &ncomp); CeedChk(ierr); 35 36 ierr = CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu); CeedChk(ierr); 37 ierr = CeedVectorGetArray(v, CEED_MEM_HOST, &vv); CeedChk(ierr); 38 // Restriction from lvector to evector 39 // Perform: v = r * u 40 if (tmode == CEED_NOTRANSPOSE) { 41 // No indicies provided, Identity Restriction 42 if (!impl->indices) { 43 for (CeedInt e = 0; e < nblk*blksize; e+=blksize) 44 for (CeedInt j = 0; j < blksize; j++) 45 for (CeedInt k = 0; k < ncomp*elemsize; k++) 46 vv[e*elemsize*ncomp + k*blksize + j] 47 = uu[CeedIntMin(e+j,nelem-1)*ncomp*elemsize + k]; 48 } else { 49 // Indicies provided, standard or blocked restriction 50 // vv has shape [elemsize, ncomp, nelem], row-major 51 // uu has shape [ndof, ncomp] 52 for (CeedInt e = 0; e < nblk*blksize; e+=blksize) 53 for (CeedInt d = 0; d < ncomp; d++) 54 for (CeedInt i = 0; i < elemsize*blksize; i++) 55 vv[i+elemsize*(d*blksize+ncomp*e)] 56 = uu[lmode == CEED_NOTRANSPOSE 57 ? impl->indices[i+elemsize*e]+ndof*d 58 : d+ncomp*impl->indices[i+elemsize*e]]; 59 } 60 } else { 61 // Restriction from evector to lvector 62 // Performing v += r^T * u 63 // No indicies provided, Identity Restriction 64 if (!impl->indices) { 65 for (CeedInt e = 0; e < nblk*blksize; e+=blksize) 66 for (CeedInt j = 0; j < CeedIntMin(blksize, nelem-e); j++) 67 for (CeedInt k = 0; k < ncomp*elemsize; k++) 68 vv[(e+j)*ncomp*elemsize + k] += uu[e*elemsize*ncomp + k*blksize + j]; 69 } else { 70 // Indicies provided, standard or blocked restriction 71 // uu has shape [elemsize, ncomp, nelem] 72 // vv has shape [ndof, ncomp] 73 for (CeedInt e = 0; e < nblk*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]+ndof*d 80 : d+ncomp*impl->indices[j+e*elemsize]] 81 += uu[j+elemsize*(d*blksize+ncomp*e)]; 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 CeedElemRestrictionDestroy_Ref(CeedElemRestriction r) { 93 int ierr; 94 CeedElemRestriction_Ref *impl; 95 ierr = CeedElemRestrictionGetData(r, (void *)&impl); CeedChk(ierr); 96 97 ierr = CeedFree(&impl->indices_allocated); CeedChk(ierr); 98 ierr = CeedFree(&impl); CeedChk(ierr); 99 return 0; 100 } 101 102 int CeedElemRestrictionCreate_Ref(CeedMemType mtype, CeedCopyMode cmode, 103 const CeedInt *indices, CeedElemRestriction r) { 104 int ierr; 105 CeedElemRestriction_Ref *impl; 106 CeedInt elemsize, nelem; 107 ierr = CeedElemRestrictionGetNumElements(r, &nelem); CeedChk(ierr); 108 ierr = CeedElemRestrictionGetElementSize(r, &elemsize); CeedChk(ierr); 109 Ceed ceed; 110 ierr = CeedElemRestrictionGetCeed(r, &ceed); CeedChk(ierr); 111 112 if (mtype != CEED_MEM_HOST) 113 return CeedError(ceed, 1, "Only MemType = HOST supported"); 114 ierr = CeedCalloc(1,&impl); CeedChk(ierr); 115 switch (cmode) { 116 case CEED_COPY_VALUES: 117 ierr = CeedMalloc(nelem*elemsize, &impl->indices_allocated); 118 CeedChk(ierr); 119 memcpy(impl->indices_allocated, indices, 120 nelem * elemsize * sizeof(indices[0])); 121 impl->indices = impl->indices_allocated; 122 break; 123 case CEED_OWN_POINTER: 124 impl->indices_allocated = (CeedInt *)indices; 125 impl->indices = impl->indices_allocated; 126 break; 127 case CEED_USE_POINTER: 128 impl->indices = indices; 129 } 130 131 ierr = CeedElemRestrictionSetData(r, (void *)&impl); CeedChk(ierr); 132 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Apply", 133 CeedElemRestrictionApply_Ref); CeedChk(ierr); 134 ierr = CeedSetBackendFunction(ceed, "ElemRestriction", r, "Destroy", 135 CeedElemRestrictionDestroy_Ref); CeedChk(ierr); 136 return 0; 137 } 138