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");
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;
128 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
129 RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(currentLevel,
"UnAmalgamationInfo");
132 bool doExperimentalWrap = pL.get<
bool>(
"lightweight wrap");
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") {
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;
183 if (A->GetFixedBlockSize() == 1 && threshold == STS::zero()) {
185 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
188 ArrayRCP<const bool > boundaryNodes;
190 graph->SetBoundaryNodeMap(boundaryNodes);
191 numTotal = A->getNodeNumEntries();
194 GO numLocalBoundaryNodes = 0;
195 GO numGlobalBoundaryNodes = 0;
196 for (LO i = 0; i < boundaryNodes.size(); ++i)
197 if (boundaryNodes[i])
198 numLocalBoundaryNodes++;
199 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
200 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
204 Set(currentLevel,
"DofsPerNode", 1);
205 Set(currentLevel,
"Graph", graph);
207 }
else if (A->GetFixedBlockSize() == 1 && threshold != STS::zero()) {
212 ArrayRCP<LO> rows (A->getNodeNumRows()+1);
213 ArrayRCP<LO> columns(A->getNodeNumEntries());
216 const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
217 const ArrayRCP<bool> boundaryNodes(A->getNodeNumRows(),
false);
222 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getNodeNumElements()); ++row) {
223 size_t nnz = A->getNumEntriesInLocalRow(row);
224 ArrayView<const LO> indices;
225 ArrayView<const SC> vals;
226 A->getLocalRowView(row, indices, vals);
234 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
235 LO col = indices[colID];
238 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
239 typename STS::magnitudeType aij = STS::magnitude(vals[colID]*vals[colID]);
241 if (aij > aiiajj || row == col) {
242 columns[realnnz++] = col;
254 boundaryNodes[row] =
true;
256 rows[row+1] = realnnz;
258 columns.resize(realnnz);
260 numTotal = A->getNodeNumEntries();
262 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, A->getRowMap(), A->getColMap(),
"thresholded graph of A"));
263 graph->SetBoundaryNodeMap(boundaryNodes);
265 GO numLocalBoundaryNodes = 0;
266 GO numGlobalBoundaryNodes = 0;
267 for (LO i = 0; i < boundaryNodes.size(); ++i)
268 if (boundaryNodes[i])
269 numLocalBoundaryNodes++;
270 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
271 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
274 Set(currentLevel,
"Graph", graph);
275 Set(currentLevel,
"DofsPerNode", 1);
277 }
else if (A->GetFixedBlockSize() > 1 && threshold == STS::zero()) {
280 const RCP<const Map> rowMap = A->getRowMap();
281 const RCP<const Map> colMap = A->getColMap();
283 graphType =
"amalgamated";
289 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
290 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
291 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
292 Array<LO> colTranslation = *(amalInfo->getColTranslation());
295 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
298 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
299 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
301 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
307 ArrayRCP<const bool > pointBoundaryNodes;
311 LO blkSize = A->GetFixedBlockSize();
313 LO blkPartSize = A->GetFixedBlockSize();
314 if (A->IsView(
"stridedMaps") ==
true) {
315 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
316 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
318 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
319 blkId = strMap->getStridedBlockId();
321 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
327 Array<LO> indicesExtra;
328 for (LO row = 0; row < numRows; row++) {
329 ArrayView<const LO> indices;
330 indicesExtra.resize(0);
338 bool isBoundary =
false;
340 for (LO j = 0; j < blkPartSize; j++) {
341 if (!pointBoundaryNodes[row*blkPartSize+j]) {
350 MergeRows(*A, row, indicesExtra, colTranslation);
352 indicesExtra.push_back(row);
353 indices = indicesExtra;
354 numTotal += indices.size();
358 LO nnz = indices.size(), rownnz = 0;
359 for (LO colID = 0; colID < nnz; colID++) {
360 LO col = indices[colID];
361 columns[realnnz++] = col;
372 amalgBoundaryNodes[row] =
true;
374 rows[row+1] = realnnz;
376 columns.resize(realnnz);
378 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
379 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
382 GO numLocalBoundaryNodes = 0;
383 GO numGlobalBoundaryNodes = 0;
385 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
386 if (amalgBoundaryNodes[i])
387 numLocalBoundaryNodes++;
389 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
390 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
392 <<
" agglomerated Dirichlet nodes" << std::endl;
395 Set(currentLevel,
"Graph", graph);
396 Set(currentLevel,
"DofsPerNode", blkSize);
398 }
else if (A->GetFixedBlockSize() > 1 && threshold != STS::zero()) {
401 const RCP<const Map> rowMap = A->getRowMap();
402 const RCP<const Map> colMap = A->getColMap();
404 graphType =
"amalgamated";
410 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
411 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
412 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
413 Array<LO> colTranslation = *(amalInfo->getColTranslation());
416 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
419 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
420 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
422 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
428 ArrayRCP<const bool > pointBoundaryNodes;
432 LO blkSize = A->GetFixedBlockSize();
434 LO blkPartSize = A->GetFixedBlockSize();
435 if (A->IsView(
"stridedMaps") ==
true) {
436 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
437 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
439 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
440 blkId = strMap->getStridedBlockId();
442 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
447 const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
452 Array<LO> indicesExtra;
453 for (LO row = 0; row < numRows; row++) {
454 ArrayView<const LO> indices;
455 indicesExtra.resize(0);
463 bool isBoundary =
false;
465 for (LO j = 0; j < blkPartSize; j++) {
466 if (!pointBoundaryNodes[row*blkPartSize+j]) {
477 indicesExtra.push_back(row);
478 indices = indicesExtra;
479 numTotal += indices.size();
483 LO nnz = indices.size(), rownnz = 0;
484 for (LO colID = 0; colID < nnz; colID++) {
485 LO col = indices[colID];
486 columns[realnnz++] = col;
497 amalgBoundaryNodes[row] =
true;
499 rows[row+1] = realnnz;
501 columns.resize(realnnz);
503 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
504 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
507 GO numLocalBoundaryNodes = 0;
508 GO numGlobalBoundaryNodes = 0;
510 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
511 if (amalgBoundaryNodes[i])
512 numLocalBoundaryNodes++;
514 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
515 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
517 <<
" agglomerated Dirichlet nodes" << std::endl;
520 Set(currentLevel,
"Graph", graph);
521 Set(currentLevel,
"DofsPerNode", blkSize);
524 }
else if (algo ==
"distance laplacian") {
525 LO blkSize = A->GetFixedBlockSize();
526 GO indexBase = A->getRowMap()->getIndexBase();
532 RCP<dxMV> Coords = Get< RCP<Xpetra::MultiVector<double,LO,GO,NO> > >(currentLevel,
"Coordinates");
538 ArrayRCP<const bool > pointBoundaryNodes;
541 if ( (blkSize == 1) && (threshold == STS::zero()) ) {
543 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
544 graph->SetBoundaryNodeMap(pointBoundaryNodes);
545 graphType=
"unamalgamated";
546 numTotal = A->getNodeNumEntries();
549 GO numLocalBoundaryNodes = 0;
550 GO numGlobalBoundaryNodes = 0;
551 for (LO i = 0; i < pointBoundaryNodes.size(); ++i)
552 if (pointBoundaryNodes[i])
553 numLocalBoundaryNodes++;
554 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
555 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
559 Set(currentLevel,
"DofsPerNode", blkSize);
560 Set(currentLevel,
"Graph", graph);
575 TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getNodeNumElements()/blkSize != Coords->getLocalLength(),
Exceptions::Incompatible,
576 "Coordinate vector length (" << Coords->getLocalLength() <<
") is incompatible with number of rows in A (" << A->getRowMap()->getNodeNumElements() <<
") by modulo block size ("<< blkSize <<
").");
578 const RCP<const Map> colMap = A->getColMap();
579 RCP<const Map> uniqueMap, nonUniqueMap;
580 Array<LO> colTranslation;
582 uniqueMap = A->getRowMap();
583 nonUniqueMap = A->getColMap();
584 graphType=
"unamalgamated";
587 uniqueMap = Coords->getMap();
589 "Different index bases for matrix and coordinates");
593 graphType =
"amalgamated";
595 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
597 RCP<dxMV> ghostedCoords;
598 RCP<Vector> ghostedLaplDiag;
599 Teuchos::ArrayRCP<SC> ghostedLaplDiagData;
600 if (threshold != STS::zero()) {
602 RCP<const Import> importer;
605 importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
607 ghostedCoords = Xpetra::MultiVectorFactory<double,LO,GO,NO>::Build(nonUniqueMap, Coords->getNumVectors());
608 ghostedCoords->doImport(*Coords, *importer, Xpetra::INSERT);
611 RCP<Vector> localLaplDiag = VectorFactory::Build(uniqueMap);
612 ArrayRCP<SC> localLaplDiagData = localLaplDiag->getDataNonConst(0);
613 Array<LO> indicesExtra;
614 for (LO row = 0; row < numRows; row++) {
615 ArrayView<const LO> indices;
616 indicesExtra.resize(0);
619 ArrayView<const SC> vals;
620 A->getLocalRowView(row, indices, vals);
624 MergeRows(*A, row, indicesExtra, colTranslation);
625 indices = indicesExtra;
628 LO nnz = indices.size();
629 for (LO colID = 0; colID < nnz; colID++) {
630 LO col = indices[colID];
636 ghostedLaplDiag = VectorFactory::Build(nonUniqueMap);
637 ghostedLaplDiag->doImport(*localLaplDiag, *importer, Xpetra::INSERT);
638 ghostedLaplDiagData = ghostedLaplDiag->getDataNonConst(0);
641 GetOStream(
Runtime0) <<
"Skipping distance laplacian construction due to 0 threshold" << std::endl;
647 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
648 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
650 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
654 Array<LO> indicesExtra;
655 for (LO row = 0; row < numRows; row++) {
656 ArrayView<const LO> indices;
657 indicesExtra.resize(0);
660 ArrayView<const SC> vals;
661 A->getLocalRowView(row, indices, vals);
665 bool isBoundary =
false;
667 for (LO j = 0; j < blkSize; j++) {
668 if (!pointBoundaryNodes[row*blkSize+j]) {
676 MergeRows(*A, row, indicesExtra, colTranslation);
678 indicesExtra.push_back(row);
679 indices = indicesExtra;
681 numTotal += indices.size();
683 LO nnz = indices.size(), rownnz = 0;
684 if (threshold != STS::zero()) {
685 for (LO colID = 0; colID < nnz; colID++) {
686 LO col = indices[colID];
689 columns[realnnz++] = col;
695 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
696 typename STS::magnitudeType aij = STS::magnitude(laplVal*laplVal);
699 columns[realnnz++] = col;
708 for (LO colID = 0; colID < nnz; colID++) {
709 LO col = indices[colID];
710 columns[realnnz++] = col;
722 amalgBoundaryNodes[row] =
true;
724 rows[row+1] = realnnz;
726 columns.resize(realnnz);
728 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
729 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
732 GO numLocalBoundaryNodes = 0;
733 GO numGlobalBoundaryNodes = 0;
735 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
736 if (amalgBoundaryNodes[i])
737 numLocalBoundaryNodes++;
739 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
740 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
742 <<
" using threshold " << dirichletThreshold << std::endl;
745 Set(currentLevel,
"Graph", graph);
746 Set(currentLevel,
"DofsPerNode", blkSize);
751 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
752 GO numGlobalTotal, numGlobalDropped;
755 GetOStream(
Statistics0) <<
"Number of dropped entries in " << graphType <<
" matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
756 if (numGlobalTotal != 0)
757 GetOStream(
Statistics0) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
765 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
767 GetOStream(
Runtime0) <<
"algorithm = \"" <<
"failsafe" <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
768 Set(currentLevel,
"Filtering", (threshold != STS::zero()));
770 RCP<const Map> rowMap = A->getRowMap();
771 RCP<const Map> colMap = A->getColMap();
774 GO indexBase = rowMap->getIndexBase();
778 if(A->IsView(
"stridedMaps") &&
779 Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
780 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
781 RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap());
782 TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
783 blockdim = strMap->getFixedBlockSize();
784 offset = strMap->getOffset();
785 oldView = A->SwitchToView(oldView);
786 GetOStream(
Statistics0) <<
"CoalesceDropFactory::Build():" <<
" found blockdim=" << blockdim <<
" from strided maps. offset=" << offset << std::endl;
787 }
else GetOStream(
Statistics0) <<
"CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
791 RCP<const Map> nodeMap = amalInfo->getNodeRowMap();
792 GetOStream(
Statistics0) <<
"CoalesceDropFactory: nodeMap " << nodeMap->getNodeNumElements() <<
"/" << nodeMap->getGlobalNumElements() <<
" elements" << std::endl;
795 RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, 10, Xpetra::DynamicProfile);
797 LO numRows = A->getRowMap()->getNodeNumElements();
798 LO numNodes = nodeMap->getNodeNumElements();
799 const ArrayRCP<bool> amalgBoundaryNodes(numNodes,
false);
800 const ArrayRCP<int> numberDirichletRowsPerNode(numNodes, 0);
801 bool bIsDiagonalEntry =
false;
806 for(LO row=0; row<numRows; row++) {
808 GO grid = rowMap->getGlobalElement(row);
811 bIsDiagonalEntry =
false;
816 size_t nnz = A->getNumEntriesInLocalRow(row);
817 Teuchos::ArrayView<const LO> indices;
818 Teuchos::ArrayView<const SC> vals;
819 A->getLocalRowView(row, indices, vals);
821 RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(
new std::vector<GO>);
823 for(LO col=0; col<Teuchos::as<LO>(nnz); col++) {
824 GO gcid = colMap->getGlobalElement(indices[col]);
828 cnodeIds->push_back(cnodeId);
830 if (grid == gcid) bIsDiagonalEntry =
true;
834 if(realnnz == 1 && bIsDiagonalEntry ==
true) {
835 LO lNodeId = nodeMap->getLocalElement(nodeId);
836 numberDirichletRowsPerNode[lNodeId] += 1;
837 if (numberDirichletRowsPerNode[lNodeId] == blockdim)
838 amalgBoundaryNodes[lNodeId] =
true;
841 Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp( cnodeIds );
843 if(arr_cnodeIds.size() > 0 )
844 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
847 crsGraph->fillComplete(nodeMap,nodeMap);
850 RCP<GraphBase> graph = rcp(
new Graph(crsGraph,
"amalgamated graph of A"));
853 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
856 GO numLocalBoundaryNodes = 0;
857 GO numGlobalBoundaryNodes = 0;
858 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
859 if (amalgBoundaryNodes[i])
860 numLocalBoundaryNodes++;
861 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
862 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
868 Set(currentLevel,
"DofsPerNode", blockdim);
869 Set(currentLevel,
"Graph", graph);
875 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
877 typedef typename ArrayView<const LO>::size_type size_type;
880 LO blkSize = A.GetFixedBlockSize();
881 if (A.IsView(
"stridedMaps") ==
true) {
882 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
883 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
885 if (strMap->getStridedBlockId() > -1)
886 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
890 size_t nnz = 0, pos = 0;
891 for (LO j = 0; j < blkSize; j++)
892 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
902 ArrayView<const LO> inds;
903 ArrayView<const SC> vals;
904 for (LO j = 0; j < blkSize; j++) {
905 A.getLocalRowView(row*blkSize+j, inds, vals);
906 size_type numIndices = inds.size();
912 cols[pos++] = translation[inds[0]];
913 for (size_type k = 1; k < numIndices; k++) {
914 LO nodeID = translation[inds[k]];
918 if (nodeID != cols[pos-1])
919 cols[pos++] = nodeID;
926 std::sort(cols.begin(), cols.end());
928 for (
size_t j = 1; j < nnz; j++)
929 if (cols[j] != cols[pos])
930 cols[++pos] = cols[j];
934 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
936 typedef typename ArrayView<const LO>::size_type size_type;
937 typedef Teuchos::ScalarTraits<SC> STS;
940 LO blkSize = A.GetFixedBlockSize();
941 if (A.IsView(
"stridedMaps") ==
true) {
942 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
943 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
945 if (strMap->getStridedBlockId() > -1)
946 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
950 size_t nnz = 0, pos = 0;
951 for (LO j = 0; j < blkSize; j++)
952 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
962 ArrayView<const LO> inds;
963 ArrayView<const SC> vals;
964 for (LO j = 0; j < blkSize; j++) {
965 A.getLocalRowView(row*blkSize+j, inds, vals);
966 size_type numIndices = inds.size();
973 for (size_type k = 0; k < numIndices; k++) {
975 LO nodeID = translation[inds[k]];
978 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold*ghostedDiagVals[dofID]*ghostedDiagVals[row*blkSize+j]);
979 typename STS::magnitudeType aij = STS::magnitude(vals[k]*vals[k]);
982 if (aij > aiiajj || (row*blkSize+j == dofID)) {
988 if (nodeID != prevNodeID) {
989 cols[pos++] = nodeID;
999 std::sort(cols.begin(), cols.end());
1001 for (
size_t j = 1; j < nnz; j++)
1002 if (cols[j] != cols[pos])
1003 cols[++pos] = cols[j];
1010 #endif // MUELU_COALESCEDROPFACTORY_DEF_HPP Important warning messages (one line)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
RCP< PreDropFunctionBaseClass > predrop_
MueLu representation of a compressed row storage graph.
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &v, LocalOrdinal i0, LocalOrdinal i1)
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.
virtual const Teuchos::ParameterList & GetParameterList() const
Exception throws to report incompatible objects (like maps).
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
translate global (row/column) id to global amalgamation block id
VerbLevel GetVerbLevel() const
Get the verbosity level.
Print statistics that do not involve significant additional computation.
CoalesceDropFactory()
Constructor.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
#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 Set(Level &level, const std::string &varName, const T &data) const
void Build(Level ¤tLevel) const
Build an object with this factory.
void DeclareInput(Level ¤tLevel) const
Input.
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
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.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.))
Exception throws to report errors in the internal logical of the program.
void Input(Level &level, const std::string &varName) const
void MergeRowsWithDropping(const Matrix &A, const LO row, const ArrayRCP< const SC > &ghostedDiagVals, SC threshold, Array< LO > &cols, const Array< LO > &translation) const