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