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>
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");
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);
97 if (!restrictionMode_)
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");
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;
137 const bool backwards = pL.
get<
bool>(
"backwards");
142 for (
int k = 0; k < numFactManagers; k++) {
143 int i = (backwards ? numFactManagers-1 - k : k);
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());
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());
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());
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(),
235 A->getRangeMap()->getComm(),
237 stridedRgFullMap->getOffset());
239 fullRangeMap = MapFactory::Build(
240 A->getRangeMap()->lib(),
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());
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(),
261 A->getDomainMap()->getComm(),
263 stridedDoFullMap->getOffset());
265 fullDomainMap = MapFactory::Build(
266 A->getDomainMap()->lib(),
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());
278 for (
size_t i = 0; i < subBlockPRangeMaps.size(); i++)
279 for (
size_t j = 0; j < subBlockPRangeMaps.size(); j++)
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.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Namespace for MueLu classes and methods.
virtual const RCP< const FactoryBase > GetFactory(const std::string &varName) const =0
Get.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
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()