47#ifndef MUELU_BLOCKEDPFACTORY_DEF_HPP_
48#define MUELU_BLOCKEDPFACTORY_DEF_HPP_
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>
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>
63#include "MueLu_TentativePFactory.hpp"
65#include "MueLu_FactoryManager.hpp"
67#include "MueLu_HierarchyUtils.hpp"
71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
73 RCP<ParameterList> validParamList = rcp(
new ParameterList());
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");
78 return validParamList;
81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
83 Input(fineLevel,
"A");
85 const ParameterList& pL = GetParameterList();
86 const bool backwards = pL.get<
bool>(
"backwards");
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];
96 if (!restrictionMode_)
97 coarseLevel.
DeclareInput(
"P", factManager->GetFactory(
"P").get(),
this);
99 coarseLevel.
DeclareInput(
"R", factManager->GetFactory(
"R").get(),
this);
103 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
107 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
111 RCP<Matrix> Ain = Get< RCP<Matrix> >(fineLevel,
"A");
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.");
116 const int numFactManagers = FactManager_.size();
120 "Number of block rows [" << A->Rows() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
122 "Number of block cols [" << A->Cols() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
125 std::vector<RCP<Matrix> > subBlockP (numFactManagers);
126 std::vector<RCP<const Map> > subBlockPRangeMaps (numFactManagers);
127 std::vector<RCP<const Map> > subBlockPDomainMaps(numFactManagers);
129 std::vector<GO> fullRangeMapVector;
130 std::vector<GO> fullDomainMapVector;
132 RCP<const MapExtractor> rangeAMapExtractor = A->getRangeMapExtractor();
133 RCP<const MapExtractor> domainAMapExtractor = A->getDomainMapExtractor();
135 const ParameterList& pL = GetParameterList();
136 const bool backwards = pL.get<
bool>(
"backwards");
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;
146 const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
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());
155 TEUCHOS_TEST_FOR_EXCEPTION(subBlockP[i]->IsView(
"stridedMaps") ==
false,
Exceptions::BadCast,
156 "subBlock P operator [" << i <<
"] has no strided map information.");
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(),
166 strPartialMap->getStridedBlockId(),
167 strPartialMap->getOffset());
171 ArrayView<const GO> nodeRangeMap = subBlockPRangeMaps[i]->getLocalElementList();
172 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
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(),
180 strPartialMap2->getStridedBlockId(),
181 strPartialMap2->getOffset());
185 ArrayView<const GO> nodeDomainMap = subBlockPDomainMaps[i]->getLocalElementList();
186 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
196 bool bDoSortRangeGIDs = rangeAMapExtractor->getThyraMode() ? false :
true;
197 bool bDoSortDomainGIDs = domainAMapExtractor->getThyraMode() ? false :
true;
199 if (bDoSortRangeGIDs) {
200 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
201 fullRangeMapVector.erase(std::unique(fullRangeMapVector.begin(), fullRangeMapVector.end()), fullRangeMapVector.end());
203 if (bDoSortDomainGIDs) {
204 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
205 fullDomainMapVector.erase(std::unique(fullDomainMapVector.begin(), fullDomainMapVector.end()), fullDomainMapVector.end());
209 GO rangeIndexBase = 0;
210 GO domainIndexBase = 0;
211 if (!restrictionMode_) {
213 rangeIndexBase = A->getRangeMap() ->getIndexBase();
214 domainIndexBase = A->getDomainMap()->getIndexBase();
218 rangeIndexBase = A->getDomainMap()->getIndexBase();
219 domainIndexBase = A->getRangeMap()->getIndexBase();
225 RCP<const StridedMap> stridedRgFullMap = rcp_dynamic_cast<const StridedMap>(rangeAMapExtractor->getFullMap());
226 RCP<const Map > fullRangeMap = Teuchos::null;
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(),
237 A->getRangeMap()->getComm(),
239 stridedRgFullMap->getOffset());
241 fullRangeMap = MapFactory::Build(
242 A->getRangeMap()->lib(),
243 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
246 A->getRangeMap()->getComm());
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) {
254 "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information!");
256 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
257 fullDomainMap = StridedMapFactory::Build(
258 A->getDomainMap()->lib(),
259 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
263 A->getDomainMap()->getComm(),
265 stridedDoFullMap->getOffset());
267 fullDomainMap = MapFactory::Build(
268 A->getDomainMap()->lib(),
269 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
272 A->getDomainMap()->getComm());
276 RCP<const MapExtractor> rangeMapExtractor = MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, rangeAMapExtractor->getThyraMode());
277 RCP<const MapExtractor> domainMapExtractor = MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, domainAMapExtractor->getThyraMode());
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++)
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);
288 P->setMatrix(i, j, Teuchos::null);
294 if (!restrictionMode_) {
296 coarseLevel.
Set(
"P", rcp_dynamic_cast<Matrix>(P),
this);
298 coarseLevel.
Set(
"CoarseMap",P->getBlockedDomainMap(),
this);
303 coarseLevel.
Set(
"R", Teuchos::rcp_dynamic_cast<Matrix>(P),
this);
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.
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)