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_SmootherFactory.hpp" 66 #include "MueLu_FactoryManager.hpp" 67 #include "MueLu_Utilities.hpp" 69 #include "MueLu_HierarchyUtils.hpp" 73 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 RCP<ParameterList> validParamList = rcp(
new ParameterList());
77 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A (block matrix)");
78 validParamList->set<
bool > (
"backwards",
false,
"Forward/backward order");
80 return validParamList;
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
85 Input(fineLevel,
"A");
87 const ParameterList& pL = GetParameterList();
88 const bool backwards = pL.get<
bool>(
"backwards");
90 const int numFactManagers = FactManager_.size();
91 for (
int k = 0; k < numFactManagers; k++) {
92 int i = (backwards ? numFactManagers-1 - k : k);
93 const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
98 if (!restrictionMode_)
99 coarseLevel.
DeclareInput(
"P", factManager->GetFactory(
"P").get(),
this);
101 coarseLevel.
DeclareInput(
"R", factManager->GetFactory(
"R").get(),
this);
105 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
109 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 RCP<Matrix> Ain = Get< RCP<Matrix> >(fineLevel,
"A");
115 RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
116 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(),
Exceptions::BadCast,
"Input matrix A is not a BlockedCrsMatrix.");
118 const int numFactManagers = FactManager_.size();
122 "Number of block rows [" << A->Rows() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
124 "Number of block cols [" << A->Cols() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
127 std::vector<RCP<Matrix> > subBlockP (numFactManagers);
128 std::vector<RCP<const Map> > subBlockPRangeMaps (numFactManagers);
129 std::vector<RCP<const Map> > subBlockPDomainMaps(numFactManagers);
131 std::vector<GO> fullRangeMapVector;
132 std::vector<GO> fullDomainMapVector;
134 RCP<const MapExtractor> rangeAMapExtractor = A->getRangeMapExtractor();
135 RCP<const MapExtractor> domainAMapExtractor = A->getDomainMapExtractor();
137 const ParameterList& pL = GetParameterList();
138 const bool backwards = pL.get<
bool>(
"backwards");
143 for (
int k = 0; k < numFactManagers; k++) {
144 int i = (backwards ? numFactManagers-1 - k : k);
145 if (restrictionMode_) GetOStream(
Runtime1) <<
"Generating R for block " << k <<
"/"<<numFactManagers <<std::endl;
146 else GetOStream(
Runtime1) <<
"Generating P for block " << k <<
"/"<<numFactManagers <<std::endl;
148 const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
153 if (!restrictionMode_) subBlockP[i] = coarseLevel.
Get<RCP<Matrix> >(
"P", factManager->GetFactory(
"P").get());
154 else subBlockP[i] = coarseLevel.
Get<RCP<Matrix> >(
"R", factManager->GetFactory(
"R").get());
157 TEUCHOS_TEST_FOR_EXCEPTION(subBlockP[i]->IsView(
"stridedMaps") ==
false,
Exceptions::BadCast,
158 "subBlock P operator [" << i <<
"] has no strided map information.");
163 RCP<const StridedMap> strPartialMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(subBlockP[i]->getRowMap(
"stridedMaps"));
164 std::vector<size_t> stridedRgData = strPartialMap->getStridingData();
165 subBlockPRangeMaps[i] = StridedMapFactory::Build(
166 subBlockP[i]->getRangeMap(),
168 strPartialMap->getStridedBlockId(),
169 strPartialMap->getOffset());
173 ArrayView<const GO> nodeRangeMap = subBlockPRangeMaps[i]->getNodeElementList();
174 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
177 RCP<const StridedMap> strPartialMap2 = Teuchos::rcp_dynamic_cast<
const StridedMap>(subBlockP[i]->getColMap(
"stridedMaps"));
178 std::vector<size_t> stridedRgData2 = strPartialMap2->getStridingData();
179 subBlockPDomainMaps[i] = StridedMapFactory::Build(
180 subBlockP[i]->getDomainMap(),
182 strPartialMap2->getStridedBlockId(),
183 strPartialMap2->getOffset());
187 ArrayView<const GO> nodeDomainMap = subBlockPDomainMaps[i]->getNodeElementList();
188 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
198 bool bDoSortRangeGIDs = rangeAMapExtractor->getThyraMode() ? false :
true;
199 bool bDoSortDomainGIDs = domainAMapExtractor->getThyraMode() ? false :
true;
201 if (bDoSortRangeGIDs) {
202 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
203 fullRangeMapVector.erase(std::unique(fullRangeMapVector.begin(), fullRangeMapVector.end()), fullRangeMapVector.end());
205 if (bDoSortDomainGIDs) {
206 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
207 fullDomainMapVector.erase(std::unique(fullDomainMapVector.begin(), fullDomainMapVector.end()), fullDomainMapVector.end());
211 GO rangeIndexBase = 0;
212 GO domainIndexBase = 0;
213 if (!restrictionMode_) {
215 rangeIndexBase = A->getRangeMap() ->getIndexBase();
216 domainIndexBase = A->getDomainMap()->getIndexBase();
220 rangeIndexBase = A->getDomainMap()->getIndexBase();
221 domainIndexBase = A->getRangeMap()->getIndexBase();
227 RCP<const StridedMap> stridedRgFullMap = rcp_dynamic_cast<
const StridedMap>(rangeAMapExtractor->getFullMap());
228 RCP<const Map > fullRangeMap = Teuchos::null;
230 ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0, fullRangeMapVector.size());
231 if (stridedRgFullMap != Teuchos::null) {
232 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
233 fullRangeMap = StridedMapFactory::Build(
234 A->getRangeMap()->lib(),
235 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
239 A->getRangeMap()->getComm(),
241 stridedRgFullMap->getOffset());
243 fullRangeMap = MapFactory::Build(
244 A->getRangeMap()->lib(),
245 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
248 A->getRangeMap()->getComm());
251 ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
252 RCP<const StridedMap> stridedDoFullMap = rcp_dynamic_cast<
const StridedMap>(domainAMapExtractor->getFullMap());
253 RCP<const Map > fullDomainMap = null;
254 if (stridedDoFullMap != null) {
256 "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
258 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
259 fullDomainMap = StridedMapFactory::Build(
260 A->getDomainMap()->lib(),
261 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
265 A->getDomainMap()->getComm(),
267 stridedDoFullMap->getOffset());
269 fullDomainMap = MapFactory::Build(
270 A->getDomainMap()->lib(),
271 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
274 A->getDomainMap()->getComm());
278 RCP<const MapExtractor> rangeMapExtractor = MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, rangeAMapExtractor->getThyraMode());
279 RCP<const MapExtractor> domainMapExtractor = MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, domainAMapExtractor->getThyraMode());
281 RCP<BlockedCrsMatrix> P = rcp(
new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 10));
282 for (
size_t i = 0; i < subBlockPRangeMaps.size(); i++)
283 for (
size_t j = 0; j < subBlockPRangeMaps.size(); j++)
285 RCP<CrsMatrixWrap> crsOpii = rcp_dynamic_cast<CrsMatrixWrap>(subBlockP[i]);
286 TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null, Xpetra::Exceptions::BadCast,
287 "Block [" << i <<
","<< j <<
"] must be of type CrsMatrixWrap.");
288 P->setMatrix(i, i, crsOpii);
290 P->setMatrix(i, j, Teuchos::null);
296 if (!restrictionMode_) {
298 coarseLevel.
Set(
"P", rcp_dynamic_cast<Matrix>(P),
this);
300 coarseLevel.
Set(
"CoarseMap",P->getBlockedDomainMap(),
this);
305 coarseLevel.
Set(
"R", Teuchos::rcp_dynamic_cast<Matrix>(P),
this);
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
Exception indicating invalid cast attempted.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
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.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
Class that holds all level-specific information.
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method 'Level::SetFactoryManager()'.
Description of what is happening (more verbose)
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()