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