ceed-ref-restriction.c (7fe495dce79ea51a5ec2e0980c38a0f44563dade) ceed-ref-restriction.c (06cdd269d0df9b74ee1fc8aae222caec83eae53c)
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.

--- 23 unchanged lines hidden (view full) ---

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++)
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.

--- 23 unchanged lines hidden (view full) ---

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];
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
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 }
44 // uu is (ndof x ncomp)
45 for (CeedInt e = 0; e < nblk*blksize; e+=blksize)
46 for (CeedInt d = 0; d < ncomp; d++)
47 for (CeedInt i = 0; i < elemsize*blksize; i++)
48 vv[i+elemsize*(d+ncomp*e)]
49 = uu[lmode == CEED_NOTRANSPOSE
50 ? impl->indices[i+elemsize*e]+r->ndof*d
51 : d+ncomp*impl->indices[i+elemsize*e]];
57 }
58 } else {
59 // Note: in transpose mode, we perform: v += r^t * u
60 if (!impl->indices) {
52 }
53 } else {
54 // Note: in transpose mode, we perform: v += r^t * u
55 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++)
56 for (CeedInt e = 0; e < nblk*blksize; e+=blksize)
57 for (CeedInt j = 0; j < CeedIntMin(blksize, r->nelem-e); j++)
65 for (CeedInt k = 0; k < ncomp*elemsize; k++)
66 vv[(e+j)*ncomp*elemsize + k] += uu[e*elemsize*ncomp + k*blksize + j];
58 for (CeedInt k = 0; k < ncomp*elemsize; k++)
59 vv[(e+j)*ncomp*elemsize + k] += uu[e*elemsize*ncomp + k*blksize + j];
67 }
68 } else {
60 } 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 }
61 // uu is (elemsize x ncomp x nelem)
62 // vv is (ndof x ncomp)
63 for (CeedInt e = 0; e < nblk*blksize; e+=blksize) {
64 for (CeedInt d = 0; d < ncomp; d++)
65 for (CeedInt i = 0; i < elemsize*blksize; i+=blksize)
66 for (CeedInt j = i; j < i+CeedIntMin(blksize, r->nelem-e); j++)
67 vv[lmode == CEED_NOTRANSPOSE
68 ? impl->indices[j+e*elemsize]+r->ndof*d
69 : d+ncomp*impl->indices[j+e*elemsize]]
70 += uu[j+elemsize*(d+ncomp*e)];
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;

--- 40 unchanged lines hidden --- 		71 }
72 }
73 }
74 ierr = CeedVectorRestoreArrayRead(u, &uu); CeedChk(ierr);
75 ierr = CeedVectorRestoreArray(v, &vv); CeedChk(ierr);
76 if (request != CEED_REQUEST_IMMEDIATE && request != CEED_REQUEST_ORDERED)
77 *request = NULL;
78 return 0;

--- 40 unchanged lines hidden ---