MueLu  Version of the Day
MueLu_SubBlockAFactory_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_SUBBLOCKAFACTORY_DEF_HPP_
48 #define MUELU_SUBBLOCKAFACTORY_DEF_HPP_
49 
50 
52 
53 #include <Xpetra_BlockedCrsMatrix.hpp>
54 #include <Xpetra_MapExtractor.hpp>
55 #include <Xpetra_Matrix.hpp>
56 #include <Xpetra_StridedMapFactory.hpp>
57 
58 #include "MueLu_Level.hpp"
59 #include "MueLu_Monitor.hpp"
60 
61 namespace MueLu {
62 
63  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
65  RCP<ParameterList> validParamList = rcp(new ParameterList());
66 
67  validParamList->set<RCP<const FactoryBase>>("A", MueLu::NoFactory::getRCP(), "Generating factory for A.");
68 
69  validParamList->set<int>("block row", 0, "Block row of subblock matrix A");
70  validParamList->set<int>("block col", 0, "Block column of subblock matrix A");
71 
72  validParamList->set<std::string >("Range map: Striding info", "{}", "Striding information for range map");
73  validParamList->set<LocalOrdinal>("Range map: Strided block id", -1, "Strided block id for range map");
74  validParamList->set<std::string >("Domain map: Striding info", "{}", "Striding information for domain map");
75  validParamList->set<LocalOrdinal>("Domain map: Strided block id", -1, "Strided block id for domain map");
76 
77  return validParamList;
78  }
79 
80  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
82  Input(currentLevel, "A");
83  }
84 
85  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
87  FactoryMonitor m(*this, "Build", currentLevel);
88 
89  const ParameterList& pL = GetParameterList();
90  size_t row = Teuchos::as<size_t>(pL.get<int>("block row"));
91  size_t col = Teuchos::as<size_t>(pL.get<int>("block col"));
92 
93  RCP<Matrix> Ain = currentLevel.Get< RCP<Matrix> >("A",this->GetFactory("A").get());
94  RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
95 
96  TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::BadCast, "Input matrix A is not a BlockedCrsMatrix.");
97  TEUCHOS_TEST_FOR_EXCEPTION(row > A->Rows(), Exceptions::RuntimeError, "row [" << row << "] > A.Rows() [" << A->Rows() << "].");
98  TEUCHOS_TEST_FOR_EXCEPTION(col > A->Cols(), Exceptions::RuntimeError, "col [" << col << "] > A.Cols() [" << A->Cols() << "].");
99 
100  // get sub-matrix
101  RCP<Matrix> Op = A->getMatrix(row, col);
102 
103  // Check if it is a BlockedCrsMatrix object
104  // If it is a BlockedCrsMatrix object (most likely a ReorderedBlockedCrsMatrix)
105  // we have to distinguish whether it is a 1x1 leaf block in the ReorderedBlockedCrsMatrix
106  // or a nxm block. If it is a 1x1 leaf block, we "unpack" it and return the underlying
107  // CrsMatrixWrap object.
108  RCP<BlockedCrsMatrix> bOp = rcp_dynamic_cast<BlockedCrsMatrix>(Op);
109  if (bOp != Teuchos::null) {
110  // check if it is a 1x1 leaf block
111  if (bOp->Rows() == 1 && bOp->Cols() == 1) {
112  // return the unwrapped CrsMatrixWrap object underneath
113  Op = bOp->getCrsMatrix();
114  TEUCHOS_TEST_FOR_EXCEPTION(rcp_dynamic_cast<CrsMatrixWrap>(Op) == Teuchos::null, Exceptions::BadCast,
115  "SubBlockAFactory::Build: sub block A[" << row << "," << col << "] must be a single block CrsMatrixWrap object!");
116  } else {
117  // If it is a regular nxm blocked operator, just return it.
118  // We do not set any kind of striding or blocking information as this
119  // usually would not make sense for general blocked operators
120  GetOStream(Statistics1) << "A(" << row << "," << col << ") is a " << bOp->Rows() << "x" << bOp->Cols() << " block matrix" << std::endl;
121  GetOStream(Statistics2) << "with altogether " << bOp->getGlobalNumRows() << "x" << bOp->getGlobalNumCols() << " rows and columns." << std::endl;
122  currentLevel.Set("A", Op, this);
123  return;
124  }
125  }
126 
127  // The sub-block is not a BlockedCrsMatrix object, that is, we expect
128  // it to be of type CrsMatrixWrap allowing direct access to the corresponding
129  // data. For a single block CrsMatrixWrap type matrix we can/should set the
130  // corresponding striding/blocking information for the algorithms to work
131  // properly.
132  //
133  // TAW: In fact, a 1x1 BlockedCrsMatrix object also allows to access the data
134  // directly, but this feature is nowhere really used in the algorithms.
135  // So let's keep checking for the CrsMatrixWrap class to avoid screwing
136  // things up
137  //
138  TEUCHOS_TEST_FOR_EXCEPTION(rcp_dynamic_cast<CrsMatrixWrap>(Op) == Teuchos::null, Exceptions::BadCast,
139  "SubBlockAFactory::Build: sub block A[" << row << "," << col << "] is NOT a BlockedCrsMatrix, but also NOT a CrsMatrixWrap object. This cannot be.");
140 
141  // strided maps for range and domain map of sub matrix
142  RCP<const StridedMap> stridedRangeMap = Teuchos::null;
143  RCP<const StridedMap> stridedDomainMap = Teuchos::null;
144 
145  // check for user-specified striding information from XML file
146  std::vector<size_t> rangeStridingInfo;
147  std::vector<size_t> domainStridingInfo;
148  LocalOrdinal rangeStridedBlockId = 0;
149  LocalOrdinal domainStridedBlockId = 0;
150  bool bRangeUserSpecified = CheckForUserSpecifiedBlockInfo(true, rangeStridingInfo, rangeStridedBlockId);
151  bool bDomainUserSpecified = CheckForUserSpecifiedBlockInfo(false, domainStridingInfo, domainStridedBlockId);
152  TEUCHOS_TEST_FOR_EXCEPTION(bRangeUserSpecified != bDomainUserSpecified, Exceptions::RuntimeError,
153  "MueLu::SubBlockAFactory[" << row << "," << col << "]: the user has to specify either both domain and range map or none.");
154 
155  // extract map information from MapExtractor
156  RCP<const MapExtractor> rangeMapExtractor = A->getRangeMapExtractor();
157  RCP<const MapExtractor> domainMapExtractor = A->getDomainMapExtractor();
158 
159  RCP<const Map> rangeMap = rangeMapExtractor ->getMap(row);
160  RCP<const Map> domainMap = domainMapExtractor->getMap(col);
161 
162  // use user-specified striding information if available. Otherwise try to use internal striding information from the submaps!
163  if(bRangeUserSpecified) stridedRangeMap = rcp(new StridedMap(rangeMap, rangeStridingInfo, rangeMap->getIndexBase(), rangeStridedBlockId, 0));
164  else stridedRangeMap = rcp_dynamic_cast<const StridedMap>(rangeMap);
165 
166  if(bDomainUserSpecified) stridedDomainMap = rcp(new StridedMap(domainMap, domainStridingInfo, domainMap->getIndexBase(), domainStridedBlockId, 0));
167  else stridedDomainMap = rcp_dynamic_cast<const StridedMap>(domainMap);
168 
169  // In case that both user-specified and internal striding information from the submaps
170  // does not contain valid striding information, try to extract it from the global maps
171  // in the map extractor.
172  if (stridedRangeMap.is_null()) {
173  RCP<const Map> fullRangeMap = rangeMapExtractor->getFullMap();
174  RCP<const StridedMap> stridedFullRangeMap = rcp_dynamic_cast<const StridedMap>(fullRangeMap);
175  TEUCHOS_TEST_FOR_EXCEPTION(stridedFullRangeMap.is_null(), Exceptions::BadCast, "Full rangeMap is not a strided map.");
176 
177  std::vector<size_t> stridedData = stridedFullRangeMap->getStridingData();
178  if (stridedData.size() == 1 && row > 0) {
179  // We have block matrices. use striding block information 0
180  stridedRangeMap = StridedMapFactory::Build(rangeMap, stridedData, 0, stridedFullRangeMap->getOffset());
181 
182  } else {
183  // We have strided matrices. use striding information of the corresponding block
184  stridedRangeMap = StridedMapFactory::Build(rangeMap, stridedData, row, stridedFullRangeMap->getOffset());
185  }
186  }
187 
188  if (stridedDomainMap.is_null()) {
189  RCP<const Map> fullDomainMap = domainMapExtractor->getFullMap();
190  RCP<const StridedMap> stridedFullDomainMap = rcp_dynamic_cast<const StridedMap>(fullDomainMap);
191  TEUCHOS_TEST_FOR_EXCEPTION(stridedFullDomainMap.is_null(), Exceptions::BadCast, "Full domainMap is not a strided map");
192 
193  std::vector<size_t> stridedData = stridedFullDomainMap->getStridingData();
194  if (stridedData.size() == 1 && col > 0) {
195  // We have block matrices. use striding block information 0
196  stridedDomainMap = StridedMapFactory::Build(domainMap, stridedData, 0, stridedFullDomainMap->getOffset());
197 
198  } else {
199  // We have strided matrices. use striding information of the corresponding block
200  stridedDomainMap = StridedMapFactory::Build(domainMap, stridedData, col, stridedFullDomainMap->getOffset());
201  }
202  }
203 
204  TEUCHOS_TEST_FOR_EXCEPTION(stridedRangeMap.is_null(), Exceptions::BadCast, "rangeMap " << row << " is not a strided map.");
205  TEUCHOS_TEST_FOR_EXCEPTION(stridedDomainMap.is_null(), Exceptions::BadCast, "domainMap " << col << " is not a strided map.");
206 
207  GetOStream(Statistics1) << "A(" << row << "," << col << ") is a single block and has strided maps:"
208  << "\n range map fixed block size = " << stridedRangeMap ->getFixedBlockSize() << ", strided block id = " << stridedRangeMap ->getStridedBlockId()
209  << "\n domain map fixed block size = " << stridedDomainMap->getFixedBlockSize() << ", strided block id = " << stridedDomainMap->getStridedBlockId() << std::endl;
210  GetOStream(Statistics2) << "A(" << row << "," << col << ") has " << Op->getGlobalNumRows() << "x" << Op->getGlobalNumCols() << " rows and columns." << std::endl;
211 
212  // TODO do we really need that? we moved the code to getMatrix...
213  if (Op->IsView("stridedMaps") == true)
214  Op->RemoveView("stridedMaps");
215  Op->CreateView("stridedMaps", stridedRangeMap, stridedDomainMap);
216 
217  TEUCHOS_TEST_FOR_EXCEPTION(Op->IsView("stridedMaps") == false, Exceptions::RuntimeError, "Failed to set \"stridedMaps\" view.");
218 
219  currentLevel.Set("A", Op, this);
220  }
221 
222  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
224  CheckForUserSpecifiedBlockInfo(bool bRange, std::vector<size_t>& stridingInfo, LocalOrdinal& stridedBlockId) const {
225  const ParameterList & pL = GetParameterList();
226 
227  if(bRange == true)
228  stridedBlockId = pL.get<LocalOrdinal>("Range map: Strided block id");
229  else
230  stridedBlockId = pL.get<LocalOrdinal>("Domain map: Strided block id");
231 
232  // read in stridingInfo from parameter list and fill the internal member variable
233  // read the data only if the parameter "Striding info" exists and is non-empty
234  std::string str;
235  if(bRange == true) str = std::string("Range map: Striding info");
236  else str = std::string("Domain map: Striding info");
237 
238  if(pL.isParameter(str)) {
239  std::string strStridingInfo = pL.get<std::string>(str);
240  if(strStridingInfo.empty() == false) {
241  Array<size_t> arrayVal = Teuchos::fromStringToArray<size_t>(strStridingInfo);
242  stridingInfo = Teuchos::createVector(arrayVal);
243  }
244  }
245 
246  if(stridingInfo.size() > 0) return true;
247  return false;
248  }
249 
250 } // namespace MueLu
251 
252 #endif /* MUELU_SUBBLOCKAFACTORY_DEF_HPP_ */
Exception indicating invalid cast attempted.
MueLu::DefaultLocalOrdinal LocalOrdinal
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.
bool CheckForUserSpecifiedBlockInfo(bool bRange, std::vector< size_t > &stridingInfo, LocalOrdinal &stridedBlockId) const
RCP< const ParameterList > GetValidParameterList() const override
Input.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print more statistics.
void Build(Level &currentLevel) const override
Build an object with this factory.
Namespace for MueLu classes and methods.
Print even more statistics.
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())
Exception throws to report errors in the internal logical of the program.
void DeclareInput(Level &currentLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
static const RCP< const NoFactory > getRCP()
Static Get() functions.