46 #ifndef MUELU_COALESCEDROPFACTORY_DEF_HPP 47 #define MUELU_COALESCEDROPFACTORY_DEF_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>
123 typedef Teuchos::ScalarTraits<SC> STS;
125 if (predrop_ != Teuchos::null)
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;
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);
201 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
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);
272 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
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);
391 GetOStream(
Statistics1) <<
"Detected " << 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]) {
475 MergeRowsWithDropping(*A, row, ghostedDiagVals, threshold, indicesExtra, colTranslation);
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);
516 GetOStream(
Statistics1) <<
"Detected " << 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);
556 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
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);
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);
741 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" agglomerated Dirichlet nodes" 742 <<
" using threshold " << dirichletThreshold << std::endl;
745 Set(currentLevel,
"Graph", graph);
746 Set(currentLevel,
"DofsPerNode", blkSize);
750 if ((GetVerbLevel() &
Statistics1) && !(A->GetFixedBlockSize() > 1 && threshold != STS::zero())) {
751 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
752 GO numGlobalTotal, numGlobalDropped;
755 GetOStream(
Statistics1) <<
"Number of dropped entries in " << graphType <<
" matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
756 if (numGlobalTotal != 0)
757 GetOStream(
Statistics1) <<
" (" << 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) {
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(
Statistics1) <<
"CoalesceDropFactory::Build():" <<
" found blockdim=" << blockdim <<
" from strided maps. offset=" << offset << std::endl;
787 }
else GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
791 RCP<const Map> nodeMap = amalInfo->getNodeRowMap();
792 GetOStream(
Statistics1) <<
"CoalesceDropFactory: nodeMap " << nodeMap->getNodeNumElements() <<
"/" << nodeMap->getGlobalNumElements() <<
" elements" << std::endl;
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);
863 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
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)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
virtual std::string description() const=0
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.
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
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.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
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.