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