48#ifndef MUELU_REORDERBLOCKAFACTORY_DEF_HPP_
49#define MUELU_REORDERBLOCKAFACTORY_DEF_HPP_
54#include <Xpetra_BlockReorderManager.hpp>
55#include <Xpetra_BlockedCrsMatrix.hpp>
56#include <Xpetra_ReorderedBlockedCrsMatrix.hpp>
57#include <Xpetra_MapExtractor.hpp>
58#include <Xpetra_MapExtractorFactory.hpp>
59#include <Xpetra_Matrix.hpp>
60#include <Xpetra_MatrixUtils.hpp>
68 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
70 RCP<ParameterList> validParamList = rcp(
new ParameterList());
74 validParamList->set< std::string > (
"Reorder Type",
"",
"String describing the reordering of blocks");
77 validParamList->set< RCP<const FactoryBase> >(
"Map1", Teuchos::null,
"Generating factory of the fine level map associated with the (1,1) block in your n x n block matrix.");
78 validParamList->set< RCP<const FactoryBase> >(
"Map2", Teuchos::null,
"Generating factory of the fine level map associated with the (2,2) block in your n x n block matrix.");
79 validParamList->set< RCP<const FactoryBase> >(
"Map3", Teuchos::null,
"Generating factory of the fine level map associated with the (3,3) block in your n x n block matrix.");
80 validParamList->set< RCP<const FactoryBase> >(
"Map4", Teuchos::null,
"Generating factory of the fine level map associated with the (4,4) block in your n x n block matrix.");
81 validParamList->set< RCP<const FactoryBase> >(
"Map5", Teuchos::null,
"Generating factory of the fine level map associated with the (5,5) block in your n x n block matrix.");
82 validParamList->set< RCP<const FactoryBase> >(
"Map6", Teuchos::null,
"Generating factory of the fine level map associated with the (6,6) block in your n x n block matrix.");
83 validParamList->set< RCP<const FactoryBase> >(
"Map7", Teuchos::null,
"Generating factory of the fine level map associated with the (7,7) block in your n x n block matrix.");
84 validParamList->set< RCP<const FactoryBase> >(
"Map8", Teuchos::null,
"Generating factory of the fine level map associated with the (8,8) block in your n x n block matrix.");
85 validParamList->set< RCP<const FactoryBase> >(
"Map9", Teuchos::null,
"Generating factory of the fine level map associated with the (9,9) block in your n x n block matrix.");
87 return validParamList;
90 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 Input(currentLevel,
"A");
95 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
99 const ParameterList& pL = GetParameterList();
100 std::string reorderStr = pL.get<std::string>(
"Reorder Type");
102 RCP<Matrix> Ain = Get<RCP<Matrix> >(currentLevel,
"A");
104 RCP<BlockedCrsMatrix> A = rcp_dynamic_cast<BlockedCrsMatrix>(Ain);
108 if (A == Teuchos::null && currentLevel.
GetLevelID() == 0) {
110 GetOStream(
Warnings0) <<
"Split input matrix (Warning: this is a rather expensive operation)" << std::endl;
112 std::vector<Teuchos::RCP<const Map> > xmaps;
114 for(
int it = 1; it < 10; it++) {
115 std::stringstream ss;
118 RCP<const Map> submap = currentLevel.
Get< RCP<const Map> >(ss.str(),
NoFactory::get());
119 GetOStream(
Runtime1) <<
"Use user-given submap #" << it <<
": length dimension=" << submap->getGlobalNumElements() << std::endl;
120 xmaps.push_back(submap);
124 bool bThyraMode =
false;
125 RCP<const MapExtractor> map_extractor = Xpetra::MapExtractorFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Build(Ain->getRowMap(),xmaps,bThyraMode);
131 Teuchos::RCP<Xpetra::BlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > bOp =
132 Xpetra::MatrixUtils<Scalar, LocalOrdinal, GlobalOrdinal, Node>::SplitMatrix(*Ain,map_extractor,map_extractor,Teuchos::null,bThyraMode);
134 TEUCHOS_TEST_FOR_EXCEPTION(Ain->getGlobalNumRows() != bOp->getGlobalNumRows(),
Exceptions::RuntimeError,
"Split operator not consistent with input operator (different number of rows).");
135 TEUCHOS_TEST_FOR_EXCEPTION(Ain->getLocalNumRows() != bOp->getLocalNumRows(),
Exceptions::RuntimeError,
"Split operator not consistent with input operator (different number of node rows).");
136 TEUCHOS_TEST_FOR_EXCEPTION(Ain->getLocalNumEntries() != bOp->getLocalNumEntries(),
Exceptions::RuntimeError,
"Split operator not consistent with input operator (different number of local entries).");
137 TEUCHOS_TEST_FOR_EXCEPTION(Ain->getGlobalNumEntries() != bOp->getGlobalNumEntries(),
Exceptions::RuntimeError,
"Split operator not consistent with input operator (different number of global entries).");
143 TEUCHOS_TEST_FOR_EXCEPTION(A.is_null(),
Exceptions::BadCast,
"Input matrix A is not a BlockedCrsMatrix.");
144 GetOStream(
Statistics1) <<
"Got a " << A->Rows() <<
"x" << A->Cols() <<
" blocked operator as input" << std::endl;
146 Teuchos::RCP<const Xpetra::BlockReorderManager> brm = Xpetra::blockedReorderFromString(reorderStr);
147 GetOStream(
Debug) <<
"Reordering A using " << brm->toString() << std::endl;
149 Teuchos::RCP<const ReorderedBlockedCrsMatrix> brop =
150 Teuchos::rcp_dynamic_cast<const ReorderedBlockedCrsMatrix>(Xpetra::buildReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>(brm, A));
152 TEUCHOS_TEST_FOR_EXCEPTION(brop.is_null(),
Exceptions::RuntimeError,
"Block reordering of " << A->Rows() <<
"x" << A->Cols() <<
" blocked operator failed.");
154 GetOStream(
Statistics1) <<
"Reordering A using " << brm->toString() <<
" block gives a " << brop->Rows() <<
"x" << brop->Cols() <<
" blocked operators" << std::endl;
155 GetOStream(
Debug) <<
"Reordered operator has " << brop->getRangeMap()->getGlobalNumElements() <<
" rows and " << brop->getDomainMap()->getGlobalNumElements() <<
" columns" << std::endl;
156 GetOStream(
Debug) <<
"Reordered operator: Use of Thyra style gids = " << brop->getRangeMapExtractor()->getThyraMode() << std::endl;
159 Teuchos::RCP<ReorderedBlockedCrsMatrix> bret =
160 Teuchos::rcp_const_cast<ReorderedBlockedCrsMatrix>(brop);
166 currentLevel.
Set(
"A", Teuchos::rcp_dynamic_cast<Matrix>(bret),
this);
Exception indicating invalid cast attempted.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
bool IsAvailable(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need's value has been saved.
int GetLevelID() const
Return level number.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static const RCP< const NoFactory > getRCP()
Static Get() functions.
static const NoFactory * get()
RCP< const ParameterList > GetValidParameterList() const
Input.
void Build(Level ¤tLevel) const
Build an object with this factory.
void DeclareInput(Level ¤tLevel) const
Specifies the data that this class needs, and the factories that generate that data.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Statistics1
Print more statistics.
@ Runtime1
Description of what is happening (more verbose)