46 #ifndef MUELU_COALESCEDROPFACTORY_DEF_HPP 47 #define MUELU_COALESCEDROPFACTORY_DEF_HPP 49 #include <Xpetra_CrsGraphFactory.hpp> 50 #include <Xpetra_CrsGraph.hpp> 51 #include <Xpetra_ImportFactory.hpp> 52 #include <Xpetra_MapFactory.hpp> 53 #include <Xpetra_Map.hpp> 54 #include <Xpetra_Matrix.hpp> 55 #include <Xpetra_MultiVectorFactory.hpp> 56 #include <Xpetra_MultiVector.hpp> 57 #include <Xpetra_StridedMap.hpp> 58 #include <Xpetra_VectorFactory.hpp> 59 #include <Xpetra_Vector.hpp> 63 #include "MueLu_AmalgamationFactory.hpp" 64 #include "MueLu_AmalgamationInfo.hpp" 67 #include "MueLu_Graph.hpp" 69 #include "MueLu_LWGraph.hpp" 72 #include "MueLu_PreDropFunctionBaseClass.hpp" 73 #include "MueLu_PreDropFunctionConstVal.hpp" 74 #include "MueLu_Utilities.hpp" 78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
80 RCP<ParameterList> validParamList = rcp(
new ParameterList());
82 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 87 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
88 validParamList->getEntry(
"aggregation: drop scheme").setValidator(
89 rcp(
new validatorType(Teuchos::tuple<std::string>(
"classical",
"distance laplacian"),
"aggregation: drop scheme")));
91 #undef SET_VALID_ENTRY 92 validParamList->set<
bool > (
"lightweight wrap",
true,
"Experimental option for lightweight graph access");
94 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
95 validParamList->set< RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory for UnAmalgamationInfo");
96 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory for Coordinates");
98 return validParamList;
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
106 Input(currentLevel,
"A");
107 Input(currentLevel,
"UnAmalgamationInfo");
109 const ParameterList& pL = GetParameterList();
110 if (pL.get<
bool>(
"lightweight wrap") ==
true) {
111 if (pL.get<std::string>(
"aggregation: drop scheme") ==
"distance laplacian")
112 Input(currentLevel,
"Coordinates");
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
121 typedef Xpetra::MultiVector<double,LO,GO,NO> dxMV;
123 typedef Teuchos::ScalarTraits<SC> STS;
125 if (predrop_ != Teuchos::null)
126 GetOStream(
Parameters0) << predrop_->description();
128 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
129 RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(currentLevel,
"UnAmalgamationInfo");
131 const ParameterList & pL = GetParameterList();
132 bool doExperimentalWrap = pL.get<
bool>(
"lightweight wrap");
134 GetOStream(
Parameters0) <<
"lightweight wrap = " << doExperimentalWrap << std::endl;
147 if (doExperimentalWrap) {
148 std::string algo = pL.get<std::string>(
"aggregation: drop scheme");
150 TEUCHOS_TEST_FOR_EXCEPTION(predrop_ != null && algo !=
"classical",
Exceptions::RuntimeError,
"Dropping function must not be provided for \"" << algo <<
"\" algorithm");
151 TEUCHOS_TEST_FOR_EXCEPTION(algo !=
"classical" && algo !=
"distance laplacian",
Exceptions::RuntimeError,
"\"algorithm\" must be one of (classical|distance laplacian)");
153 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
154 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
155 Set(currentLevel,
"Filtering", (threshold != STS::zero()));
157 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<
double>(
"aggregation: Dirichlet threshold")));
159 GO numDropped = 0, numTotal = 0;
160 std::string graphType =
"unamalgamated";
161 if (algo ==
"classical") {
162 if (predrop_ == null) {
167 if (predrop_ != null) {
170 "MueLu::CoalesceFactory::Build: cast to PreDropFunctionConstVal failed.");
172 SC newt = predropConstVal->GetThreshold();
173 if (newt != threshold) {
174 GetOStream(
Warnings0) <<
"switching threshold parameter from " << threshold <<
" (list) to " << newt <<
" (user function" << std::endl;
181 if (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && A->hasCrsGraph()) {
183 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
185 ArrayRCP<const bool > boundaryNodes;
187 graph->SetBoundaryNodeMap(boundaryNodes);
188 numTotal = A->getNodeNumEntries();
191 GO numLocalBoundaryNodes = 0;
192 GO numGlobalBoundaryNodes = 0;
193 for (LO i = 0; i < boundaryNodes.size(); ++i)
194 if (boundaryNodes[i])
195 numLocalBoundaryNodes++;
196 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
197 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
198 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
201 Set(currentLevel,
"DofsPerNode", 1);
202 Set(currentLevel,
"Graph", graph);
204 }
else if ( (A->GetFixedBlockSize() == 1 && threshold != STS::zero()) ||
205 (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && !A->hasCrsGraph())) {
211 ArrayRCP<LO> rows (A->getNodeNumRows()+1);
212 ArrayRCP<LO> columns(A->getNodeNumEntries());
215 const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
216 const ArrayRCP<bool> boundaryNodes(A->getNodeNumRows(),
false);
221 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getNodeNumElements()); ++row) {
222 size_t nnz = A->getNumEntriesInLocalRow(row);
223 ArrayView<const LO> indices;
224 ArrayView<const SC> vals;
225 A->getLocalRowView(row, indices, vals);
233 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
234 LO col = indices[colID];
237 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
238 typename STS::magnitudeType aij = STS::magnitude(vals[colID]*vals[colID]);
240 if (aij > aiiajj || row == col) {
241 columns[realnnz++] = col;
253 boundaryNodes[row] =
true;
255 rows[row+1] = realnnz;
257 columns.resize(realnnz);
258 numTotal = A->getNodeNumEntries();
260 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, A->getRowMap(), A->getColMap(),
"thresholded graph of A"));
261 graph->SetBoundaryNodeMap(boundaryNodes);
263 GO numLocalBoundaryNodes = 0;
264 GO numGlobalBoundaryNodes = 0;
265 for (LO i = 0; i < boundaryNodes.size(); ++i)
266 if (boundaryNodes[i])
267 numLocalBoundaryNodes++;
268 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
269 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
270 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
272 Set(currentLevel,
"Graph", graph);
273 Set(currentLevel,
"DofsPerNode", 1);
275 }
else if (A->GetFixedBlockSize() > 1 && threshold == STS::zero()) {
277 const RCP<const Map> rowMap = A->getRowMap();
278 const RCP<const Map> colMap = A->getColMap();
280 graphType =
"amalgamated";
286 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
287 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
288 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
289 Array<LO> colTranslation = *(amalInfo->getColTranslation());
292 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
295 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
296 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
298 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
304 ArrayRCP<const bool > pointBoundaryNodes;
308 LO blkSize = A->GetFixedBlockSize();
310 LO blkPartSize = A->GetFixedBlockSize();
311 if (A->IsView(
"stridedMaps") ==
true) {
312 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
313 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
315 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
316 blkId = strMap->getStridedBlockId();
318 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
324 Array<LO> indicesExtra;
325 for (LO row = 0; row < numRows; row++) {
326 ArrayView<const LO> indices;
327 indicesExtra.resize(0);
335 bool isBoundary =
false;
337 for (LO j = 0; j < blkPartSize; j++) {
338 if (!pointBoundaryNodes[row*blkPartSize+j]) {
347 MergeRows(*A, row, indicesExtra, colTranslation);
349 indicesExtra.push_back(row);
350 indices = indicesExtra;
351 numTotal += indices.size();
355 LO nnz = indices.size(), rownnz = 0;
356 for (LO colID = 0; colID < nnz; colID++) {
357 LO col = indices[colID];
358 columns[realnnz++] = col;
369 amalgBoundaryNodes[row] =
true;
371 rows[row+1] = realnnz;
373 columns.resize(realnnz);
375 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
376 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
379 GO numLocalBoundaryNodes = 0;
380 GO numGlobalBoundaryNodes = 0;
382 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
383 if (amalgBoundaryNodes[i])
384 numLocalBoundaryNodes++;
386 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
387 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
388 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
389 <<
" agglomerated Dirichlet nodes" << std::endl;
392 Set(currentLevel,
"Graph", graph);
393 Set(currentLevel,
"DofsPerNode", blkSize);
395 }
else if (A->GetFixedBlockSize() > 1 && threshold != STS::zero()) {
397 const RCP<const Map> rowMap = A->getRowMap();
398 const RCP<const Map> colMap = A->getColMap();
400 graphType =
"amalgamated";
406 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
407 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
408 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
409 Array<LO> colTranslation = *(amalInfo->getColTranslation());
412 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
415 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
416 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
418 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
424 ArrayRCP<const bool > pointBoundaryNodes;
428 LO blkSize = A->GetFixedBlockSize();
430 LO blkPartSize = A->GetFixedBlockSize();
431 if (A->IsView(
"stridedMaps") ==
true) {
432 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
433 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
435 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
436 blkId = strMap->getStridedBlockId();
438 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
443 const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
448 Array<LO> indicesExtra;
449 for (LO row = 0; row < numRows; row++) {
450 ArrayView<const LO> indices;
451 indicesExtra.resize(0);
459 bool isBoundary =
false;
461 for (LO j = 0; j < blkPartSize; j++) {
462 if (!pointBoundaryNodes[row*blkPartSize+j]) {
471 MergeRowsWithDropping(*A, row, ghostedDiagVals, threshold, indicesExtra, colTranslation);
473 indicesExtra.push_back(row);
474 indices = indicesExtra;
475 numTotal += indices.size();
479 LO nnz = indices.size(), rownnz = 0;
480 for (LO colID = 0; colID < nnz; colID++) {
481 LO col = indices[colID];
482 columns[realnnz++] = col;
493 amalgBoundaryNodes[row] =
true;
495 rows[row+1] = realnnz;
497 columns.resize(realnnz);
499 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
500 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
503 GO numLocalBoundaryNodes = 0;
504 GO numGlobalBoundaryNodes = 0;
506 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
507 if (amalgBoundaryNodes[i])
508 numLocalBoundaryNodes++;
510 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
511 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
512 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
513 <<
" agglomerated Dirichlet nodes" << std::endl;
516 Set(currentLevel,
"Graph", graph);
517 Set(currentLevel,
"DofsPerNode", blkSize);
520 }
else if (algo ==
"distance laplacian") {
521 LO blkSize = A->GetFixedBlockSize();
522 GO indexBase = A->getRowMap()->getIndexBase();
528 RCP<dxMV> Coords = Get< RCP<Xpetra::MultiVector<double,LO,GO,NO> > >(currentLevel,
"Coordinates");
534 ArrayRCP<const bool > pointBoundaryNodes;
537 if ( (blkSize == 1) && (threshold == STS::zero()) ) {
539 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
540 graph->SetBoundaryNodeMap(pointBoundaryNodes);
541 graphType=
"unamalgamated";
542 numTotal = A->getNodeNumEntries();
545 GO numLocalBoundaryNodes = 0;
546 GO numGlobalBoundaryNodes = 0;
547 for (LO i = 0; i < pointBoundaryNodes.size(); ++i)
548 if (pointBoundaryNodes[i])
549 numLocalBoundaryNodes++;
550 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
551 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
552 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
555 Set(currentLevel,
"DofsPerNode", blkSize);
556 Set(currentLevel,
"Graph", graph);
571 TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getNodeNumElements()/blkSize != Coords->getLocalLength(),
Exceptions::Incompatible,
572 "Coordinate vector length (" << Coords->getLocalLength() <<
") is incompatible with number of rows in A (" << A->getRowMap()->getNodeNumElements() <<
") by modulo block size ("<< blkSize <<
").");
574 const RCP<const Map> colMap = A->getColMap();
575 RCP<const Map> uniqueMap, nonUniqueMap;
576 Array<LO> colTranslation;
578 uniqueMap = A->getRowMap();
579 nonUniqueMap = A->getColMap();
580 graphType=
"unamalgamated";
583 uniqueMap = Coords->getMap();
585 "Different index bases for matrix and coordinates");
589 graphType =
"amalgamated";
591 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
593 RCP<dxMV> ghostedCoords;
594 RCP<Vector> ghostedLaplDiag;
595 Teuchos::ArrayRCP<SC> ghostedLaplDiagData;
596 if (threshold != STS::zero()) {
598 RCP<const Import> importer;
601 if (blkSize == 1 && A->getCrsGraph()->getImporter() != Teuchos::null) {
602 GetOStream(
Warnings1) <<
"Using existing importer from matrix graph" << std::endl;
603 importer = A->getCrsGraph()->getImporter();
605 GetOStream(
Warnings0) <<
"Constructing new importer instance" << std::endl;
606 importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
609 ghostedCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(nonUniqueMap, Coords->getNumVectors());
612 ghostedCoords->doImport(*Coords, *importer, Xpetra::INSERT);
616 RCP<Vector> localLaplDiag = VectorFactory::Build(uniqueMap);
617 ArrayRCP<SC> localLaplDiagData = localLaplDiag->getDataNonConst(0);
618 Array<LO> indicesExtra;
619 Teuchos::Array<Teuchos::ArrayRCP<const double>> coordData;
620 if (threshold != STS::zero()) {
621 const size_t numVectors = ghostedCoords->getNumVectors();
622 coordData.reserve(numVectors);
623 for (
size_t j = 0; j < numVectors; j++) {
624 Teuchos::ArrayRCP<const double> tmpData=ghostedCoords->getData(j);
625 coordData.push_back(tmpData);
630 for (LO row = 0; row < numRows; row++) {
631 ArrayView<const LO> indices;
634 ArrayView<const SC> vals;
635 A->getLocalRowView(row, indices, vals);
639 indicesExtra.resize(0);
640 MergeRows(*A, row, indicesExtra, colTranslation);
641 indices = indicesExtra;
644 LO nnz = indices.size();
645 for (LO colID = 0; colID < nnz; colID++) {
646 const LO col = indices[colID];
656 ghostedLaplDiag = VectorFactory::Build(nonUniqueMap);
657 ghostedLaplDiag->doImport(*localLaplDiag, *importer, Xpetra::INSERT);
658 ghostedLaplDiagData = ghostedLaplDiag->getDataNonConst(0);
662 GetOStream(
Runtime0) <<
"Skipping distance laplacian construction due to 0 threshold" << std::endl;
668 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
669 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
671 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
675 Array<LO> indicesExtra;
678 Teuchos::Array<Teuchos::ArrayRCP<const double>> coordData;
679 if (threshold != STS::zero()) {
680 const size_t numVectors = ghostedCoords->getNumVectors();
681 coordData.reserve(numVectors);
682 for (
size_t j = 0; j < numVectors; j++) {
683 Teuchos::ArrayRCP<const double> tmpData=ghostedCoords->getData(j);
684 coordData.push_back(tmpData);
687 for (LO row = 0; row < numRows; row++) {
688 ArrayView<const LO> indices;
689 indicesExtra.resize(0);
692 ArrayView<const SC> vals;
693 A->getLocalRowView(row, indices, vals);
697 bool isBoundary =
false;
699 for (LO j = 0; j < blkSize; j++) {
700 if (!pointBoundaryNodes[row*blkSize+j]) {
708 MergeRows(*A, row, indicesExtra, colTranslation);
710 indicesExtra.push_back(row);
711 indices = indicesExtra;
713 numTotal += indices.size();
715 LO nnz = indices.size(), rownnz = 0;
716 if (threshold != STS::zero()) {
717 for (LO colID = 0; colID < nnz; colID++) {
718 LO col = indices[colID];
721 columns[realnnz++] = col;
727 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
728 typename STS::magnitudeType aij = STS::magnitude(laplVal*laplVal);
731 columns[realnnz++] = col;
740 for (LO colID = 0; colID < nnz; colID++) {
741 LO col = indices[colID];
742 columns[realnnz++] = col;
754 amalgBoundaryNodes[row] =
true;
756 rows[row+1] = realnnz;
759 columns.resize(realnnz);
761 RCP<GraphBase> graph;
764 graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
765 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
769 GO numLocalBoundaryNodes = 0;
770 GO numGlobalBoundaryNodes = 0;
772 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
773 if (amalgBoundaryNodes[i])
774 numLocalBoundaryNodes++;
776 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
777 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
778 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" agglomerated Dirichlet nodes" 779 <<
" using threshold " << dirichletThreshold << std::endl;
782 Set(currentLevel,
"Graph", graph);
783 Set(currentLevel,
"DofsPerNode", blkSize);
787 if ((GetVerbLevel() &
Statistics1) && !(A->GetFixedBlockSize() > 1 && threshold != STS::zero())) {
788 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
789 GO numGlobalTotal, numGlobalDropped;
792 GetOStream(
Statistics1) <<
"Number of dropped entries in " << graphType <<
" matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
793 if (numGlobalTotal != 0)
794 GetOStream(
Statistics1) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
801 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
803 GetOStream(
Runtime0) <<
"algorithm = \"" <<
"failsafe" <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
804 Set(currentLevel,
"Filtering", (threshold != STS::zero()));
806 RCP<const Map> rowMap = A->getRowMap();
807 RCP<const Map> colMap = A->getColMap();
810 GO indexBase = rowMap->getIndexBase();
814 if(A->IsView(
"stridedMaps") &&
815 Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
816 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
817 RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap());
818 TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
819 blockdim = strMap->getFixedBlockSize();
820 offset = strMap->getOffset();
821 oldView = A->SwitchToView(oldView);
822 GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build():" <<
" found blockdim=" << blockdim <<
" from strided maps. offset=" << offset << std::endl;
823 }
else GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
827 RCP<const Map> nodeMap = amalInfo->getNodeRowMap();
828 GetOStream(
Statistics1) <<
"CoalesceDropFactory: nodeMap " << nodeMap->getNodeNumElements() <<
"/" << nodeMap->getGlobalNumElements() <<
" elements" << std::endl;
831 RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, 10, Xpetra::DynamicProfile);
833 LO numRows = A->getRowMap()->getNodeNumElements();
834 LO numNodes = nodeMap->getNodeNumElements();
835 const ArrayRCP<bool> amalgBoundaryNodes(numNodes,
false);
836 const ArrayRCP<int> numberDirichletRowsPerNode(numNodes, 0);
837 bool bIsDiagonalEntry =
false;
842 for(LO row=0; row<numRows; row++) {
844 GO grid = rowMap->getGlobalElement(row);
847 bIsDiagonalEntry =
false;
852 size_t nnz = A->getNumEntriesInLocalRow(row);
853 Teuchos::ArrayView<const LO> indices;
854 Teuchos::ArrayView<const SC> vals;
855 A->getLocalRowView(row, indices, vals);
857 RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(
new std::vector<GO>);
859 for(LO col=0; col<Teuchos::as<LO>(nnz); col++) {
860 GO gcid = colMap->getGlobalElement(indices[col]);
864 cnodeIds->push_back(cnodeId);
866 if (grid == gcid) bIsDiagonalEntry =
true;
870 if(realnnz == 1 && bIsDiagonalEntry ==
true) {
871 LO lNodeId = nodeMap->getLocalElement(nodeId);
872 numberDirichletRowsPerNode[lNodeId] += 1;
873 if (numberDirichletRowsPerNode[lNodeId] == blockdim)
874 amalgBoundaryNodes[lNodeId] =
true;
877 Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp( cnodeIds );
879 if(arr_cnodeIds.size() > 0 )
880 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
883 crsGraph->fillComplete(nodeMap,nodeMap);
886 RCP<GraphBase> graph = rcp(
new Graph(crsGraph,
"amalgamated graph of A"));
889 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
892 GO numLocalBoundaryNodes = 0;
893 GO numGlobalBoundaryNodes = 0;
894 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
895 if (amalgBoundaryNodes[i])
896 numLocalBoundaryNodes++;
897 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
898 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
899 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
904 Set(currentLevel,
"DofsPerNode", blockdim);
905 Set(currentLevel,
"Graph", graph);
911 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
913 typedef typename ArrayView<const LO>::size_type size_type;
916 LO blkSize = A.GetFixedBlockSize();
917 if (A.IsView(
"stridedMaps") ==
true) {
918 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
919 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
921 if (strMap->getStridedBlockId() > -1)
922 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
926 size_t nnz = 0, pos = 0;
927 for (LO j = 0; j < blkSize; j++)
928 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
938 ArrayView<const LO> inds;
939 ArrayView<const SC> vals;
940 for (LO j = 0; j < blkSize; j++) {
941 A.getLocalRowView(row*blkSize+j, inds, vals);
942 size_type numIndices = inds.size();
948 cols[pos++] = translation[inds[0]];
949 for (size_type k = 1; k < numIndices; k++) {
950 LO nodeID = translation[inds[k]];
954 if (nodeID != cols[pos-1])
955 cols[pos++] = nodeID;
962 std::sort(cols.begin(), cols.end());
964 for (
size_t j = 1; j < nnz; j++)
965 if (cols[j] != cols[pos])
966 cols[++pos] = cols[j];
970 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
972 typedef typename ArrayView<const LO>::size_type size_type;
973 typedef Teuchos::ScalarTraits<SC> STS;
976 LO blkSize = A.GetFixedBlockSize();
977 if (A.IsView(
"stridedMaps") ==
true) {
978 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
979 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
981 if (strMap->getStridedBlockId() > -1)
982 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
986 size_t nnz = 0, pos = 0;
987 for (LO j = 0; j < blkSize; j++)
988 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
998 ArrayView<const LO> inds;
999 ArrayView<const SC> vals;
1000 for (LO j = 0; j < blkSize; j++) {
1001 A.getLocalRowView(row*blkSize+j, inds, vals);
1002 size_type numIndices = inds.size();
1004 if (numIndices == 0)
1009 for (size_type k = 0; k < numIndices; k++) {
1011 LO nodeID = translation[inds[k]];
1014 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold*ghostedDiagVals[dofID]*ghostedDiagVals[row*blkSize+j]);
1015 typename STS::magnitudeType aij = STS::magnitude(vals[k]*vals[k]);
1018 if (aij > aiiajj || (row*blkSize+j == dofID)) {
1024 if (nodeID != prevNodeID) {
1025 cols[pos++] = nodeID;
1026 prevNodeID = nodeID;
1035 std::sort(cols.begin(), cols.end());
1037 for (
size_t j = 1; j < nnz; j++)
1038 if (cols[j] != cols[pos])
1039 cols[++pos] = cols[j];
1046 #endif // MUELU_COALESCEDROPFACTORY_DEF_HPP Important warning messages (one line)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
MueLu representation of a compressed row storage graph.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
Timer to be used in factories. Similar to Monitor but with additional timers.
One-liner description of what is happening.
Namespace for MueLu classes and methods.
Exception throws to report incompatible objects (like maps).
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
translate global (row/column) id to global amalgamation block id
CoalesceDropFactory()
Constructor.
void DeclareInput(Level ¤tLevel) const
Input.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
#define SET_VALID_ENTRY(name)
static void AmalgamateMap(const Map &sourceMap, const Matrix &A, RCP< const Map > &amalgamatedMap, Array< LO > &translation)
Method to create merged map for systems of PDEs.
void MergeRowsWithDropping(const Matrix &A, const LO row, const ArrayRCP< const SC > &ghostedDiagVals, SC threshold, Array< LO > &cols, const Array< LO > &translation) const
void Build(Level ¤tLevel) const
Build an object with this factory.
void MergeRows(const Matrix &A, const LO row, Array< LO > &cols, const Array< LO > &translation) const
Method to merge rows of matrix for systems of PDEs.
static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Lightweight MueLu representation of a compressed row storage graph.
Exception throws to report errors in the internal logical of the program.