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);