47 #ifndef MUELU_BLOCKEDPFACTORY_DEF_HPP_ 48 #define MUELU_BLOCKEDPFACTORY_DEF_HPP_ 50 #include <Xpetra_VectorFactory.hpp> 51 #include <Xpetra_ImportFactory.hpp> 52 #include <Xpetra_ExportFactory.hpp> 53 #include <Xpetra_CrsMatrixWrap.hpp> 55 #include <Xpetra_BlockedCrsMatrix.hpp> 56 #include <Xpetra_Map.hpp> 57 #include <Xpetra_MapFactory.hpp> 58 #include <Xpetra_MapExtractor.hpp> 59 #include <Xpetra_MapExtractorFactory.hpp> 62 #include "MueLu_TentativePFactory.hpp" 64 #include "MueLu_SmootherFactory.hpp" 65 #include "MueLu_FactoryManager.hpp" 66 #include "MueLu_Utilities.hpp" 68 #include "MueLu_HierarchyUtils.hpp" 72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
74 RCP<ParameterList> validParamList = rcp(
new ParameterList());
76 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A (block matrix)");
77 validParamList->set<
bool > (
"backwards",
false,
"Forward/backward order");
79 return validParamList;
82 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
84 Input(fineLevel,
"A");
86 const ParameterList& pL = GetParameterList();
87 const bool backwards = pL.get<
bool>(
"backwards");
89 const int numFactManagers = FactManager_.size();
90 for (
int k = 0; k < numFactManagers; k++) {
91 int i = (backwards ? numFactManagers-1 - k : k);
92 const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
97 if (!restrictionMode_)
98 coarseLevel.
DeclareInput(
"P", factManager->GetFactory(
"P").get(),
this);
100 coarseLevel.
DeclareInput(
"R", factManager->GetFactory(
"R").get(),
this);
104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
108 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
112 RCP<Matrix> Ain = Get< RCP<Matrix> >(fineLevel,
"A");
114 RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
115 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(),
Exceptions::BadCast,
"Input matrix A is not a BlockedCrsMatrix.");
117 const int numFactManagers = FactManager_.size();
121 "Number of block rows [" << A->Rows() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
123 "Number of block cols [" << A->Cols() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
126 std::vector<RCP<Matrix> > subBlockP (numFactManagers);
127 std::vector<RCP<const Map> > subBlockPRangeMaps (numFactManagers);
128 std::vector<RCP<const Map> > subBlockPDomainMaps(numFactManagers);
130 std::vector<GO> fullRangeMapVector;
131 std::vector<GO> fullDomainMapVector;
133 RCP<const MapExtractor> rangeAMapExtractor = A->getRangeMapExtractor();
134 RCP<const MapExtractor> domainAMapExtractor = A->getDomainMapExtractor();
136 const ParameterList& pL = GetParameterList();
137 const bool backwards = pL.get<
bool>(
"backwards");
142 for (
int k = 0; k < numFactManagers; k++) {
143 int i = (backwards ? numFactManagers-1 - k : k);
144 const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
149 if (!restrictionMode_) subBlockP[i] = coarseLevel.
Get<RCP<Matrix> >(
"P", factManager->GetFactory(
"P").get());
150 else subBlockP[i] = coarseLevel.
Get<RCP<Matrix> >(
"R", factManager->GetFactory(
"R").get());
153 TEUCHOS_TEST_FOR_EXCEPTION(subBlockP[i]->IsView(
"stridedMaps") ==
false,
Exceptions::BadCast,
154 "subBlock P operator [" << i <<
"] has no strided map information.");
159 RCP<const StridedMap> strPartialMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(subBlockP[i]->getRowMap(
"stridedMaps"));
160 std::vector<size_t> stridedRgData = strPartialMap->getStridingData();
161 subBlockPRangeMaps[i] = StridedMapFactory::Build(
162 subBlockP[i]->getRangeMap(),
164 strPartialMap->getStridedBlockId(),
165 strPartialMap->getOffset());
169 ArrayView<const GO> nodeRangeMap = subBlockPRangeMaps[i]->getNodeElementList();
170 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
173 RCP<const StridedMap> strPartialMap2 = Teuchos::rcp_dynamic_cast<
const StridedMap>(subBlockP[i]->getColMap(
"stridedMaps"));
174 std::vector<size_t> stridedRgData2 = strPartialMap2->getStridingData();
175 subBlockPDomainMaps[i] = StridedMapFactory::Build(
176 subBlockP[i]->getDomainMap(),
178 strPartialMap2->getStridedBlockId(),
179 strPartialMap2->getOffset());
183 ArrayView<const GO> nodeDomainMap = subBlockPDomainMaps[i]->getNodeElementList();
184 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
194 bool bDoSortRangeGIDs = rangeAMapExtractor->getThyraMode() ? false :
true;
195 bool bDoSortDomainGIDs = domainAMapExtractor->getThyraMode() ? false :
true;
197 if (bDoSortRangeGIDs) {
198 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
199 fullRangeMapVector.erase(std::unique(fullRangeMapVector.begin(), fullRangeMapVector.end()), fullRangeMapVector.end());
201 if (bDoSortDomainGIDs) {
202 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
203 fullDomainMapVector.erase(std::unique(fullDomainMapVector.begin(), fullDomainMapVector.end()), fullDomainMapVector.end());
207 GO rangeIndexBase = 0;
208 GO domainIndexBase = 0;
209 if (!restrictionMode_) {
211 rangeIndexBase = A->getRangeMap() ->getIndexBase();
212 domainIndexBase = A->getDomainMap()->getIndexBase();
216 rangeIndexBase = A->getDomainMap()->getIndexBase();
217 domainIndexBase = A->getRangeMap()->getIndexBase();
223 RCP<const StridedMap> stridedRgFullMap = rcp_dynamic_cast<
const StridedMap>(rangeAMapExtractor->getFullMap());
224 RCP<const Map > fullRangeMap = Teuchos::null;
226 ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0, fullRangeMapVector.size());
227 if (stridedRgFullMap != Teuchos::null) {
228 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
229 fullRangeMap = StridedMapFactory::Build(
230 A->getRangeMap()->lib(),
231 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
235 A->getRangeMap()->getComm(),
237 stridedRgFullMap->getOffset());
239 fullRangeMap = MapFactory::Build(
240 A->getRangeMap()->lib(),
241 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
244 A->getRangeMap()->getComm());
247 ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
248 RCP<const StridedMap> stridedDoFullMap = rcp_dynamic_cast<
const StridedMap>(domainAMapExtractor->getFullMap());
249 RCP<const Map > fullDomainMap = null;
250 if (stridedDoFullMap != null) {
252 "MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
254 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
255 fullDomainMap = StridedMapFactory::Build(
256 A->getDomainMap()->lib(),
257 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
261 A->getDomainMap()->getComm(),
263 stridedDoFullMap->getOffset());
265 fullDomainMap = MapFactory::Build(
266 A->getDomainMap()->lib(),
267 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
270 A->getDomainMap()->getComm());
274 RCP<const MapExtractor> rangeMapExtractor = MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, rangeAMapExtractor->getThyraMode());
275 RCP<const MapExtractor> domainMapExtractor = MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, domainAMapExtractor->getThyraMode());
277 RCP<BlockedCrsMatrix> P = rcp(
new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 10));
278 for (
size_t i = 0; i < subBlockPRangeMaps.size(); i++)
279 for (
size_t j = 0; j < subBlockPRangeMaps.size(); j++)
281 RCP<CrsMatrixWrap> crsOpii = rcp_dynamic_cast<CrsMatrixWrap>(subBlockP[i]);
282 TEUCHOS_TEST_FOR_EXCEPTION(crsOpii == Teuchos::null, Xpetra::Exceptions::BadCast,
283 "Block [" << i <<
","<< j <<
"] must be of type CrsMatrixWrap.");
284 P->setMatrix(i, i, crsOpii);
286 P->setMatrix(i, j, Teuchos::null);
292 if (!restrictionMode_) {
294 coarseLevel.
Set(
"P", rcp_dynamic_cast<Matrix>(P),
this);
300 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()'.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()