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>
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>
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");
103 if (fineLevel.GetLevelID() == 0) {
104 dofStatus = Get<Teuchos::Array<char> >(fineLevel,
"DofStatus");
109 bool bHasZeroDiagonal =
false;
113 for (decltype(dirOrNot.
size()) i = 0; i < dirOrNot.
size(); ++i) {
114 if (dirOrNot[i] ==
true) dofStatus[i] =
'p';
127 amalgPcrs->getAllValues(amalgRowPtr, amalgCols, amalgVals);
130 size_t paddedNrows = amalgP->getRowMap()->getLocalNumElements() * Teuchos::as<size_t>(maxDofPerNode);
138 if (fineIsPadded ==
true || fineLevel.GetLevelID() > 0) {
146 for (decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
148 size_t rowLength = amalgRowPtr[i + 1] - amalgRowPtr[i];
151 for (
int j = 0; j < maxDofPerNode; j++) {
152 newPRowPtr[i * maxDofPerNode + j] = cnt;
153 if (dofStatus[i * maxDofPerNode + j] ==
's') {
155 for (
size_t k = 0; k < rowLength; k++) {
156 newPCols[cnt] = amalgCols[k + amalgRowPtr[i]] * maxDofPerNode + j;
157 newPVals[cnt++] = amalgVals[k + amalgRowPtr[i]];
163 newPRowPtr[paddedNrows] = cnt;
164 rowCount = paddedNrows;
172 for (decltype(amalgRowPtr.size()) i = 0; i < amalgRowPtr.size() - 1; i++) {
174 size_t rowLength = amalgRowPtr[i + 1] - amalgRowPtr[i];
177 for (
int j = 0; j < maxDofPerNode; j++) {
180 if (dofStatus[i * maxDofPerNode + j] ==
's') {
181 newPRowPtr[rowCount++] = cnt;
183 for (
size_t k = 0; k < rowLength; k++) {
184 newPCols[cnt] = amalgCols[k + amalgRowPtr[i]] * maxDofPerNode + j;
185 newPVals[cnt++] = amalgVals[k + amalgRowPtr[i]];
188 if (dofStatus[i * maxDofPerNode + j] ==
'd') {
189 newPRowPtr[rowCount++] = cnt;
193 newPRowPtr[rowCount] = cnt;
199 std::vector<size_t> stridingInfo(1, maxDofPerNode);
201 GlobalOrdinal nCoarseDofs = amalgP->getDomainMap()->getLocalNumElements() * maxDofPerNode;
202 GlobalOrdinal indexBase = amalgP->getDomainMap()->getIndexBase();
203 RCP<const Map> coarseDomainMap = StridedMapFactory::Build(amalgP->getDomainMap()->lib(),
208 amalgP->getDomainMap()->getComm(),
212 size_t nColCoarseDofs = Teuchos::as<size_t>(amalgP->getColMap()->getLocalNumElements() * maxDofPerNode);
214 for (
size_t c = 0; c < amalgP->getColMap()->getLocalNumElements(); ++c) {
215 GlobalOrdinal gid = (amalgP->getColMap()->getGlobalElement(c) - indexBase) * maxDofPerNode + indexBase;
217 for (
int i = 0; i < maxDofPerNode; ++i) {
218 unsmooshColMapGIDs[c * maxDofPerNode + i] = gid + i;
221 Teuchos::RCP<Map> coarseColMap = MapFactory::Build(amalgP->getDomainMap()->lib(),
223 unsmooshColMapGIDs(),
225 amalgP->getDomainMap()->getComm());
230 maxDofPerNode * amalgP->getLocalMaxNumRowEntries());
231 for (
size_t i = 0; i < rowCount; i++) {
232 unamalgPCrs->insertLocalValues(i,
233 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());
240 Set(coarseLevel,
"P", unamalgP);
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)
void Build(Level &fineLevel, Level &coarseLevel) const
Build an object with this factory.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultGlobalOrdinal GlobalOrdinal
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
Class that holds all level-specific information.
UnsmooshFactory()
Constructor.
static Teuchos::ArrayRCP< const bool > DetectDirichletRowsExt(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, bool &bHasZeroDiagonal, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::zero())
Detect Dirichlet rows (extended version)
int GetLevelID() const
Return level number.
Exception throws to report errors in the internal logical of the program.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.