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" 87 template<
class real_type,
class LO>
97 DropTol(real_type val_, real_type diag_, LO col_,
bool drop_)
101 real_type
val {Teuchos::ScalarTraits<real_type>::zero()};
102 real_type
diag {Teuchos::ScalarTraits<real_type>::zero()};
103 LO
col {Teuchos::OrdinalTraits<LO>::invalid()};
109 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
111 RCP<ParameterList> validParamList = rcp(
new ParameterList());
113 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 118 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
119 validParamList->getEntry(
"aggregation: drop scheme").setValidator(
120 rcp(
new validatorType(Teuchos::tuple<std::string>(
"classical",
"distance laplacian"),
"aggregation: drop scheme")));
124 #undef SET_VALID_ENTRY 125 validParamList->set<
bool > (
"lightweight wrap",
true,
"Experimental option for lightweight graph access");
127 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
128 validParamList->set< RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory for UnAmalgamationInfo");
129 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory for Coordinates");
131 return validParamList;
134 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
137 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
139 Input(currentLevel,
"A");
140 Input(currentLevel,
"UnAmalgamationInfo");
142 const ParameterList& pL = GetParameterList();
143 if (pL.get<
bool>(
"lightweight wrap") ==
true) {
144 if (pL.get<std::string>(
"aggregation: drop scheme") ==
"distance laplacian")
145 Input(currentLevel,
"Coordinates");
150 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
155 typedef Teuchos::ScalarTraits<SC> STS;
156 typedef typename STS::magnitudeType real_type;
157 typedef Xpetra::MultiVector<real_type,LO,GO,NO> RealValuedMultiVector;
159 if (predrop_ != Teuchos::null)
160 GetOStream(
Parameters0) << predrop_->description();
162 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
163 RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(currentLevel,
"UnAmalgamationInfo");
165 const ParameterList & pL = GetParameterList();
166 bool doExperimentalWrap = pL.get<
bool>(
"lightweight wrap");
168 GetOStream(
Parameters0) <<
"lightweight wrap = " << doExperimentalWrap << std::endl;
181 if (doExperimentalWrap) {
182 std::string algo = pL.get<std::string>(
"aggregation: drop scheme");
184 TEUCHOS_TEST_FOR_EXCEPTION(predrop_ != null && algo !=
"classical",
Exceptions::RuntimeError,
"Dropping function must not be provided for \"" << algo <<
"\" algorithm");
185 TEUCHOS_TEST_FOR_EXCEPTION(algo !=
"classical" && algo !=
"distance laplacian",
Exceptions::RuntimeError,
"\"algorithm\" must be one of (classical|distance laplacian)");
187 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
188 std::string distanceLaplacianAlgoStr = pL.get<std::string>(
"aggregation: distance laplacian algo");
189 std::string classicalAlgoStr = pL.get<std::string>(
"aggregation: classical algo");
190 real_type realThreshold = STS::magnitude(threshold);
193 #ifdef HAVE_MUELU_DEBUG 194 int distanceLaplacianCutVerbose = 0;
196 #ifdef DJS_READ_ENV_VARIABLES 197 if (getenv(
"MUELU_DROP_TOLERANCE_MODE")) {
198 distanceLaplacianAlgoStr = std::string(getenv(
"MUELU_DROP_TOLERANCE_MODE"));
201 if (getenv(
"MUELU_DROP_TOLERANCE_THRESHOLD")) {
202 auto tmp = atoi(getenv(
"MUELU_DROP_TOLERANCE_THRESHOLD"));
203 realThreshold = 1e-4*tmp;
206 # ifdef HAVE_MUELU_DEBUG 207 if (getenv(
"MUELU_DROP_TOLERANCE_VERBOSE")) {
208 distanceLaplacianCutVerbose = atoi(getenv(
"MUELU_DROP_TOLERANCE_VERBOSE"));
214 enum decisionAlgoType {defaultAlgo, unscaled_cut, scaled_cut};
216 decisionAlgoType distanceLaplacianAlgo = defaultAlgo;
217 decisionAlgoType classicalAlgo = defaultAlgo;
218 if (algo ==
"distance laplacian") {
219 if (distanceLaplacianAlgoStr ==
"default")
220 distanceLaplacianAlgo = defaultAlgo;
221 else if (distanceLaplacianAlgoStr ==
"unscaled cut")
222 distanceLaplacianAlgo = unscaled_cut;
223 else if (distanceLaplacianAlgoStr ==
"scaled cut")
224 distanceLaplacianAlgo = scaled_cut;
226 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut), not \"" << distanceLaplacianAlgoStr <<
"\"");
227 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\" distance laplacian algorithm = \"" << distanceLaplacianAlgoStr <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
228 }
else if (algo ==
"classical") {
229 if (classicalAlgoStr ==
"default")
230 classicalAlgo = defaultAlgo;
231 else if (classicalAlgoStr ==
"unscaled cut")
232 classicalAlgo = unscaled_cut;
233 else if (classicalAlgoStr ==
"scaled cut")
234 classicalAlgo = scaled_cut;
236 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"\"aggregation: classical algo\" must be one of (default|unscaled cut|scaled cut), not \"" << classicalAlgoStr <<
"\"");
237 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\" classical algorithm = \"" << classicalAlgoStr <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
240 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
241 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
243 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<
double>(
"aggregation: Dirichlet threshold")));
245 GO numDropped = 0, numTotal = 0;
246 std::string graphType =
"unamalgamated";
247 if (algo ==
"classical") {
248 if (predrop_ == null) {
253 if (predrop_ != null) {
256 "MueLu::CoalesceFactory::Build: cast to PreDropFunctionConstVal failed.");
258 SC newt = predropConstVal->GetThreshold();
259 if (newt != threshold) {
260 GetOStream(
Warnings0) <<
"switching threshold parameter from " << threshold <<
" (list) to " << newt <<
" (user function" << std::endl;
267 if (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && A->hasCrsGraph()) {
269 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
271 ArrayRCP<const bool > boundaryNodes;
273 graph->SetBoundaryNodeMap(boundaryNodes);
274 numTotal = A->getNodeNumEntries();
277 GO numLocalBoundaryNodes = 0;
278 GO numGlobalBoundaryNodes = 0;
279 for (LO i = 0; i < boundaryNodes.size(); ++i)
280 if (boundaryNodes[i])
281 numLocalBoundaryNodes++;
282 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
283 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
284 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
287 Set(currentLevel,
"DofsPerNode", 1);
288 Set(currentLevel,
"Graph", graph);
290 }
else if ( (A->GetFixedBlockSize() == 1 && threshold != STS::zero()) ||
291 (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && !A->hasCrsGraph())) {
297 ArrayRCP<LO> rows (A->getNodeNumRows()+1);
298 ArrayRCP<LO> columns(A->getNodeNumEntries());
301 const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
306 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getNodeNumElements()); ++row) {
307 size_t nnz = A->getNumEntriesInLocalRow(row);
308 ArrayView<const LO> indices;
309 ArrayView<const SC> vals;
310 A->getLocalRowView(row, indices, vals);
312 if(classicalAlgo == defaultAlgo) {
319 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
320 LO col = indices[colID];
323 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
324 typename STS::magnitudeType aij = STS::magnitude(vals[colID]*vals[colID]);
326 if (aij > aiiajj || row == col) {
327 columns[realnnz++] = col;
332 rows[row+1] = realnnz;
338 std::vector<DropTol> drop_vec;
339 drop_vec.reserve(nnz);
340 const real_type zero = Teuchos::ScalarTraits<real_type>::zero();
341 const real_type one = Teuchos::ScalarTraits<real_type>::one();
345 for (LO colID = 0; colID < (LO)nnz; colID++) {
346 LO col = indices[colID];
348 drop_vec.emplace_back( zero, one, colID,
false);
353 if(boundaryNodes[colID])
continue;
354 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
355 typename STS::magnitudeType aij = STS::magnitude(vals[colID]*vals[colID]);
356 drop_vec.emplace_back(aij, aiiajj, colID,
false);
359 const size_t n = drop_vec.size();
361 if (classicalAlgo == unscaled_cut) {
362 std::sort( drop_vec.begin(), drop_vec.end()
363 , [](DropTol
const& a, DropTol
const& b) {
364 return a.val > b.val;
368 for (
size_t i=1; i<n; ++i) {
370 auto const& x = drop_vec[i-1];
371 auto const& y = drop_vec[i];
374 if (a > realThreshold*b) {
376 #ifdef HAVE_MUELU_DEBUG 377 if (distanceLaplacianCutVerbose) {
378 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
383 drop_vec[i].drop = drop;
385 }
else if (classicalAlgo == scaled_cut) {
386 std::sort( drop_vec.begin(), drop_vec.end()
387 , [](DropTol
const& a, DropTol
const& b) {
388 return a.val/a.diag > b.val/b.diag;
393 for (
size_t i=1; i<n; ++i) {
395 auto const& x = drop_vec[i-1];
396 auto const& y = drop_vec[i];
397 auto a = x.val/x.diag;
398 auto b = y.val/y.diag;
399 if (a > realThreshold*b) {
402 #ifdef HAVE_MUELU_DEBUG 403 if (distanceLaplacianCutVerbose) {
404 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
411 drop_vec[i].drop = drop;
415 std::sort( drop_vec.begin(), drop_vec.end()
416 , [](DropTol
const& a, DropTol
const& b) {
417 return a.col < b.col;
421 for (LO idxID =0; idxID<(LO)drop_vec.size(); idxID++) {
422 LO col = indices[drop_vec[idxID].col];
425 columns[realnnz++] = col;
430 if (!drop_vec[idxID].drop) {
431 columns[realnnz++] = col;
438 rows[row+1] = realnnz;
443 columns.resize(realnnz);
444 numTotal = A->getNodeNumEntries();
445 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, A->getRowMap(), A->getColMap(),
"thresholded graph of A"));
446 graph->SetBoundaryNodeMap(boundaryNodes);
448 GO numLocalBoundaryNodes = 0;
449 GO numGlobalBoundaryNodes = 0;
450 for (LO i = 0; i < boundaryNodes.size(); ++i)
451 if (boundaryNodes[i])
452 numLocalBoundaryNodes++;
453 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
454 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
455 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
457 Set(currentLevel,
"Graph", graph);
458 Set(currentLevel,
"DofsPerNode", 1);
460 }
else if (A->GetFixedBlockSize() > 1 && threshold == STS::zero()) {
462 const RCP<const Map> rowMap = A->getRowMap();
463 const RCP<const Map> colMap = A->getColMap();
465 graphType =
"amalgamated";
471 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
472 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
473 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
474 Array<LO> colTranslation = *(amalInfo->getColTranslation());
477 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
480 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
481 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
483 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
489 ArrayRCP<const bool > pointBoundaryNodes;
493 LO blkSize = A->GetFixedBlockSize();
495 LO blkPartSize = A->GetFixedBlockSize();
496 if (A->IsView(
"stridedMaps") ==
true) {
497 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
498 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
500 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
501 blkId = strMap->getStridedBlockId();
503 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
509 Array<LO> indicesExtra;
510 for (LO row = 0; row < numRows; row++) {
511 ArrayView<const LO> indices;
512 indicesExtra.resize(0);
520 bool isBoundary =
false;
522 for (LO j = 0; j < blkPartSize; j++) {
523 if (!pointBoundaryNodes[row*blkPartSize+j]) {
532 MergeRows(*A, row, indicesExtra, colTranslation);
534 indicesExtra.push_back(row);
535 indices = indicesExtra;
536 numTotal += indices.size();
540 LO nnz = indices.size(), rownnz = 0;
541 for (LO colID = 0; colID < nnz; colID++) {
542 LO col = indices[colID];
543 columns[realnnz++] = col;
554 amalgBoundaryNodes[row] =
true;
556 rows[row+1] = realnnz;
558 columns.resize(realnnz);
560 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
561 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
564 GO numLocalBoundaryNodes = 0;
565 GO numGlobalBoundaryNodes = 0;
567 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
568 if (amalgBoundaryNodes[i])
569 numLocalBoundaryNodes++;
571 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
572 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
573 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
574 <<
" agglomerated Dirichlet nodes" << std::endl;
577 Set(currentLevel,
"Graph", graph);
578 Set(currentLevel,
"DofsPerNode", blkSize);
580 }
else if (A->GetFixedBlockSize() > 1 && threshold != STS::zero()) {
582 const RCP<const Map> rowMap = A->getRowMap();
583 const RCP<const Map> colMap = A->getColMap();
585 graphType =
"amalgamated";
591 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
592 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
593 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
594 Array<LO> colTranslation = *(amalInfo->getColTranslation());
597 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
600 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
601 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
603 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
609 ArrayRCP<const bool > pointBoundaryNodes;
613 LO blkSize = A->GetFixedBlockSize();
615 LO blkPartSize = A->GetFixedBlockSize();
616 if (A->IsView(
"stridedMaps") ==
true) {
617 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
618 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
620 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
621 blkId = strMap->getStridedBlockId();
623 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
628 const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
633 Array<LO> indicesExtra;
634 for (LO row = 0; row < numRows; row++) {
635 ArrayView<const LO> indices;
636 indicesExtra.resize(0);
644 bool isBoundary =
false;
646 for (LO j = 0; j < blkPartSize; j++) {
647 if (!pointBoundaryNodes[row*blkPartSize+j]) {
656 MergeRowsWithDropping(*A, row, ghostedDiagVals, threshold, indicesExtra, colTranslation);
658 indicesExtra.push_back(row);
659 indices = indicesExtra;
660 numTotal += indices.size();
664 LO nnz = indices.size(), rownnz = 0;
665 for (LO colID = 0; colID < nnz; colID++) {
666 LO col = indices[colID];
667 columns[realnnz++] = col;
678 amalgBoundaryNodes[row] =
true;
680 rows[row+1] = realnnz;
682 columns.resize(realnnz);
684 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
685 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
688 GO numLocalBoundaryNodes = 0;
689 GO numGlobalBoundaryNodes = 0;
691 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
692 if (amalgBoundaryNodes[i])
693 numLocalBoundaryNodes++;
695 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
696 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
697 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
698 <<
" agglomerated Dirichlet nodes" << std::endl;
701 Set(currentLevel,
"Graph", graph);
702 Set(currentLevel,
"DofsPerNode", blkSize);
705 }
else if (algo ==
"distance laplacian") {
706 LO blkSize = A->GetFixedBlockSize();
707 GO indexBase = A->getRowMap()->getIndexBase();
713 RCP<RealValuedMultiVector> Coords = Get< RCP<RealValuedMultiVector > >(currentLevel,
"Coordinates");
719 ArrayRCP<const bool > pointBoundaryNodes;
722 if ( (blkSize == 1) && (threshold == STS::zero()) ) {
724 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
725 graph->SetBoundaryNodeMap(pointBoundaryNodes);
726 graphType=
"unamalgamated";
727 numTotal = A->getNodeNumEntries();
730 GO numLocalBoundaryNodes = 0;
731 GO numGlobalBoundaryNodes = 0;
732 for (LO i = 0; i < pointBoundaryNodes.size(); ++i)
733 if (pointBoundaryNodes[i])
734 numLocalBoundaryNodes++;
735 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
736 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
737 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
740 Set(currentLevel,
"DofsPerNode", blkSize);
741 Set(currentLevel,
"Graph", graph);
756 TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getNodeNumElements()/blkSize != Coords->getLocalLength(),
Exceptions::Incompatible,
757 "Coordinate vector length (" << Coords->getLocalLength() <<
") is incompatible with number of rows in A (" << A->getRowMap()->getNodeNumElements() <<
") by modulo block size ("<< blkSize <<
").");
759 const RCP<const Map> colMap = A->getColMap();
760 RCP<const Map> uniqueMap, nonUniqueMap;
761 Array<LO> colTranslation;
763 uniqueMap = A->getRowMap();
764 nonUniqueMap = A->getColMap();
765 graphType=
"unamalgamated";
768 uniqueMap = Coords->getMap();
770 "Different index bases for matrix and coordinates");
774 graphType =
"amalgamated";
776 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
778 RCP<RealValuedMultiVector> ghostedCoords;
779 RCP<Vector> ghostedLaplDiag;
780 Teuchos::ArrayRCP<SC> ghostedLaplDiagData;
781 if (threshold != STS::zero()) {
783 RCP<const Import> importer;
786 if (blkSize == 1 && A->getCrsGraph()->getImporter() != Teuchos::null) {
787 GetOStream(
Warnings1) <<
"Using existing importer from matrix graph" << std::endl;
788 importer = A->getCrsGraph()->getImporter();
790 GetOStream(
Warnings0) <<
"Constructing new importer instance" << std::endl;
791 importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
794 ghostedCoords = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(nonUniqueMap, Coords->getNumVectors());
797 ghostedCoords->doImport(*Coords, *importer, Xpetra::INSERT);
801 RCP<Vector> localLaplDiag = VectorFactory::Build(uniqueMap);
802 ArrayRCP<SC> localLaplDiagData = localLaplDiag->getDataNonConst(0);
803 Array<LO> indicesExtra;
804 Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
805 if (threshold != STS::zero()) {
806 const size_t numVectors = ghostedCoords->getNumVectors();
807 coordData.reserve(numVectors);
808 for (
size_t j = 0; j < numVectors; j++) {
809 Teuchos::ArrayRCP<const real_type> tmpData=ghostedCoords->getData(j);
810 coordData.push_back(tmpData);
815 for (LO row = 0; row < numRows; row++) {
816 ArrayView<const LO> indices;
819 ArrayView<const SC> vals;
820 A->getLocalRowView(row, indices, vals);
824 indicesExtra.resize(0);
825 MergeRows(*A, row, indicesExtra, colTranslation);
826 indices = indicesExtra;
829 LO nnz = indices.size();
830 for (LO colID = 0; colID < nnz; colID++) {
831 const LO col = indices[colID];
841 ghostedLaplDiag = VectorFactory::Build(nonUniqueMap);
842 ghostedLaplDiag->doImport(*localLaplDiag, *importer, Xpetra::INSERT);
843 ghostedLaplDiagData = ghostedLaplDiag->getDataNonConst(0);
847 GetOStream(
Runtime0) <<
"Skipping distance laplacian construction due to 0 threshold" << std::endl;
853 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
854 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
856 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
860 Array<LO> indicesExtra;
863 Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
864 if (threshold != STS::zero()) {
865 const size_t numVectors = ghostedCoords->getNumVectors();
866 coordData.reserve(numVectors);
867 for (
size_t j = 0; j < numVectors; j++) {
868 Teuchos::ArrayRCP<const real_type> tmpData=ghostedCoords->getData(j);
869 coordData.push_back(tmpData);
873 for (LO row = 0; row < numRows; row++) {
874 ArrayView<const LO> indices;
875 indicesExtra.resize(0);
876 bool isBoundary =
false;
879 ArrayView<const SC> vals;
880 A->getLocalRowView(row, indices, vals);
881 isBoundary = pointBoundaryNodes[row];
884 for (LO j = 0; j < blkSize; j++) {
885 if (!pointBoundaryNodes[row*blkSize+j]) {
893 MergeRows(*A, row, indicesExtra, colTranslation);
895 indicesExtra.push_back(row);
896 indices = indicesExtra;
898 numTotal += indices.size();
900 LO nnz = indices.size(), rownnz = 0;
902 if (threshold != STS::zero()) {
904 if (distanceLaplacianAlgo == defaultAlgo) {
906 for (LO colID = 0; colID < nnz; colID++) {
908 LO col = indices[colID];
911 columns[realnnz++] = col;
917 if(isBoundary)
continue;
921 real_type aiiajj = STS::magnitude(realThreshold*realThreshold * ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
922 real_type aij = STS::magnitude(laplVal*laplVal);
925 columns[realnnz++] = col;
934 std::vector<DropTol> drop_vec;
935 drop_vec.reserve(nnz);
936 const real_type zero = Teuchos::ScalarTraits<real_type>::zero();
937 const real_type one = Teuchos::ScalarTraits<real_type>::one();
940 for (LO colID = 0; colID < nnz; colID++) {
942 LO col = indices[colID];
945 drop_vec.emplace_back( zero, one, colID,
false);
949 if(isBoundary)
continue;
952 real_type aiiajj = STS::magnitude(ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
953 real_type aij = STS::magnitude(laplVal*laplVal);
955 drop_vec.emplace_back(aij, aiiajj, colID,
false);
958 const size_t n = drop_vec.size();
960 if (distanceLaplacianAlgo == unscaled_cut) {
962 std::sort( drop_vec.begin(), drop_vec.end()
963 , [](DropTol
const& a, DropTol
const& b) {
964 return a.val > b.val;
969 for (
size_t i=1; i<n; ++i) {
971 auto const& x = drop_vec[i-1];
972 auto const& y = drop_vec[i];
975 if (a > realThreshold*b) {
977 #ifdef HAVE_MUELU_DEBUG 978 if (distanceLaplacianCutVerbose) {
979 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
984 drop_vec[i].drop = drop;
987 else if (distanceLaplacianAlgo == scaled_cut) {
989 std::sort( drop_vec.begin(), drop_vec.end()
990 , [](DropTol
const& a, DropTol
const& b) {
991 return a.val/a.diag > b.val/b.diag;
996 for (
size_t i=1; i<n; ++i) {
998 auto const& x = drop_vec[i-1];
999 auto const& y = drop_vec[i];
1000 auto a = x.val/x.diag;
1001 auto b = y.val/y.diag;
1002 if (a > realThreshold*b) {
1004 #ifdef HAVE_MUELU_DEBUG 1005 if (distanceLaplacianCutVerbose) {
1006 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
1011 drop_vec[i].drop = drop;
1015 std::sort( drop_vec.begin(), drop_vec.end()
1016 , [](DropTol
const& a, DropTol
const& b) {
1017 return a.col < b.col;
1021 for (LO idxID =0; idxID<(LO)drop_vec.size(); idxID++) {
1022 LO col = indices[drop_vec[idxID].col];
1027 columns[realnnz++] = col;
1032 if (!drop_vec[idxID].drop) {
1033 columns[realnnz++] = col;
1042 for (LO colID = 0; colID < nnz; colID++) {
1043 LO col = indices[colID];
1044 columns[realnnz++] = col;
1056 amalgBoundaryNodes[row] =
true;
1058 rows[row+1] = realnnz;
1062 columns.resize(realnnz);
1064 RCP<GraphBase> graph;
1067 graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
1068 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
1072 GO numLocalBoundaryNodes = 0;
1073 GO numGlobalBoundaryNodes = 0;
1075 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
1076 if (amalgBoundaryNodes[i])
1077 numLocalBoundaryNodes++;
1079 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1080 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1081 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" agglomerated Dirichlet nodes" 1082 <<
" using threshold " << dirichletThreshold << std::endl;
1085 Set(currentLevel,
"Graph", graph);
1086 Set(currentLevel,
"DofsPerNode", blkSize);
1090 if ((GetVerbLevel() &
Statistics1) && !(A->GetFixedBlockSize() > 1 && threshold != STS::zero())) {
1091 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1092 GO numGlobalTotal, numGlobalDropped;
1095 GetOStream(
Statistics1) <<
"Number of dropped entries in " << graphType <<
" matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1096 if (numGlobalTotal != 0)
1097 GetOStream(
Statistics1) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
1104 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
1106 GetOStream(
Runtime0) <<
"algorithm = \"" <<
"failsafe" <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
1107 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
1109 RCP<const Map> rowMap = A->getRowMap();
1110 RCP<const Map> colMap = A->getColMap();
1113 GO indexBase = rowMap->getIndexBase();
1117 if(A->IsView(
"stridedMaps") &&
1118 Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1119 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
1120 RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap());
1121 TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
1122 blockdim = strMap->getFixedBlockSize();
1123 offset = strMap->getOffset();
1124 oldView = A->SwitchToView(oldView);
1125 GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build():" <<
" found blockdim=" << blockdim <<
" from strided maps. offset=" << offset << std::endl;
1126 }
else GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
1130 RCP<const Map> nodeMap = amalInfo->getNodeRowMap();
1131 GetOStream(
Statistics1) <<
"CoalesceDropFactory: nodeMap " << nodeMap->getNodeNumElements() <<
"/" << nodeMap->getGlobalNumElements() <<
" elements" << std::endl;
1134 RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, A->getNodeMaxNumRowEntries()*blockdim);
1136 LO numRows = A->getRowMap()->getNodeNumElements();
1137 LO numNodes = nodeMap->getNodeNumElements();
1138 const ArrayRCP<bool> amalgBoundaryNodes(numNodes,
false);
1139 const ArrayRCP<int> numberDirichletRowsPerNode(numNodes, 0);
1140 bool bIsDiagonalEntry =
false;
1145 for(LO row=0; row<numRows; row++) {
1147 GO grid = rowMap->getGlobalElement(row);
1150 bIsDiagonalEntry =
false;
1155 size_t nnz = A->getNumEntriesInLocalRow(row);
1156 Teuchos::ArrayView<const LO> indices;
1157 Teuchos::ArrayView<const SC> vals;
1158 A->getLocalRowView(row, indices, vals);
1160 RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(
new std::vector<GO>);
1162 for(LO col=0; col<Teuchos::as<LO>(nnz); col++) {
1163 GO gcid = colMap->getGlobalElement(indices[col]);
1165 if(vals[col]!=STS::zero()) {
1167 cnodeIds->push_back(cnodeId);
1169 if (grid == gcid) bIsDiagonalEntry =
true;
1173 if(realnnz == 1 && bIsDiagonalEntry ==
true) {
1174 LO lNodeId = nodeMap->getLocalElement(nodeId);
1175 numberDirichletRowsPerNode[lNodeId] += 1;
1176 if (numberDirichletRowsPerNode[lNodeId] == blockdim)
1177 amalgBoundaryNodes[lNodeId] =
true;
1180 Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp( cnodeIds );
1182 if(arr_cnodeIds.size() > 0 )
1183 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
1186 crsGraph->fillComplete(nodeMap,nodeMap);
1189 RCP<GraphBase> graph = rcp(
new Graph(crsGraph,
"amalgamated graph of A"));
1192 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
1195 GO numLocalBoundaryNodes = 0;
1196 GO numGlobalBoundaryNodes = 0;
1197 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
1198 if (amalgBoundaryNodes[i])
1199 numLocalBoundaryNodes++;
1200 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1201 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1202 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
1207 Set(currentLevel,
"DofsPerNode", blockdim);
1208 Set(currentLevel,
"Graph", graph);
1214 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1216 typedef typename ArrayView<const LO>::size_type size_type;
1219 LO blkSize = A.GetFixedBlockSize();
1220 if (A.IsView(
"stridedMaps") ==
true) {
1221 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
1222 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
1224 if (strMap->getStridedBlockId() > -1)
1225 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1229 size_t nnz = 0, pos = 0;
1230 for (LO j = 0; j < blkSize; j++)
1231 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
1241 ArrayView<const LO> inds;
1242 ArrayView<const SC> vals;
1243 for (LO j = 0; j < blkSize; j++) {
1244 A.getLocalRowView(row*blkSize+j, inds, vals);
1245 size_type numIndices = inds.size();
1247 if (numIndices == 0)
1251 cols[pos++] = translation[inds[0]];
1252 for (size_type k = 1; k < numIndices; k++) {
1253 LO nodeID = translation[inds[k]];
1257 if (nodeID != cols[pos-1])
1258 cols[pos++] = nodeID;
1265 std::sort(cols.begin(), cols.end());
1267 for (
size_t j = 1; j < nnz; j++)
1268 if (cols[j] != cols[pos])
1269 cols[++pos] = cols[j];
1273 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1275 typedef typename ArrayView<const LO>::size_type size_type;
1276 typedef Teuchos::ScalarTraits<SC> STS;
1279 LO blkSize = A.GetFixedBlockSize();
1280 if (A.IsView(
"stridedMaps") ==
true) {
1281 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
1282 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
1284 if (strMap->getStridedBlockId() > -1)
1285 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1289 size_t nnz = 0, pos = 0;
1290 for (LO j = 0; j < blkSize; j++)
1291 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
1301 ArrayView<const LO> inds;
1302 ArrayView<const SC> vals;
1303 for (LO j = 0; j < blkSize; j++) {
1304 A.getLocalRowView(row*blkSize+j, inds, vals);
1305 size_type numIndices = inds.size();
1307 if (numIndices == 0)
1312 for (size_type k = 0; k < numIndices; k++) {
1314 LO nodeID = translation[inds[k]];
1317 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold*ghostedDiagVals[dofID]*ghostedDiagVals[row*blkSize+j]);
1318 typename STS::magnitudeType aij = STS::magnitude(vals[k]*vals[k]);
1321 if (aij > aiiajj || (row*blkSize+j == dofID)) {
1327 if (nodeID != prevNodeID) {
1328 cols[pos++] = nodeID;
1329 prevNodeID = nodeID;
1338 std::sort(cols.begin(), cols.end());
1340 for (
size_t j = 1; j < nnz; j++)
1341 if (cols[j] != cols[pos])
1342 cols[++pos] = cols[j];
1349 #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.
DropTol & operator=(DropTol const &)=default
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.
DropTol(real_type val_, real_type diag_, LO col_, bool drop_)
Exception throws to report errors in the internal logical of the program.