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