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