MueLu  Version of the Day
MueLu_RegionRFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 //
3 // ***********************************************************************
4 //
5 // MueLu: A package for multigrid based preconditioning
6 // Copyright 2012 Sandia Corporation
7 //
8 // Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9 // the U.S. Government retains certain rights in this software.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact
39 // Jonathan Hu (jhu@sandia.gov)
40 // Andrey Prokopenko (aprokop@sandia.gov)
41 // Ray Tuminaro (rstumin@sandia.gov)
42 //
43 // ***********************************************************************
44 //
45 // @HEADER
46 #ifndef MUELU_REGIONRFACTORY_DEF_HPP
47 #define MUELU_REGIONRFACTORY_DEF_HPP
48 
49 #include <Xpetra_Matrix.hpp>
50 #include <Xpetra_CrsGraphFactory.hpp>
51 
52 // #include "MueLu_PFactory.hpp"
53 // #include "MueLu_FactoryManagerBase.hpp"
54 #include "MueLu_Monitor.hpp"
55 #include "MueLu_MasterList.hpp"
56 
57 #include "MueLu_RegionRFactory_decl.hpp"
58 
59 namespace MueLu {
60 
61  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
62  RCP<const ParameterList> RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
63  GetValidParameterList() const {
64  RCP<ParameterList> validParamList = rcp(new ParameterList());
65 
66  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null,
67  "Generating factory of the matrix A");
68  validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
69  "Number of spatial dimensions in the problem.");
70  validParamList->set<RCP<const FactoryBase> >("lNodesPerDim", Teuchos::null,
71  "Number of local nodes per spatial dimension on the fine grid.");
72  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null,
73  "Fine level nullspace used to construct the coarse level nullspace.");
74  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null,
75  "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
76  validParamList->set<bool> ("keep coarse coords", false, "Flag to keep coordinates for special coarse grid solve");
77 
78  return validParamList;
79  } // GetValidParameterList()
80 
81  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
82  void RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
83  DeclareInput(Level& fineLevel, Level& /* coarseLevel */) const {
84 
85  Input(fineLevel, "A");
86  Input(fineLevel, "numDimensions");
87  Input(fineLevel, "lNodesPerDim");
88  Input(fineLevel, "Nullspace");
89  Input(fineLevel, "Coordinates");
90 
91  } // DeclareInput()
92 
93  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
94  void RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
95  Build(Level& fineLevel, Level& coarseLevel) const {
96 
97  // Set debug outputs based on environment variable
98  RCP<Teuchos::FancyOStream> out;
99  if(const char* dbg = std::getenv("MUELU_REGIONRFACTORY_DEBUG")) {
100  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
101  out->setShowAllFrontMatter(false).setShowProcRank(true);
102  } else {
103  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
104  }
105 
106  *out << "Starting RegionRFactory::Build." << std::endl;
107 
108  // First get the inputs from the fineLevel
109  const int numDimensions = Get<int>(fineLevel, "numDimensions");
110  Array<LO> lFineNodesPerDim(3, Teuchos::OrdinalTraits<LO>::one());
111  {
112  Array<LO> lNodesPerDim = Get<Array<LO> >(fineLevel, "lNodesPerDim");
113  for(int dim = 0; dim < numDimensions; ++dim) {
114  lFineNodesPerDim[dim] = lNodesPerDim[dim];
115  }
116  }
117  *out << "numDimensions " << numDimensions << " and lFineNodesPerDim: " << lFineNodesPerDim
118  << std::endl;
119 
120  // Let us check that the inputs verify our assumptions
121  if(numDimensions < 1 || numDimensions > 3) {
122  throw std::runtime_error("numDimensions must be 1, 2 or 3!");
123  }
124  for(int dim = 0; dim < numDimensions; ++dim) {
125  if(lFineNodesPerDim[dim] % 3 != 1) {
126  throw std::runtime_error("The number of fine node in each direction need to be 3n+1");
127  }
128  }
129  Array<LO> lCoarseNodesPerDim(3, Teuchos::OrdinalTraits<LO>::one());
130 
131  const RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
132 
133  RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
134  fineCoordinates = Get< RCP<realvaluedmultivector_type> >(fineLevel, "Coordinates");
135  if(static_cast<int>(fineCoordinates->getNumVectors()) != numDimensions) {
136  throw std::runtime_error("The number of vectors in the coordinates is not equal to numDimensions!");
137  }
138 
139  // Let us create R and pass it down to the
140  // appropriate specialization and see what we
141  // get back!
142  RCP<Matrix> R;
143 
144  if(numDimensions == 1) {
145  throw std::runtime_error("RegionRFactory no implemented for 1D case yet.");
146  } else if(numDimensions == 2) {
147  throw std::runtime_error("RegionRFactory no implemented for 2D case yet.");
148  } else if(numDimensions == 3) {
149  Build3D(numDimensions, lFineNodesPerDim, A, fineCoordinates,
150  R, coarseCoordinates, lCoarseNodesPerDim);
151  }
152 
153  const Teuchos::ParameterList& pL = GetParameterList();
154 
155  // Reuse pattern if available (multiple solve)
156  RCP<ParameterList> Tparams;
157  if(pL.isSublist("matrixmatrix: kernel params"))
158  Tparams=rcp(new ParameterList(pL.sublist("matrixmatrix: kernel params")));
159  else
160  Tparams= rcp(new ParameterList);
161 
162  // R->describe(*out, Teuchos::VERB_EXTREME);
163  *out << "Compute P=R^t" << std::endl;
164  // By default, we don't need global constants for transpose
165  Tparams->set("compute global constants: temporaries",Tparams->get("compute global constants: temporaries", false));
166  Tparams->set("compute global constants", Tparams->get("compute global constants",false));
167  std::string label = "MueLu::RegionR-transR" + Teuchos::toString(coarseLevel.GetLevelID());
168  RCP<Matrix> P = Utilities::Transpose(*R, true, label, Tparams);
169 
170  *out << "Compute coarse nullspace" << std::endl;
171  RCP<MultiVector> fineNullspace = Get<RCP<MultiVector> >(fineLevel, "Nullspace");
172  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(R->getRowMap(),
173  fineNullspace->getNumVectors());
174  R->apply(*fineNullspace, *coarseNullspace, Teuchos::NO_TRANS, Teuchos::ScalarTraits<SC>::one(),
175  Teuchos::ScalarTraits<SC>::zero());
176 
177  *out << "Set data on coarse level" << std::endl;
178  Set(coarseLevel, "numDimensions", numDimensions);
179  Set(coarseLevel, "lNodesPerDim", lCoarseNodesPerDim);
180  Set(coarseLevel, "Nullspace", coarseNullspace);
181  Set(coarseLevel, "Coordinates", coarseCoordinates);
182  if(pL.get<bool>("keep coarse coords")) {
183  coarseLevel.Set<RCP<realvaluedmultivector_type> >("Coordinates2", coarseCoordinates, NoFactory::get());
184  }
185 
186  R->SetFixedBlockSize(A->GetFixedBlockSize());
187  P->SetFixedBlockSize(A->GetFixedBlockSize());
188 
189  Set(coarseLevel, "R", R);
190  Set(coarseLevel, "P", P);
191 
192  } // Build()
193 
194  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
195  void RegionRFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
196  Build3D(const int numDimensions,
197  Teuchos::Array<LocalOrdinal>& lFineNodesPerDim,
198  const RCP<Matrix>& A,
199  const RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> >& fineCoordinates,
200  RCP<Matrix>& R,
201  RCP<Xpetra::MultiVector<typename Teuchos::ScalarTraits<Scalar>::coordinateType, LocalOrdinal, GlobalOrdinal, Node> >& coarseCoordinates,
202  Teuchos::Array<LocalOrdinal>& lCoarseNodesPerDim) const {
203  using local_matrix_type = typename CrsMatrix::local_matrix_type;
204  using local_graph_type = typename local_matrix_type::staticcrsgraph_type;
205  using row_map_type = typename local_matrix_type::row_map_type::non_const_type;
206  using entries_type = typename local_matrix_type::index_type::non_const_type;
207  using values_type = typename local_matrix_type::values_type::non_const_type;
208 
209  // Set debug outputs based on environment variable
210  RCP<Teuchos::FancyOStream> out;
211  if(const char* dbg = std::getenv("MUELU_REGIONRFACTORY_DEBUG")) {
212  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
213  out->setShowAllFrontMatter(false).setShowProcRank(true);
214  } else {
215  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
216  }
217 
218 
219  // Now compute number of coarse grid points
220  for(int dim = 0; dim < numDimensions; ++dim) {
221  lCoarseNodesPerDim[dim] = lFineNodesPerDim[dim] / 3 + 1;
222  }
223  *out << "lCoarseNodesPerDim " << lCoarseNodesPerDim << std::endl;
224 
225  // Grab the block size here and multiply all existing offsets by it
226  const LO blkSize = A->GetFixedBlockSize();
227  *out << "blkSize " << blkSize << std::endl;
228 
229  // Based on lCoarseNodesPerDim and lFineNodesPerDim
230  // we can compute numRows, numCols and NNZ for R
231  const LO numRows = blkSize*lCoarseNodesPerDim[0]*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[2];
232  const LO numCols = blkSize*lFineNodesPerDim[0]*lFineNodesPerDim[1]*lFineNodesPerDim[2];
233 
234  // Create the coarse coordinates multivector
235  // so we can fill it on the fly while computing
236  // the restriction operator
237  RCP<Map> rowMap = MapFactory::Build(A->getRowMap()->lib(),
238  Teuchos::OrdinalTraits<GO>::invalid(),
239  numRows,
240  A->getRowMap()->getIndexBase(),
241  A->getRowMap()->getComm());
242 
243  RCP<Map> coordRowMap = MapFactory::Build(A->getRowMap()->lib(),
244  Teuchos::OrdinalTraits<GO>::invalid(),
245  lCoarseNodesPerDim[0]*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[2],
246  A->getRowMap()->getIndexBase(),
247  A->getRowMap()->getComm());
248 
249  coarseCoordinates = Xpetra::MultiVectorFactory<real_type, LO, GO, NO>::Build(coordRowMap,
250  numDimensions);
251  Array<ArrayRCP<const real_type> > fineCoordData(numDimensions);
252  Array<ArrayRCP<real_type> > coarseCoordData(numDimensions);
253  for(int dim = 0; dim < numDimensions; ++dim) {
254  fineCoordData[dim] = fineCoordinates->getData(dim);
255  coarseCoordData[dim] = coarseCoordinates->getDataNonConst(dim);
256  }
257 
258  // Let us set some parameter that will be useful
259  // while constructing R
260 
261  // Length of interpolation stencils based on geometry
262  const LO cornerStencilLength = 27;
263  const LO edgeStencilLength = 45;
264  const LO faceStencilLength = 75;
265  const LO interiorStencilLength = 125;
266 
267  // Number of corner, edge, face and interior nodes
268  const LO numCorners = 8;
269  const LO numEdges = 4*(lCoarseNodesPerDim[0] - 2)
270  + 4*(lCoarseNodesPerDim[1] - 2)
271  + 4*(lCoarseNodesPerDim[2] - 2);
272  const LO numFaces = 2*(lCoarseNodesPerDim[0] - 2)*(lCoarseNodesPerDim[1] - 2)
273  + 2*(lCoarseNodesPerDim[0] - 2)*(lCoarseNodesPerDim[2] - 2)
274  + 2*(lCoarseNodesPerDim[1] - 2)*(lCoarseNodesPerDim[2] - 2);
275  const LO numInteriors = (lCoarseNodesPerDim[0] - 2)*(lCoarseNodesPerDim[1] - 2)
276  *(lCoarseNodesPerDim[2] - 2);
277 
278  const LO nnz = (numCorners*cornerStencilLength + numEdges*edgeStencilLength
279  + numFaces*faceStencilLength + numInteriors*interiorStencilLength)*blkSize;
280 
281  // Having the number of rows and columns we can genrate
282  // the appropriate maps for our operator.
283 
284  *out << "R statistics:" << std::endl
285  << " -numRows= " << numRows << std::endl
286  << " -numCols= " << numCols << std::endl
287  << " -nnz= " << nnz << std::endl;
288 
289  row_map_type row_map(Kokkos::ViewAllocateWithoutInitializing("row_map"), numRows + 1);
290  typename row_map_type::HostMirror row_map_h = Kokkos::create_mirror_view(row_map);
291 
292  entries_type entries(Kokkos::ViewAllocateWithoutInitializing("entries"), nnz);
293  typename entries_type::HostMirror entries_h = Kokkos::create_mirror_view(entries);
294 
295  values_type values(Kokkos::ViewAllocateWithoutInitializing("values"), nnz);
296  typename values_type::HostMirror values_h = Kokkos::create_mirror_view(values);
297 
298  // Compute the basic interpolation
299  // coefficients for 1D rate of 3
300  // coarsening.
301  Array<SC> coeffs({1.0/3.0, 2.0/3.0, 1.0, 2.0/3.0, 1.0/3.0});
302  row_map_h(0) = 0;
303 
304  // Define some offsets that
305  // will be needed often later on
306  const LO edgeLineOffset = 2*cornerStencilLength + (lCoarseNodesPerDim[0] - 2)*edgeStencilLength;
307  const LO faceLineOffset = 2*edgeStencilLength + (lCoarseNodesPerDim[0] - 2)*faceStencilLength;
308  const LO interiorLineOffset = 2*faceStencilLength
309  + (lCoarseNodesPerDim[0] - 2)*interiorStencilLength;
310 
311  const LO facePlaneOffset = 2*edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset;
312  const LO interiorPlaneOffset = 2*faceLineOffset + (lCoarseNodesPerDim[1] - 2)*interiorLineOffset;
313 
314  // Let us take care of the corners
315  // first since we always have
316  // corners to deal with!
317  {
318  // Corner 1
319  LO coordRowIdx = 0, rowIdx = 0, coordColumnOffset = 0, columnOffset = 0, entryOffset = 0;
320  for(LO l = 0; l < blkSize; ++l) {
321  for(LO k = 0; k < 3; ++k) {
322  for(LO j = 0; j < 3; ++j) {
323  for(LO i = 0; i < 3; ++i) {
324  entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
325  + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i)*blkSize + l;
326  values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k + 2]*coeffs[j + 2]*coeffs[i + 2];
327  }
328  }
329  }
330  }
331  for(LO l = 0; l < blkSize; ++l) {
332  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
333  }
334  for(int dim = 0; dim <numDimensions; ++dim) {
335  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
336  }
337 
338  // Corner 5
339  coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
340  rowIdx = coordRowIdx*blkSize;
341  coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
342  columnOffset = coordColumnOffset*blkSize;
343  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
344  for(LO l = 0; l < blkSize; ++l) {
345  for(LO k = 0; k < 3; ++k) {
346  for(LO j = 0; j < 3; ++j) {
347  for(LO i = 0; i < 3; ++i) {
348  entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
349  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i)*blkSize + l;
350  values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i + 2];
351  }
352  }
353  }
354  }
355  for(LO l = 0; l < blkSize; ++l) {
356  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
357  }
358  for(int dim = 0; dim <numDimensions; ++dim) {
359  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
360  }
361 
362  // Corner 2
363  coordRowIdx = (lCoarseNodesPerDim[0] - 1);
364  rowIdx = coordRowIdx*blkSize;
365  coordColumnOffset = (lFineNodesPerDim[0] - 1);
366  columnOffset = coordColumnOffset*blkSize;
367  entryOffset = (cornerStencilLength + (lCoarseNodesPerDim[0] - 2)*edgeStencilLength)*blkSize;
368  for(LO l = 0; l < blkSize; ++l) {
369  for(LO k = 0; k < 3; ++k) {
370  for(LO j = 0; j < 3; ++j) {
371  for(LO i = 0; i < 3; ++i) {
372  entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
373  + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0]
374  + j*lFineNodesPerDim[0] + (i - 2))*blkSize + l;
375  values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k + 2]*coeffs[j + 2]*coeffs[i];
376  }
377  }
378  }
379  }
380  for(LO l = 0; l < blkSize; ++l) {
381  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
382  }
383  for(int dim = 0; dim <numDimensions; ++dim) {
384  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
385  }
386 
387  // Corner 6
388  coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
389  rowIdx = coordRowIdx*blkSize;
390  coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
391  columnOffset = coordColumnOffset*blkSize;
392  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
393  for(LO l = 0; l < blkSize; ++l) {
394  for(LO k = 0; k < 3; ++k) {
395  for(LO j = 0; j < 3; ++j) {
396  for(LO i = 0; i < 3; ++i) {
397  entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
398  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
399  values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i];
400  }
401  }
402  }
403  }
404  for(LO l = 0; l < blkSize; ++l) {
405  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
406  }
407  for(int dim = 0; dim <numDimensions; ++dim) {
408  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
409  }
410 
411  // Corner 3
412  coordRowIdx = (lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0];
413  rowIdx = coordRowIdx*blkSize;
414  coordColumnOffset = (lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0];
415  columnOffset = coordColumnOffset*blkSize;
416  entryOffset = (edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset)*blkSize;
417  for(LO l = 0; l < blkSize; ++l) {
418  for(LO k = 0; k < 3; ++k) {
419  for(LO j = 0; j < 3; ++j) {
420  for(LO i = 0; i < 3; ++i) {
421  entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
422  + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0]
423  + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
424  values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i + 2];
425  }
426  }
427  }
428  }
429  for(LO l = 0; l < blkSize; ++l) {
430  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
431  }
432  for(int dim = 0; dim <numDimensions; ++dim) {
433  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
434  }
435 
436  // Corner 7
437  coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
438  rowIdx = coordRowIdx*blkSize;
439  coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
440  columnOffset = coordColumnOffset*blkSize;
441  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
442  for(LO l = 0; l < blkSize; ++l) {
443  for(LO k = 0; k < 3; ++k) {
444  for(LO j = 0; j < 3; ++j) {
445  for(LO i = 0; i < 3; ++i) {
446  entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
447  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
448  values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i + 2];
449  }
450  }
451  }
452  }
453  for(LO l = 0; l < blkSize; ++l) {
454  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
455  }
456  for(int dim = 0; dim <numDimensions; ++dim) {
457  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
458  }
459 
460  // Corner 4
461  coordRowIdx = (lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0] - 1);
462  rowIdx = coordRowIdx*blkSize;
463  coordColumnOffset = (lFineNodesPerDim[1]*lFineNodesPerDim[0] - 1);
464  columnOffset = coordColumnOffset*blkSize;
465  entryOffset = (edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset +
466  cornerStencilLength + (lCoarseNodesPerDim[0] - 2)*edgeStencilLength)*blkSize;
467  for(LO l = 0; l < blkSize; ++l) {
468  for(LO k = 0; k < 3; ++k) {
469  for(LO j = 0; j < 3; ++j) {
470  for(LO i = 0; i < 3; ++i) {
471  entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
472  + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0]
473  + (j - 2)*lFineNodesPerDim[0] + (i - 2))*blkSize + l;
474  values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i];
475  }
476  }
477  }
478  }
479  for(LO l = 0; l < blkSize; ++l) {
480  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
481  }
482  for(int dim = 0; dim <numDimensions; ++dim) {
483  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
484  }
485 
486  // Corner 8
487  coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
488  rowIdx = coordRowIdx*blkSize;
489  coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
490  columnOffset = coordColumnOffset*blkSize;
491  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
492  for(LO l = 0; l < blkSize; ++l) {
493  for(LO k = 0; k < 3; ++k) {
494  for(LO j = 0; j < 3; ++j) {
495  for(LO i = 0; i < 3; ++i) {
496  entries_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = columnOffset
497  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + (i - 2))*blkSize + l;
498  values_h(entryOffset + k*9 + j*3 + i + cornerStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
499  }
500  }
501  }
502  }
503  for(LO l = 0; l < blkSize; ++l) {
504  row_map_h(rowIdx + 1 + l) = entryOffset + cornerStencilLength*(l+1);
505  }
506  for(int dim = 0; dim <numDimensions; ++dim) {
507  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
508  }
509  } // Corners are done!
510 
511  // Edges along 0 direction
512  if(lCoarseNodesPerDim[0] - 2 > 0) {
513 
514  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
515  for(LO edgeIdx = 0; edgeIdx < lCoarseNodesPerDim[0] - 2; ++edgeIdx) {
516 
517  // Edge 0
518  coordRowIdx = (edgeIdx + 1);
519  rowIdx = coordRowIdx*blkSize;
520  coordColumnOffset = (edgeIdx + 1)*3;
521  columnOffset = coordColumnOffset*blkSize;
522  entryOffset = (cornerStencilLength + edgeIdx*edgeStencilLength)*blkSize;
523  for(LO l = 0; l < blkSize; ++l) {
524  for(LO k = 0; k < 3; ++k) {
525  for(LO j = 0; j < 3; ++j) {
526  for(LO i = 0; i < 5; ++i) {
527  entries_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = columnOffset
528  + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
529  values_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = coeffs[k + 2]*coeffs[j + 2]*coeffs[i];
530  }
531  }
532  }
533  }
534  for(LO l = 0; l < blkSize; ++l) {
535  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
536  }
537  for(int dim = 0; dim <numDimensions; ++dim) {
538  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
539  }
540 
541  // Edge 1
542  coordRowIdx = ((lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0] + edgeIdx + 1);
543  rowIdx = coordRowIdx*blkSize;
544  coordColumnOffset = ((lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0] + (edgeIdx + 1)*3);
545  columnOffset = coordColumnOffset*blkSize;
546  entryOffset = (edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset
547  + cornerStencilLength + edgeIdx*edgeStencilLength)*blkSize;
548  for(LO l = 0; l < blkSize; ++l) {
549  for(LO k = 0; k < 3; ++k) {
550  for(LO j = 0; j < 3; ++j) {
551  for(LO i = 0; i < 5; ++i) {
552  entries_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = columnOffset
553  + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
554  values_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i];
555  }
556  }
557  }
558  }
559  for(LO l = 0; l < blkSize; ++l) {
560  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
561  }
562  for(int dim = 0; dim <numDimensions; ++dim) {
563  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
564  }
565 
566  // Edge 2
567  coordRowIdx = ((lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
568  + edgeIdx + 1);
569  rowIdx = coordRowIdx*blkSize;
570  coordColumnOffset = ((lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
571  + (edgeIdx + 1)*3);
572  columnOffset = coordColumnOffset*blkSize;
573  entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset
574  + cornerStencilLength + edgeIdx*edgeStencilLength)*blkSize;
575  for(LO l = 0; l < blkSize; ++l) {
576  for(LO k = 0; k < 3; ++k) {
577  for(LO j = 0; j < 3; ++j) {
578  for(LO i = 0; i < 5; ++i) {
579  entries_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = columnOffset
580  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
581  values_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i];
582  }
583  }
584  }
585  }
586  for(LO l = 0; l < blkSize; ++l) {
587  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
588  }
589  for(int dim = 0; dim <numDimensions; ++dim) {
590  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
591  }
592 
593  // Edge 3
594  coordRowIdx = ((lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
595  + (lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0] + edgeIdx + 1);
596  rowIdx = coordRowIdx*blkSize;
597  coordColumnOffset = ((lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
598  + (lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0] + (edgeIdx + 1)*3);
599  columnOffset = coordColumnOffset*blkSize;
600  entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset
601  + edgeLineOffset + (lCoarseNodesPerDim[1] - 2)*faceLineOffset
602  + cornerStencilLength + edgeIdx*edgeStencilLength)*blkSize;
603  for(LO l = 0; l < blkSize; ++l) {
604  for(LO k = 0; k < 3; ++k) {
605  for(LO j = 0; j < 3; ++j) {
606  for(LO i = 0; i < 5; ++i) {
607  entries_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = columnOffset
608  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
609  + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
610  values_h(entryOffset + k*15 + j*5 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
611  }
612  }
613  }
614  }
615  for(LO l = 0; l < blkSize; ++l) {
616  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
617  }
618  for(int dim = 0; dim <numDimensions; ++dim) {
619  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
620  }
621  }
622  }
623 
624  // Edges along 1 direction
625  if(lCoarseNodesPerDim[1] - 2 > 0) {
626 
627  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
628  for(LO edgeIdx = 0; edgeIdx < lCoarseNodesPerDim[1] - 2; ++edgeIdx) {
629 
630  // Edge 0
631  coordRowIdx = (edgeIdx + 1)*lCoarseNodesPerDim[0];
632  rowIdx = coordRowIdx*blkSize;
633  coordColumnOffset = (edgeIdx + 1)*3*lFineNodesPerDim[0];
634  columnOffset = coordColumnOffset*blkSize;
635  entryOffset = (edgeLineOffset + edgeIdx*faceLineOffset)*blkSize;
636  for(LO l = 0; l < blkSize; ++l) {
637  for(LO k = 0; k < 3; ++k) {
638  for(LO j = 0; j < 5; ++j) {
639  for(LO i = 0; i < 3; ++i) {
640  entries_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = columnOffset
641  + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
642  values_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i + 2];
643  }
644  }
645  }
646  }
647  for(LO l = 0; l < blkSize; ++l) {
648  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
649  }
650  for(int dim = 0; dim <numDimensions; ++dim) {
651  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
652  }
653 
654  // Edge 1
655  coordRowIdx = ((edgeIdx + 1)*lCoarseNodesPerDim[0] + lCoarseNodesPerDim[0] - 1);
656  rowIdx = coordRowIdx*blkSize;
657  coordColumnOffset = ((edgeIdx + 1)*3*lFineNodesPerDim[0] + lFineNodesPerDim[0] - 1);
658  columnOffset = coordColumnOffset*blkSize;
659  entryOffset = (edgeLineOffset + edgeIdx*faceLineOffset
660  + edgeStencilLength + (lCoarseNodesPerDim[0] - 2)*faceStencilLength)*blkSize;
661  for(LO l = 0; l < blkSize; ++l) {
662  for(LO k = 0; k < 3; ++k) {
663  for(LO j = 0; j < 5; ++j) {
664  for(LO i = 0; i < 3; ++i) {
665  entries_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = columnOffset
666  + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
667  values_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i];
668  }
669  }
670  }
671  }
672  for(LO l = 0; l < blkSize; ++l) {
673  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
674  }
675  for(int dim = 0; dim <numDimensions; ++dim) {
676  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
677  }
678 
679  // Edge 2
680  coordRowIdx = ((lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
681  + (edgeIdx + 1)*lCoarseNodesPerDim[0]);
682  rowIdx = coordRowIdx*blkSize;
683  coordColumnOffset = ((lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
684  + (edgeIdx + 1)*3*lFineNodesPerDim[0]);
685  columnOffset = coordColumnOffset*blkSize;
686  entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset
687  + edgeLineOffset + edgeIdx*faceLineOffset)*blkSize;
688  for(LO l = 0; l < blkSize; ++l) {
689  for(LO k = 0; k < 3; ++k) {
690  for(LO j = 0; j < 5; ++j) {
691  for(LO i = 0; i < 3; ++i) {
692  entries_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = columnOffset
693  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
694  values_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i + 2];
695  }
696  }
697  }
698  }
699  for(LO l = 0; l < blkSize; ++l) {
700  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
701  }
702  for(int dim = 0; dim <numDimensions; ++dim) {
703  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
704  }
705 
706  // Edge 3
707  coordRowIdx = ((lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
708  + (edgeIdx + 1)*lCoarseNodesPerDim[0] + lCoarseNodesPerDim[0] - 1);
709  rowIdx = coordRowIdx*blkSize;
710  coordColumnOffset = ((lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
711  + (edgeIdx + 1)*3*lFineNodesPerDim[0] + lFineNodesPerDim[0] - 1);
712  columnOffset = coordColumnOffset*blkSize;
713  entryOffset = (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset
714  + edgeLineOffset + edgeIdx*faceLineOffset
715  + edgeStencilLength + (lCoarseNodesPerDim[0] - 2)*faceStencilLength)*blkSize;
716  for(LO l = 0; l < blkSize; ++l) {
717  for(LO k = 0; k < 3; ++k) {
718  for(LO j = 0; j < 5; ++j) {
719  for(LO i = 0; i < 3; ++i) {
720  entries_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = columnOffset
721  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
722  + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
723  values_h(entryOffset + k*15 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
724  }
725  }
726  }
727  }
728  for(LO l = 0; l < blkSize; ++l) {
729  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
730  }
731  for(int dim = 0; dim <numDimensions; ++dim) {
732  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
733  }
734  }
735  }
736 
737  // Edges along 2 direction
738  if(lCoarseNodesPerDim[2] - 2 > 0) {
739 
740  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
741  for(LO edgeIdx = 0; edgeIdx < lCoarseNodesPerDim[2] - 2; ++edgeIdx) {
742 
743  // Edge 0
744  coordRowIdx = (edgeIdx + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
745  rowIdx = coordRowIdx*blkSize;
746  coordColumnOffset = (edgeIdx + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0];
747  columnOffset = coordColumnOffset*blkSize;
748  entryOffset = (facePlaneOffset + edgeIdx*interiorPlaneOffset)*blkSize;
749  for(LO l = 0; l < blkSize; ++l) {
750  for(LO k = 0; k < 5; ++k) {
751  for(LO j = 0; j < 3; ++j) {
752  for(LO i = 0; i < 3; ++i) {
753  entries_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = columnOffset
754  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i)*blkSize + l;
755  values_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i + 2];
756  }
757  }
758  }
759  }
760  for(LO l = 0; l < blkSize; ++l) {
761  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
762  }
763  for(int dim = 0; dim <numDimensions; ++dim) {
764  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
765  }
766 
767  // Edge 1
768  coordRowIdx = ((edgeIdx + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
769  + lCoarseNodesPerDim[0] - 1);
770  rowIdx = coordRowIdx*blkSize;
771  coordColumnOffset = ((edgeIdx + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0]
772  + lFineNodesPerDim[0] - 1);
773  columnOffset = coordColumnOffset*blkSize;
774  entryOffset = (facePlaneOffset + faceLineOffset - edgeStencilLength
775  + edgeIdx*interiorPlaneOffset)*blkSize;
776  for(LO l = 0; l < blkSize; ++l) {
777  for(LO k = 0; k < 5; ++k) {
778  for(LO j = 0; j < 3; ++j) {
779  for(LO i = 0; i < 3; ++i) {
780  entries_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = columnOffset
781  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
782  values_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i];
783  }
784  }
785  }
786  }
787  for(LO l = 0; l < blkSize; ++l) {
788  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
789  }
790  for(int dim = 0; dim <numDimensions; ++dim) {
791  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
792  }
793 
794  // Edge 2
795  coordRowIdx = ((edgeIdx + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
796  + (lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0]);
797  rowIdx = coordRowIdx*blkSize;
798  coordColumnOffset = ((edgeIdx + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0]
799  + (lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0]);
800  columnOffset = coordColumnOffset*blkSize;
801  entryOffset = (facePlaneOffset + edgeIdx*interiorPlaneOffset + faceLineOffset
802  + (lCoarseNodesPerDim[1] - 2)*interiorLineOffset)*blkSize;
803  for(LO l = 0; l < blkSize; ++l) {
804  for(LO k = 0; k < 5; ++k) {
805  for(LO j = 0; j < 3; ++j) {
806  for(LO i = 0; i < 3; ++i) {
807  entries_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = columnOffset
808  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
809  values_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i + 2];
810  }
811  }
812  }
813  }
814  for(LO l = 0; l < blkSize; ++l) {
815  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
816  }
817  for(int dim = 0; dim <numDimensions; ++dim) {
818  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
819  }
820 
821  // Edge 3
822  coordRowIdx = ((edgeIdx + 2)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0] - 1);
823  rowIdx = coordRowIdx*blkSize;
824  coordColumnOffset = ((edgeIdx + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0]
825  + lFineNodesPerDim[1]*lFineNodesPerDim[0] - 1);
826  columnOffset = coordColumnOffset*blkSize;
827  entryOffset = (facePlaneOffset + (edgeIdx + 1)*interiorPlaneOffset - edgeStencilLength)*blkSize;
828  for(LO l = 0; l < blkSize; ++l) {
829  for(LO k = 0; k < 5; ++k) {
830  for(LO j = 0; j < 3; ++j) {
831  for(LO i = 0; i < 3; ++i) {
832  entries_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = columnOffset
833  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
834  + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
835  values_h(entryOffset + k*9 + j*3 + i + edgeStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
836  }
837  }
838  }
839  }
840  for(LO l = 0; l < blkSize; ++l) {
841  row_map_h(rowIdx + 1 + l) = entryOffset + edgeStencilLength*(l+1);
842  }
843  for(int dim = 0; dim <numDimensions; ++dim) {
844  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
845  }
846  }
847  }
848 
849  // Faces in 0-1 plane
850  if((lCoarseNodesPerDim[0] - 2 > 0) && (lCoarseNodesPerDim[1] - 2 > 0)) {
851 
852  Array<LO> gridIdx(3);
853  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
854  for(LO faceIdx=0; faceIdx < (lCoarseNodesPerDim[1]-2)*(lCoarseNodesPerDim[0]-2); ++faceIdx) {
855 
856  // Face 0
857  coordRowIdx = ((gridIdx[1] + 1)*lCoarseNodesPerDim[0] + gridIdx[0] + 1);
858  rowIdx = coordRowIdx*blkSize;
859  coordColumnOffset = 3*((gridIdx[1] + 1)*lFineNodesPerDim[0] + gridIdx[0] + 1);
860  columnOffset = coordColumnOffset*blkSize;
861  entryOffset = (edgeLineOffset + edgeStencilLength
862  + gridIdx[1]*faceLineOffset + gridIdx[0]*faceStencilLength)*blkSize;
863  for(LO l = 0; l < blkSize; ++l) {
864  for(LO k = 0; k < 3; ++k) {
865  for(LO j = 0; j < 5; ++j) {
866  for(LO i = 0; i < 5; ++i) {
867  entries_h(entryOffset + k*25 + j*5 + i + faceStencilLength*l) = columnOffset
868  + (k*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
869  values_h(entryOffset + k*25 + j*5 + i + faceStencilLength*l) = coeffs[k + 2]*coeffs[j]*coeffs[i];
870  }
871  }
872  }
873  }
874  for(LO l = 0; l < blkSize; ++l) {
875  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
876  }
877  for(int dim = 0; dim <numDimensions; ++dim) {
878  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
879  }
880 
881  // Face 1
882  coordRowIdx += (lCoarseNodesPerDim[2] - 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0];
883  rowIdx = coordRowIdx*blkSize;
884  coordColumnOffset += (lFineNodesPerDim[2] - 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0];
885  columnOffset = coordColumnOffset*blkSize;
886  entryOffset += (facePlaneOffset + (lCoarseNodesPerDim[2] - 2)*interiorPlaneOffset)*blkSize;
887  for(LO l = 0; l < blkSize; ++l) {
888  for(LO k = 0; k < 3; ++k) {
889  for(LO j = 0; j < 5; ++j) {
890  for(LO i = 0; i < 5; ++i) {
891  entries_h(entryOffset + k*25 + j*5 + i + faceStencilLength*l) = columnOffset
892  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
893  + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
894  values_h(entryOffset + k*25 + j*5 + i + faceStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
895  }
896  }
897  }
898  }
899  for(LO l = 0; l < blkSize; ++l) {
900  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
901  }
902  for(int dim = 0; dim <numDimensions; ++dim) {
903  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
904  }
905 
906  // Last step in the loop
907  // update the grid indices
908  // for next grid point
909  ++gridIdx[0];
910  if(gridIdx[0] == lCoarseNodesPerDim[0] - 2) {
911  gridIdx[0] = 0;
912  ++gridIdx[1];
913  }
914  }
915  }
916 
917  // Faces in 0-2 plane
918  if((lCoarseNodesPerDim[0] - 2 > 0) && (lCoarseNodesPerDim[2] - 2 > 0)) {
919 
920  Array<LO> gridIdx(3);
921  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
922  for(LO faceIdx=0; faceIdx < (lCoarseNodesPerDim[2]-2)*(lCoarseNodesPerDim[0]-2); ++faceIdx) {
923 
924  // Face 0
925  coordRowIdx = ((gridIdx[2] + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0] + (gridIdx[0] + 1));
926  rowIdx = coordRowIdx*blkSize;
927  coordColumnOffset = ((gridIdx[2] + 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
928  + gridIdx[0] + 1)*3;
929  columnOffset = coordColumnOffset*blkSize;
930  entryOffset = (facePlaneOffset + gridIdx[2]*interiorPlaneOffset + edgeStencilLength
931  + gridIdx[0]*faceStencilLength)*blkSize;
932  for(LO l = 0; l < blkSize; ++l) {
933  for(LO k = 0; k < 5; ++k) {
934  for(LO j = 0; j < 3; ++j) {
935  for(LO i = 0; i < 5; ++i) {
936  entries_h(entryOffset + k*15 + j*5 + i + faceStencilLength*l) = columnOffset
937  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + j*lFineNodesPerDim[0] + i - 2)*blkSize + l;
938  values_h(entryOffset + k*15 + j*5 + i + faceStencilLength*l) = coeffs[k]*coeffs[j + 2]*coeffs[i];
939  }
940  }
941  }
942  }
943  for(LO l = 0; l < blkSize; ++l) {
944  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
945  }
946  for(int dim = 0; dim <numDimensions; ++dim) {
947  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
948  }
949 
950  // Face 1
951  coordRowIdx += (lCoarseNodesPerDim[1] - 1)*lCoarseNodesPerDim[0];
952  rowIdx = coordRowIdx*blkSize;
953  coordColumnOffset += (lFineNodesPerDim[1] - 1)*lFineNodesPerDim[0];
954  columnOffset = coordColumnOffset*blkSize;
955  entryOffset += (faceLineOffset + (lCoarseNodesPerDim[1] - 2)*interiorLineOffset)*blkSize;
956  for(LO l = 0; l < blkSize; ++l) {
957  for(LO k = 0; k < 5; ++k) {
958  for(LO j = 0; j < 3; ++j) {
959  for(LO i = 0; i < 5; ++i) {
960  entries_h(entryOffset + k*15 + j*5 + i + faceStencilLength*l) = columnOffset
961  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
962  + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
963  values_h(entryOffset + k*15 + j*5 + i + faceStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
964  }
965  }
966  }
967  }
968  for(LO l = 0; l < blkSize; ++l) {
969  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
970  }
971  for(int dim = 0; dim <numDimensions; ++dim) {
972  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
973  }
974 
975  // Last step in the loop
976  // update the grid indices
977  // for next grid point
978  ++gridIdx[0];
979  if(gridIdx[0] == lCoarseNodesPerDim[0] - 2) {
980  gridIdx[0] = 0;
981  ++gridIdx[2];
982  }
983  }
984  }
985 
986  // Faces in 1-2 plane
987  if((lCoarseNodesPerDim[1] - 2 > 0) && (lCoarseNodesPerDim[2] - 2 > 0)) {
988 
989  Array<LO> gridIdx(3);
990  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
991  for(LO faceIdx=0; faceIdx < (lCoarseNodesPerDim[2]-2)*(lCoarseNodesPerDim[1]-2); ++faceIdx) {
992 
993  // Face 0
994  coordRowIdx = ((gridIdx[2] + 1)*lCoarseNodesPerDim[1]*lCoarseNodesPerDim[0]
995  + (gridIdx[1] + 1)*lCoarseNodesPerDim[0]);
996  rowIdx = coordRowIdx*blkSize;
997  coordColumnOffset = ((gridIdx[2] + 1)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
998  + (gridIdx[1] + 1)*lFineNodesPerDim[0])*3;
999  columnOffset = coordColumnOffset*blkSize;
1000  entryOffset = (facePlaneOffset + gridIdx[2]*interiorPlaneOffset + faceLineOffset
1001  + gridIdx[1]*interiorLineOffset)*blkSize;
1002  for(LO l = 0; l < blkSize; ++l) {
1003  for(LO k = 0; k < 5; ++k) {
1004  for(LO j = 0; j < 5; ++j) {
1005  for(LO i = 0; i < 3; ++i) {
1006  entries_h(entryOffset + k*15 + j*3 + i + faceStencilLength*l) = columnOffset
1007  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0] + (j - 2)*lFineNodesPerDim[0] + i)*blkSize + l;
1008  values_h(entryOffset + k*15 + j*3 + i + faceStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i + 2];
1009  }
1010  }
1011  }
1012  }
1013  for(LO l = 0; l < blkSize; ++l) {
1014  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
1015  }
1016  for(int dim = 0; dim <numDimensions; ++dim) {
1017  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
1018  }
1019 
1020  // Face 1
1021  coordRowIdx += (lCoarseNodesPerDim[0] - 1);
1022  rowIdx = coordRowIdx*blkSize;
1023  coordColumnOffset += (lFineNodesPerDim[0] - 1);
1024  columnOffset = coordColumnOffset*blkSize;
1025  entryOffset += (faceStencilLength + (lCoarseNodesPerDim[0] - 2)*interiorStencilLength)*blkSize;
1026  for(LO l = 0; l < blkSize; ++l) {
1027  for(LO k = 0; k < 5; ++k) {
1028  for(LO j = 0; j < 5; ++j) {
1029  for(LO i = 0; i < 3; ++i) {
1030  entries_h(entryOffset + k*15 + j*3 + i + faceStencilLength*l) = columnOffset
1031  + ((k - 2)*lFineNodesPerDim[1]*lFineNodesPerDim[0]
1032  + (j - 2)*lFineNodesPerDim[0] + i - 2)*blkSize + l;
1033  values_h(entryOffset + k*15 + j*3 + i + faceStencilLength*l) = coeffs[k]*coeffs[j]*coeffs[i];
1034  }
1035  }
1036  }
1037  }
1038  for(LO l = 0; l < blkSize; ++l) {
1039  row_map_h(rowIdx + 1 + l) = entryOffset + faceStencilLength*(l+1);
1040  }
1041  for(int dim = 0; dim <numDimensions; ++dim) {
1042  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
1043  }
1044 
1045  // Last step in the loop
1046  // update the grid indices
1047  // for next grid point
1048  ++gridIdx[1];
1049  if(gridIdx[1] == lCoarseNodesPerDim[1] - 2) {
1050  gridIdx[1] = 0;
1051  ++gridIdx[2];
1052  }
1053  }
1054  }
1055 
1056  if(numInteriors > 0) {
1057  // Allocate and compute arrays
1058  // containing column offsets
1059  // and values associated with
1060  // interior points
1061  LO countRowEntries = 0;
1062  Array<LO> coordColumnOffsets(125);
1063  for(LO k = -2; k < 3; ++k) {
1064  for(LO j = -2; j < 3; ++j) {
1065  for(LO i = -2; i < 3; ++i) {
1066  coordColumnOffsets[countRowEntries] = k*lFineNodesPerDim[1]*lFineNodesPerDim[0]
1067  + j*lFineNodesPerDim[0] + i;
1068  ++countRowEntries;
1069  }
1070  }
1071  }
1072 
1073  LO countValues = 0;
1074  Array<SC> interiorValues(125);
1075  for(LO k = 0; k < 5; ++k) {
1076  for(LO j = 0; j < 5; ++j) {
1077  for(LO i = 0; i < 5; ++i) {
1078  interiorValues[countValues] = coeffs[k]*coeffs[j]*coeffs[i];
1079  ++countValues;
1080  }
1081  }
1082  }
1083 
1084  LO coordRowIdx, rowIdx, coordColumnOffset, columnOffset, entryOffset;
1085  Array<LO> gridIdx(3);
1086  for(LO interiorIdx = 0; interiorIdx < numInteriors; ++interiorIdx) {
1087  coordRowIdx = ((gridIdx[2] + 1)*lCoarseNodesPerDim[0]*lCoarseNodesPerDim[1]
1088  + (gridIdx[1] + 1)*lCoarseNodesPerDim[0]
1089  + gridIdx[0] + 1);
1090  rowIdx = coordRowIdx*blkSize;
1091  coordColumnOffset = ((gridIdx[2] + 1)*3*lFineNodesPerDim[1]*lFineNodesPerDim[0]
1092  + (gridIdx[1] + 1)*3*lFineNodesPerDim[0] + (gridIdx[0] + 1)*3);
1093  columnOffset = coordColumnOffset*blkSize;
1094  entryOffset = (facePlaneOffset + faceLineOffset + faceStencilLength
1095  + gridIdx[2]*interiorPlaneOffset + gridIdx[1]*interiorLineOffset
1096  + gridIdx[0]*interiorStencilLength)*blkSize;
1097  for(LO l = 0; l < blkSize; ++l) {
1098  row_map_h(rowIdx + 1 + l) = entryOffset + interiorStencilLength*(l+1);
1099  }
1100  // Fill the column indices
1101  // and values in the approproate
1102  // views.
1103  for(LO l = 0; l < blkSize; ++l) {
1104  for(LO entryIdx = 0; entryIdx < interiorStencilLength; ++entryIdx) {
1105  entries_h(entryOffset + entryIdx + interiorStencilLength*l) = columnOffset + coordColumnOffsets[entryIdx]*blkSize + l;
1106  values_h(entryOffset + entryIdx + interiorStencilLength*l) = interiorValues[entryIdx];
1107  }
1108  }
1109  for(int dim = 0; dim <numDimensions; ++dim) {
1110  coarseCoordData[dim][coordRowIdx] = fineCoordData[dim][coordColumnOffset];
1111  }
1112 
1113  // Last step in the loop
1114  // update the grid indices
1115  // for next grid point
1116  ++gridIdx[0];
1117  if(gridIdx[0] == lCoarseNodesPerDim[0] - 2) {
1118  gridIdx[0] = 0;
1119  ++gridIdx[1];
1120  if(gridIdx[1] == lCoarseNodesPerDim[1] - 2) {
1121  gridIdx[1] = 0;
1122  ++gridIdx[2];
1123  }
1124  }
1125  }
1126  }
1127 
1128  Kokkos::deep_copy(row_map, row_map_h);
1129  Kokkos::deep_copy(entries, entries_h);
1130  Kokkos::deep_copy(values, values_h);
1131 
1132  local_graph_type localGraph(entries, row_map);
1133  local_matrix_type localR("R", numCols, values, localGraph);
1134 
1135  R = MatrixFactory::Build(localR, // the local data
1136  rowMap, // rowMap
1137  A->getRowMap(), // colMap
1138  A->getRowMap(), // domainMap == colMap
1139  rowMap, // rangeMap == rowMap
1140  Teuchos::null); // params for optimized construction
1141 
1142  } // Build3D()
1143 
1144 } //namespace MueLu
1145 
1146 #define MUELU_REGIONRFACTORY_SHORT
1147 #endif // MUELU_REGIONRFACTORY_DEF_HPP
LocalOrdinal
MueLu::DefaultLocalOrdinal LocalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:50
MueLu::toString
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Definition: MueLu_Utilities_decl.hpp:1306
MueLu::RegionRFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_RegionRFactory_def.hpp:83
MueLu
Namespace for MueLu classes and methods.
Definition: MueLu_BrickAggregationFactory_decl.hpp:78
MueLu::RegionRFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_RegionRFactory_def.hpp:95
Node
MueLu::DefaultNode Node
Definition: MueLu_UseDefaultTypes.hpp:52
MueLu_RegionRFactory_decl.hpp
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:59
MueLu_MasterList.hpp
MueLu::Level::GetLevelID
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
GlobalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Definition: MueLu_UseDefaultTypes.hpp:51
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
MueLu::Utilities::Transpose
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > &params=Teuchos::null)
Definition: MueLu_Utilities_def.hpp:542
MueLu::Level::Set
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Definition: MueLu_Level.hpp:157
MueLu_Monitor.hpp
MueLu::RegionRFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Input.
Definition: MueLu_RegionRFactory_def.hpp:63
MueLu::RegionRFactory::Build3D
void Build3D(const int numDimensions, Array< LO > &lFineNodesPerDim, const RCP< Matrix > &A, const RCP< realvaluedmultivector_type > &fineCoordinates, RCP< Matrix > &R, RCP< realvaluedmultivector_type > &coarseCoordinates, Array< LO > &lCoarseNodesPerDim) const
Definition: MueLu_RegionRFactory_def.hpp:196