MueLu Version of the Day
Loading...
Searching...
No Matches
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
58
59namespace MueLu {
60
61 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
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>
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>
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>
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
MueLu::DefaultLocalOrdinal LocalOrdinal
MueLu::DefaultGlobalOrdinal GlobalOrdinal
MueLu::DefaultNode Node
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
int GetLevelID() const
Return level number.
Definition: MueLu_Level.cpp:76
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const NoFactory * get()
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
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const
Input.
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)
Namespace for MueLu classes and methods.