MueLu  Version of the Day
MueLu_BlockedPFactory_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 
47 #ifndef MUELU_BLOCKEDPFACTORY_DEF_HPP_
48 #define MUELU_BLOCKEDPFACTORY_DEF_HPP_
49 
50 #include <Xpetra_MultiVectorFactory.hpp>
51 #include <Xpetra_VectorFactory.hpp>
52 #include <Xpetra_ImportFactory.hpp>
53 #include <Xpetra_ExportFactory.hpp>
54 #include <Xpetra_CrsMatrixWrap.hpp>
55 
56 #include <Xpetra_BlockedCrsMatrix.hpp>
57 #include <Xpetra_Map.hpp>
58 #include <Xpetra_MapFactory.hpp>
59 #include <Xpetra_MapExtractor.hpp>
60 #include <Xpetra_MapExtractorFactory.hpp>
61 
63 #include "MueLu_TentativePFactory.hpp"
64 #include "MueLu_FactoryBase.hpp"
65 #include "MueLu_FactoryManager.hpp"
66 #include "MueLu_Monitor.hpp"
67 #include "MueLu_HierarchyUtils.hpp"
68 
69 namespace MueLu {
70 
71  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
73  RCP<ParameterList> validParamList = rcp(new ParameterList());
74 
75  validParamList->set< RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A (block matrix)");
76  validParamList->set< bool > ("backwards", false, "Forward/backward order");
77 
78  return validParamList;
79  }
80 
81  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
83  Input(fineLevel, "A");
84 
85  const ParameterList& pL = GetParameterList();
86  const bool backwards = pL.get<bool>("backwards");
87 
88  const int numFactManagers = FactManager_.size();
89  for (int k = 0; k < numFactManagers; k++) {
90  int i = (backwards ? numFactManagers-1 - k : k);
91  const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
92 
93  SetFactoryManager fineSFM (rcpFromRef(fineLevel), factManager);
94  SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), factManager);
95 
96  if (!restrictionMode_)
97  coarseLevel.DeclareInput("P", factManager->GetFactory("P").get(), this);
98  else
99  coarseLevel.DeclareInput("R", factManager->GetFactory("R").get(), this);
100  }
101  }
102 
103  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
104  void BlockedPFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::BuildP(Level& /* fineLevel */, Level& /* coarseLevel */) const
105  { }
106 
107  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
109  FactoryMonitor m(*this, "Build", coarseLevel);
110 
111  RCP<Matrix> Ain = Get< RCP<Matrix> >(fineLevel, "A");
112 
113  RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
114  TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::BadCast, "Input matrix A is not a BlockedCrsMatrix.");
115 
116  const int numFactManagers = FactManager_.size();
117 
118  // Plausibility check
119  TEUCHOS_TEST_FOR_EXCEPTION(A->Rows() != as<size_t>(numFactManagers), Exceptions::RuntimeError,
120  "Number of block rows [" << A->Rows() << "] does not match the number of SubFactorManagers [" << numFactManagers << "]");
121  TEUCHOS_TEST_FOR_EXCEPTION(A->Cols() != as<size_t>(numFactManagers), Exceptions::RuntimeError,
122  "Number of block cols [" << A->Cols() << "] does not match the number of SubFactorManagers [" << numFactManagers << "]");
123 
124  // Build blocked prolongator
125  std::vector<RCP<Matrix> > subBlockP (numFactManagers);
126  std::vector<RCP<const Map> > subBlockPRangeMaps (numFactManagers);
127  std::vector<RCP<const Map> > subBlockPDomainMaps(numFactManagers);
128 
129  std::vector<GO> fullRangeMapVector;
130  std::vector<GO> fullDomainMapVector;
131 
132  RCP<const MapExtractor> rangeAMapExtractor = A->getRangeMapExtractor();
133  RCP<const MapExtractor> domainAMapExtractor = A->getDomainMapExtractor();
134 
135  const ParameterList& pL = GetParameterList();
136  const bool backwards = pL.get<bool>("backwards");
137 
138  // Build and store the subblocks and the corresponding range and domain
139  // maps. Since we put together the full range and domain map from the
140  // submaps, we do not have to use the maps from blocked A
141  for (int k = 0; k < numFactManagers; k++) {
142  int i = (backwards ? numFactManagers-1 - k : k);
143  if (restrictionMode_) GetOStream(Runtime1) << "Generating R for block " << k <<"/"<<numFactManagers <<std::endl;
144  else GetOStream(Runtime1) << "Generating P for block " << k <<"/"<<numFactManagers <<std::endl;
145 
146  const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
147 
148  SetFactoryManager fineSFM (rcpFromRef(fineLevel), factManager);
149  SetFactoryManager coarseSFM(rcpFromRef(coarseLevel), factManager);
150 
151  if (!restrictionMode_) subBlockP[i] = coarseLevel.Get<RCP<Matrix> >("P", factManager->GetFactory("P").get());
152  else subBlockP[i] = coarseLevel.Get<RCP<Matrix> >("R", factManager->GetFactory("R").get());
153 
154  // Check if prolongator/restrictor operators have strided maps
155  TEUCHOS_TEST_FOR_EXCEPTION(subBlockP[i]->IsView("stridedMaps") == false, Exceptions::BadCast,
156  "subBlock P operator [" << i << "] has no strided map information.");
157 
158  // Append strided row map (= range map) to list of range maps.
159  // TAW the row map GIDs extracted here represent the actual row map GIDs.
160  // No support for nested operators! (at least if Thyra style gids are used)
161  RCP<const StridedMap> strPartialMap = Teuchos::rcp_dynamic_cast<const StridedMap>(subBlockP[i]->getRowMap("stridedMaps"));
162  std::vector<size_t> stridedRgData = strPartialMap->getStridingData();
163  subBlockPRangeMaps[i] = StridedMapFactory::Build(
164  subBlockP[i]->getRangeMap(), // actual global IDs (Thyra or Xpetra)
165  stridedRgData,
166  strPartialMap->getStridedBlockId(),
167  strPartialMap->getOffset());
168  //subBlockPRangeMaps[i] = subBlockP[i]->getRowMap("stridedMaps");
169 
170  // Use plain range map to determine the DOF ids
171  ArrayView<const GO> nodeRangeMap = subBlockPRangeMaps[i]->getNodeElementList();
172  fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
173 
174  // Append strided col map (= domain map) to list of range maps.
175  RCP<const StridedMap> strPartialMap2 = Teuchos::rcp_dynamic_cast<const StridedMap>(subBlockP[i]->getColMap("stridedMaps"));
176  std::vector<size_t> stridedRgData2 = strPartialMap2->getStridingData();
177  subBlockPDomainMaps[i] = StridedMapFactory::Build(
178  subBlockP[i]->getDomainMap(), // actual global IDs (Thyra or Xpetra)
179  stridedRgData2,
180  strPartialMap2->getStridedBlockId(),
181  strPartialMap2->getOffset());
182  //subBlockPDomainMaps[i] = subBlockP[i]->getColMap("stridedMaps");
183 
184  // Use plain domain map to determine the DOF ids
185  ArrayView<const GO> nodeDomainMap = subBlockPDomainMaps[i]->getNodeElementList();
186  fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
187  }
188 
189  // check if GIDs for full maps have to be sorted:
190  // For the Thyra mode ordering they do not have to be sorted since the GIDs are
191  // numbered as 0...n1,0...,n2 (starting with zero for each subblock). The MapExtractor
192  // generates unique GIDs during the construction.
193  // For Xpetra style, the GIDs have to be reordered. Such that one obtains a ordered
194  // list of GIDs in an increasing ordering. In Xpetra, the GIDs are all unique through
195  // out all submaps.
196  bool bDoSortRangeGIDs = rangeAMapExtractor->getThyraMode() ? false : true;
197  bool bDoSortDomainGIDs = domainAMapExtractor->getThyraMode() ? false : true;
198 
199  if (bDoSortRangeGIDs) {
200  sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
201  fullRangeMapVector.erase(std::unique(fullRangeMapVector.begin(), fullRangeMapVector.end()), fullRangeMapVector.end());
202  }
203  if (bDoSortDomainGIDs) {
204  sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
205  fullDomainMapVector.erase(std::unique(fullDomainMapVector.begin(), fullDomainMapVector.end()), fullDomainMapVector.end());
206  }
207 
208  // extract map index base from maps of blocked A
209  GO rangeIndexBase = 0;
210  GO domainIndexBase = 0;
211  if (!restrictionMode_) {
212  // Prolongation mode: just use index base of range and domain map of A
213  rangeIndexBase = A->getRangeMap() ->getIndexBase();
214  domainIndexBase = A->getDomainMap()->getIndexBase();
215 
216  } else {
217  // Restriction mode: switch range and domain map for blocked restriction operator
218  rangeIndexBase = A->getDomainMap()->getIndexBase();
219  domainIndexBase = A->getRangeMap()->getIndexBase();
220  }
221 
222  // Build full range map.
223  // If original range map has striding information, then transfer it to the
224  // new range map
225  RCP<const StridedMap> stridedRgFullMap = rcp_dynamic_cast<const StridedMap>(rangeAMapExtractor->getFullMap());
226  RCP<const Map > fullRangeMap = Teuchos::null;
227 
228  ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0, fullRangeMapVector.size());
229  if (stridedRgFullMap != Teuchos::null) {
230  std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
231  fullRangeMap = StridedMapFactory::Build(
232  A->getRangeMap()->lib(),
233  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
234  fullRangeMapGIDs,
235  rangeIndexBase,
236  stridedData,
237  A->getRangeMap()->getComm(),
238  -1, /* the full map vector should always have strided block id -1! */
239  stridedRgFullMap->getOffset());
240  } else {
241  fullRangeMap = MapFactory::Build(
242  A->getRangeMap()->lib(),
243  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
244  fullRangeMapGIDs,
245  rangeIndexBase,
246  A->getRangeMap()->getComm());
247  }
248 
249  ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
250  RCP<const StridedMap> stridedDoFullMap = rcp_dynamic_cast<const StridedMap>(domainAMapExtractor->getFullMap());
251  RCP<const Map > fullDomainMap = null;
252  if (stridedDoFullMap != null) {
253  TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap == Teuchos::null, Exceptions::BadCast,
254  "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information!");
255 
256  std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
257  fullDomainMap = StridedMapFactory::Build(
258  A->getDomainMap()->lib(),
259  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
260  fullDomainMapGIDs,
261  domainIndexBase,
262  stridedData2,
263  A->getDomainMap()->getComm(),
264  -1, /* the full map vector should always have strided block id -1! */
265  stridedDoFullMap->getOffset());
266  } else {
267  fullDomainMap = MapFactory::Build(
268  A->getDomainMap()->lib(),
269  Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
270  fullDomainMapGIDs,
271  domainIndexBase,
272  A->getDomainMap()->getComm());
273  }
274 
275  // Build map extractors
276  RCP<const MapExtractor> rangeMapExtractor = MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, rangeAMapExtractor->getThyraMode());
277  RCP<const MapExtractor> domainMapExtractor = MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, domainAMapExtractor->getThyraMode());
278 
279  RCP<BlockedCrsMatrix> P = rcp(new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 10));
280  for (size_t i = 0; i < subBlockPRangeMaps.size(); i++)
281  for (size_t j = 0; j < subBlockPRangeMaps.size(); j++)
282  if (i == j) {
283  RCP<CrsMatrixWrap> crsOpii = rcp_dynamic_cast<CrsMatrixWrap>(subBlockP[i]);
284  TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null, Xpetra::Exceptions::BadCast,
285  "Block [" << i << ","<< j << "] must be of type CrsMatrixWrap.");
286  P->setMatrix(i, i, crsOpii);
287  } else {
288  P->setMatrix(i, j, Teuchos::null);
289  }
290 
291  P->fillComplete();
292 
293  // Level Set
294  if (!restrictionMode_) {
295  // Prolongation mode
296  coarseLevel.Set("P", rcp_dynamic_cast<Matrix>(P), this);
297  // Stick the CoarseMap on the level if somebody wants it (useful for repartitioning)
298  coarseLevel.Set("CoarseMap",P->getBlockedDomainMap(),this);
299  } else {
300  // Restriction mode
301  // We do not have to transpose the blocked R operator since the subblocks
302  // on the diagonal are already valid R subblocks
303  coarseLevel.Set("R", Teuchos::rcp_dynamic_cast<Matrix>(P), this);
304  }
305 
306  }
307 
308 } // namespace MueLu
309 
310 #endif /* MUELU_BLOCKEDPFACTORY_DEF_HPP_ */
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Exception indicating invalid cast attempted.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method &#39;Level::SetFactoryManager()&#39;.
Description of what is happening (more verbose)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()