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