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  * MueLu_SubBlockAFactory_def.hpp
48  *
49  * Created on: 02.01.2012
50  * Author: tobias
51  */
52 
53 #ifndef MUELU_SUBBLOCKAFACTORY_DEF_HPP_
54 #define MUELU_SUBBLOCKAFACTORY_DEF_HPP_
55 
56 
58 
59 #include <Xpetra_BlockedCrsMatrix.hpp>
60 #include <Xpetra_CrsMatrix.hpp>
61 #include <Xpetra_CrsMatrixWrap.hpp>
62 #include <Xpetra_MapExtractor.hpp>
63 #include <Xpetra_Matrix.hpp>
64 #include <Xpetra_StridedMapFactory.hpp>
65 
66 #include "MueLu_Level.hpp"
67 #include "MueLu_Monitor.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", MueLu::NoFactory::getRCP(), "Generating factory for A.");
76  validParamList->set< int > ("block row", 0, "Block row of subblock matrix A");
77  validParamList->set< int > ("block col", 0, "Block column of subblock matrix A");
78 
79  validParamList->set< std::string > ("Range map: Striding info", "{}", "Striding information for range map");
80  validParamList->set< LocalOrdinal > ("Range map: Strided block id", -1, "Strided block id for range map" );
81  validParamList->set< std::string > ("Domain map: Striding info", "{}", "Striding information for domain map");
82  validParamList->set< LocalOrdinal > ("Domain map: Strided block id", -1, "Strided block id for domain map" );
83 
84  return validParamList;
85  }
86 
87  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
89  Input(currentLevel, "A");
90  }
91 
92  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
94  const ParameterList& pL = GetParameterList();
95  size_t row = Teuchos::as<size_t>(pL.get<int>("block row"));
96  size_t col = Teuchos::as<size_t>(pL.get<int>("block col"));
97 
98  RCP<Matrix> Ain = Get<RCP<Matrix> >(currentLevel, "A");
99  RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
100 
101  TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(), Exceptions::BadCast, "Input matrix A is not a BlockedCrsMatrix.");
102  TEUCHOS_TEST_FOR_EXCEPTION(row > A->Rows(), Exceptions::RuntimeError, "row [" << row << "] > A.Rows() [" << A->Rows() << "].");
103  TEUCHOS_TEST_FOR_EXCEPTION(col > A->Cols(), Exceptions::RuntimeError, "col [" << col << "] > A.Cols() [" << A->Cols() << "].");
104 
105  RCP<CrsMatrixWrap> Op = Teuchos::rcp(new CrsMatrixWrap(A->getMatrix(row, col)));
106 
107  // strided maps for range and domain map of sub matrix
108  RCP<const StridedMap> srangeMap = Teuchos::null;
109  RCP<const StridedMap> sdomainMap = Teuchos::null;
110 
111  // check for user-specified striding information from XML file
112 
113  std::vector<size_t> rangeStridingInfo;
114  std::vector<size_t> domainStridingInfo;
115  LocalOrdinal rangeStridedBlockId = 0;
116  LocalOrdinal domainStridedBlockId = 0;
117  bool bRangeUserSpecified = CheckForUserSpecifiedBlockInfo(true, rangeStridingInfo, rangeStridedBlockId);
118  bool bDomainUserSpecified = CheckForUserSpecifiedBlockInfo(false, domainStridingInfo, domainStridedBlockId);
119  TEUCHOS_TEST_FOR_EXCEPTION(bRangeUserSpecified != bDomainUserSpecified, Exceptions::RuntimeError, "MueLu::SubBlockAFactory[" << row << "," << col << "]: the user has to specify either both domain and range map or none.");
120 
121  // extract map information from MapExtractor
122  RCP<const MapExtractor> rangeMapExtractor = A->getRangeMapExtractor();
123  RCP<const MapExtractor> domainMapExtractor = A->getDomainMapExtractor();
124 
125  RCP<const Map> rangeMap = rangeMapExtractor ->getMap(row);
126  RCP<const Map> domainMap = domainMapExtractor->getMap(col);
127 
128  // use user-specified striding information if available. Otherwise try to use internal striding information from the submaps!
129  if(bRangeUserSpecified) srangeMap = Teuchos::rcp(new StridedMap(rangeMap,rangeStridingInfo,rangeMap->getIndexBase(),rangeStridedBlockId,0));
130  else srangeMap = rcp_dynamic_cast<const StridedMap>(rangeMap);
131 
132  if(bDomainUserSpecified) sdomainMap = Teuchos::rcp(new StridedMap(domainMap,domainStridingInfo,domainMap->getIndexBase(),domainStridedBlockId,0));
133  else sdomainMap = rcp_dynamic_cast<const StridedMap>(domainMap);
134 
135  // In case that both user-specified and internal striding information from the submaps
136  // does not contain valid striding information, try to extract it from the global maps
137  // in the map extractor.
138  if (srangeMap.is_null()) {
139  RCP<const Map> fullRangeMap = rangeMapExtractor->getFullMap();
140  RCP<const StridedMap> sFullRangeMap = rcp_dynamic_cast<const StridedMap>(fullRangeMap);
141  TEUCHOS_TEST_FOR_EXCEPTION(sFullRangeMap.is_null(), Exceptions::BadCast, "Full rangeMap is not a strided map.");
142 
143  std::vector<size_t> stridedData = sFullRangeMap->getStridingData();
144  if (stridedData.size() == 1 && row > 0) {
145  // We have block matrices. use striding block information 0
146  srangeMap = StridedMapFactory::Build(rangeMap, stridedData, 0, sFullRangeMap->getOffset());
147 
148  } else {
149  // We have strided matrices. use striding information of the corresponding block
150  srangeMap = StridedMapFactory::Build(rangeMap, stridedData, row, sFullRangeMap->getOffset());
151  }
152  }
153 
154  if (sdomainMap.is_null()) {
155  RCP<const Map> fullDomainMap = domainMapExtractor->getFullMap();
156  RCP<const StridedMap> sFullDomainMap = rcp_dynamic_cast<const StridedMap>(fullDomainMap);
157  TEUCHOS_TEST_FOR_EXCEPTION(sFullDomainMap.is_null(), Exceptions::BadCast, "Full domainMap is not a strided map");
158 
159  std::vector<size_t> stridedData = sFullDomainMap->getStridingData();
160  if (stridedData.size() == 1 && col > 0) {
161  // We have block matrices. use striding block information 0
162  sdomainMap = StridedMapFactory::Build(domainMap, stridedData, 0, sFullDomainMap->getOffset());
163 
164  } else {
165  // We have strided matrices. use striding information of the corresponding block
166  sdomainMap = StridedMapFactory::Build(domainMap, stridedData, col, sFullDomainMap->getOffset());
167  }
168  }
169 
170  TEUCHOS_TEST_FOR_EXCEPTION(srangeMap.is_null(), Exceptions::BadCast, "rangeMap " << row << " is not a strided map.");
171  TEUCHOS_TEST_FOR_EXCEPTION(sdomainMap.is_null(), Exceptions::BadCast, "domainMap " << col << " is not a strided map.");
172 
173  GetOStream(Statistics1) << "A(" << row << "," << col << ") has strided maps:"
174  << "\n range map fixed block size = " << srangeMap ->getFixedBlockSize() << ", strided block id = " << srangeMap ->getStridedBlockId()
175  << "\n domain map fixed block size = " << sdomainMap->getFixedBlockSize() << ", strided block id = " << sdomainMap->getStridedBlockId() << std::endl;
176 
177  if (Op->IsView("stridedMaps") == true)
178  Op->RemoveView("stridedMaps");
179  Op->CreateView("stridedMaps", srangeMap, sdomainMap);
180 
181  TEUCHOS_TEST_FOR_EXCEPTION(Op->IsView("stridedMaps") == false, Exceptions::RuntimeError, "Failed to set \"stridedMaps\" view.");
182 
183  currentLevel.Set("A", rcp_dynamic_cast<Matrix>(Op), this);
184  }
185 
186  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
187  bool SubBlockAFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::CheckForUserSpecifiedBlockInfo(bool bRange, std::vector<size_t>& stridingInfo, LocalOrdinal& stridedBlockId) const {
188  const ParameterList & pL = GetParameterList();
189 
190  if(bRange == true)
191  stridedBlockId = pL.get<LocalOrdinal>("Range map: Strided block id");
192  else
193  stridedBlockId = pL.get<LocalOrdinal>("Domain map: Strided block id");
194 
195  // read in stridingInfo from parameter list and fill the internal member variable
196  // read the data only if the parameter "Striding info" exists and is non-empty
197  std::string str;
198  if(bRange == true) str = std::string("Range map: Striding info");
199  else str = std::string("Domain map: Striding info");
200  if(pL.isParameter(str)) {
201  std::string strStridingInfo = pL.get<std::string>(str);
202  if(strStridingInfo.empty() == false) {
203  Teuchos::Array<size_t> arrayVal = Teuchos::fromStringToArray<size_t>(strStridingInfo);
204  stridingInfo = Teuchos::createVector(arrayVal);
205  }
206  }
207 
208  if(stridingInfo.size() > 0) return true;
209  return false;
210  }
211 
212 } // namespace MueLu
213 
214 #endif /* MUELU_SUBBLOCKAFACTORY_DEF_HPP_ */
Exception indicating invalid cast attempted.
Print more statistics.
bool CheckForUserSpecifiedBlockInfo(bool bRange, std::vector< size_t > &stridingInfo, LocalOrdinal &stridedBlockId) const
Namespace for MueLu classes and methods.
void DeclareInput(Level &currentLevel) const
Specifies the data that this class needs, and the factories that generate that data.
void Build(Level &currentLevel) const
Build an object with this factory.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
RCP< const ParameterList > GetValidParameterList() const
Input.
Exception throws to report errors in the internal logical of the program.
static const RCP< const NoFactory > getRCP()
Static Get() functions.