46 #ifndef MUELU_COALESCEDROPFACTORY_DEF_HPP
47 #define MUELU_COALESCEDROPFACTORY_DEF_HPP
53 #include <Xpetra_MapFactory.hpp>
54 #include <Xpetra_Map.hpp>
56 #include <Xpetra_MultiVectorFactory.hpp>
58 #include <Xpetra_StridedMap.hpp>
66 #include "MueLu_AmalgamationFactory.hpp"
67 #include "MueLu_AmalgamationInfo.hpp"
69 #include "MueLu_LWGraph.hpp"
72 #include "MueLu_LWGraph.hpp"
75 #include "MueLu_PreDropFunctionConstVal.hpp"
76 #include "MueLu_Utilities.hpp"
78 #ifdef HAVE_XPETRA_TPETRA
79 #include "Tpetra_CrsGraphTransposer.hpp"
93 template <
class real_type,
class LO>
102 DropTol(real_type val_, real_type diag_,
LO col_,
bool drop_)
118 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
122 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name))
130 SET_VALID_ENTRY(
"aggregation: distance laplacian directional weights");
136 validParamList->
getEntry(
"aggregation: drop scheme").
setValidator(
rcp(
new validatorType(Teuchos::tuple<std::string>(
"signed classical sa",
"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");
145 validParamList->
set<
RCP<const FactoryBase>>(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory for UnAmalgamationInfo");
149 return validParamList;
152 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
154 : predrop_(Teuchos::null) {}
156 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
158 Input(currentLevel,
"A");
159 Input(currentLevel,
"UnAmalgamationInfo");
162 if (pL.
get<
bool>(
"lightweight wrap") ==
true) {
163 std::string algo = pL.
get<std::string>(
"aggregation: drop scheme");
164 if (algo ==
"distance laplacian" || algo ==
"block diagonal distance laplacian") {
165 Input(currentLevel,
"Coordinates");
167 if (algo ==
"signed classical sa")
169 else if (algo.find(
"block diagonal") != std::string::npos || algo.find(
"signed classical") != std::string::npos) {
170 Input(currentLevel,
"BlockNumber");
175 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
180 typedef typename STS::magnitudeType real_type;
184 if (predrop_ != Teuchos::null)
185 GetOStream(
Parameters0) << predrop_->description();
187 RCP<Matrix> realA = Get<RCP<Matrix>>(currentLevel,
"A");
190 bool doExperimentalWrap = pL.get<
bool>(
"lightweight wrap");
192 GetOStream(
Parameters0) <<
"lightweight wrap = " << doExperimentalWrap << std::endl;
193 std::string algo = pL.get<std::string>(
"aggregation: drop scheme");
194 const bool aggregationMayCreateDirichlet = pL.get<
bool>(
"aggregation: dropping may create Dirichlet");
199 bool use_block_algorithm =
false;
200 LO interleaved_blocksize = as<LO>(pL.get<
int>(
"aggregation: block diagonal: interleaved blocksize"));
201 bool useSignedClassicalRS =
false;
202 bool useSignedClassicalSA =
false;
203 bool generateColoringGraph =
false;
207 typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<
double>(
"aggregation: row sum drop tol"));
212 if (algo ==
"distance laplacian") {
214 Coords = Get<RCP<RealValuedMultiVector>>(currentLevel,
"Coordinates");
216 }
else if (algo ==
"signed classical sa") {
217 useSignedClassicalSA =
true;
220 }
else if (algo ==
"signed classical" || algo ==
"block diagonal colored signed classical" || algo ==
"block diagonal signed classical") {
221 useSignedClassicalRS =
true;
226 if (!importer.is_null()) {
229 ghostedBlockNumber->doImport(*BlockNumber, *importer,
Xpetra::INSERT);
231 ghostedBlockNumber = BlockNumber;
233 g_block_id = ghostedBlockNumber->getData(0);
235 if (algo ==
"block diagonal colored signed classical")
236 generateColoringGraph =
true;
240 }
else if (algo ==
"block diagonal") {
242 BlockDiagonalize(currentLevel, realA,
false);
244 }
else if (algo ==
"block diagonal classical" || algo ==
"block diagonal distance laplacian") {
246 use_block_algorithm =
true;
247 RCP<Matrix> filteredMatrix = BlockDiagonalize(currentLevel, realA,
true);
248 if (algo ==
"block diagonal distance laplacian") {
251 if (OldCoords->getLocalLength() != realA->getLocalNumRows()) {
252 LO dim = (
LO)OldCoords->getNumVectors();
253 Coords = RealValuedMultiVectorFactory::Build(realA->getRowMap(), dim);
254 for (
LO k = 0; k < dim; k++) {
257 for (
LO i = 0; i < (
LO)OldCoords->getLocalLength(); i++) {
258 LO new_base = i * dim;
259 for (
LO j = 0; j < interleaved_blocksize; j++)
260 new_vec[new_base + j] = old_vec[i];
266 algo =
"distance laplacian";
267 }
else if (algo ==
"block diagonal classical") {
279 enum { NO_WEIGHTS = 0,
282 int use_dlap_weights = NO_WEIGHTS;
283 if (algo ==
"distance laplacian") {
284 LO dim = (
LO)Coords->getNumVectors();
286 bool non_unity =
false;
287 for (
LO i = 0; !non_unity && i < (
LO)dlap_weights.
size(); i++) {
288 if (dlap_weights[i] != 1.0) {
293 LO blocksize = use_block_algorithm ? as<LO>(pL.get<
int>(
"aggregation: block diagonal: interleaved blocksize")) : 1;
294 if ((
LO)dlap_weights.
size() == dim)
295 use_dlap_weights = SINGLE_WEIGHTS;
296 else if ((
LO)dlap_weights.
size() == blocksize * dim)
297 use_dlap_weights = BLOCK_WEIGHTS;
300 "length of 'aggregation: distance laplacian directional weights' must equal the coordinate dimension OR the coordinate dimension times the blocksize");
303 GetOStream(
Statistics1) <<
"Using distance laplacian weights: " << dlap_weights << std::endl;
318 if (doExperimentalWrap) {
324 if (pL.get<
bool>(
"aggregation: use ml scaling of drop tol"))
325 threshold = pL.get<
double>(
"aggregation: drop tol") / pow(2.0, currentLevel.GetLevelID());
327 threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
329 std::string distanceLaplacianAlgoStr = pL.get<std::string>(
"aggregation: distance laplacian algo");
330 std::string classicalAlgoStr = pL.get<std::string>(
"aggregation: classical algo");
331 real_type realThreshold = STS::magnitude(threshold);
335 #ifdef HAVE_MUELU_DEBUG
336 int distanceLaplacianCutVerbose = 0;
338 #ifdef DJS_READ_ENV_VARIABLES
339 if (getenv(
"MUELU_DROP_TOLERANCE_MODE")) {
340 distanceLaplacianAlgoStr = std::string(getenv(
"MUELU_DROP_TOLERANCE_MODE"));
343 if (getenv(
"MUELU_DROP_TOLERANCE_THRESHOLD")) {
344 auto tmp = atoi(getenv(
"MUELU_DROP_TOLERANCE_THRESHOLD"));
345 realThreshold = 1e-4 * tmp;
348 #ifdef HAVE_MUELU_DEBUG
349 if (getenv(
"MUELU_DROP_TOLERANCE_VERBOSE")) {
350 distanceLaplacianCutVerbose = atoi(getenv(
"MUELU_DROP_TOLERANCE_VERBOSE"));
356 enum decisionAlgoType { defaultAlgo,
359 scaled_cut_symmetric };
361 decisionAlgoType distanceLaplacianAlgo = defaultAlgo;
362 decisionAlgoType classicalAlgo = defaultAlgo;
363 if (algo ==
"distance laplacian") {
364 if (distanceLaplacianAlgoStr ==
"default")
365 distanceLaplacianAlgo = defaultAlgo;
366 else if (distanceLaplacianAlgoStr ==
"unscaled cut")
367 distanceLaplacianAlgo = unscaled_cut;
368 else if (distanceLaplacianAlgoStr ==
"scaled cut")
369 distanceLaplacianAlgo = scaled_cut;
370 else if (distanceLaplacianAlgoStr ==
"scaled cut symmetric")
371 distanceLaplacianAlgo = scaled_cut_symmetric;
374 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\" distance laplacian algorithm = \"" << distanceLaplacianAlgoStr <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
375 }
else if (algo ==
"classical") {
376 if (classicalAlgoStr ==
"default")
377 classicalAlgo = defaultAlgo;
378 else if (classicalAlgoStr ==
"unscaled cut")
379 classicalAlgo = unscaled_cut;
380 else if (classicalAlgoStr ==
"scaled cut")
381 classicalAlgo = scaled_cut;
384 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\" classical algorithm = \"" << classicalAlgoStr <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
387 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
388 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
390 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<
double>(
"aggregation: Dirichlet threshold")));
396 GO numDropped = 0, numTotal = 0;
397 std::string graphType =
"unamalgamated";
416 const LO BlockSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
419 if (algo ==
"classical") {
420 if (predrop_ == null) {
425 if (predrop_ != null) {
428 "MueLu::CoalesceFactory::Build: cast to PreDropFunctionConstVal failed.");
431 if (newt != threshold) {
432 GetOStream(
Warnings0) <<
"switching threshold parameter from " << threshold <<
" (list) to " << newt <<
" (user function" << std::endl;
439 if (BlockSize == 1 && threshold == STS::zero() && !useSignedClassicalRS && !useSignedClassicalSA && A->hasCrsGraph()) {
448 numTotal = A->getLocalNumEntries();
451 GO numLocalBoundaryNodes = 0;
452 GO numGlobalBoundaryNodes = 0;
453 for (
size_t i = 0; i < boundaryNodes.size(); ++i)
454 if (boundaryNodes[i])
455 numLocalBoundaryNodes++;
457 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
458 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
461 Set(currentLevel,
"DofsPerNode", 1);
462 Set(currentLevel,
"Graph", graph);
464 }
else if ((BlockSize == 1 && threshold != STS::zero()) ||
465 (BlockSize == 1 && threshold == STS::zero() && !A->hasCrsGraph()) ||
466 (BlockSize == 1 && useSignedClassicalRS) ||
467 (BlockSize == 1 && useSignedClassicalSA)) {
473 typename LWGraph::row_type::non_const_type
rows(
"rows", A->getLocalNumRows() + 1);
474 typename LWGraph::entries_type::non_const_type columns(
"columns", A->getLocalNumEntries());
476 using MT =
typename STS::magnitudeType;
481 if (useSignedClassicalRS) {
482 if (ghostedBlockNumber.
is_null()) {
485 GetOStream(
Statistics1) <<
"Calculated max point off-diagonal" << std::endl;
489 GetOStream(
Statistics1) <<
"Calculating max block off-diagonal" << std::endl;
493 ghostedDiagVals = ghostedDiag->getData(0);
496 if (rowSumTol > 0.) {
497 if (ghostedBlockNumber.
is_null()) {
499 GetOStream(
Statistics1) <<
"Applying point row sum criterion." << std::endl;
503 GetOStream(
Statistics1) <<
"Applying block row sum criterion." << std::endl;
510 for (
LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); ++row) {
511 size_t nnz = A->getNumEntriesInLocalRow(row);
512 bool rowIsDirichlet = boundaryNodes[row];
515 A->getLocalRowView(row, indices, vals);
517 if (classicalAlgo == defaultAlgo) {
524 if (useSignedClassicalRS) {
526 for (
LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
527 LO col = indices[colID];
528 MT max_neg_aik = realThreshold * STS::real(negMaxOffDiagonal[row]);
529 MT neg_aij = -STS::real(vals[colID]);
534 if ((!rowIsDirichlet && (g_block_id.
is_null() || g_block_id[row] == g_block_id[col]) && neg_aij > max_neg_aik) || row == col) {
535 columns[realnnz++] = col;
540 rows(row + 1) = realnnz;
541 }
else if (useSignedClassicalSA) {
543 for (
LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
544 LO col = indices[colID];
546 bool is_nonpositive = STS::real(vals[colID]) <= 0;
547 MT aiiajj = STS::magnitude(threshold * threshold * ghostedDiagVals[col] * ghostedDiagVals[row]);
548 MT aij = is_nonpositive ? STS::magnitude(vals[colID] * vals[colID]) : (-STS::magnitude(vals[colID] * vals[colID]));
554 if ((!rowIsDirichlet && aij > aiiajj) || row == col) {
555 columns(realnnz++) = col;
560 rows[row + 1] = realnnz;
563 for (
LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
564 LO col = indices[colID];
565 MT aiiajj = STS::magnitude(threshold * threshold * ghostedDiagVals[col] * ghostedDiagVals[row]);
566 MT aij = STS::magnitude(vals[colID] * vals[colID]);
568 if ((!rowIsDirichlet && aij > aiiajj) || row == col) {
569 columns(realnnz++) = col;
574 rows(row + 1) = realnnz;
580 std::vector<DropTol> drop_vec;
581 drop_vec.reserve(nnz);
588 for (
LO colID = 0; colID < (
LO)nnz; colID++) {
589 LO col = indices[colID];
591 drop_vec.emplace_back(zero, one, colID,
false);
596 if (boundaryNodes[colID])
continue;
597 typename STS::magnitudeType aiiajj = STS::magnitude(threshold * threshold * ghostedDiagVals[col] * ghostedDiagVals[row]);
598 typename STS::magnitudeType aij = STS::magnitude(vals[colID] * vals[colID]);
599 drop_vec.emplace_back(aij, aiiajj, colID,
false);
602 const size_t n = drop_vec.size();
604 if (classicalAlgo == unscaled_cut) {
605 std::sort(drop_vec.begin(), drop_vec.end(), [](DropTol
const& a, DropTol
const& b) {
606 return a.val > b.val;
610 for (
size_t i = 1; i < n; ++i) {
612 auto const& x = drop_vec[i - 1];
613 auto const& y = drop_vec[i];
616 if (a > realThreshold * b) {
618 #ifdef HAVE_MUELU_DEBUG
619 if (distanceLaplacianCutVerbose) {
620 std::cout <<
"DJS: KEEP, N, ROW: " << i + 1 <<
", " << n <<
", " << row << std::endl;
625 drop_vec[i].drop = drop;
627 }
else if (classicalAlgo == scaled_cut) {
628 std::sort(drop_vec.begin(), drop_vec.end(), [](DropTol
const& a, DropTol
const& b) {
629 return a.val / a.diag > b.val / b.diag;
634 for (
size_t i = 1; i < n; ++i) {
636 auto const& x = drop_vec[i - 1];
637 auto const& y = drop_vec[i];
638 auto a = x.val / x.diag;
639 auto b = y.val / y.diag;
640 if (a > realThreshold * b) {
643 #ifdef HAVE_MUELU_DEBUG
644 if (distanceLaplacianCutVerbose) {
645 std::cout <<
"DJS: KEEP, N, ROW: " << i + 1 <<
", " << n <<
", " << row << std::endl;
651 drop_vec[i].drop = drop;
655 std::sort(drop_vec.begin(), drop_vec.end(), [](DropTol
const& a, DropTol
const& b) {
656 return a.col < b.col;
659 for (
LO idxID = 0; idxID < (
LO)drop_vec.size(); idxID++) {
660 LO col = indices[drop_vec[idxID].col];
663 columns[realnnz++] = col;
668 if (!drop_vec[idxID].drop) {
669 columns[realnnz++] = col;
676 rows[row + 1] = realnnz;
680 numTotal = A->getLocalNumEntries();
682 if (aggregationMayCreateDirichlet) {
684 for (
LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); ++row) {
686 boundaryNodes[row] =
true;
690 RCP<LWGraph> graph =
rcp(
new LWGraph(
rows, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), A->getRowMap(), A->getColMap(),
"thresholded graph of A"));
693 GO numLocalBoundaryNodes = 0;
694 GO numGlobalBoundaryNodes = 0;
695 for (
size_t i = 0; i < boundaryNodes.size(); ++i)
696 if (boundaryNodes(i))
697 numLocalBoundaryNodes++;
699 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
700 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
702 Set(currentLevel,
"Graph", graph);
703 Set(currentLevel,
"DofsPerNode", 1);
706 if (generateColoringGraph) {
709 BlockDiagonalizeGraph(graph, ghostedBlockNumber, colorGraph, importer);
710 Set(currentLevel,
"Coloring Graph", colorGraph);
730 }
else if (BlockSize > 1 && threshold == STS::zero()) {
735 graphType =
"amalgamated";
743 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
744 Array<LO> colTranslation = *(amalInfo->getColTranslation());
747 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
750 typename LWGraph::row_type::non_const_type
rows(
"rows", numRows + 1);
751 typename LWGraph::entries_type::non_const_type columns(
"columns", A->getLocalNumEntries());
754 Kokkos::deep_copy(amalgBoundaryNodes,
false);
766 LO blkSize = A->GetFixedBlockSize();
768 LO blkPartSize = A->GetFixedBlockSize();
769 if (A->IsView(
"stridedMaps") ==
true) {
773 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
774 blkId = strMap->getStridedBlockId();
776 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
783 for (
LO row = 0; row < numRows; row++) {
793 bool isBoundary =
false;
794 if (pL.get<
bool>(
"aggregation: greedy Dirichlet") ==
true) {
795 for (
LO j = 0; j < blkPartSize; j++) {
796 if (pointBoundaryNodes[row * blkPartSize + j]) {
803 for (
LO j = 0; j < blkPartSize; j++) {
804 if (!pointBoundaryNodes[row * blkPartSize + j]) {
814 MergeRows(*A, row, indicesExtra, colTranslation);
817 indices = indicesExtra;
818 numTotal += indices.
size();
822 LO nnz = indices.
size(), rownnz = 0;
823 for (
LO colID = 0; colID < nnz; colID++) {
824 LO col = indices[colID];
825 columns(realnnz++) = col;
836 amalgBoundaryNodes[row] =
true;
838 rows(row + 1) = realnnz;
841 RCP<LWGraph> graph =
rcp(
new LWGraph(rows, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
845 GO numLocalBoundaryNodes = 0;
846 GO numGlobalBoundaryNodes = 0;
848 for (
size_t i = 0; i < amalgBoundaryNodes.size(); ++i)
849 if (amalgBoundaryNodes(i))
850 numLocalBoundaryNodes++;
853 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
854 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
855 <<
" agglomerated Dirichlet nodes" << std::endl;
858 Set(currentLevel,
"Graph", graph);
859 Set(currentLevel,
"DofsPerNode", blkSize);
861 }
else if (BlockSize > 1 && threshold != STS::zero()) {
865 graphType =
"amalgamated";
873 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
874 Array<LO> colTranslation = *(amalInfo->getColTranslation());
877 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
880 typename LWGraph::row_type::non_const_type
rows(
"rows", numRows + 1);
881 typename LWGraph::entries_type::non_const_type columns(
"columns", A->getLocalNumEntries());
884 Kokkos::deep_copy(amalgBoundaryNodes,
false);
895 LO blkSize = A->GetFixedBlockSize();
897 LO blkPartSize = A->GetFixedBlockSize();
898 if (A->IsView(
"stridedMaps") ==
true) {
902 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
903 blkId = strMap->getStridedBlockId();
905 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
916 for (
LO row = 0; row < numRows; row++) {
926 bool isBoundary =
false;
927 if (pL.get<
bool>(
"aggregation: greedy Dirichlet") ==
true) {
928 for (
LO j = 0; j < blkPartSize; j++) {
929 if (pointBoundaryNodes[row * blkPartSize + j]) {
936 for (
LO j = 0; j < blkPartSize; j++) {
937 if (!pointBoundaryNodes[row * blkPartSize + j]) {
947 MergeRowsWithDropping(*A, row, ghostedDiagVals, threshold, indicesExtra, colTranslation);
950 indices = indicesExtra;
951 numTotal += indices.
size();
955 LO nnz = indices.
size(), rownnz = 0;
956 for (
LO colID = 0; colID < nnz; colID++) {
957 LO col = indices[colID];
958 columns[realnnz++] = col;
969 amalgBoundaryNodes[row] =
true;
971 rows[row + 1] = realnnz;
975 RCP<LWGraph> graph =
rcp(
new LWGraph(rows, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
979 GO numLocalBoundaryNodes = 0;
980 GO numGlobalBoundaryNodes = 0;
982 for (
size_t i = 0; i < amalgBoundaryNodes.size(); ++i)
983 if (amalgBoundaryNodes(i))
984 numLocalBoundaryNodes++;
987 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
988 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
989 <<
" agglomerated Dirichlet nodes" << std::endl;
992 Set(currentLevel,
"Graph", graph);
993 Set(currentLevel,
"DofsPerNode", blkSize);
996 }
else if (algo ==
"distance laplacian") {
997 LO blkSize = A->GetFixedBlockSize();
998 GO indexBase = A->getRowMap()->getIndexBase();
1013 if ((blkSize == 1) && (threshold == STS::zero())) {
1017 graphType =
"unamalgamated";
1018 numTotal = A->getLocalNumEntries();
1021 GO numLocalBoundaryNodes = 0;
1022 GO numGlobalBoundaryNodes = 0;
1023 for (
size_t i = 0; i < pointBoundaryNodes.size(); ++i)
1024 if (pointBoundaryNodes(i))
1025 numLocalBoundaryNodes++;
1027 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1028 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
1031 Set(currentLevel,
"DofsPerNode", blkSize);
1032 Set(currentLevel,
"Graph", graph);
1048 "Coordinate vector length (" << Coords->getLocalLength() <<
") is incompatible with number of rows in A (" << A->getRowMap()->getLocalNumElements() <<
") by modulo block size (" << blkSize <<
").");
1054 uniqueMap = A->getRowMap();
1055 nonUniqueMap = A->getColMap();
1056 graphType =
"unamalgamated";
1059 uniqueMap = Coords->getMap();
1061 "Different index bases for matrix and coordinates");
1065 graphType =
"amalgamated";
1067 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
1072 if (threshold != STS::zero()) {
1077 if (blkSize == 1 && realA->getCrsGraph()->getImporter() != Teuchos::null) {
1078 GetOStream(
Warnings1) <<
"Using existing importer from matrix graph" << std::endl;
1079 importer = realA->getCrsGraph()->getImporter();
1081 GetOStream(
Warnings0) <<
"Constructing new importer instance" << std::endl;
1082 importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
1092 RCP<Vector> localLaplDiag = VectorFactory::Build(uniqueMap);
1095 if (threshold != STS::zero()) {
1096 const size_t numVectors = ghostedCoords->getNumVectors();
1097 coordData.
reserve(numVectors);
1098 for (
size_t j = 0; j < numVectors; j++) {
1105 ArrayRCP<SC> localLaplDiagData = localLaplDiag->getDataNonConst(0);
1106 for (
LO row = 0; row < numRows; row++) {
1111 A->getLocalRowView(row, indices, vals);
1116 MergeRows(*A, row, indicesExtra, colTranslation);
1117 indices = indicesExtra;
1121 bool haveAddedToDiag =
false;
1122 for (
LO colID = 0; colID < nnz; colID++) {
1123 const LO col = indices[colID];
1126 if (use_dlap_weights == SINGLE_WEIGHTS) {
1131 }
else if (use_dlap_weights == BLOCK_WEIGHTS) {
1132 int block_id = row % interleaved_blocksize;
1133 int block_start = block_id * interleaved_blocksize;
1139 haveAddedToDiag =
true;
1144 if (!haveAddedToDiag)
1145 localLaplDiagData[row] = STS::rmax();
1150 ghostedLaplDiag = VectorFactory::Build(nonUniqueMap);
1151 ghostedLaplDiag->doImport(*localLaplDiag, *importer,
Xpetra::INSERT);
1152 ghostedLaplDiagData = ghostedLaplDiag->getDataNonConst(0);
1156 GetOStream(
Runtime0) <<
"Skipping distance laplacian construction due to 0 threshold" << std::endl;
1162 typename LWGraph::row_type::non_const_type
rows(
"rows", numRows + 1);
1163 typename LWGraph::entries_type::non_const_type columns(
"columns", A->getLocalNumEntries());
1165 #ifdef HAVE_MUELU_DEBUG
1167 for (
LO i = 0; i < (
LO)columns.size(); i++) columns[i] = -666;
1172 bool use_stop_array = threshold != STS::zero() && distanceLaplacianAlgo == scaled_cut_symmetric;
1175 rows_stop.
resize(numRows);
1178 Kokkos::deep_copy(amalgBoundaryNodes,
false);
1187 if (threshold != STS::zero()) {
1188 const size_t numVectors = ghostedCoords->getNumVectors();
1189 coordData.
reserve(numVectors);
1190 for (
size_t j = 0; j < numVectors; j++) {
1197 for (
LO row = 0; row < numRows; row++) {
1200 bool isBoundary =
false;
1204 A->getLocalRowView(row, indices, vals);
1205 isBoundary = pointBoundaryNodes[row];
1208 for (
LO j = 0; j < blkSize; j++) {
1209 if (!pointBoundaryNodes[row * blkSize + j]) {
1217 MergeRows(*A, row, indicesExtra, colTranslation);
1220 indices = indicesExtra;
1222 numTotal += indices.
size();
1224 LO nnz = indices.
size(), rownnz = 0;
1226 if (use_stop_array) {
1228 realnnz =
rows(row);
1231 if (threshold != STS::zero()) {
1233 if (distanceLaplacianAlgo == defaultAlgo) {
1235 for (
LO colID = 0; colID < nnz; colID++) {
1236 LO col = indices[colID];
1239 columns(realnnz++) = col;
1245 if (isBoundary)
continue;
1248 if (use_dlap_weights == SINGLE_WEIGHTS) {
1250 }
else if (use_dlap_weights == BLOCK_WEIGHTS) {
1251 int block_id = row % interleaved_blocksize;
1252 int block_start = block_id * interleaved_blocksize;
1257 real_type aiiajj = STS::magnitude(realThreshold * realThreshold * ghostedLaplDiagData[row] * ghostedLaplDiagData[col]);
1258 real_type aij = STS::magnitude(laplVal * laplVal);
1261 columns(realnnz++) = col;
1270 std::vector<DropTol> drop_vec;
1271 drop_vec.reserve(nnz);
1276 for (
LO colID = 0; colID < nnz; colID++) {
1277 LO col = indices[colID];
1280 drop_vec.emplace_back(zero, one, colID,
false);
1284 if (isBoundary)
continue;
1287 if (use_dlap_weights == SINGLE_WEIGHTS) {
1289 }
else if (use_dlap_weights == BLOCK_WEIGHTS) {
1290 int block_id = row % interleaved_blocksize;
1291 int block_start = block_id * interleaved_blocksize;
1297 real_type aiiajj = STS::magnitude(ghostedLaplDiagData[row] * ghostedLaplDiagData[col]);
1298 real_type aij = STS::magnitude(laplVal * laplVal);
1300 drop_vec.emplace_back(aij, aiiajj, colID,
false);
1303 const size_t n = drop_vec.size();
1305 if (distanceLaplacianAlgo == unscaled_cut) {
1306 std::sort(drop_vec.begin(), drop_vec.end(), [](DropTol
const& a, DropTol
const& b) {
1307 return a.val > b.val;
1311 for (
size_t i = 1; i < n; ++i) {
1313 auto const& x = drop_vec[i - 1];
1314 auto const& y = drop_vec[i];
1317 if (a > realThreshold * b) {
1319 #ifdef HAVE_MUELU_DEBUG
1320 if (distanceLaplacianCutVerbose) {
1321 std::cout <<
"DJS: KEEP, N, ROW: " << i + 1 <<
", " << n <<
", " << row << std::endl;
1326 drop_vec[i].drop = drop;
1328 }
else if (distanceLaplacianAlgo == scaled_cut || distanceLaplacianAlgo == scaled_cut_symmetric) {
1329 std::sort(drop_vec.begin(), drop_vec.end(), [](DropTol
const& a, DropTol
const& b) {
1330 return a.val / a.diag > b.val / b.diag;
1334 for (
size_t i = 1; i < n; ++i) {
1336 auto const& x = drop_vec[i - 1];
1337 auto const& y = drop_vec[i];
1338 auto a = x.val / x.diag;
1339 auto b = y.val / y.diag;
1340 if (a > realThreshold * b) {
1342 #ifdef HAVE_MUELU_DEBUG
1343 if (distanceLaplacianCutVerbose) {
1344 std::cout <<
"DJS: KEEP, N, ROW: " << i + 1 <<
", " << n <<
", " << row << std::endl;
1349 drop_vec[i].drop = drop;
1353 std::sort(drop_vec.begin(), drop_vec.end(), [](DropTol
const& a, DropTol
const& b) {
1354 return a.col < b.col;
1357 for (
LO idxID = 0; idxID < (
LO)drop_vec.size(); idxID++) {
1358 LO col = indices[drop_vec[idxID].col];
1362 columns(realnnz++) = col;
1368 if (!drop_vec[idxID].drop) {
1369 columns(realnnz++) = col;
1380 for (
LO colID = 0; colID < nnz; colID++) {
1381 LO col = indices[colID];
1382 columns(realnnz++) = col;
1394 amalgBoundaryNodes[row] =
true;
1398 rows_stop[row] = rownnz + rows[row];
1400 rows[row + 1] = realnnz;
1405 if (use_stop_array) {
1408 for (
LO row = 0; row < numRows; row++) {
1409 for (
LO colidx = rows[row]; colidx < rows_stop[row]; colidx++) {
1410 LO col = columns[colidx];
1411 if (col >= numRows)
continue;
1414 for (
LO t_col =
rows(col); !found && t_col < rows_stop[col]; t_col++) {
1415 if (columns[t_col] == row)
1420 if (!found && !pointBoundaryNodes[col] && Teuchos::as<typename LWGraph::row_type::value_type>(rows_stop[col]) < rows[col + 1]) {
1421 LO new_idx = rows_stop[col];
1423 columns[new_idx] = row;
1431 LO current_start = 0;
1432 for (
LO row = 0; row < numRows; row++) {
1433 LO old_start = current_start;
1434 for (
LO col =
rows(row); col < rows_stop[row]; col++) {
1435 if (current_start != col) {
1436 columns(current_start) = columns(col);
1440 rows[row] = old_start;
1442 rows(numRows) = realnnz = current_start;
1448 graph =
rcp(
new LWGraph(rows, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
1453 GO numLocalBoundaryNodes = 0;
1454 GO numGlobalBoundaryNodes = 0;
1456 for (
size_t i = 0; i < amalgBoundaryNodes.size(); ++i)
1457 if (amalgBoundaryNodes(i))
1458 numLocalBoundaryNodes++;
1461 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1462 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" agglomerated Dirichlet nodes"
1463 <<
" using threshold " << dirichletThreshold << std::endl;
1466 Set(currentLevel,
"Graph", graph);
1467 Set(currentLevel,
"DofsPerNode", blkSize);
1471 if ((GetVerbLevel() &
Statistics1) && !(A->GetFixedBlockSize() > 1 && threshold != STS::zero())) {
1473 GO numGlobalTotal, numGlobalDropped;
1476 GetOStream(
Statistics1) <<
"Number of dropped entries in " << graphType <<
" matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1477 if (numGlobalTotal != 0)
1478 GetOStream(
Statistics1) <<
" (" << 100 * Teuchos::as<double>(numGlobalDropped) / Teuchos::as<double>(numGlobalTotal) <<
"%)";
1485 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
1487 GetOStream(
Runtime0) <<
"algorithm = \""
1489 <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
1490 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
1496 GO indexBase = rowMap->getIndexBase();
1500 if (A->IsView(
"stridedMaps") &&
1501 Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1505 blockdim = strMap->getFixedBlockSize();
1506 offset = strMap->getOffset();
1507 oldView = A->SwitchToView(oldView);
1508 GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build():"
1509 <<
" found blockdim=" << blockdim <<
" from strided maps. offset=" << offset << std::endl;
1511 GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
1516 GetOStream(
Statistics1) <<
"CoalesceDropFactory: nodeMap " << nodeMap->getLocalNumElements() <<
"/" << nodeMap->getGlobalNumElements() <<
" elements" << std::endl;
1519 RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, A->getLocalMaxNumRowEntries() * blockdim);
1521 LO numRows = A->getRowMap()->getLocalNumElements();
1522 LO numNodes = nodeMap->getLocalNumElements();
1524 Kokkos::deep_copy(amalgBoundaryNodes,
false);
1525 const ArrayRCP<int> numberDirichletRowsPerNode(numNodes, 0);
1526 bool bIsDiagonalEntry =
false;
1531 for (
LO row = 0; row < numRows; row++) {
1533 GO grid = rowMap->getGlobalElement(row);
1536 bIsDiagonalEntry =
false;
1541 size_t nnz = A->getNumEntriesInLocalRow(row);
1544 A->getLocalRowView(row, indices, vals);
1548 for (
LO col = 0; col < Teuchos::as<LO>(nnz); col++) {
1549 GO gcid = colMap->getGlobalElement(indices[col]);
1551 if (vals[col] != STS::zero()) {
1553 cnodeIds->push_back(cnodeId);
1555 if (grid == gcid) bIsDiagonalEntry =
true;
1559 if (realnnz == 1 && bIsDiagonalEntry ==
true) {
1560 LO lNodeId = nodeMap->getLocalElement(nodeId);
1561 numberDirichletRowsPerNode[lNodeId] += 1;
1562 if (numberDirichletRowsPerNode[lNodeId] == blockdim)
1563 amalgBoundaryNodes[lNodeId] =
true;
1568 if (arr_cnodeIds.
size() > 0)
1569 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
1572 crsGraph->fillComplete(nodeMap, nodeMap);
1581 GO numLocalBoundaryNodes = 0;
1582 GO numGlobalBoundaryNodes = 0;
1583 for (
size_t i = 0; i < amalgBoundaryNodes.size(); ++i)
1584 if (amalgBoundaryNodes(i))
1585 numLocalBoundaryNodes++;
1587 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1588 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
1593 Set(currentLevel,
"DofsPerNode", blockdim);
1594 Set(currentLevel,
"Graph", graph);
1600 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1605 LO blkSize = A.GetFixedBlockSize();
1606 if (A.IsView(
"stridedMaps") ==
true) {
1610 if (strMap->getStridedBlockId() > -1)
1611 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1615 size_t nnz = 0, pos = 0;
1616 for (
LO j = 0; j < blkSize; j++)
1617 nnz += A.getNumEntriesInLocalRow(row * blkSize + j);
1629 for (
LO j = 0; j < blkSize; j++) {
1630 A.getLocalRowView(row * blkSize + j, inds, vals);
1631 size_type numIndices = inds.
size();
1633 if (numIndices == 0)
1637 cols[pos++] = translation[inds[0]];
1638 for (size_type k = 1; k < numIndices; k++) {
1639 LO nodeID = translation[inds[k]];
1643 if (nodeID != cols[pos - 1])
1644 cols[pos++] = nodeID;
1651 std::sort(cols.
begin(), cols.
end());
1653 for (
size_t j = 1; j < nnz; j++)
1654 if (cols[j] != cols[pos])
1655 cols[++pos] = cols[j];
1659 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1665 LO blkSize = A.GetFixedBlockSize();
1666 if (A.IsView(
"stridedMaps") ==
true) {
1670 if (strMap->getStridedBlockId() > -1)
1671 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1675 size_t nnz = 0, pos = 0;
1676 for (
LO j = 0; j < blkSize; j++)
1677 nnz += A.getNumEntriesInLocalRow(row * blkSize + j);
1689 for (
LO j = 0; j < blkSize; j++) {
1690 A.getLocalRowView(row * blkSize + j, inds, vals);
1691 size_type numIndices = inds.
size();
1693 if (numIndices == 0)
1698 for (size_type k = 0; k < numIndices; k++) {
1700 LO nodeID = translation[inds[k]];
1703 typename STS::magnitudeType aiiajj = STS::magnitude(threshold * threshold * ghostedDiagVals[dofID] * ghostedDiagVals[row * blkSize + j]);
1704 typename STS::magnitudeType aij = STS::magnitude(vals[k] * vals[k]);
1707 if (aij > aiiajj || (row * blkSize + j == dofID)) {
1713 if (nodeID != prevNodeID) {
1714 cols[pos++] = nodeID;
1715 prevNodeID = nodeID;
1724 std::sort(cols.
begin(), cols.
end());
1726 for (
size_t j = 1; j < nnz; j++)
1727 if (cols[j] != cols[pos])
1728 cols[++pos] = cols[j];
1734 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1739 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.
get<
double>(
"aggregation: Dirichlet threshold")));
1740 const typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.
get<
double>(
"aggregation: row sum drop tol"));
1744 GetOStream(
Statistics1) <<
"Using BlockDiagonal Graph before dropping (with provided blocking)" << std::endl;
1751 ghostedBlockNumber->doImport(*BlockNumber, *importer,
Xpetra::INSERT);
1753 ghostedBlockNumber = BlockNumber;
1761 typename CrsMatrix::local_matrix_type::row_map_type::HostMirror::non_const_type rows_mat;
1762 typename LWGraph::row_type::non_const_type rows_graph;
1763 typename LWGraph::entries_type::non_const_type columns;
1764 typename CrsMatrix::local_matrix_type::values_type::HostMirror::non_const_type values;
1767 if (generate_matrix) {
1768 crs_matrix_wrap =
rcp(
new CrsMatrixWrap(A->getRowMap(), A->getColMap(), 0));
1769 rows_mat =
typename CrsMatrix::local_matrix_type::row_map_type::HostMirror::non_const_type(
"rows_mat", A->getLocalNumRows() + 1);
1771 rows_graph =
typename LWGraph::row_type::non_const_type(
"rows_graph", A->getLocalNumRows() + 1);
1773 columns =
typename LWGraph::entries_type::non_const_type(
"columns", A->getLocalNumEntries());
1774 values =
typename CrsMatrix::local_matrix_type::values_type::HostMirror::non_const_type(
"values", A->getLocalNumEntries());
1777 GO numDropped = 0, numTotal = 0;
1778 for (
LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); ++row) {
1779 LO row_block = row_block_number[row];
1780 size_t nnz = A->getNumEntriesInLocalRow(row);
1783 A->getLocalRowView(row, indices, vals);
1786 for (
LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
1787 LO col = indices[colID];
1788 LO col_block = col_block_number[col];
1790 if (row_block == col_block) {
1791 if (generate_matrix) values[realnnz] = vals[colID];
1792 columns[realnnz++] = col;
1797 if (generate_matrix)
1798 rows_mat[row + 1] = realnnz;
1800 rows_graph[row + 1] = realnnz;
1807 numTotal = A->getLocalNumEntries();
1810 GO numLocalBoundaryNodes = 0;
1811 GO numGlobalBoundaryNodes = 0;
1812 for (
size_t i = 0; i < boundaryNodes.size(); ++i)
1813 if (boundaryNodes(i))
1814 numLocalBoundaryNodes++;
1816 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1817 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
1819 GO numGlobalTotal, numGlobalDropped;
1822 GetOStream(
Statistics1) <<
"Number of dropped entries in block-diagonalized matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1823 if (numGlobalTotal != 0)
1824 GetOStream(
Statistics1) <<
" (" << 100 * Teuchos::as<double>(numGlobalDropped) / Teuchos::as<double>(numGlobalTotal) <<
"%)";
1828 Set(currentLevel,
"Filtering",
true);
1830 if (generate_matrix) {
1835 typename CrsMatrix::local_matrix_type::row_map_type>::value) {
1836 crs_matrix_wrap->getCrsMatrix()->setAllValues(rows_mat, columns, values);
1838 auto rows_mat2 =
typename CrsMatrix::local_matrix_type::row_map_type::non_const_type(
"rows_mat2", rows_mat.extent(0));
1839 Kokkos::deep_copy(rows_mat2, rows_mat);
1840 auto columns2 =
typename CrsMatrix::local_graph_type::entries_type::non_const_type(
"columns2", columns.extent(0));
1841 Kokkos::deep_copy(columns2, columns);
1842 auto values2 =
typename CrsMatrix::local_matrix_type::values_type::non_const_type(
"values2", values.extent(0));
1843 Kokkos::deep_copy(values2, values);
1844 crs_matrix_wrap->getCrsMatrix()->setAllValues(rows_mat2, columns2, values2);
1846 crs_matrix_wrap->getCrsMatrix()->expertStaticFillComplete(A->getColMap(), A->getRowMap());
1848 RCP<LWGraph> graph =
rcp(
new LWGraph(rows_graph, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), A->getRowMap(), A->getColMap(),
"block-diagonalized graph of A"));
1850 Set(currentLevel,
"Graph", graph);
1853 Set(currentLevel,
"DofsPerNode", 1);
1854 return crs_matrix_wrap;
1857 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1862 const bool localizeColoringGraph = pL.get<
bool>(
"aggregation: coloring: localize color graph");
1864 GetOStream(
Statistics1) <<
"Using BlockDiagonal Graph after Dropping (with provided blocking)";
1865 if (localizeColoringGraph)
1866 GetOStream(
Statistics1) <<
", with localization" << std::endl;
1868 GetOStream(
Statistics1) <<
", without localization" << std::endl;
1876 typename LWGraph::row_type::non_const_type rows_graph(
"rows_graph", inputGraph->
GetNodeNumVertices() + 1);
1877 typename LWGraph::entries_type::non_const_type columns(
"columns", inputGraph->
GetNodeNumEdges());
1880 GO numDropped = 0, numTotal = 0;
1881 const LO numRows = Teuchos::as<LO>(inputGraph->
GetDomainMap()->getLocalNumElements());
1882 if (localizeColoringGraph) {
1883 for (
LO row = 0; row < numRows; ++row) {
1884 LO row_block = row_block_number[row];
1888 for (
LO colID = 0; colID < Teuchos::as<LO>(indices.length); colID++) {
1889 LO col = indices(colID);
1890 LO col_block = col_block_number[col];
1892 if ((row_block == col_block) && (col < numRows)) {
1893 columns(realnnz++) = col;
1898 rows_graph(row + 1) = realnnz;
1905 boundaryNodesVector->getDataNonConst(0)[i] = boundaryNodes[i];
1908 boundaryColumnVector->doImport(*boundaryNodesVector, *importer,
Xpetra::INSERT);
1909 auto boundaryColumn = boundaryColumnVector->getData(0);
1911 for (
LO row = 0; row < numRows; ++row) {
1912 LO row_block = row_block_number[row];
1916 for (
LO colID = 0; colID < Teuchos::as<LO>(indices.length); colID++) {
1917 LO col = indices(colID);
1918 LO col_block = col_block_number[col];
1920 if ((row_block == col_block) && ((row == col) || (boundaryColumn[col] == 0))) {
1921 columns(realnnz++) = col;
1926 rows_graph(row + 1) = realnnz;
1934 GO numGlobalTotal, numGlobalDropped;
1937 GetOStream(
Statistics1) <<
"Number of dropped entries in block-diagonalized matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1938 if (numGlobalTotal != 0)
1939 GetOStream(
Statistics1) <<
" (" << 100 * Teuchos::as<double>(numGlobalDropped) / Teuchos::as<double>(numGlobalTotal) <<
"%)";
1943 if (localizeColoringGraph) {
1944 outputGraph =
rcp(
new LWGraph(rows_graph, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), inputGraph->
GetDomainMap(), inputGraph->
GetImportMap(),
"block-diagonalized graph of A"));
1948 #ifdef HAVE_XPETRA_TPETRA
1949 auto outputGraph2 =
rcp(
new LWGraph(rows_graph, Kokkos::subview(columns, Kokkos::make_pair(0, realnnz)), inputGraph->
GetDomainMap(), inputGraph->
GetImportMap(),
"block-diagonalized graph of A"));
1951 auto tpGraph =
Xpetra::toTpetra(rcp_const_cast<const CrsGraph>(outputGraph2->GetCrsGraph()));
1953 auto tpGraphSym = sym->symmetrize();
1954 auto lclGraphSym = tpGraphSym->getLocalGraphHost();
1955 auto colIndsSym = lclGraphSym.entries;
1957 auto rowsSym = tpGraphSym->getLocalRowPtrsHost();
1958 typename LWGraph::row_type::non_const_type rows_graphSym(
"rows_graphSym", rowsSym.size());
1959 for (
size_t row = 0; row < rowsSym.size(); row++)
1960 rows_graphSym(row) = rowsSym(row);
1969 #endif // MUELU_COALESCEDROPFACTORY_DEF_HPP
Important warning messages (one line)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
typename local_graph_type::row_map_type row_type
Exception indicating invalid cast attempted.
void reserve(size_type n)
void BlockDiagonalizeGraph(const RCP< LWGraph > &inputGraph, const RCP< LocalOrdinalVector > &ghostedBlockNumber, RCP< LWGraph > &outputGraph, RCP< const Import > &importer) const
#define MueLu_sumAll(rcpComm, in, out)
static void ApplyRowSumCriterionHost(const Matrix &A, const typename Teuchos::ScalarTraits< Scalar >::magnitudeType rowSumTol, Kokkos::View< bool *, Kokkos::HostSpace > &dirichletRows)
KOKKOS_INLINE_FUNCTION void SetBoundaryNodeMap(const boundary_nodes_type bndry)
Set boolean array indicating which rows correspond to Dirichlet boundaries.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
void setValidator(RCP< const ParameterEntryValidator > const &validator)
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
Apply Rowsum Criterion.
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static void Write(const std::string &fileName, const Xpetra::Map< LocalOrdinal, GlobalOrdinal, Node > &M)
DropTol & operator=(DropTol const &)=default
Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > BlockDiagonalize(Level ¤tLevel, const RCP< Matrix > &A, bool generate_matrix) const
One-liner description of what is happening.
KOKKOS_INLINE_FUNCTION size_type GetNodeNumVertices() const
Return number of graph vertices.
KOKKOS_INLINE_FUNCTION size_type GetNodeNumEdges() const
Return number of graph edges.
KOKKOS_INLINE_FUNCTION const boundary_nodes_type GetBoundaryNodeMap() const
Returns map with global ids of boundary nodes.
Exception throws to report incompatible objects (like maps).
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Return vector containing: max_{i=k}(-a_ik), for each for i in the matrix.
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
Translate global (row/column) id to global amalgamation block id.
Scalar GetThreshold() const
Return threshold value.
void resize(const size_type n, const T &val=T())
CoalesceDropFactory()
Constructor.
Kokkos::View< bool *, memory_space > boundary_nodes_type
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Magnitude >::zero(), bool count_twos_as_dirichlet=false)
Detect Dirichlet rows.
#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.
static RCP< Vector > GetMatrixOverlappedDiagonal(const Matrix &A)
Extract Overlapped Matrix Diagonal.
void resize(size_type new_size, const value_type &x=value_type())
void Build(Level ¤tLevel) const
Build an object with this factory.
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
void DeclareInput(Level ¤tLevel) const
Input.
KOKKOS_INLINE_FUNCTION neighbor_vertices_type getNeighborVertices(LO i) const
Return the list of vertices adjacent to the vertex 'v'.
const RCP< const Map > GetImportMap() const
Return overlapping import map (nodes).
RCP< const Tpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > toTpetra(const RCP< const CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > &graph)
void push_back(const value_type &x)
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.
RCP< const CrsGraph< int, GlobalOrdinal, Node > > toXpetra(const Epetra_CrsGraph &g)
const RCP< const Map > GetDomainMap() const
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.
#define TEUCHOS_ASSERT(assertion_test)
static Kokkos::View< bool *, typename Kokkos::HostSpace > DetectDirichletRows_kokkos_host(const Matrix &A, const Magnitude &tol=Teuchos::ScalarTraits< typename Teuchos::ScalarTraits< SC >::magnitudeType >::zero(), const bool count_twos_as_dirichlet=false)
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar >> &v, LocalOrdinal i0, LocalOrdinal i1)
Squared distance between two rows in a multivector.
ParameterEntry & getEntry(const std::string &name)
void MergeRowsWithDropping(const Matrix &A, const LO row, const ArrayRCP< const SC > &ghostedDiagVals, SC threshold, Array< LO > &cols, const Array< LO > &translation) const