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