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