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