MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_GeometricInterpolationPFactory_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // MueLu: A package for multigrid based preconditioning
4 //
5 // Copyright 2012 NTESS and the MueLu contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
11 #define MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
12 
13 #include "Xpetra_CrsGraph.hpp"
14 #include "Xpetra_CrsMatrixUtils.hpp"
15 
16 #include "MueLu_MasterList.hpp"
17 #include "MueLu_Monitor.hpp"
18 #include "MueLu_Aggregates.hpp"
19 #include "Xpetra_TpetraCrsMatrix.hpp"
20 
21 // Including this one last ensure that the short names of the above headers are defined properly
22 #include "MueLu_GeometricInterpolationPFactory_decl.hpp"
23 
24 namespace MueLu {
25 
26 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
27 RCP<const ParameterList> GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetValidParameterList() const {
28  RCP<ParameterList> validParamList = rcp(new ParameterList());
29 
30 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
31  SET_VALID_ENTRY("interp: build coarse coordinates");
32 #undef SET_VALID_ENTRY
33 
34  // general variables needed in GeometricInterpolationPFactory
35  validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null,
36  "Generating factory of the matrix A");
37  validParamList->set<RCP<const FactoryBase> >("Aggregates", Teuchos::null,
38  "Aggregates generated by StructuredAggregationFactory used to construct a piece-constant prolongator.");
39  validParamList->set<RCP<const FactoryBase> >("prolongatorGraph", Teuchos::null,
40  "Graph generated by StructuredAggregationFactory used to construct a piece-linear prolongator.");
41  validParamList->set<RCP<const FactoryBase> >("Coordinates", Teuchos::null,
42  "Fine level coordinates used to construct piece-wise linear prolongator and coarse level coordinates.");
43  validParamList->set<RCP<const FactoryBase> >("coarseCoordinatesFineMap", Teuchos::null,
44  "map of the coarse coordinates' GIDs as indexed on the fine mesh.");
45  validParamList->set<RCP<const FactoryBase> >("coarseCoordinatesMap", Teuchos::null,
46  "map of the coarse coordinates' GIDs as indexed on the coarse mesh.");
47  validParamList->set<RCP<const FactoryBase> >("Nullspace", Teuchos::null,
48  "Fine level nullspace used to construct the coarse level nullspace.");
49  validParamList->set<RCP<const FactoryBase> >("numDimensions", Teuchos::null,
50  "Number of spacial dimensions in the problem.");
51  validParamList->set<RCP<const FactoryBase> >("lCoarseNodesPerDim", Teuchos::null,
52  "Number of nodes per spatial dimension on the coarse grid.");
53  validParamList->set<RCP<const FactoryBase> >("structuredInterpolationOrder", Teuchos::null,
54  "Interpolation order for constructing the prolongator.");
55  validParamList->set<bool>("keep coarse coords", false, "Flag to keep coordinates for special coarse grid solve");
56  validParamList->set<bool>("interp: remove small entries", true, "Remove small interpolation coeficient from prolongator to reduce fill-in on coarse level");
57 
58  return validParamList;
59 }
60 
61 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
62 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
63  DeclareInput(Level& fineLevel, Level& /* coarseLevel */) const {
64  const ParameterList& pL = GetParameterList();
65 
66  Input(fineLevel, "A");
67  Input(fineLevel, "Nullspace");
68  Input(fineLevel, "numDimensions");
69  Input(fineLevel, "prolongatorGraph");
70  Input(fineLevel, "lCoarseNodesPerDim");
71  Input(fineLevel, "structuredInterpolationOrder");
72 
73  if (pL.get<bool>("interp: build coarse coordinates") ||
74  Get<int>(fineLevel, "structuredInterpolationOrder") == 1) {
75  Input(fineLevel, "Coordinates");
76  Input(fineLevel, "coarseCoordinatesFineMap");
77  Input(fineLevel, "coarseCoordinatesMap");
78  }
79 }
80 
81 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
82 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
83  Build(Level& fineLevel, Level& coarseLevel) const {
84  return BuildP(fineLevel, coarseLevel);
85 }
86 
87 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
88 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
89  BuildP(Level& fineLevel, Level& coarseLevel) const {
90  FactoryMonitor m(*this, "BuildP", coarseLevel);
91 
92  // Set debug outputs based on environment variable
93  RCP<Teuchos::FancyOStream> out;
94  if (const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
95  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
96  out->setShowAllFrontMatter(false).setShowProcRank(true);
97  } else {
98  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
99  }
100 
101  *out << "Starting GeometricInterpolationPFactory::BuildP." << std::endl;
102 
103  // Get inputs from the parameter list
104  const ParameterList& pL = GetParameterList();
105  const bool removeSmallEntries = pL.get<bool>("interp: remove small entries");
106  const bool buildCoarseCoordinates = pL.get<bool>("interp: build coarse coordinates");
107  const int interpolationOrder = Get<int>(fineLevel, "structuredInterpolationOrder");
108  const int numDimensions = Get<int>(fineLevel, "numDimensions");
109 
110  // Declared main input/outputs to be retrieved and placed on the fine resp. coarse level
111  RCP<Matrix> A = Get<RCP<Matrix> >(fineLevel, "A");
112  RCP<const CrsGraph> prolongatorGraph = Get<RCP<CrsGraph> >(fineLevel, "prolongatorGraph");
113  RCP<realvaluedmultivector_type> fineCoordinates, coarseCoordinates;
114  RCP<Matrix> P;
115 
116  // Check if we need to build coarse coordinates as they are used if we construct
117  // a linear interpolation prolongator
118  if (buildCoarseCoordinates || (interpolationOrder == 1)) {
119  SubFactoryMonitor sfm(*this, "BuildCoordinates", coarseLevel);
120  RCP<const Map> coarseCoordsFineMap = Get<RCP<const Map> >(fineLevel, "coarseCoordinatesFineMap");
121  RCP<const Map> coarseCoordsMap = Get<RCP<const Map> >(fineLevel, "coarseCoordinatesMap");
122 
123  fineCoordinates = Get<RCP<realvaluedmultivector_type> >(fineLevel, "Coordinates");
124  coarseCoordinates = Xpetra::MultiVectorFactory<real_type, LO, GO, Node>::Build(coarseCoordsFineMap,
125  fineCoordinates->getNumVectors());
126  RCP<const Import> coordsImporter = ImportFactory::Build(fineCoordinates->getMap(),
127  coarseCoordsFineMap);
128  coarseCoordinates->doImport(*fineCoordinates, *coordsImporter, Xpetra::INSERT);
129  coarseCoordinates->replaceMap(coarseCoordsMap);
130 
131  Set(coarseLevel, "Coordinates", coarseCoordinates);
132 
133  if (pL.get<bool>("keep coarse coords")) {
134  coarseLevel.Set<RCP<realvaluedmultivector_type> >("Coordinates2", coarseCoordinates, NoFactory::get());
135  }
136  }
137 
138  *out << "Fine and coarse coordinates have been loaded from the fine level and set on the coarse level." << std::endl;
139 
140  if (interpolationOrder == 0) {
141  SubFactoryMonitor sfm(*this, "BuildConstantP", coarseLevel);
142  // Compute the prolongator using piece-wise constant interpolation
143  BuildConstantP(P, prolongatorGraph, A);
144  } else if (interpolationOrder == 1) {
145  // Compute the prolongator using piece-wise linear interpolation
146  // First get all the required coordinates to compute the local part of P
147  RCP<const Map> prolongatorColMap = prolongatorGraph->getColMap();
148 
149  const size_t dofsPerNode = static_cast<size_t>(A->GetFixedBlockSize());
150  const size_t numColIndices = prolongatorColMap->getLocalNumElements();
151  TEUCHOS_TEST_FOR_EXCEPTION((numColIndices % dofsPerNode) != 0,
152  Exceptions::RuntimeError,
153  "Something went wrong, the number of columns in the prolongator is not a multiple of dofsPerNode!");
154  const size_t numGhostCoords = numColIndices / dofsPerNode;
155  const GO indexBase = prolongatorColMap->getIndexBase();
156  const GO coordIndexBase = fineCoordinates->getMap()->getIndexBase();
157 
158  ArrayView<const GO> prolongatorColIndices = prolongatorColMap->getLocalElementList();
159  Array<GO> ghostCoordIndices(numGhostCoords);
160  for (size_t ghostCoordIdx = 0; ghostCoordIdx < numGhostCoords; ++ghostCoordIdx) {
161  ghostCoordIndices[ghostCoordIdx] = (prolongatorColIndices[ghostCoordIdx * dofsPerNode] - indexBase) / dofsPerNode + coordIndexBase;
162  }
163  RCP<Map> ghostCoordMap = MapFactory::Build(fineCoordinates->getMap()->lib(),
164  prolongatorColMap->getGlobalNumElements() / dofsPerNode,
165  ghostCoordIndices(),
166  coordIndexBase,
167  fineCoordinates->getMap()->getComm());
168 
169  RCP<realvaluedmultivector_type> ghostCoordinates = Xpetra::MultiVectorFactory<real_type, LO, GO, NO>::Build(ghostCoordMap,
170  fineCoordinates->getNumVectors());
171  RCP<const Import> ghostImporter = ImportFactory::Build(coarseCoordinates->getMap(),
172  ghostCoordMap);
173  ghostCoordinates->doImport(*coarseCoordinates, *ghostImporter, Xpetra::INSERT);
174 
175  BuildLinearP(coarseLevel, A, prolongatorGraph, fineCoordinates,
176  ghostCoordinates, numDimensions, removeSmallEntries, P);
177  }
178 
179  *out << "The prolongator matrix has been built." << std::endl;
180 
181  {
182  SubFactoryMonitor sfm(*this, "BuildNullspace", coarseLevel);
183  // Build the coarse nullspace
184  RCP<MultiVector> fineNullspace = Get<RCP<MultiVector> >(fineLevel, "Nullspace");
185  RCP<MultiVector> coarseNullspace = MultiVectorFactory::Build(P->getDomainMap(),
186  fineNullspace->getNumVectors());
187 
188  using helpers = Xpetra::Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
189  if (helpers::isTpetraBlockCrs(A)) {
190  // FIXME: BlockCrs doesn't currently support transpose apply, so we have to do this the hard way
191  RCP<Matrix> Ptrans = Utilities::Transpose(*P);
192  Ptrans->apply(*fineNullspace, *coarseNullspace, Teuchos::NO_TRANS, Teuchos::ScalarTraits<SC>::one(),
193  Teuchos::ScalarTraits<SC>::zero());
194  } else {
195  P->apply(*fineNullspace, *coarseNullspace, Teuchos::TRANS, Teuchos::ScalarTraits<SC>::one(),
196  Teuchos::ScalarTraits<SC>::zero());
197  }
198 
199  Set(coarseLevel, "Nullspace", coarseNullspace);
200  }
201 
202  *out << "The coarse nullspace is constructed and set on the coarse level." << std::endl;
203 
204  Set(coarseLevel, "P", P);
205 
206  *out << "GeometricInterpolationPFactory::BuildP has completed." << std::endl;
207 
208 } // BuildP
209 
210 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
211 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
212  BuildConstantP(RCP<Matrix>& P, RCP<const CrsGraph>& prolongatorGraph, RCP<Matrix>& A) const {
213  // Set debug outputs based on environment variable
214  RCP<Teuchos::FancyOStream> out;
215  if (const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
216  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
217  out->setShowAllFrontMatter(false).setShowProcRank(true);
218  } else {
219  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
220  }
221 
222  *out << "BuildConstantP" << std::endl;
223 
224  std::vector<size_t> strideInfo(1);
225  strideInfo[0] = A->GetFixedBlockSize();
226  RCP<const StridedMap> stridedDomainMap =
227  StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
228 
229  *out << "Call prolongator constructor" << std::endl;
230 
231  using helpers = Xpetra::Helpers<Scalar, LocalOrdinal, GlobalOrdinal, Node>;
232  if (helpers::isTpetraBlockCrs(A)) {
233  SC one = Teuchos::ScalarTraits<SC>::one();
234  SC zero = Teuchos::ScalarTraits<SC>::zero();
235  LO NSDim = A->GetStorageBlockSize();
236 
237  // Build the exploded Map
238  RCP<const Map> BlockMap = prolongatorGraph->getDomainMap();
239  Teuchos::ArrayView<const GO> block_dofs = BlockMap->getLocalElementList();
240  Teuchos::Array<GO> point_dofs(block_dofs.size() * NSDim);
241  for (LO i = 0, ct = 0; i < block_dofs.size(); i++) {
242  for (LO j = 0; j < NSDim; j++) {
243  point_dofs[ct] = block_dofs[i] * NSDim + j;
244  ct++;
245  }
246  }
247 
248  RCP<const Map> PointMap = MapFactory::Build(BlockMap->lib(),
249  BlockMap->getGlobalNumElements() * NSDim,
250  point_dofs(),
251  BlockMap->getIndexBase(),
252  BlockMap->getComm());
253  strideInfo[0] = A->GetFixedBlockSize();
254  RCP<const StridedMap> stridedPointMap = StridedMapFactory::Build(PointMap, strideInfo);
255 
256  RCP<Xpetra::CrsMatrix<SC, LO, GO, NO> > P_xpetra = Xpetra::CrsMatrixFactory<SC, LO, GO, NO>::BuildBlock(prolongatorGraph, PointMap, A->getRangeMap(), NSDim);
257  RCP<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> > P_tpetra = rcp_dynamic_cast<Xpetra::TpetraBlockCrsMatrix<SC, LO, GO, NO> >(P_xpetra);
258  if (P_tpetra.is_null()) throw std::runtime_error("BuildConstantP Matrix factory did not return a Tpetra::BlockCrsMatrix");
259  RCP<CrsMatrixWrap> P_wrap = rcp(new CrsMatrixWrap(P_xpetra));
260 
261  // NOTE: Assumes block-diagonal prolongation
262  Teuchos::Array<LO> temp(1);
263  Teuchos::ArrayView<const LO> indices;
264  Teuchos::Array<Scalar> block(NSDim * NSDim, zero);
265  for (LO i = 0; i < NSDim; i++)
266  block[i * NSDim + i] = one;
267  for (LO i = 0; i < (int)prolongatorGraph->getLocalNumRows(); i++) {
268  prolongatorGraph->getLocalRowView(i, indices);
269  for (LO j = 0; j < (LO)indices.size(); j++) {
270  temp[0] = indices[j];
271  P_tpetra->replaceLocalValues(i, temp(), block());
272  }
273  }
274 
275  P = P_wrap;
276  if (A->IsView("stridedMaps") == true) {
277  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedPointMap);
278  } else {
279  P->CreateView("stridedMaps", P->getRangeMap(), PointMap);
280  }
281  } else {
282  // Create the prolongator matrix and its associated objects
283  RCP<ParameterList> dummyList = rcp(new ParameterList());
284  P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
285  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
286  PCrs->setAllToScalar(1.0);
287  PCrs->fillComplete();
288 
289  // set StridingInformation of P
290  if (A->IsView("stridedMaps") == true)
291  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
292  else
293  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
294  }
295 
296 } // BuildConstantP
297 
298 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
299 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
300  BuildLinearP(Level& coarseLevel,
301  RCP<Matrix>& A, RCP<const CrsGraph>& prolongatorGraph,
302  RCP<realvaluedmultivector_type>& fineCoordinates,
303  RCP<realvaluedmultivector_type>& ghostCoordinates,
304  const int numDimensions, const bool removeSmallEntries,
305  RCP<Matrix>& P) const {
306  // Set debug outputs based on environment variable
307  RCP<Teuchos::FancyOStream> out;
308  if (const char* dbg = std::getenv("MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEBUG")) {
309  out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
310  out->setShowAllFrontMatter(false).setShowProcRank(true);
311  } else {
312  out = Teuchos::getFancyOStream(rcp(new Teuchos::oblackholestream()));
313  }
314 
315  // Compute 2^numDimensions using bit logic to avoid round-off errors
316  const int numInterpolationPoints = 1 << numDimensions;
317  const int dofsPerNode = A->GetFixedBlockSize() / A->GetStorageBlockSize();
318  ;
319 
320  RCP<ParameterList> dummyList = rcp(new ParameterList());
321  P = rcp(new CrsMatrixWrap(prolongatorGraph, dummyList));
322  RCP<CrsMatrix> PCrs = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
323  PCrs->resumeFill(); // The Epetra matrix is considered filled at this point.
324 
325  {
326  *out << "Entering BuildLinearP" << std::endl;
327  SubFactoryMonitor sfm(*this, "BuildLinearP", coarseLevel);
328 
329  // Extract coordinates for interpolation stencil calculations
330  const LO numFineNodes = fineCoordinates->getLocalLength();
331  const LO numGhostNodes = ghostCoordinates->getLocalLength();
332  Array<ArrayRCP<const real_type> > fineCoords(3);
333  Array<ArrayRCP<const real_type> > ghostCoords(3);
334  const real_type realZero = Teuchos::as<real_type>(0.0);
335  ArrayRCP<real_type> fineZero(numFineNodes, realZero);
336  ArrayRCP<real_type> ghostZero(numGhostNodes, realZero);
337  for (int dim = 0; dim < 3; ++dim) {
338  if (dim < numDimensions) {
339  fineCoords[dim] = fineCoordinates->getData(dim);
340  ghostCoords[dim] = ghostCoordinates->getData(dim);
341  } else {
342  fineCoords[dim] = fineZero;
343  ghostCoords[dim] = ghostZero;
344  }
345  }
346 
347  *out << "Coordinates extracted from the multivectors!" << std::endl;
348 
349  { // Construct the linear interpolation prolongator
350  LO interpolationNodeIdx = 0, rowIdx = 0;
351  ArrayView<const LO> colIndices;
352  Array<SC> values;
353  Array<Array<real_type> > coords(numInterpolationPoints + 1);
354  Array<real_type> stencil(numInterpolationPoints);
355  for (LO nodeIdx = 0; nodeIdx < numFineNodes; ++nodeIdx) {
356  if (PCrs->getNumEntriesInLocalRow(nodeIdx * dofsPerNode) == 1) {
357  values.resize(1);
358  values[0] = 1.0;
359  for (LO dof = 0; dof < dofsPerNode; ++dof) {
360  rowIdx = nodeIdx * dofsPerNode + dof;
361  prolongatorGraph->getLocalRowView(rowIdx, colIndices);
362  PCrs->replaceLocalValues(rowIdx, colIndices, values());
363  }
364  } else {
365  // Extract the coordinates associated with the current node
366  // and the neighboring coarse nodes
367  coords[0].resize(3);
368  for (int dim = 0; dim < 3; ++dim) {
369  coords[0][dim] = fineCoords[dim][nodeIdx];
370  }
371  prolongatorGraph->getLocalRowView(nodeIdx * dofsPerNode, colIndices);
372  for (int interpolationIdx = 0; interpolationIdx < numInterpolationPoints; ++interpolationIdx) {
373  coords[interpolationIdx + 1].resize(3);
374  interpolationNodeIdx = colIndices[interpolationIdx] / dofsPerNode;
375  for (int dim = 0; dim < 3; ++dim) {
376  coords[interpolationIdx + 1][dim] = ghostCoords[dim][interpolationNodeIdx];
377  }
378  }
379  RCP<Teuchos::TimeMonitor> tm = rcp(new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer("Compute Linear Interpolation")));
380  ComputeLinearInterpolationStencil(numDimensions, numInterpolationPoints, coords, stencil);
381  tm = Teuchos::null;
382  values.resize(numInterpolationPoints);
383  for (LO valueIdx = 0; valueIdx < numInterpolationPoints; ++valueIdx) {
384  values[valueIdx] = Teuchos::as<SC>(stencil[valueIdx]);
385  }
386 
387  // Set values in all the rows corresponding to nodeIdx
388  for (LO dof = 0; dof < dofsPerNode; ++dof) {
389  rowIdx = nodeIdx * dofsPerNode + dof;
390  prolongatorGraph->getLocalRowView(rowIdx, colIndices);
391  PCrs->replaceLocalValues(rowIdx, colIndices, values());
392  }
393  } // Check for Coarse vs. Fine point
394  } // Loop over fine nodes
395  } // Construct the linear interpolation prolongator
396 
397  *out << "The calculation of the interpolation stencils has completed." << std::endl;
398 
399  PCrs->fillComplete();
400  }
401 
402  *out << "All values in P have been set and fillComplete has been performed." << std::endl;
403 
404  // Note lbv Jan 29 2019: this should be handle at aggregation level
405  // if the user really does not want potential d2 neighbors on coarse grid
406  // that way we would avoid a new graph construction...
407 
408  // Check if we want to remove small entries from P
409  // to reduce stencil growth on next level.
410  if (removeSmallEntries) {
411  *out << "Entering remove small entries" << std::endl;
412  SubFactoryMonitor sfm(*this, "remove small entries", coarseLevel);
413 
414  ArrayRCP<const size_t> rowptrOrig;
415  ArrayRCP<const LO> colindOrig;
416  ArrayRCP<const Scalar> valuesOrig;
417  PCrs->getAllValues(rowptrOrig, colindOrig, valuesOrig);
418 
419  const size_t numRows = static_cast<size_t>(rowptrOrig.size() - 1);
420  ArrayRCP<size_t> rowPtr(numRows + 1);
421  ArrayRCP<size_t> nnzOnRows(numRows);
422  rowPtr[0] = 0;
423  size_t countRemovedEntries = 0;
424  for (size_t rowIdx = 0; rowIdx < numRows; ++rowIdx) {
425  for (size_t entryIdx = rowptrOrig[rowIdx]; entryIdx < rowptrOrig[rowIdx + 1]; ++entryIdx) {
426  if (Teuchos::ScalarTraits<Scalar>::magnitude(valuesOrig[entryIdx]) < 1e-6) {
427  ++countRemovedEntries;
428  }
429  }
430  rowPtr[rowIdx + 1] = rowptrOrig[rowIdx + 1] - countRemovedEntries;
431  nnzOnRows[rowIdx] = rowPtr[rowIdx + 1] - rowPtr[rowIdx];
432  }
433  GetOStream(Statistics1) << "interp: number of small entries removed= " << countRemovedEntries << " / " << rowptrOrig[numRows] << std::endl;
434 
435  size_t countKeptEntries = 0;
436  ArrayRCP<LO> colInd(rowPtr[numRows]);
437  ArrayRCP<SC> values(rowPtr[numRows]);
438  for (size_t entryIdx = 0; entryIdx < rowptrOrig[numRows]; ++entryIdx) {
439  if (Teuchos::ScalarTraits<Scalar>::magnitude(valuesOrig[entryIdx]) > 1e-6) {
440  colInd[countKeptEntries] = colindOrig[entryIdx];
441  values[countKeptEntries] = valuesOrig[entryIdx];
442  ++countKeptEntries;
443  }
444  }
445 
446  P = rcp(new CrsMatrixWrap(prolongatorGraph->getRowMap(),
447  prolongatorGraph->getColMap(),
448  nnzOnRows));
449  RCP<CrsMatrix> PCrsSqueezed = rcp_dynamic_cast<CrsMatrixWrap>(P)->getCrsMatrix();
450  PCrsSqueezed->resumeFill(); // The Epetra matrix is considered filled at this point.
451  PCrsSqueezed->setAllValues(rowPtr, colInd, values);
452  PCrsSqueezed->expertStaticFillComplete(prolongatorGraph->getDomainMap(),
453  prolongatorGraph->getRangeMap());
454  }
455 
456  std::vector<size_t> strideInfo(1);
457  strideInfo[0] = dofsPerNode;
458  RCP<const StridedMap> stridedDomainMap =
459  StridedMapFactory::Build(prolongatorGraph->getDomainMap(), strideInfo);
460 
461  *out << "The strided maps of P have been computed" << std::endl;
462 
463  // set StridingInformation of P
464  if (A->IsView("stridedMaps") == true) {
465  P->CreateView("stridedMaps", A->getRowMap("stridedMaps"), stridedDomainMap);
466  } else {
467  P->CreateView("stridedMaps", P->getRangeMap(), stridedDomainMap);
468  }
469 
470 } // BuildLinearP
471 
472 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
473 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
474  ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints,
475  const Array<Array<real_type> > coord,
476  Array<real_type>& stencil) const {
477  // 7 8 Find xi, eta and zeta such that
478  // x---------x
479  // /| /| Rx = x_p - sum N_i(xi,eta,zeta)x_i = 0
480  // 5/ | 6/ | Ry = y_p - sum N_i(xi,eta,zeta)y_i = 0
481  // x---------x | Rz = z_p - sum N_i(xi,eta,zeta)z_i = 0
482  // | | *P | |
483  // | x------|--x We can do this with a Newton solver:
484  // | /3 | /4 We will start with initial guess (xi,eta,zeta) = (0,0,0)
485  // |/ |/ We compute the Jacobian and iterate until convergence...
486  // z y x---------x
487  // | / 1 2 Once we have (xi,eta,zeta), we can evaluate all N_i which
488  // |/ give us the weights for the interpolation stencil!
489  // o---x
490  //
491 
492  Teuchos::SerialDenseMatrix<LO, real_type> Jacobian(numDimensions, numDimensions);
493  Teuchos::SerialDenseVector<LO, real_type> residual(numDimensions);
494  Teuchos::SerialDenseVector<LO, real_type> solutionDirection(numDimensions);
495  Teuchos::SerialDenseVector<LO, real_type> paramCoords(numDimensions);
496  Teuchos::SerialDenseSolver<LO, real_type> problem;
497  int iter = 0, max_iter = 5;
498  real_type functions[4][8], norm_ref = 1.0, norm2 = 1.0, tol = 1.0e-5;
499  paramCoords.size(numDimensions);
500 
501  while ((iter < max_iter) && (norm2 > tol * norm_ref)) {
502  ++iter;
503  norm2 = 0.0;
504  solutionDirection.size(numDimensions);
505  residual.size(numDimensions);
506  Jacobian = 0.0;
507 
508  // Compute Jacobian and Residual
509  GetInterpolationFunctions(numDimensions, paramCoords, functions);
510  for (LO i = 0; i < numDimensions; ++i) {
511  residual(i) = coord[0][i]; // Add coordinates from point of interest
512  for (LO k = 0; k < numInterpolationPoints; ++k) {
513  residual(i) -= functions[0][k] * coord[k + 1][i]; // Remove contribution from all coarse points
514  }
515  if (iter == 1) {
516  norm_ref += residual(i) * residual(i);
517  if (i == numDimensions - 1) {
518  norm_ref = std::sqrt(norm_ref);
519  }
520  }
521 
522  for (LO j = 0; j < numDimensions; ++j) {
523  for (LO k = 0; k < numInterpolationPoints; ++k) {
524  Jacobian(i, j) += functions[j + 1][k] * coord[k + 1][i];
525  }
526  }
527  }
528 
529  // Set Jacobian, Vectors and solve problem
530  problem.setMatrix(Teuchos::rcp(&Jacobian, false));
531  problem.setVectors(Teuchos::rcp(&solutionDirection, false), Teuchos::rcp(&residual, false));
532  if (problem.shouldEquilibrate()) {
533  problem.factorWithEquilibration(true);
534  }
535  problem.solve();
536 
537  for (LO i = 0; i < numDimensions; ++i) {
538  paramCoords(i) = paramCoords(i) + solutionDirection(i);
539  }
540 
541  // Recompute Residual norm
542  GetInterpolationFunctions(numDimensions, paramCoords, functions);
543  for (LO i = 0; i < numDimensions; ++i) {
544  real_type tmp = coord[0][i];
545  for (LO k = 0; k < numInterpolationPoints; ++k) {
546  tmp -= functions[0][k] * coord[k + 1][i];
547  }
548  norm2 += tmp * tmp;
549  tmp = 0.0;
550  }
551  norm2 = std::sqrt(norm2);
552  }
553 
554  // Load the interpolation values onto the stencil.
555  for (LO i = 0; i < numInterpolationPoints; ++i) {
556  stencil[i] = functions[0][i];
557  }
558 
559 } // End ComputeLinearInterpolationStencil
560 
561 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
562 void GeometricInterpolationPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::
563  GetInterpolationFunctions(const LO numDimensions,
564  const Teuchos::SerialDenseVector<LO, real_type> parametricCoordinates,
565  real_type functions[4][8]) const {
566  real_type xi = 0.0, eta = 0.0, zeta = 0.0, denominator = 0.0;
567  if (numDimensions == 1) {
568  xi = parametricCoordinates[0];
569  denominator = 2.0;
570  } else if (numDimensions == 2) {
571  xi = parametricCoordinates[0];
572  eta = parametricCoordinates[1];
573  denominator = 4.0;
574  } else if (numDimensions == 3) {
575  xi = parametricCoordinates[0];
576  eta = parametricCoordinates[1];
577  zeta = parametricCoordinates[2];
578  denominator = 8.0;
579  }
580 
581  functions[0][0] = (1.0 - xi) * (1.0 - eta) * (1.0 - zeta) / denominator;
582  functions[0][1] = (1.0 + xi) * (1.0 - eta) * (1.0 - zeta) / denominator;
583  functions[0][2] = (1.0 - xi) * (1.0 + eta) * (1.0 - zeta) / denominator;
584  functions[0][3] = (1.0 + xi) * (1.0 + eta) * (1.0 - zeta) / denominator;
585  functions[0][4] = (1.0 - xi) * (1.0 - eta) * (1.0 + zeta) / denominator;
586  functions[0][5] = (1.0 + xi) * (1.0 - eta) * (1.0 + zeta) / denominator;
587  functions[0][6] = (1.0 - xi) * (1.0 + eta) * (1.0 + zeta) / denominator;
588  functions[0][7] = (1.0 + xi) * (1.0 + eta) * (1.0 + zeta) / denominator;
589 
590  functions[1][0] = -(1.0 - eta) * (1.0 - zeta) / denominator;
591  functions[1][1] = (1.0 - eta) * (1.0 - zeta) / denominator;
592  functions[1][2] = -(1.0 + eta) * (1.0 - zeta) / denominator;
593  functions[1][3] = (1.0 + eta) * (1.0 - zeta) / denominator;
594  functions[1][4] = -(1.0 - eta) * (1.0 + zeta) / denominator;
595  functions[1][5] = (1.0 - eta) * (1.0 + zeta) / denominator;
596  functions[1][6] = -(1.0 + eta) * (1.0 + zeta) / denominator;
597  functions[1][7] = (1.0 + eta) * (1.0 + zeta) / denominator;
598 
599  functions[2][0] = -(1.0 - xi) * (1.0 - zeta) / denominator;
600  functions[2][1] = -(1.0 + xi) * (1.0 - zeta) / denominator;
601  functions[2][2] = (1.0 - xi) * (1.0 - zeta) / denominator;
602  functions[2][3] = (1.0 + xi) * (1.0 - zeta) / denominator;
603  functions[2][4] = -(1.0 - xi) * (1.0 + zeta) / denominator;
604  functions[2][5] = -(1.0 + xi) * (1.0 + zeta) / denominator;
605  functions[2][6] = (1.0 - xi) * (1.0 + zeta) / denominator;
606  functions[2][7] = (1.0 + xi) * (1.0 + zeta) / denominator;
607 
608  functions[3][0] = -(1.0 - xi) * (1.0 - eta) / denominator;
609  functions[3][1] = -(1.0 + xi) * (1.0 - eta) / denominator;
610  functions[3][2] = -(1.0 - xi) * (1.0 + eta) / denominator;
611  functions[3][3] = -(1.0 + xi) * (1.0 + eta) / denominator;
612  functions[3][4] = (1.0 - xi) * (1.0 - eta) / denominator;
613  functions[3][5] = (1.0 + xi) * (1.0 - eta) / denominator;
614  functions[3][6] = (1.0 - xi) * (1.0 + eta) / denominator;
615  functions[3][7] = (1.0 + xi) * (1.0 + eta) / denominator;
616 
617 } // End GetInterpolationFunctions
618 
619 } // namespace MueLu
620 
621 #endif // MUELU_GEOMETRICINTERPOLATIONPFACTORY_DEF_HPP
Xpetra::MultiVectorFactory::Build
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
MueLu::removeSmallEntries
Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > removeSmallEntries(Teuchos::RCP< Xpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType threshold, const bool keepDiagonal)
Definition: MueLu_UtilitiesBase_def.hpp:56
Teuchos::basic_FancyOStream::setShowProcRank
basic_FancyOStream & setShowProcRank(const bool showProcRank)
GO
GlobalOrdinal GO
Teuchos::ParameterList::get
T & get(const std::string &name, T def_value)
MueLu::FactoryMonitor
Timer to be used in factories. Similar to Monitor but with additional timers.
Definition: MueLu_Monitor.hpp:165
TEUCHOS_TEST_FOR_EXCEPTION
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Teuchos::basic_oblackholestream
MueLu::Statistics1
Print more statistics.
Definition: MueLu_VerbosityLevel.hpp:34
Teuchos::ArrayView::size
size_type size() const
MueLu::GeometricInterpolationPFactory::BuildLinearP
void BuildLinearP(Level &coarseLevel, RCP< Matrix > &A, RCP< const CrsGraph > &prolongatorGraph, RCP< realvaluedmultivector_type > &fineCoordinates, RCP< realvaluedmultivector_type > &ghostCoordinates, const int numDimensions, const bool keepD2, RCP< Matrix > &P) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:300
LO
LocalOrdinal LO
MueLu_GeometricInterpolationPFactory_decl.hpp
Teuchos::ParameterList::set
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Teuchos::ArrayRCP::size
size_type size() const
MueLu::NoFactory::get
static const NoFactory * get()
Definition: MueLu_NoFactory.cpp:34
MueLu::GeometricInterpolationPFactory::GetInterpolationFunctions
void GetInterpolationFunctions(const LO numDimensions, const Teuchos::SerialDenseVector< LO, real_type > parametricCoordinates, real_type functions[4][8]) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:563
Xpetra_CrsGraph.hpp
MueLu_MasterList.hpp
Teuchos::TimeMonitor::getNewTimer
static RCP< Time > getNewTimer(const std::string &name)
rcp
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::GeometricInterpolationPFactory::Build
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:83
MueLu::Level
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:63
MueLu::GeometricInterpolationPFactory::ComputeLinearInterpolationStencil
void ComputeLinearInterpolationStencil(const int numDimensions, const int numInterpolationPoints, const Array< Array< real_type > > coord, Array< real_type > &stencil) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:474
MueLu::SubFactoryMonitor
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
Definition: MueLu_Monitor.hpp:218
Teuchos::SerialDenseSolver::factorWithEquilibration
void factorWithEquilibration(bool flag)
SET_VALID_ENTRY
#define SET_VALID_ENTRY(name)
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:496
Teuchos::basic_FancyOStream::setShowAllFrontMatter
basic_FancyOStream & setShowAllFrontMatter(const bool showAllFrontMatter)
Teuchos::Array::resize
void resize(size_type new_size, const value_type &x=value_type())
MueLu::Level::Set
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Definition: MueLu_Level.hpp:120
MueLu::GeometricInterpolationPFactory::DeclareInput
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:63
iter
TransListIter iter
Teuchos::SerialDenseVector::size
int size(OrdinalType length_in)
Teuchos::SerialDenseSolver::setVectors
int setVectors(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &X, const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &B)
MueLu::GeometricInterpolationPFactory::GetValidParameterList
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:27
SC
Scalar SC
MueLu_Monitor.hpp
residual
void residual(const Operator< SC, LO, GO, NO > &Aop, const MultiVector< SC, LO, GO, NO > &X_in, const MultiVector< SC, LO, GO, NO > &B_in, MultiVector< SC, LO, GO, NO > &R_in)
MueLu::Exceptions::RuntimeError
Exception throws to report errors in the internal logical of the program.
Definition: MueLu_Exceptions.hpp:36
MueLu::GeometricInterpolationPFactory::BuildConstantP
void BuildConstantP(RCP< Matrix > &P, RCP< const CrsGraph > &prolongatorGraph, RCP< Matrix > &A) const
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:212
MueLu::GeometricInterpolationPFactory::real_type
typename Teuchos::ScalarTraits< SC >::coordinateType real_type
Definition: MueLu_GeometricInterpolationPFactory_decl.hpp:33
Xpetra_CrsMatrixUtils.hpp
MueLu::GeometricInterpolationPFactory::BuildP
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Definition: MueLu_GeometricInterpolationPFactory_def.hpp:89
Teuchos::SerialDenseSolver::setMatrix
int setMatrix(const RCP< SerialDenseMatrix< OrdinalType, ScalarType > > &A)