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 CeedElemRestriction_Ref *impl = r->data; 26 int ierr; 27 const CeedScalar *uu; 28 CeedScalar *vv; 29 CeedInt nblk = r->nblk, blksize = r->blksize, elemsize = r->elemsize, 30 ncomp=r->ncomp; 31 32 ierr = CeedVectorGetArrayRead(u, CEED_MEM_HOST, &uu); CeedChk(ierr); 33 ierr = CeedVectorGetArray(v, CEED_MEM_HOST, &vv); CeedChk(ierr); 34 if (tmode == CEED_NOTRANSPOSE) { 35 // Perform: v = r * u 36 if (!impl->indices) { 37 for (CeedInt e = 0; e < nblk*blksize; e+=blksize) 38 for (CeedInt j = 0; j < blksize; j++) 39 for (CeedInt k = 0; k < ncomp*elemsize; k++) 40 vv[e*elemsize*ncomp + k*blksize + j] 41 = uu[CeedIntMin(e+j,r->nelem-1)*ncomp*elemsize + k]; 42 } else { 43 // vv is (elemsize x ncomp x nelem), column-major 44 // uu is (ndof x ncomp) 45 for (CeedInt e = 0; e < nblk*blksize; e+=blksize) 46 for (CeedInt d = 0; d < ncomp; d++) 47 for (CeedInt i = 0; i < elemsize*blksize; i++) 48 vv[i+elemsize*(d+ncomp*e)] 49 = uu[lmode == CEED_NOTRANSPOSE 50 ? impl->indices[i+elemsize*e]+r->ndof*d 51 : d+ncomp*impl->indices[i+elemsize*e]]; 52 } 53 } else { 54 // Note: in transpose mode, we perform: v += r^t * u 55 if (!impl->indices) { 56 for (CeedInt e = 0; e < nblk*blksize; e+=blksize) 57 for (CeedInt j = 0; j < CeedIntMin(blksize, r->nelem-e); j++) 58 for (CeedInt k = 0; k < ncomp*elemsize; k++) 59 vv[(e+j)*ncomp*elemsize + k] += uu[e*elemsize*ncomp + k*blksize + j]; 60 } else { 61 // uu is (elemsize x ncomp x nelem) 62 // vv is (ndof x ncomp) 63 for (CeedInt e = 0; e < nblk*blksize; e+=blksize) { 64 for (CeedInt d = 0; d < ncomp; d++) 65 for (CeedInt i = 0; i < elemsize*blksize; i+=blksize) 66 for (CeedInt j = i; j < i+CeedIntMin(blksize, r->nelem-e); j++) 67 vv[lmode == CEED_NOTRANSPOSE 68 ? impl->indices[j+e*elemsize]+r->ndof*d 69 : d+ncomp*impl->indices[j+e*elemsize]] 70 += uu[j+elemsize*(d+ncomp*e)]; 71 } 72 } 73 } 74 ierr = CeedVectorRestoreArrayRead(u, &uu); CeedChk(ierr); 75 ierr = CeedVectorRestoreArray(v, &vv); CeedChk(ierr); 76 if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED) 77 *request = NULL; 78 return 0; 79 } 80 81 static int CeedElemRestrictionDestroy_Ref(CeedElemRestriction r) { 82 CeedElemRestriction_Ref *impl = r->data; 83 int ierr; 84 85 ierr = CeedFree(&impl->indices_allocated); CeedChk(ierr); 86 ierr = CeedFree(&r->data); CeedChk(ierr); 87 return 0; 88 } 89 90 int CeedElemRestrictionCreate_Ref(CeedElemRestriction r, 91 CeedMemType mtype, 92 CeedCopyMode cmode, const CeedInt *indices) { 93 int ierr; 94 CeedElemRestriction_Ref *impl; 95 96 if (mtype != CEED_MEM_HOST) 97 return CeedError(r->ceed, 1, "Only MemType = HOST supported"); 98 ierr = CeedCalloc(1,&impl); CeedChk(ierr); 99 switch (cmode) { 100 case CEED_COPY_VALUES: 101 ierr = CeedMalloc(r->nelem*r->elemsize, &impl->indices_allocated); 102 CeedChk(ierr); 103 memcpy(impl->indices_allocated, indices, 104 r->nelem * r->elemsize * sizeof(indices[0])); 105 impl->indices = impl->indices_allocated; 106 break; 107 case CEED_OWN_POINTER: 108 impl->indices_allocated = (CeedInt *)indices; 109 impl->indices = impl->indices_allocated; 110 break; 111 case CEED_USE_POINTER: 112 impl->indices = indices; 113 } 114 r->data = impl; 115 r->Apply = CeedElemRestrictionApply_Ref; 116 r->Destroy = CeedElemRestrictionDestroy_Ref; 117 return 0; 118 } 119