47 #ifndef MUELU_BLOCKEDPFACTORY_DEF_HPP_
48 #define MUELU_BLOCKEDPFACTORY_DEF_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);
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.
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)
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 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()'.
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()
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.