47#ifndef PACKAGES_MUELU_SRC_GRAPH_MUELU_UNSMOOSHFACTORY_DEF_HPP_
48#define PACKAGES_MUELU_SRC_GRAPH_MUELU_UNSMOOSHFACTORY_DEF_HPP_
56 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
59 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
61 RCP<ParameterList> validParamList = rcp(
new ParameterList());
62 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory for unamalgamated matrix. Row map of (unamalgamted) output prolongation operator should match row map of this A.");
63 validParamList->set< RCP<const FactoryBase> >(
"P", Teuchos::null,
"Generating factory of the (amalgamated) prolongator P");
64 validParamList->set< RCP<const FactoryBase> >(
"DofStatus", Teuchos::null,
"Generating factory for dofStatus array (usually the VariableDofLaplacdianFactory)");
66 validParamList->set<
int > (
"maxDofPerNode", 1,
"Maximum number of DOFs per node");
67 validParamList->set<
bool > (
"fineIsPadded" ,
false,
"true if finest level input matrix is padded");
69 return validParamList;
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 Input(fineLevel,
"A");
76 Input(coarseLevel,
"P");
81 Input(fineLevel,
"DofStatus");
84 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
87 typedef Teuchos::ScalarTraits<SC> STS;
89 const ParameterList & pL = GetParameterList();
92 RCP<Matrix> unamalgA = Get< RCP<Matrix> >(fineLevel,
"A");
93 RCP<Matrix> amalgP = Get< RCP<Matrix> >(coarseLevel,
"P");
96 int maxDofPerNode = pL.get<
int> (
"maxDofPerNode");
97 bool fineIsPadded = pL.get<
bool>(
"fineIsPadded");
102 Teuchos::Array<char> dofStatus;
104 dofStatus = Get<Teuchos::Array<char> >(fineLevel,
"DofStatus");
107 dofStatus = Teuchos::Array<char>(unamalgA->getRowMap()->getLocalNumElements() ,
's');
109 bool bHasZeroDiagonal =
false;
112 TEUCHOS_TEST_FOR_EXCEPTION(dirOrNot.size() != dofStatus.size(),
MueLu::Exceptions::RuntimeError,
"MueLu::UnsmooshFactory::Build: inconsistent number of coarse DBC array and dofStatus array. dirOrNot.size() = " << dirOrNot.size() <<
" dofStatus.size() = " << dofStatus.size());
113 for(
decltype(dirOrNot.size()) i = 0; i < dirOrNot.size(); ++i) {
114 if(dirOrNot[i] ==
true) dofStatus[i] =
'p';
122 Teuchos::ArrayRCP<const size_t> amalgRowPtr(amalgP->getLocalNumRows());
123 Teuchos::ArrayRCP<const LocalOrdinal> amalgCols(amalgP->getLocalNumEntries());
124 Teuchos::ArrayRCP<const Scalar> amalgVals(amalgP->getLocalNumEntries());
125 Teuchos::RCP<CrsMatrixWrap> amalgPwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(amalgP);
126 Teuchos::RCP<CrsMatrix> amalgPcrs = amalgPwrap->getCrsMatrix();
127 amalgPcrs->getAllValues(amalgRowPtr, amalgCols, amalgVals);
130 size_t paddedNrows = amalgP->getRowMap()->getLocalNumElements() * Teuchos::as<size_t>(maxDofPerNode);
133 Teuchos::ArrayRCP<size_t> newPRowPtr(paddedNrows+1);
134 Teuchos::ArrayRCP<LocalOrdinal> newPCols(amalgP->getLocalNumEntries() * maxDofPerNode);
135 Teuchos::ArrayRCP<Scalar> newPVals(amalgP->getLocalNumEntries() * maxDofPerNode);
138 if(fineIsPadded ==
true || fineLevel.
GetLevelID() > 0) {
147 for (
decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
149 size_t rowLength = amalgRowPtr[i+1] - amalgRowPtr[i];
152 for(
int j = 0; j < maxDofPerNode; j++) {
153 newPRowPtr[i*maxDofPerNode+j] = cnt;
154 if (dofStatus[i*maxDofPerNode+j] ==
's') {
156 for (
size_t k = 0; k < rowLength; k++) {
157 newPCols[cnt ] = amalgCols[k+amalgRowPtr[i]] * maxDofPerNode + j;
158 newPVals[cnt++] = amalgVals[k+amalgRowPtr[i]];
165 newPRowPtr[paddedNrows] = cnt;
166 rowCount = paddedNrows;
174 for (
decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
176 size_t rowLength = amalgRowPtr[i+1] - amalgRowPtr[i];
179 for(
int j = 0; j < maxDofPerNode; j++) {
182 if (dofStatus[i*maxDofPerNode+j] ==
's') {
183 newPRowPtr[rowCount++] = cnt;
185 for (
size_t k = 0; k < rowLength; k++) {
186 newPCols[cnt ] = amalgCols[k+amalgRowPtr[i]] * maxDofPerNode + j;
187 newPVals[cnt++] = amalgVals[k+amalgRowPtr[i]];
191 if (dofStatus[i*maxDofPerNode+j] ==
'd') {
192 newPRowPtr[rowCount++] = cnt;
196 newPRowPtr[rowCount] = cnt;
202 std::vector<size_t> stridingInfo(1, maxDofPerNode);
204 GlobalOrdinal nCoarseDofs = amalgP->getDomainMap()->getLocalNumElements() * maxDofPerNode;
205 GlobalOrdinal indexBase = amalgP->getDomainMap()->getIndexBase();
206 RCP<const Map> coarseDomainMap = StridedMapFactory::Build(amalgP->getDomainMap()->lib(),
207 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
211 amalgP->getDomainMap()->getComm(),
215 size_t nColCoarseDofs = Teuchos::as<size_t>(amalgP->getColMap()->getLocalNumElements() * maxDofPerNode);
216 Teuchos::Array<GlobalOrdinal> unsmooshColMapGIDs(nColCoarseDofs);
217 for(
size_t c = 0; c < amalgP->getColMap()->getLocalNumElements(); ++c) {
218 GlobalOrdinal gid = (amalgP->getColMap()->getGlobalElement(c)-indexBase) * maxDofPerNode + indexBase;
220 for(
int i = 0; i < maxDofPerNode; ++i) {
221 unsmooshColMapGIDs[c * maxDofPerNode + i] = gid + i;
224 Teuchos::RCP<Map> coarseColMap = MapFactory::Build(amalgP->getDomainMap()->lib(),
225 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
226 unsmooshColMapGIDs(),
228 amalgP->getDomainMap()->getComm());
231 Teuchos::RCP<CrsMatrix> unamalgPCrs = CrsMatrixFactory::Build(unamalgA->getRowMap(),
233 maxDofPerNode*amalgP->getLocalMaxNumRowEntries());
234 for (
size_t i = 0; i < rowCount; i++) {
235 unamalgPCrs->insertLocalValues(i,
236 newPCols.view(newPRowPtr[i], newPRowPtr[i+1] - newPRowPtr[i]),
237 newPVals.view(newPRowPtr[i], newPRowPtr[i+1] - newPRowPtr[i]));
239 unamalgPCrs->fillComplete(coarseDomainMap, unamalgA->getRowMap());
241 Teuchos::RCP<Matrix> unamalgP = Teuchos::rcp(
new CrsMatrixWrap(unamalgPCrs));
243 Set(coarseLevel,
"P",unamalgP);
MueLu::DefaultGlobalOrdinal GlobalOrdinal
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.
int GetLevelID() const
Return level number.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
UnsmooshFactory()
Constructor.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Namespace for MueLu classes and methods.