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;
103 if(fineLevel.GetLevelID() == 0) {
104 dofStatus = Get<Teuchos::Array<char> >(fineLevel,
"DofStatus");
105 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(dofStatus.size()) == Teuchos::as<size_t>(unamalgA->getRowMap()->getNodeNumElements()),
MueLu::Exceptions::RuntimeError,
"MueLu::UnsmooshFactory::Build: User provided dofStatus on level 0 does not fit to size of unamalgamted A");
108 dofStatus = Teuchos::Array<char>(unamalgA->getRowMap()->getNodeNumElements() ,
's');
110 bool bHasZeroDiagonal =
false;
113 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());
114 for(decltype(dirOrNot.size()) i = 0; i < dirOrNot.size(); ++i) {
115 if(dirOrNot[i] ==
true) dofStatus[i] =
'p';
123 Teuchos::ArrayRCP<const size_t> amalgRowPtr(amalgP->getNodeNumRows());
124 Teuchos::ArrayRCP<const LocalOrdinal> amalgCols(amalgP->getNodeNumEntries());
125 Teuchos::ArrayRCP<const Scalar> amalgVals(amalgP->getNodeNumEntries());
126 Teuchos::RCP<CrsMatrixWrap> amalgPwrap = Teuchos::rcp_dynamic_cast<CrsMatrixWrap>(amalgP);
127 Teuchos::RCP<CrsMatrix> amalgPcrs = amalgPwrap->getCrsMatrix();
128 amalgPcrs->getAllValues(amalgRowPtr, amalgCols, amalgVals);
131 size_t paddedNrows = amalgP->getRowMap()->getNodeNumElements() * Teuchos::as<size_t>(maxDofPerNode);
134 Teuchos::ArrayRCP<size_t> newPRowPtr(paddedNrows+1);
135 Teuchos::ArrayRCP<LocalOrdinal> newPCols(amalgP->getNodeNumEntries() * maxDofPerNode);
136 Teuchos::ArrayRCP<Scalar> newPVals(amalgP->getNodeNumEntries() * maxDofPerNode);
139 if(fineIsPadded ==
true || fineLevel.GetLevelID() > 0) {
148 for (decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
150 size_t rowLength = amalgRowPtr[i+1] - amalgRowPtr[i];
153 for(
int j = 0; j < maxDofPerNode; j++) {
154 newPRowPtr[i*maxDofPerNode+j] = cnt;
155 if (dofStatus[i*maxDofPerNode+j] ==
's') {
157 for (
size_t k = 0; k < rowLength; k++) {
158 newPCols[cnt ] = amalgCols[k+amalgRowPtr[i]] * maxDofPerNode + j;
159 newPVals[cnt++] = amalgVals[k+amalgRowPtr[i]];
166 newPRowPtr[paddedNrows] = cnt;
167 rowCount = paddedNrows;
175 for (decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
177 size_t rowLength = amalgRowPtr[i+1] - amalgRowPtr[i];
180 for(
int j = 0; j < maxDofPerNode; j++) {
183 if (dofStatus[i*maxDofPerNode+j] ==
's') {
184 newPRowPtr[rowCount++] = cnt;
186 for (
size_t k = 0; k < rowLength; k++) {
187 newPCols[cnt ] = amalgCols[k+amalgRowPtr[i]] * maxDofPerNode + j;
188 newPVals[cnt++] = amalgVals[k+amalgRowPtr[i]];
192 if (dofStatus[i*maxDofPerNode+j] ==
'd') {
193 newPRowPtr[rowCount++] = cnt;
197 newPRowPtr[rowCount] = cnt;
203 std::vector<size_t> stridingInfo(1,maxDofPerNode);
205 GlobalOrdinal nCoarseDofs = amalgP->getDomainMap()->getNodeNumElements() * maxDofPerNode;
206 GlobalOrdinal indexBase = amalgP->getDomainMap()->getIndexBase();
207 RCP<const Map> coarseDomainMap = StridedMapFactory::Build(amalgP->getDomainMap()->lib(),
208 Teuchos::OrdinalTraits<Xpetra::global_size_t>::invalid(),
212 amalgP->getDomainMap()->getComm(),
216 size_t nColCoarseDofs = Teuchos::as<size_t>(amalgP->getColMap()->getNodeNumElements() * maxDofPerNode);
217 Teuchos::Array<GlobalOrdinal> unsmooshColMapGIDs(nColCoarseDofs);
218 for(
size_t c = 0; c < amalgP->getColMap()->getNodeNumElements(); ++c) {
219 GlobalOrdinal gid = amalgP->getColMap()->getGlobalElement(c) * maxDofPerNode;
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(),coarseColMap, 3);
232 for (decltype(rowCount) i = 0; i < rowCount; i++) {
233 unamalgPCrs->insertLocalValues(i, newPCols.view(newPRowPtr[i],newPRowPtr[i+1]-newPRowPtr[i]),
234 newPVals.view(newPRowPtr[i],newPRowPtr[i+1]-newPRowPtr[i]));
236 unamalgPCrs->fillComplete(coarseDomainMap,unamalgA->getRowMap());
238 Teuchos::RCP<Matrix> unamalgP = Teuchos::rcp(
new CrsMatrixWrap(unamalgPCrs));
240 Set(coarseLevel,
"P",unamalgP);
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
int GetLevelID() const
Return level number.
Class that holds all level-specific information.
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
UnsmooshFactory()
Constructor.
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Exception throws to report errors in the internal logical of the program.