53 #ifndef MUELU_BLOCKEDPFACTORY_DEF_HPP_ 54 #define MUELU_BLOCKEDPFACTORY_DEF_HPP_ 56 #include <Xpetra_VectorFactory.hpp> 57 #include <Xpetra_ImportFactory.hpp> 58 #include <Xpetra_ExportFactory.hpp> 59 #include <Xpetra_CrsMatrixWrap.hpp> 61 #include <Xpetra_BlockedCrsMatrix.hpp> 62 #include <Xpetra_Map.hpp> 63 #include <Xpetra_MapFactory.hpp> 64 #include <Xpetra_MapExtractor.hpp> 65 #include <Xpetra_MapExtractorFactory.hpp> 68 #include "MueLu_TentativePFactory.hpp" 70 #include "MueLu_SmootherFactory.hpp" 71 #include "MueLu_FactoryManager.hpp" 72 #include "MueLu_Utilities.hpp" 74 #include "MueLu_HierarchyHelpers.hpp" 78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 RCP<ParameterList> validParamList = rcp(
new ParameterList());
82 validParamList->set< RCP<const FactoryBase> >(
"A", null,
"Generating factory of the matrix A (block matrix)");
83 validParamList->set<
bool > (
"backwards",
false,
"Forward/backward order");
85 return validParamList;
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
90 Input(fineLevel,
"A");
92 const ParameterList& pL = GetParameterList();
93 const bool backwards = pL.get<
bool>(
"backwards");
95 const int numFactManagers = FactManager_.size();
96 for (
int k = 0; k < numFactManagers; k++) {
97 int i = (backwards ? numFactManagers-1 - k : k);
98 const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
103 if (!restrictionMode_)
104 coarseLevel.
DeclareInput(
"P", factManager->GetFactory(
"P").get(),
this);
106 coarseLevel.
DeclareInput(
"R", factManager->GetFactory(
"R").get(),
this);
111 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
117 RCP<Matrix> Ain = Get< RCP<Matrix> >(fineLevel,
"A");
119 RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
120 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(),
Exceptions::BadCast,
"Input matrix A is not a BlockedCrsMatrix.");
122 const int numFactManagers = FactManager_.size();
126 "Number of block rows [" << A->Rows() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
128 "Number of block cols [" << A->Cols() <<
"] does not match the number of SubFactorManagers [" << numFactManagers <<
"]");
132 std::vector<RCP<Matrix> > subBlockP (numFactManagers);
133 std::vector<RCP<const Map> > subBlockPRangeMaps (numFactManagers);
134 std::vector<RCP<const Map> > subBlockPDomainMaps(numFactManagers);
136 std::vector<GO> fullRangeMapVector;
137 std::vector<GO> fullDomainMapVector;
139 const ParameterList& pL = GetParameterList();
140 const bool backwards = pL.get<
bool>(
"backwards");
145 for (
int k = 0; k < numFactManagers; k++) {
146 int i = (backwards ? numFactManagers-1 - k : k);
147 const RCP<const FactoryManagerBase>& factManager = FactManager_[i];
152 if (!restrictionMode_) subBlockP[i] = coarseLevel.
Get<RCP<Matrix> >(
"P", factManager->GetFactory(
"P").get());
153 else subBlockP[i] = coarseLevel.
Get<RCP<Matrix> >(
"R", factManager->GetFactory(
"R").get());
156 TEUCHOS_TEST_FOR_EXCEPTION(subBlockP[i]->IsView(
"stridedMaps") ==
false,
Exceptions::BadCast,
157 "subBlock P operator [" << i <<
"] has no strided map information.");
160 subBlockPRangeMaps[i] = subBlockP[i]->getRowMap(
"stridedMaps");
163 ArrayView<const GO> nodeRangeMap = subBlockPRangeMaps[i]->getNodeElementList();
164 fullRangeMapVector.insert(fullRangeMapVector.end(), nodeRangeMap.begin(), nodeRangeMap.end());
167 subBlockPDomainMaps[i] = subBlockP[i]->getColMap(
"stridedMaps");
170 ArrayView<const GO> nodeDomainMap = subBlockPDomainMaps[i]->getNodeElementList();
171 fullDomainMapVector.insert(fullDomainMapVector.end(), nodeDomainMap.begin(), nodeDomainMap.end());
178 RCP<const StridedMap> stridedSubBlockPRangeMap = rcp_dynamic_cast<
const StridedMap>(subBlockPRangeMaps[0]);
179 if(stridedSubBlockPRangeMap != Teuchos::null && stridedSubBlockPRangeMap->isStrided() ==
true) {
180 sort(fullRangeMapVector.begin(), fullRangeMapVector.end());
181 fullRangeMapVector.erase(std::unique(fullRangeMapVector.begin(), fullRangeMapVector.end()), fullRangeMapVector.end());
183 RCP<const StridedMap> stridedSubBlockPDomainMap = rcp_dynamic_cast<
const StridedMap>(subBlockPDomainMaps[0]);
184 if(stridedSubBlockPDomainMap != Teuchos::null && stridedSubBlockPDomainMap->isStrided() ==
true) {
185 sort(fullDomainMapVector.begin(), fullDomainMapVector.end());
186 fullDomainMapVector.erase(std::unique(fullDomainMapVector.begin(), fullDomainMapVector.end()), fullDomainMapVector.end());
190 GO rangeIndexBase = 0;
191 GO domainIndexBase = 0;
192 if (!restrictionMode_) {
194 rangeIndexBase = A->getRangeMap() ->getIndexBase();
195 domainIndexBase = A->getDomainMap()->getIndexBase();
199 rangeIndexBase = A->getDomainMap()->getIndexBase();
200 domainIndexBase = A->getRangeMap()->getIndexBase();
206 RCP<const MapExtractor> rangeAMapExtractor = A->getRangeMapExtractor();
207 RCP<const StridedMap> stridedRgFullMap = rcp_dynamic_cast<
const StridedMap>(rangeAMapExtractor->getFullMap());
208 RCP<const Map > fullRangeMap = Teuchos::null;
210 ArrayView<GO> fullRangeMapGIDs(fullRangeMapVector.size() ? &fullRangeMapVector[0] : 0, fullRangeMapVector.size());
211 if (stridedRgFullMap != Teuchos::null) {
212 std::vector<size_t> stridedData = stridedRgFullMap->getStridingData();
213 fullRangeMap = StridedMapFactory::Build(
214 A->getRangeMap()->lib(),
215 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
219 A->getRangeMap()->getComm(),
222 stridedRgFullMap->getOffset());
224 fullRangeMap = MapFactory::Build(
225 A->getRangeMap()->lib(),
226 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
229 A->getRangeMap()->getComm());
232 RCP<const MapExtractor> domainAMapExtractor = A->getDomainMapExtractor();
233 Teuchos::ArrayView<GO> fullDomainMapGIDs(fullDomainMapVector.size() ? &fullDomainMapVector[0] : 0,fullDomainMapVector.size());
234 Teuchos::RCP<const StridedMap> stridedDoFullMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(domainAMapExtractor->getFullMap());
235 Teuchos::RCP<const Map > fullDomainMap = Teuchos::null;
236 if (stridedDoFullMap != Teuchos::null) {
237 TEUCHOS_TEST_FOR_EXCEPTION(stridedDoFullMap==Teuchos::null,
Exceptions::BadCast,
"MueLu::BlockedPFactory::Build: full map in domain map extractor has no striding information! error.");
238 std::vector<size_t> stridedData2 = stridedDoFullMap->getStridingData();
239 fullDomainMap = StridedMapFactory::Build(
240 A->getDomainMap()->lib(),
241 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
245 A->getDomainMap()->getComm(),
248 stridedDoFullMap->getOffset());
250 fullDomainMap = MapFactory::Build(
251 A->getDomainMap()->lib(),
252 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
255 A->getDomainMap()->getComm());
272 bool bRangeUseThyraStyleNumbering = !this->areGidsUnique(fullRangeMapVector);
273 bool bDomainUseThyraStyleNumbering = !this->areGidsUnique(fullDomainMapVector);
274 Teuchos::RCP<const Teuchos::Comm<int> > comm = Ain->getRowMap()->getComm();
275 int maxRangeStyles = 0;
276 int maxDomainStyles = 0;
277 MueLu_maxAll(comm, bRangeUseThyraStyleNumbering ==
true ? 1 : 0, maxRangeStyles);
278 MueLu_maxAll(comm, bDomainUseThyraStyleNumbering ==
true ? 1 : 0, maxDomainStyles);
279 bRangeUseThyraStyleNumbering = (maxRangeStyles == 1) ?
true :
false;
280 bDomainUseThyraStyleNumbering = (maxDomainStyles == 1) ?
true :
false;
283 RCP<const MapExtractor> rangeMapExtractor = MapExtractorFactory::Build(fullRangeMap, subBlockPRangeMaps, bRangeUseThyraStyleNumbering);
284 RCP<const MapExtractor> domainMapExtractor = MapExtractorFactory::Build(fullDomainMap, subBlockPDomainMaps, bDomainUseThyraStyleNumbering);
286 RCP<BlockedCrsMatrix> P = rcp(
new BlockedCrsMatrix(rangeMapExtractor, domainMapExtractor, 10));
287 for (
size_t i = 0; i < subBlockPRangeMaps.size(); i++)
288 for (
size_t j = 0; j < subBlockPRangeMaps.size(); j++)
290 RCP<CrsMatrixWrap> crsOpii = rcp_dynamic_cast<CrsMatrixWrap>(subBlockP[i]);
291 RCP<CrsMatrix> crsMatii = crsOpii->getCrsMatrix();
292 P->setMatrix(i, i, crsMatii);
294 P->setMatrix(i, j, Teuchos::null);
300 if (!restrictionMode_) {
302 coarseLevel.
Set(
"P", rcp_dynamic_cast<Matrix>(P),
this);
308 coarseLevel.
Set(
"R", Teuchos::rcp_dynamic_cast<Matrix>(P),
this);
Exception indicating invalid cast attempted.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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.
#define MueLu_maxAll(rcpComm, in, out)
Namespace for MueLu classes and methods.
Class that holds all level-specific information.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
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()
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.