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