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_ExportFactory.hpp> 53 #include <Xpetra_MapFactory.hpp> 54 #include <Xpetra_Map.hpp> 55 #include <Xpetra_Matrix.hpp> 56 #include <Xpetra_MultiVectorFactory.hpp> 57 #include <Xpetra_MultiVector.hpp> 58 #include <Xpetra_StridedMap.hpp> 59 #include <Xpetra_VectorFactory.hpp> 60 #include <Xpetra_Vector.hpp> 62 #include <Xpetra_IO.hpp> 66 #include "MueLu_AmalgamationFactory.hpp" 67 #include "MueLu_AmalgamationInfo.hpp" 70 #include "MueLu_Graph.hpp" 72 #include "MueLu_LWGraph.hpp" 75 #include "MueLu_PreDropFunctionBaseClass.hpp" 76 #include "MueLu_PreDropFunctionConstVal.hpp" 77 #include "MueLu_Utilities.hpp" 79 #ifdef HAVE_XPETRA_TPETRA 80 #include "Tpetra_CrsGraphTransposer.hpp" 95 template<
class real_type,
class LO>
105 DropTol(real_type val_, real_type diag_, LO col_,
bool drop_)
109 real_type
val {Teuchos::ScalarTraits<real_type>::zero()};
110 real_type
diag {Teuchos::ScalarTraits<real_type>::zero()};
111 LO
col {Teuchos::OrdinalTraits<LO>::invalid()};
120 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
122 RCP<ParameterList> validParamList = rcp(
new ParameterList());
124 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 131 SET_VALID_ENTRY(
"aggregation: distance laplacian directional weights");
134 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
135 validParamList->getEntry(
"aggregation: drop scheme").setValidator(rcp(
new validatorType(Teuchos::tuple<std::string>(
"classical",
"distance laplacian",
"signed classical",
"block diagonal",
"block diagonal classical",
"block diagonal distance laplacian",
"block diagonal signed classical",
"block diagonal colored signed classical"),
"aggregation: drop scheme")));
141 #undef SET_VALID_ENTRY 142 validParamList->set<
bool > (
"lightweight wrap",
true,
"Experimental option for lightweight graph access");
144 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
145 validParamList->set< RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory for UnAmalgamationInfo");
146 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory for Coordinates");
147 validParamList->set< RCP<const FactoryBase> >(
"BlockNumber", Teuchos::null,
"Generating factory for BlockNUmber");
149 return validParamList;
152 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
155 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
157 Input(currentLevel,
"A");
158 Input(currentLevel,
"UnAmalgamationInfo");
160 const ParameterList& pL = GetParameterList();
161 if (pL.get<
bool>(
"lightweight wrap") ==
true) {
162 std::string algo = pL.get<std::string>(
"aggregation: drop scheme");
163 if (algo ==
"distance laplacian" || algo ==
"block diagonal distance laplacian") {
164 Input(currentLevel,
"Coordinates");
166 if (algo.find(
"block diagonal") != std::string::npos || algo.find(
"signed classical") != std::string::npos) {
167 Input(currentLevel,
"BlockNumber");
173 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
178 typedef Teuchos::ScalarTraits<SC> STS;
179 typedef typename STS::magnitudeType real_type;
180 typedef Xpetra::MultiVector<real_type,LO,GO,NO> RealValuedMultiVector;
181 typedef Xpetra::MultiVectorFactory<real_type,LO,GO,NO> RealValuedMultiVectorFactory;
183 if (predrop_ != Teuchos::null)
184 GetOStream(
Parameters0) << predrop_->description();
186 RCP<Matrix> realA = Get< RCP<Matrix> >(currentLevel,
"A");
187 RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(currentLevel,
"UnAmalgamationInfo");
188 const ParameterList & pL = GetParameterList();
189 bool doExperimentalWrap = pL.get<
bool>(
"lightweight wrap");
191 GetOStream(
Parameters0) <<
"lightweight wrap = " << doExperimentalWrap << std::endl;
192 std::string algo = pL.get<std::string>(
"aggregation: drop scheme");
194 RCP<RealValuedMultiVector> Coords;
197 bool use_block_algorithm=
false;
198 LO interleaved_blocksize = as<LO>(pL.get<
int>(
"aggregation: block diagonal: interleaved blocksize"));
199 bool useSignedClassical =
false;
200 bool generateColoringGraph =
false;
204 typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<
double>(
"aggregation: row sum drop tol"));
205 RCP<LocalOrdinalVector> ghostedBlockNumber;
206 ArrayRCP<const LO> g_block_id;
208 if(algo ==
"distance laplacian" ) {
210 Coords = Get< RCP<RealValuedMultiVector > >(currentLevel,
"Coordinates");
213 else if(algo ==
"signed classical" || algo ==
"block diagonal colored signed classical" || algo ==
"block diagonal signed classical") {
214 useSignedClassical =
true;
216 RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(currentLevel,
"BlockNumber");
218 RCP<const Import> importer = realA->getCrsGraph()->getImporter();
219 if(!importer.is_null()) {
221 ghostedBlockNumber= Xpetra::VectorFactory<LO,LO,GO,NO>::Build(importer->getTargetMap());
222 ghostedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
225 ghostedBlockNumber = BlockNumber;
227 g_block_id = ghostedBlockNumber->getData(0);
229 if(algo ==
"block diagonal colored signed classical")
230 generateColoringGraph=
true;
235 else if(algo ==
"block diagonal") {
237 BlockDiagonalize(currentLevel,realA,
false);
240 else if (algo ==
"block diagonal classical" || algo ==
"block diagonal distance laplacian") {
242 use_block_algorithm =
true;
243 RCP<Matrix> filteredMatrix = BlockDiagonalize(currentLevel,realA,
true);
244 if(algo ==
"block diagonal distance laplacian") {
246 RCP<RealValuedMultiVector> OldCoords = Get< RCP<RealValuedMultiVector > >(currentLevel,
"Coordinates");
247 if (OldCoords->getLocalLength() != realA->getNodeNumRows()) {
248 LO dim = (LO) OldCoords->getNumVectors();
249 Coords = RealValuedMultiVectorFactory::Build(realA->getRowMap(),dim);
250 for(LO k=0; k<dim; k++){
251 ArrayRCP<const real_type> old_vec = OldCoords->getData(k);
252 ArrayRCP<real_type> new_vec = Coords->getDataNonConst(k);
253 for(LO i=0; i <(LO)OldCoords->getLocalLength(); i++) {
255 for(LO j=0; j<interleaved_blocksize; j++)
256 new_vec[new_base + j] = old_vec[i];
263 algo =
"distance laplacian";
265 else if(algo ==
"block diagonal classical") {
277 Array<double> dlap_weights = pL.get<Array<double> >(
"aggregation: distance laplacian directional weights");
278 enum {NO_WEIGHTS=0, SINGLE_WEIGHTS, BLOCK_WEIGHTS};
279 int use_dlap_weights = NO_WEIGHTS;
280 if(algo ==
"distance laplacian") {
281 LO dim = (LO) Coords->getNumVectors();
283 bool non_unity =
false;
284 for (LO i=0; !non_unity && i<(LO)dlap_weights.size(); i++) {
285 if(dlap_weights[i] != 1.0) {
290 LO blocksize = use_block_algorithm ? as<LO>(pL.get<
int>(
"aggregation: block diagonal: interleaved blocksize")) : 1;
291 if((LO)dlap_weights.size() == dim)
292 use_dlap_weights = SINGLE_WEIGHTS;
293 else if((LO)dlap_weights.size() == blocksize * dim)
294 use_dlap_weights = BLOCK_WEIGHTS;
297 "length of 'aggregation: distance laplacian directional weights' must equal the coordinate dimension OR the coordinate dimension times the blocksize");
300 GetOStream(
Statistics1) <<
"Using distance laplacian weights: "<<dlap_weights<<std::endl;
315 if (doExperimentalWrap) {
316 TEUCHOS_TEST_FOR_EXCEPTION(predrop_ != null && algo !=
"classical",
Exceptions::RuntimeError,
"Dropping function must not be provided for \"" << algo <<
"\" algorithm");
317 TEUCHOS_TEST_FOR_EXCEPTION(algo !=
"classical" && algo !=
"distance laplacian" && algo !=
"signed classical",
Exceptions::RuntimeError,
"\"algorithm\" must be one of (classical|distance laplacian|signed classical)");
319 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
320 std::string distanceLaplacianAlgoStr = pL.get<std::string>(
"aggregation: distance laplacian algo");
321 std::string classicalAlgoStr = pL.get<std::string>(
"aggregation: classical algo");
322 real_type realThreshold = STS::magnitude(threshold);
326 #ifdef HAVE_MUELU_DEBUG 327 int distanceLaplacianCutVerbose = 0;
329 #ifdef DJS_READ_ENV_VARIABLES 330 if (getenv(
"MUELU_DROP_TOLERANCE_MODE")) {
331 distanceLaplacianAlgoStr = std::string(getenv(
"MUELU_DROP_TOLERANCE_MODE"));
334 if (getenv(
"MUELU_DROP_TOLERANCE_THRESHOLD")) {
335 auto tmp = atoi(getenv(
"MUELU_DROP_TOLERANCE_THRESHOLD"));
336 realThreshold = 1e-4*tmp;
339 # ifdef HAVE_MUELU_DEBUG 340 if (getenv(
"MUELU_DROP_TOLERANCE_VERBOSE")) {
341 distanceLaplacianCutVerbose = atoi(getenv(
"MUELU_DROP_TOLERANCE_VERBOSE"));
347 enum decisionAlgoType {defaultAlgo, unscaled_cut, scaled_cut, scaled_cut_symmetric};
349 decisionAlgoType distanceLaplacianAlgo = defaultAlgo;
350 decisionAlgoType classicalAlgo = defaultAlgo;
351 if (algo ==
"distance laplacian") {
352 if (distanceLaplacianAlgoStr ==
"default")
353 distanceLaplacianAlgo = defaultAlgo;
354 else if (distanceLaplacianAlgoStr ==
"unscaled cut")
355 distanceLaplacianAlgo = unscaled_cut;
356 else if (distanceLaplacianAlgoStr ==
"scaled cut")
357 distanceLaplacianAlgo = scaled_cut;
358 else if (distanceLaplacianAlgoStr ==
"scaled cut symmetric")
359 distanceLaplacianAlgo = scaled_cut_symmetric;
361 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut), not \"" << distanceLaplacianAlgoStr <<
"\"");
362 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\" distance laplacian algorithm = \"" << distanceLaplacianAlgoStr <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
363 }
else if (algo ==
"classical") {
364 if (classicalAlgoStr ==
"default")
365 classicalAlgo = defaultAlgo;
366 else if (classicalAlgoStr ==
"unscaled cut")
367 classicalAlgo = unscaled_cut;
368 else if (classicalAlgoStr ==
"scaled cut")
369 classicalAlgo = scaled_cut;
371 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"\"aggregation: classical algo\" must be one of (default|unscaled cut|scaled cut), not \"" << classicalAlgoStr <<
"\"");
372 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\" classical algorithm = \"" << classicalAlgoStr <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
375 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
376 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
378 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<
double>(
"aggregation: Dirichlet threshold")));
382 TEUCHOS_TEST_FOR_EXCEPTION(useSignedClassical && classicalAlgo != defaultAlgo,
Exceptions::RuntimeError,
"\"aggregation: classical algo\" != default is not supported for scalled classical aggregation");
385 GO numDropped = 0, numTotal = 0;
386 std::string graphType =
"unamalgamated";
387 if (algo ==
"classical") {
388 if (predrop_ == null) {
393 if (predrop_ != null) {
396 "MueLu::CoalesceFactory::Build: cast to PreDropFunctionConstVal failed.");
398 SC newt = predropConstVal->GetThreshold();
399 if (newt != threshold) {
400 GetOStream(
Warnings0) <<
"switching threshold parameter from " << threshold <<
" (list) to " << newt <<
" (user function" << std::endl;
407 if (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && !useSignedClassical && A->hasCrsGraph()) {
409 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
415 graph->SetBoundaryNodeMap(boundaryNodes);
416 numTotal = A->getNodeNumEntries();
419 GO numLocalBoundaryNodes = 0;
420 GO numGlobalBoundaryNodes = 0;
421 for (LO i = 0; i < boundaryNodes.size(); ++i)
422 if (boundaryNodes[i])
423 numLocalBoundaryNodes++;
424 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
425 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
426 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
429 Set(currentLevel,
"DofsPerNode", 1);
430 Set(currentLevel,
"Graph", graph);
432 }
else if ( (A->GetFixedBlockSize() == 1 && threshold != STS::zero()) ||
433 (A->GetFixedBlockSize() == 1 && threshold == STS::zero() && !A->hasCrsGraph()) ||
434 (A->GetFixedBlockSize() == 1 && useSignedClassical) ) {
440 ArrayRCP<LO> rows (A->getNodeNumRows()+1);
441 ArrayRCP<LO> columns(A->getNodeNumEntries());
443 using MT =
typename STS::magnitudeType;
444 RCP<Vector> ghostedDiag;
445 ArrayRCP<const SC> ghostedDiagVals;
446 ArrayRCP<const MT> negMaxOffDiagonal;
447 if(useSignedClassical) {
448 if(ghostedBlockNumber.is_null()) {
451 GetOStream(
Statistics1) <<
"Calculated max point off-diagonal" << std::endl;
456 GetOStream(
Statistics1) <<
"Calculating max block off-diagonal" << std::endl;
461 ghostedDiagVals = ghostedDiag->getData(0);
464 if (rowSumTol > 0.) {
465 if(ghostedBlockNumber.is_null()) {
467 GetOStream(
Statistics1) <<
"Applying point row sum criterion." << std::endl;
471 GetOStream(
Statistics1) <<
"Applying block row sum criterion." << std::endl;
478 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getNodeNumElements()); ++row) {
479 size_t nnz = A->getNumEntriesInLocalRow(row);
480 bool rowIsDirichlet = boundaryNodes[row];
481 ArrayView<const LO> indices;
482 ArrayView<const SC> vals;
483 A->getLocalRowView(row, indices, vals);
485 if(classicalAlgo == defaultAlgo) {
492 if(useSignedClassical) {
494 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
495 LO col = indices[colID];
496 MT max_neg_aik = realThreshold * STS::real(negMaxOffDiagonal[row]);
497 MT neg_aij = - STS::real(vals[colID]);
502 if ((!rowIsDirichlet && (g_block_id.is_null() || g_block_id[row] == g_block_id[col]) && neg_aij > max_neg_aik) || row == col) {
503 columns[realnnz++] = col;
508 rows[row+1] = realnnz;
512 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
513 LO col = indices[colID];
514 MT aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
515 MT aij = STS::magnitude(vals[colID]*vals[colID]);
517 if ((!rowIsDirichlet && aij > aiiajj) || row == col) {
518 columns[realnnz++] = col;
523 rows[row+1] = realnnz;
530 std::vector<DropTol> drop_vec;
531 drop_vec.reserve(nnz);
532 const real_type zero = Teuchos::ScalarTraits<real_type>::zero();
533 const real_type one = Teuchos::ScalarTraits<real_type>::one();
538 for (LO colID = 0; colID < (LO)nnz; colID++) {
539 LO col = indices[colID];
541 drop_vec.emplace_back( zero, one, colID,
false);
546 if(boundaryNodes[colID])
continue;
547 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
548 typename STS::magnitudeType aij = STS::magnitude(vals[colID]*vals[colID]);
549 drop_vec.emplace_back(aij, aiiajj, colID,
false);
552 const size_t n = drop_vec.size();
554 if (classicalAlgo == unscaled_cut) {
555 std::sort( drop_vec.begin(), drop_vec.end()
556 , [](DropTol
const& a, DropTol
const& b) {
557 return a.val > b.val;
561 for (
size_t i=1; i<n; ++i) {
563 auto const& x = drop_vec[i-1];
564 auto const& y = drop_vec[i];
567 if (a > realThreshold*b) {
569 #ifdef HAVE_MUELU_DEBUG 570 if (distanceLaplacianCutVerbose) {
571 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
576 drop_vec[i].drop = drop;
578 }
else if (classicalAlgo == scaled_cut) {
579 std::sort( drop_vec.begin(), drop_vec.end()
580 , [](DropTol
const& a, DropTol
const& b) {
581 return a.val/a.diag > b.val/b.diag;
586 for (
size_t i=1; i<n; ++i) {
588 auto const& x = drop_vec[i-1];
589 auto const& y = drop_vec[i];
590 auto a = x.val/x.diag;
591 auto b = y.val/y.diag;
592 if (a > realThreshold*b) {
595 #ifdef HAVE_MUELU_DEBUG 596 if (distanceLaplacianCutVerbose) {
597 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
604 drop_vec[i].drop = drop;
608 std::sort( drop_vec.begin(), drop_vec.end()
609 , [](DropTol
const& a, DropTol
const& b) {
610 return a.col < b.col;
614 for (LO idxID =0; idxID<(LO)drop_vec.size(); idxID++) {
615 LO col = indices[drop_vec[idxID].col];
618 columns[realnnz++] = col;
623 if (!drop_vec[idxID].drop) {
624 columns[realnnz++] = col;
631 rows[row+1] = realnnz;
636 columns.resize(realnnz);
637 numTotal = A->getNodeNumEntries();
638 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, A->getRowMap(), A->getColMap(),
"thresholded graph of A"));
639 graph->SetBoundaryNodeMap(boundaryNodes);
641 GO numLocalBoundaryNodes = 0;
642 GO numGlobalBoundaryNodes = 0;
643 for (LO i = 0; i < boundaryNodes.size(); ++i)
644 if (boundaryNodes[i])
645 numLocalBoundaryNodes++;
646 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
647 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
648 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
650 Set(currentLevel,
"Graph", graph);
651 Set(currentLevel,
"DofsPerNode", 1);
654 if(generateColoringGraph) {
655 RCP<GraphBase> colorGraph;
656 RCP<const Import> importer = A->getCrsGraph()->getImporter();
657 BlockDiagonalizeGraph(graph,ghostedBlockNumber,colorGraph,importer);
658 Set(currentLevel,
"Coloring Graph",colorGraph);
662 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"m_regular_graph."+std::to_string(currentLevel.GetLevelID()), *rcp_dynamic_cast<LWGraph>(graph)->GetCrsGraph());
663 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"m_color_graph."+std::to_string(currentLevel.GetLevelID()), *rcp_dynamic_cast<LWGraph>(colorGraph)->GetCrsGraph());
679 }
else if (A->GetFixedBlockSize() > 1 && threshold == STS::zero()) {
681 const RCP<const Map> rowMap = A->getRowMap();
682 const RCP<const Map> colMap = A->getColMap();
684 graphType =
"amalgamated";
690 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
691 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
692 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
693 Array<LO> colTranslation = *(amalInfo->getColTranslation());
696 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
699 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
700 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
702 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
708 ArrayRCP<bool > pointBoundaryNodes;
715 LO blkSize = A->GetFixedBlockSize();
717 LO blkPartSize = A->GetFixedBlockSize();
718 if (A->IsView(
"stridedMaps") ==
true) {
719 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
720 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
722 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
723 blkId = strMap->getStridedBlockId();
725 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
731 Array<LO> indicesExtra;
732 for (LO row = 0; row < numRows; row++) {
733 ArrayView<const LO> indices;
734 indicesExtra.resize(0);
742 bool isBoundary =
false;
743 if (pL.get<
bool>(
"aggregation: greedy Dirichlet") ==
true) {
744 for (LO j = 0; j < blkPartSize; j++) {
745 if (pointBoundaryNodes[row*blkPartSize+j]) {
752 for (LO j = 0; j < blkPartSize; j++) {
753 if (!pointBoundaryNodes[row*blkPartSize+j]) {
763 MergeRows(*A, row, indicesExtra, colTranslation);
765 indicesExtra.push_back(row);
766 indices = indicesExtra;
767 numTotal += indices.size();
771 LO nnz = indices.size(), rownnz = 0;
772 for (LO colID = 0; colID < nnz; colID++) {
773 LO col = indices[colID];
774 columns[realnnz++] = col;
785 amalgBoundaryNodes[row] =
true;
787 rows[row+1] = realnnz;
789 columns.resize(realnnz);
791 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
792 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
795 GO numLocalBoundaryNodes = 0;
796 GO numGlobalBoundaryNodes = 0;
798 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
799 if (amalgBoundaryNodes[i])
800 numLocalBoundaryNodes++;
802 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
803 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
804 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
805 <<
" agglomerated Dirichlet nodes" << std::endl;
808 Set(currentLevel,
"Graph", graph);
809 Set(currentLevel,
"DofsPerNode", blkSize);
811 }
else if (A->GetFixedBlockSize() > 1 && threshold != STS::zero()) {
813 const RCP<const Map> rowMap = A->getRowMap();
814 const RCP<const Map> colMap = A->getColMap();
815 graphType =
"amalgamated";
821 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
822 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
823 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
824 Array<LO> colTranslation = *(amalInfo->getColTranslation());
827 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
830 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
831 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
833 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
839 ArrayRCP<bool > pointBoundaryNodes;
846 LO blkSize = A->GetFixedBlockSize();
848 LO blkPartSize = A->GetFixedBlockSize();
849 if (A->IsView(
"stridedMaps") ==
true) {
850 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
851 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
853 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
854 blkId = strMap->getStridedBlockId();
856 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
861 const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
866 Array<LO> indicesExtra;
867 for (LO row = 0; row < numRows; row++) {
868 ArrayView<const LO> indices;
869 indicesExtra.resize(0);
877 bool isBoundary =
false;
878 if (pL.get<
bool>(
"aggregation: greedy Dirichlet") ==
true) {
879 for (LO j = 0; j < blkPartSize; j++) {
880 if (pointBoundaryNodes[row*blkPartSize+j]) {
887 for (LO j = 0; j < blkPartSize; j++) {
888 if (!pointBoundaryNodes[row*blkPartSize+j]) {
898 MergeRowsWithDropping(*A, row, ghostedDiagVals, threshold, indicesExtra, colTranslation);
900 indicesExtra.push_back(row);
901 indices = indicesExtra;
902 numTotal += indices.size();
906 LO nnz = indices.size(), rownnz = 0;
907 for (LO colID = 0; colID < nnz; colID++) {
908 LO col = indices[colID];
909 columns[realnnz++] = col;
920 amalgBoundaryNodes[row] =
true;
922 rows[row+1] = realnnz;
924 columns.resize(realnnz);
926 RCP<GraphBase> graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
927 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
930 GO numLocalBoundaryNodes = 0;
931 GO numGlobalBoundaryNodes = 0;
933 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
934 if (amalgBoundaryNodes[i])
935 numLocalBoundaryNodes++;
937 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
938 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
939 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
940 <<
" agglomerated Dirichlet nodes" << std::endl;
943 Set(currentLevel,
"Graph", graph);
944 Set(currentLevel,
"DofsPerNode", blkSize);
947 }
else if (algo ==
"distance laplacian") {
948 LO blkSize = A->GetFixedBlockSize();
949 GO indexBase = A->getRowMap()->getIndexBase();
960 ArrayRCP<bool > pointBoundaryNodes;
965 if ( (blkSize == 1) && (threshold == STS::zero()) ) {
967 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
968 graph->SetBoundaryNodeMap(pointBoundaryNodes);
969 graphType=
"unamalgamated";
970 numTotal = A->getNodeNumEntries();
973 GO numLocalBoundaryNodes = 0;
974 GO numGlobalBoundaryNodes = 0;
975 for (LO i = 0; i < pointBoundaryNodes.size(); ++i)
976 if (pointBoundaryNodes[i])
977 numLocalBoundaryNodes++;
978 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
979 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
980 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
983 Set(currentLevel,
"DofsPerNode", blkSize);
984 Set(currentLevel,
"Graph", graph);
999 TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getNodeNumElements()/blkSize != Coords->getLocalLength(),
Exceptions::Incompatible,
1000 "Coordinate vector length (" << Coords->getLocalLength() <<
") is incompatible with number of rows in A (" << A->getRowMap()->getNodeNumElements() <<
") by modulo block size ("<< blkSize <<
").");
1002 const RCP<const Map> colMap = A->getColMap();
1003 RCP<const Map> uniqueMap, nonUniqueMap;
1004 Array<LO> colTranslation;
1006 uniqueMap = A->getRowMap();
1007 nonUniqueMap = A->getColMap();
1008 graphType=
"unamalgamated";
1011 uniqueMap = Coords->getMap();
1013 "Different index bases for matrix and coordinates");
1017 graphType =
"amalgamated";
1019 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
1021 RCP<RealValuedMultiVector> ghostedCoords;
1022 RCP<Vector> ghostedLaplDiag;
1023 Teuchos::ArrayRCP<SC> ghostedLaplDiagData;
1024 if (threshold != STS::zero()) {
1026 RCP<const Import> importer;
1029 if (blkSize == 1 && realA->getCrsGraph()->getImporter() != Teuchos::null) {
1030 GetOStream(
Warnings1) <<
"Using existing importer from matrix graph" << std::endl;
1031 importer = realA->getCrsGraph()->getImporter();
1033 GetOStream(
Warnings0) <<
"Constructing new importer instance" << std::endl;
1034 importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
1037 ghostedCoords = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(nonUniqueMap, Coords->getNumVectors());
1040 ghostedCoords->doImport(*Coords, *importer, Xpetra::INSERT);
1044 RCP<Vector> localLaplDiag = VectorFactory::Build(uniqueMap);
1045 ArrayRCP<SC> localLaplDiagData = localLaplDiag->getDataNonConst(0);
1046 Array<LO> indicesExtra;
1047 Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
1048 if (threshold != STS::zero()) {
1049 const size_t numVectors = ghostedCoords->getNumVectors();
1050 coordData.reserve(numVectors);
1051 for (
size_t j = 0; j < numVectors; j++) {
1052 Teuchos::ArrayRCP<const real_type> tmpData=ghostedCoords->getData(j);
1053 coordData.push_back(tmpData);
1058 for (LO row = 0; row < numRows; row++) {
1059 ArrayView<const LO> indices;
1062 ArrayView<const SC> vals;
1063 A->getLocalRowView(row, indices, vals);
1067 indicesExtra.resize(0);
1068 MergeRows(*A, row, indicesExtra, colTranslation);
1069 indices = indicesExtra;
1072 LO nnz = indices.size();
1073 bool haveAddedToDiag =
false;
1074 for (LO colID = 0; colID < nnz; colID++) {
1075 const LO col = indices[colID];
1078 if(use_dlap_weights == SINGLE_WEIGHTS) {
1084 else if(use_dlap_weights == BLOCK_WEIGHTS) {
1085 int block_id = row % interleaved_blocksize;
1086 int block_start = block_id * interleaved_blocksize;
1093 haveAddedToDiag =
true;
1098 if (!haveAddedToDiag)
1099 localLaplDiagData[row] = STS::rmax();
1104 ghostedLaplDiag = VectorFactory::Build(nonUniqueMap);
1105 ghostedLaplDiag->doImport(*localLaplDiag, *importer, Xpetra::INSERT);
1106 ghostedLaplDiagData = ghostedLaplDiag->getDataNonConst(0);
1110 GetOStream(
Runtime0) <<
"Skipping distance laplacian construction due to 0 threshold" << std::endl;
1116 ArrayRCP<LO> rows = ArrayRCP<LO>(numRows+1);
1117 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getNodeNumEntries());
1119 #ifdef HAVE_MUELU_DEBUG 1121 for(LO i=0; i<(LO)columns.size(); i++) columns[i]=-666;
1125 ArrayRCP<LO> rows_stop;
1126 bool use_stop_array = threshold != STS::zero() && distanceLaplacianAlgo == scaled_cut_symmetric;
1128 rows_stop.resize(numRows);
1130 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
1135 Array<LO> indicesExtra;
1138 Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
1139 if (threshold != STS::zero()) {
1140 const size_t numVectors = ghostedCoords->getNumVectors();
1141 coordData.reserve(numVectors);
1142 for (
size_t j = 0; j < numVectors; j++) {
1143 Teuchos::ArrayRCP<const real_type> tmpData=ghostedCoords->getData(j);
1144 coordData.push_back(tmpData);
1148 ArrayView<const SC> vals;
1149 for (LO row = 0; row < numRows; row++) {
1150 ArrayView<const LO> indices;
1151 indicesExtra.resize(0);
1152 bool isBoundary =
false;
1156 A->getLocalRowView(row, indices, vals);
1157 isBoundary = pointBoundaryNodes[row];
1160 for (LO j = 0; j < blkSize; j++) {
1161 if (!pointBoundaryNodes[row*blkSize+j]) {
1169 MergeRows(*A, row, indicesExtra, colTranslation);
1171 indicesExtra.push_back(row);
1172 indices = indicesExtra;
1174 numTotal += indices.size();
1176 LO nnz = indices.size(), rownnz = 0;
1178 if(use_stop_array) {
1179 rows[row+1] = rows[row]+nnz;
1180 realnnz = rows[row];
1183 if (threshold != STS::zero()) {
1185 if (distanceLaplacianAlgo == defaultAlgo) {
1187 for (LO colID = 0; colID < nnz; colID++) {
1189 LO col = indices[colID];
1192 columns[realnnz++] = col;
1198 if(isBoundary)
continue;
1201 if(use_dlap_weights == SINGLE_WEIGHTS) {
1204 else if(use_dlap_weights == BLOCK_WEIGHTS) {
1205 int block_id = row % interleaved_blocksize;
1206 int block_start = block_id * interleaved_blocksize;
1212 real_type aiiajj = STS::magnitude(realThreshold*realThreshold * ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
1213 real_type aij = STS::magnitude(laplVal*laplVal);
1216 columns[realnnz++] = col;
1225 std::vector<DropTol> drop_vec;
1226 drop_vec.reserve(nnz);
1227 const real_type zero = Teuchos::ScalarTraits<real_type>::zero();
1228 const real_type one = Teuchos::ScalarTraits<real_type>::one();
1231 for (LO colID = 0; colID < nnz; colID++) {
1233 LO col = indices[colID];
1236 drop_vec.emplace_back( zero, one, colID,
false);
1240 if(isBoundary)
continue;
1243 if(use_dlap_weights == SINGLE_WEIGHTS) {
1246 else if(use_dlap_weights == BLOCK_WEIGHTS) {
1247 int block_id = row % interleaved_blocksize;
1248 int block_start = block_id * interleaved_blocksize;
1255 real_type aiiajj = STS::magnitude(ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
1256 real_type aij = STS::magnitude(laplVal*laplVal);
1258 drop_vec.emplace_back(aij, aiiajj, colID,
false);
1261 const size_t n = drop_vec.size();
1263 if (distanceLaplacianAlgo == unscaled_cut) {
1265 std::sort( drop_vec.begin(), drop_vec.end()
1266 , [](DropTol
const& a, DropTol
const& b) {
1267 return a.val > b.val;
1272 for (
size_t i=1; i<n; ++i) {
1274 auto const& x = drop_vec[i-1];
1275 auto const& y = drop_vec[i];
1278 if (a > realThreshold*b) {
1280 #ifdef HAVE_MUELU_DEBUG 1281 if (distanceLaplacianCutVerbose) {
1282 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
1287 drop_vec[i].drop = drop;
1290 else if (distanceLaplacianAlgo == scaled_cut || distanceLaplacianAlgo == scaled_cut_symmetric) {
1292 std::sort( drop_vec.begin(), drop_vec.end()
1293 , [](DropTol
const& a, DropTol
const& b) {
1294 return a.val/a.diag > b.val/b.diag;
1299 for (
size_t i=1; i<n; ++i) {
1301 auto const& x = drop_vec[i-1];
1302 auto const& y = drop_vec[i];
1303 auto a = x.val/x.diag;
1304 auto b = y.val/y.diag;
1305 if (a > realThreshold*b) {
1307 #ifdef HAVE_MUELU_DEBUG 1308 if (distanceLaplacianCutVerbose) {
1309 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
1314 drop_vec[i].drop = drop;
1318 std::sort( drop_vec.begin(), drop_vec.end()
1319 , [](DropTol
const& a, DropTol
const& b) {
1320 return a.col < b.col;
1324 for (LO idxID =0; idxID<(LO)drop_vec.size(); idxID++) {
1325 LO col = indices[drop_vec[idxID].col];
1330 columns[realnnz++] = col;
1336 if (!drop_vec[idxID].drop) {
1337 columns[realnnz++] = col;
1348 for (LO colID = 0; colID < nnz; colID++) {
1349 LO col = indices[colID];
1350 columns[realnnz++] = col;
1362 amalgBoundaryNodes[row] =
true;
1366 rows_stop[row] = rownnz + rows[row];
1368 rows[row+1] = realnnz;
1373 if (use_stop_array) {
1376 for (LO row = 0; row < numRows; row++) {
1377 for (LO colidx = rows[row]; colidx < rows_stop[row]; colidx++) {
1378 LO col = columns[colidx];
1379 if(col >= numRows)
continue;
1382 for(LO t_col = rows[col] ; !found && t_col < rows_stop[col]; t_col++) {
1383 if (columns[t_col] == row)
1388 if(!found && !pointBoundaryNodes[col] && rows_stop[col] < rows[col+1]) {
1389 LO new_idx = rows_stop[col];
1391 columns[new_idx] = row;
1400 for (LO row = 0; row < numRows; row++) {
1401 LO old_start = current_start;
1402 for (LO col = rows[row]; col < rows_stop[row]; col++) {
1403 if(current_start != col) {
1404 columns[current_start] = columns[col];
1408 rows[row] = old_start;
1410 rows[numRows] = realnnz = current_start;
1414 columns.resize(realnnz);
1416 RCP<GraphBase> graph;
1419 graph = rcp(
new LWGraph(rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
1420 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
1424 GO numLocalBoundaryNodes = 0;
1425 GO numGlobalBoundaryNodes = 0;
1427 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
1428 if (amalgBoundaryNodes[i])
1429 numLocalBoundaryNodes++;
1431 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1432 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1433 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" agglomerated Dirichlet nodes" 1434 <<
" using threshold " << dirichletThreshold << std::endl;
1437 Set(currentLevel,
"Graph", graph);
1438 Set(currentLevel,
"DofsPerNode", blkSize);
1442 if ((GetVerbLevel() &
Statistics1) && !(A->GetFixedBlockSize() > 1 && threshold != STS::zero())) {
1443 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1444 GO numGlobalTotal, numGlobalDropped;
1447 GetOStream(
Statistics1) <<
"Number of dropped entries in " << graphType <<
" matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1448 if (numGlobalTotal != 0)
1449 GetOStream(
Statistics1) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
1456 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
1458 GetOStream(
Runtime0) <<
"algorithm = \"" <<
"failsafe" <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
1459 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
1461 RCP<const Map> rowMap = A->getRowMap();
1462 RCP<const Map> colMap = A->getColMap();
1465 GO indexBase = rowMap->getIndexBase();
1469 if(A->IsView(
"stridedMaps") &&
1470 Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1471 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
1472 RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap());
1473 TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
1474 blockdim = strMap->getFixedBlockSize();
1475 offset = strMap->getOffset();
1476 oldView = A->SwitchToView(oldView);
1477 GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build():" <<
" found blockdim=" << blockdim <<
" from strided maps. offset=" << offset << std::endl;
1478 }
else GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
1482 RCP<const Map> nodeMap = amalInfo->getNodeRowMap();
1483 GetOStream(
Statistics1) <<
"CoalesceDropFactory: nodeMap " << nodeMap->getNodeNumElements() <<
"/" << nodeMap->getGlobalNumElements() <<
" elements" << std::endl;
1486 RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, A->getNodeMaxNumRowEntries()*blockdim);
1488 LO numRows = A->getRowMap()->getNodeNumElements();
1489 LO numNodes = nodeMap->getNodeNumElements();
1490 const ArrayRCP<bool> amalgBoundaryNodes(numNodes,
false);
1491 const ArrayRCP<int> numberDirichletRowsPerNode(numNodes, 0);
1492 bool bIsDiagonalEntry =
false;
1497 for(LO row=0; row<numRows; row++) {
1499 GO grid = rowMap->getGlobalElement(row);
1502 bIsDiagonalEntry =
false;
1507 size_t nnz = A->getNumEntriesInLocalRow(row);
1508 Teuchos::ArrayView<const LO> indices;
1509 Teuchos::ArrayView<const SC> vals;
1510 A->getLocalRowView(row, indices, vals);
1512 RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(
new std::vector<GO>);
1514 for(LO col=0; col<Teuchos::as<LO>(nnz); col++) {
1515 GO gcid = colMap->getGlobalElement(indices[col]);
1517 if(vals[col]!=STS::zero()) {
1519 cnodeIds->push_back(cnodeId);
1521 if (grid == gcid) bIsDiagonalEntry =
true;
1525 if(realnnz == 1 && bIsDiagonalEntry ==
true) {
1526 LO lNodeId = nodeMap->getLocalElement(nodeId);
1527 numberDirichletRowsPerNode[lNodeId] += 1;
1528 if (numberDirichletRowsPerNode[lNodeId] == blockdim)
1529 amalgBoundaryNodes[lNodeId] =
true;
1532 Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp( cnodeIds );
1534 if(arr_cnodeIds.size() > 0 )
1535 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
1538 crsGraph->fillComplete(nodeMap,nodeMap);
1541 RCP<GraphBase> graph = rcp(
new Graph(crsGraph,
"amalgamated graph of A"));
1544 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
1547 GO numLocalBoundaryNodes = 0;
1548 GO numGlobalBoundaryNodes = 0;
1549 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
1550 if (amalgBoundaryNodes[i])
1551 numLocalBoundaryNodes++;
1552 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1553 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1554 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
1559 Set(currentLevel,
"DofsPerNode", blockdim);
1560 Set(currentLevel,
"Graph", graph);
1567 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1569 typedef typename ArrayView<const LO>::size_type size_type;
1572 LO blkSize = A.GetFixedBlockSize();
1573 if (A.IsView(
"stridedMaps") ==
true) {
1574 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
1575 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
1577 if (strMap->getStridedBlockId() > -1)
1578 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1582 size_t nnz = 0, pos = 0;
1583 for (LO j = 0; j < blkSize; j++)
1584 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
1594 ArrayView<const LO> inds;
1595 ArrayView<const SC> vals;
1596 for (LO j = 0; j < blkSize; j++) {
1597 A.getLocalRowView(row*blkSize+j, inds, vals);
1598 size_type numIndices = inds.size();
1600 if (numIndices == 0)
1604 cols[pos++] = translation[inds[0]];
1605 for (size_type k = 1; k < numIndices; k++) {
1606 LO nodeID = translation[inds[k]];
1610 if (nodeID != cols[pos-1])
1611 cols[pos++] = nodeID;
1618 std::sort(cols.begin(), cols.end());
1620 for (
size_t j = 1; j < nnz; j++)
1621 if (cols[j] != cols[pos])
1622 cols[++pos] = cols[j];
1626 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1628 typedef typename ArrayView<const LO>::size_type size_type;
1629 typedef Teuchos::ScalarTraits<SC> STS;
1632 LO blkSize = A.GetFixedBlockSize();
1633 if (A.IsView(
"stridedMaps") ==
true) {
1634 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
1635 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<
const StridedMap>(myMap);
1637 if (strMap->getStridedBlockId() > -1)
1638 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1642 size_t nnz = 0, pos = 0;
1643 for (LO j = 0; j < blkSize; j++)
1644 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
1654 ArrayView<const LO> inds;
1655 ArrayView<const SC> vals;
1656 for (LO j = 0; j < blkSize; j++) {
1657 A.getLocalRowView(row*blkSize+j, inds, vals);
1658 size_type numIndices = inds.size();
1660 if (numIndices == 0)
1665 for (size_type k = 0; k < numIndices; k++) {
1667 LO nodeID = translation[inds[k]];
1670 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold*ghostedDiagVals[dofID]*ghostedDiagVals[row*blkSize+j]);
1671 typename STS::magnitudeType aij = STS::magnitude(vals[k]*vals[k]);
1674 if (aij > aiiajj || (row*blkSize+j == dofID)) {
1680 if (nodeID != prevNodeID) {
1681 cols[pos++] = nodeID;
1682 prevNodeID = nodeID;
1691 std::sort(cols.begin(), cols.end());
1693 for (
size_t j = 1; j < nnz; j++)
1694 if (cols[j] != cols[pos])
1695 cols[++pos] = cols[j];
1703 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1705 typedef Teuchos::ScalarTraits<SC> STS;
1707 const ParameterList & pL = GetParameterList();
1708 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<
double>(
"aggregation: Dirichlet threshold")));
1709 const typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<
double>(
"aggregation: row sum drop tol"));
1711 RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(currentLevel,
"BlockNumber");
1712 RCP<LocalOrdinalVector> ghostedBlockNumber;
1713 GetOStream(
Statistics1) <<
"Using BlockDiagonal Graph before dropping (with provided blocking)"<<std::endl;
1716 RCP<const Import> importer = A->getCrsGraph()->getImporter();
1717 if(!importer.is_null()) {
1719 ghostedBlockNumber= Xpetra::VectorFactory<LO,LO,GO,NO>::Build(importer->getTargetMap());
1720 ghostedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
1723 ghostedBlockNumber = BlockNumber;
1727 Teuchos::ArrayRCP<const LO> row_block_number = BlockNumber->getData(0);
1728 Teuchos::ArrayRCP<const LO> col_block_number = ghostedBlockNumber->getData(0);
1731 ArrayRCP<size_t> rows_mat;
1732 ArrayRCP<LO> rows_graph,columns;
1733 ArrayRCP<SC> values;
1734 RCP<CrsMatrixWrap> crs_matrix_wrap;
1736 if(generate_matrix) {
1737 crs_matrix_wrap = rcp(
new CrsMatrixWrap(A->getRowMap(), A->getColMap(), 0));
1738 crs_matrix_wrap->getCrsMatrix()->allocateAllValues(A->getNodeNumEntries(), rows_mat, columns, values);
1741 rows_graph.resize(A->getNodeNumRows()+1);
1742 columns.resize(A->getNodeNumEntries());
1743 values.resize(A->getNodeNumEntries());
1747 GO numDropped = 0, numTotal = 0;
1748 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getNodeNumElements()); ++row) {
1749 LO row_block = row_block_number[row];
1750 size_t nnz = A->getNumEntriesInLocalRow(row);
1751 ArrayView<const LO> indices;
1752 ArrayView<const SC> vals;
1753 A->getLocalRowView(row, indices, vals);
1756 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
1757 LO col = indices[colID];
1758 LO col_block = col_block_number[col];
1760 if(row_block == col_block) {
1761 if(generate_matrix) values[realnnz] = vals[colID];
1762 columns[realnnz++] = col;
1767 if(generate_matrix) rows_mat[row+1] = realnnz;
1768 else rows_graph[row+1] = realnnz;
1776 if(!generate_matrix) {
1778 values.resize(realnnz);
1779 columns.resize(realnnz);
1781 numTotal = A->getNodeNumEntries();
1784 GO numLocalBoundaryNodes = 0;
1785 GO numGlobalBoundaryNodes = 0;
1786 for (LO i = 0; i < boundaryNodes.size(); ++i)
1787 if (boundaryNodes[i])
1788 numLocalBoundaryNodes++;
1789 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1790 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1791 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
1793 GO numGlobalTotal, numGlobalDropped;
1796 GetOStream(
Statistics1) <<
"Number of dropped entries in block-diagonalized matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1797 if (numGlobalTotal != 0)
1798 GetOStream(
Statistics1) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
1802 Set(currentLevel,
"Filtering",
true);
1804 if(generate_matrix) {
1808 crs_matrix_wrap->getCrsMatrix()->setAllValues(rows_mat,columns,values);
1809 crs_matrix_wrap->getCrsMatrix()->expertStaticFillComplete(A->getColMap(), A->getRowMap());
1812 RCP<GraphBase> graph = rcp(
new LWGraph(rows_graph, columns, A->getRowMap(), A->getColMap(),
"block-diagonalized graph of A"));
1813 graph->SetBoundaryNodeMap(boundaryNodes);
1814 Set(currentLevel,
"Graph", graph);
1818 Set(currentLevel,
"DofsPerNode", 1);
1819 return crs_matrix_wrap;
1823 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1825 typedef Teuchos::ScalarTraits<SC> STS;
1827 TEUCHOS_TEST_FOR_EXCEPTION(ghostedBlockNumber.is_null(),
Exceptions::RuntimeError,
"BlockDiagonalizeGraph(): ghostedBlockNumber is null.");
1828 const ParameterList & pL = GetParameterList();
1830 const bool localizeColoringGraph = pL.get<
bool>(
"aggregation: coloring: localize color graph");
1832 GetOStream(
Statistics1) <<
"Using BlockDiagonal Graph after Dropping (with provided blocking)";
1833 if (localizeColoringGraph)
1834 GetOStream(
Statistics1) <<
", with localization" <<std::endl;
1836 GetOStream(
Statistics1) <<
", without localization" <<std::endl;
1839 Teuchos::ArrayRCP<const LO> row_block_number = ghostedBlockNumber->getData(0);
1840 Teuchos::ArrayRCP<const LO> col_block_number = ghostedBlockNumber->getData(0);
1843 ArrayRCP<size_t> rows_mat;
1844 ArrayRCP<LO> rows_graph,columns;
1846 rows_graph.resize(inputGraph->GetNodeNumVertices()+1);
1847 columns.resize(inputGraph->GetNodeNumEdges());
1850 GO numDropped = 0, numTotal = 0;
1851 const LO numRows = Teuchos::as<LO>(inputGraph->GetDomainMap()->getNodeNumElements());
1852 if (localizeColoringGraph) {
1854 for (LO row = 0; row < numRows; ++row) {
1855 LO row_block = row_block_number[row];
1856 ArrayView<const LO> indices = inputGraph->getNeighborVertices(row);
1859 for (LO colID = 0; colID < Teuchos::as<LO>(indices.size()); colID++) {
1860 LO col = indices[colID];
1861 LO col_block = col_block_number[col];
1863 if((row_block == col_block) && (col < numRows)) {
1864 columns[realnnz++] = col;
1869 rows_graph[row+1] = realnnz;
1873 Teuchos::ArrayRCP<const bool> boundaryNodes = inputGraph->GetBoundaryNodeMap();
1874 auto boundaryNodesVector = Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(inputGraph->GetDomainMap());
1875 for (
size_t i=0; i<inputGraph->GetNodeNumVertices(); i++)
1876 boundaryNodesVector->getDataNonConst(0)[i] = boundaryNodes[i];
1878 auto boundaryColumnVector = Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(inputGraph->GetImportMap());
1879 boundaryColumnVector->doImport(*boundaryNodesVector,*importer, Xpetra::INSERT);
1880 auto boundaryColumn = boundaryColumnVector->getData(0);
1882 for (LO row = 0; row < numRows; ++row) {
1883 LO row_block = row_block_number[row];
1884 ArrayView<const LO> indices = inputGraph->getNeighborVertices(row);
1887 for (LO colID = 0; colID < Teuchos::as<LO>(indices.size()); colID++) {
1888 LO col = indices[colID];
1889 LO col_block = col_block_number[col];
1891 if((row_block == col_block) && ((row == col) || (boundaryColumn[col] == 0))) {
1892 columns[realnnz++] = col;
1897 rows_graph[row+1] = realnnz;
1901 columns.resize(realnnz);
1902 numTotal = inputGraph->GetNodeNumEdges();
1905 RCP<const Teuchos::Comm<int> > comm = inputGraph->GetDomainMap()->getComm();
1906 GO numGlobalTotal, numGlobalDropped;
1909 GetOStream(
Statistics1) <<
"Number of dropped entries in block-diagonalized matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1910 if (numGlobalTotal != 0)
1911 GetOStream(
Statistics1) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
1915 if (localizeColoringGraph) {
1916 outputGraph = rcp(
new LWGraph(rows_graph, columns, inputGraph->GetDomainMap(), inputGraph->GetImportMap(),
"block-diagonalized graph of A"));
1917 outputGraph->SetBoundaryNodeMap(inputGraph->GetBoundaryNodeMap());
1919 TEUCHOS_ASSERT(inputGraph->GetDomainMap()->lib() == Xpetra::UseTpetra);
1920 #ifdef HAVE_XPETRA_TPETRA 1921 auto outputGraph2 = rcp(
new LWGraph(rows_graph, columns, inputGraph->GetDomainMap(), inputGraph->GetImportMap(),
"block-diagonalized graph of A"));
1923 auto tpGraph = Xpetra::toTpetra(rcp_const_cast<const CrsGraph>(outputGraph2->GetCrsGraph()));
1924 auto sym = rcp(
new Tpetra::CrsGraphTransposer<LocalOrdinal,GlobalOrdinal,Node>(tpGraph));
1925 auto tpGraphSym = sym->symmetrize();
1927 auto rowsSym = tpGraphSym->getNodeRowPtrs();
1928 ArrayRCP<LO> rows_graphSym;
1929 rows_graphSym.resize(rowsSym.size());
1930 for (LO row = 0; row < rowsSym.size(); row++)
1931 rows_graphSym[row] = rowsSym[row];
1932 outputGraph = rcp(
new LWGraph(rows_graphSym, tpGraphSym->getNodePackedIndices(), inputGraph->GetDomainMap(), Xpetra::toXpetra(tpGraphSym->getColMap()),
"block-diagonalized graph of A"));
1933 outputGraph->SetBoundaryNodeMap(inputGraph->GetBoundaryNodeMap());
1943 #endif // MUELU_COALESCEDROPFACTORY_DEF_HPP Important warning messages (one line)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
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.
Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > BlockDiagonalize(Level ¤tLevel, const RCP< Matrix > &A, bool generate_matrix) const
static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(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.
void BlockDiagonalizeGraph(const RCP< GraphBase > &inputGraph, const RCP< LocalOrdinalVector > &ghostedBlockNumber, RCP< GraphBase > &outputGraph, RCP< const Import > &importer) const