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>
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");
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);
98 if (!restrictionMode_)
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");
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;
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;
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());
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());
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());
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(),
239 A->getRangeMap()->getComm(),
241 stridedRgFullMap->getOffset());
243 fullRangeMap = MapFactory::Build(
244 A->getRangeMap()->lib(),
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());
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(),
265 A->getDomainMap()->getComm(),
267 stridedDoFullMap->getOffset());
269 fullDomainMap = MapFactory::Build(
270 A->getDomainMap()->lib(),
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());
282 for (
size_t i = 0; i < subBlockPRangeMaps.size(); i++)
283 for (
size_t j = 0; j < subBlockPRangeMaps.size(); j++)
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.
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()'.
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()