MueLu Version of the Day
Loading...
Searching...
No Matches
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
61namespace 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_ */
MueLu::DefaultLocalOrdinal LocalOrdinal
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const RCP< const NoFactory > getRCP()
Static Get() functions.
void DeclareInput(Level &currentLevel) const override
Specifies the data that this class needs, and the factories that generate that data.
bool CheckForUserSpecifiedBlockInfo(bool bRange, std::vector< size_t > &stridingInfo, LocalOrdinal &stridedBlockId) const
void Build(Level &currentLevel) const override
Build an object with this factory.
RCP< const ParameterList > GetValidParameterList() const override
Input.
Namespace for MueLu classes and methods.
@ Statistics2
Print even more statistics.
@ Statistics1
Print more statistics.