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");
135 typedef Teuchos::StringToIntegralParameterEntryValidator<int> validatorType;
137 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")));
143#undef SET_VALID_ENTRY
144 validParamList->set<
bool > (
"lightweight wrap",
true,
"Experimental option for lightweight graph access");
146 validParamList->set< RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A");
147 validParamList->set< RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory for UnAmalgamationInfo");
148 validParamList->set< RCP<const FactoryBase> >(
"Coordinates", Teuchos::null,
"Generating factory for Coordinates");
149 validParamList->set< RCP<const FactoryBase> >(
"BlockNumber", Teuchos::null,
"Generating factory for BlockNUmber");
151 return validParamList;
154 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
157 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
159 Input(currentLevel,
"A");
160 Input(currentLevel,
"UnAmalgamationInfo");
162 const ParameterList& pL = GetParameterList();
163 if (pL.get<
bool>(
"lightweight wrap") ==
true) {
164 std::string algo = pL.get<std::string>(
"aggregation: drop scheme");
165 if (algo ==
"distance laplacian" || algo ==
"block diagonal distance laplacian") {
166 Input(currentLevel,
"Coordinates");
168 if(algo ==
"signed classical sa")
170 else if (algo.find(
"block diagonal") != std::string::npos || algo.find(
"signed classical") != std::string::npos) {
171 Input(currentLevel,
"BlockNumber");
177 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
182 typedef Teuchos::ScalarTraits<SC> STS;
183 typedef typename STS::magnitudeType real_type;
184 typedef Xpetra::MultiVector<real_type,LO,GO,NO> RealValuedMultiVector;
185 typedef Xpetra::MultiVectorFactory<real_type,LO,GO,NO> RealValuedMultiVectorFactory;
187 if (predrop_ != Teuchos::null)
188 GetOStream(
Parameters0) << predrop_->description();
190 RCP<Matrix> realA = Get< RCP<Matrix> >(currentLevel,
"A");
191 RCP<AmalgamationInfo> amalInfo = Get< RCP<AmalgamationInfo> >(currentLevel,
"UnAmalgamationInfo");
192 const ParameterList & pL = GetParameterList();
193 bool doExperimentalWrap = pL.get<
bool>(
"lightweight wrap");
195 GetOStream(
Parameters0) <<
"lightweight wrap = " << doExperimentalWrap << std::endl;
196 std::string algo = pL.get<std::string>(
"aggregation: drop scheme");
197 const bool aggregationMayCreateDirichlet = pL.get<
bool>(
"aggregation: dropping may create Dirichlet");
199 RCP<RealValuedMultiVector> Coords;
202 bool use_block_algorithm=
false;
203 LO interleaved_blocksize = as<LO>(pL.get<
int>(
"aggregation: block diagonal: interleaved blocksize"));
204 bool useSignedClassicalRS =
false;
205 bool useSignedClassicalSA =
false;
206 bool generateColoringGraph =
false;
210 typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<
double>(
"aggregation: row sum drop tol"));
211 RCP<LocalOrdinalVector> ghostedBlockNumber;
212 ArrayRCP<const LO> g_block_id;
214 if(algo ==
"distance laplacian" ) {
216 Coords = Get< RCP<RealValuedMultiVector > >(currentLevel,
"Coordinates");
219 else if(algo ==
"signed classical sa") {
220 useSignedClassicalSA =
true;
224 else if(algo ==
"signed classical" || algo ==
"block diagonal colored signed classical" || algo ==
"block diagonal signed classical") {
225 useSignedClassicalRS =
true;
227 RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(currentLevel,
"BlockNumber");
229 RCP<const Import> importer = realA->getCrsGraph()->getImporter();
230 if(!importer.is_null()) {
232 ghostedBlockNumber= Xpetra::VectorFactory<LO,LO,GO,NO>::Build(importer->getTargetMap());
233 ghostedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
236 ghostedBlockNumber = BlockNumber;
238 g_block_id = ghostedBlockNumber->getData(0);
240 if(algo ==
"block diagonal colored signed classical")
241 generateColoringGraph=
true;
246 else if(algo ==
"block diagonal") {
248 BlockDiagonalize(currentLevel,realA,
false);
251 else if (algo ==
"block diagonal classical" || algo ==
"block diagonal distance laplacian") {
253 use_block_algorithm =
true;
254 RCP<Matrix> filteredMatrix = BlockDiagonalize(currentLevel,realA,
true);
255 if(algo ==
"block diagonal distance laplacian") {
257 RCP<RealValuedMultiVector> OldCoords = Get< RCP<RealValuedMultiVector > >(currentLevel,
"Coordinates");
258 if (OldCoords->getLocalLength() != realA->getLocalNumRows()) {
259 LO dim = (LO) OldCoords->getNumVectors();
260 Coords = RealValuedMultiVectorFactory::Build(realA->getRowMap(),dim);
261 for(LO k=0; k<dim; k++){
262 ArrayRCP<const real_type> old_vec = OldCoords->getData(k);
263 ArrayRCP<real_type> new_vec = Coords->getDataNonConst(k);
264 for(LO i=0; i <(LO)OldCoords->getLocalLength(); i++) {
266 for(LO j=0; j<interleaved_blocksize; j++)
267 new_vec[new_base + j] = old_vec[i];
274 algo =
"distance laplacian";
276 else if(algo ==
"block diagonal classical") {
288 Array<double> dlap_weights = pL.get<Array<double> >(
"aggregation: distance laplacian directional weights");
289 enum {NO_WEIGHTS=0, SINGLE_WEIGHTS, BLOCK_WEIGHTS};
290 int use_dlap_weights = NO_WEIGHTS;
291 if(algo ==
"distance laplacian") {
292 LO dim = (LO) Coords->getNumVectors();
294 bool non_unity =
false;
295 for (LO i=0; !non_unity && i<(LO)dlap_weights.size(); i++) {
296 if(dlap_weights[i] != 1.0) {
301 LO blocksize = use_block_algorithm ? as<LO>(pL.get<
int>(
"aggregation: block diagonal: interleaved blocksize")) : 1;
302 if((LO)dlap_weights.size() == dim)
303 use_dlap_weights = SINGLE_WEIGHTS;
304 else if((LO)dlap_weights.size() == blocksize * dim)
305 use_dlap_weights = BLOCK_WEIGHTS;
308 "length of 'aggregation: distance laplacian directional weights' must equal the coordinate dimension OR the coordinate dimension times the blocksize");
311 GetOStream(
Statistics1) <<
"Using distance laplacian weights: "<<dlap_weights<<std::endl;
326 if (doExperimentalWrap) {
327 TEUCHOS_TEST_FOR_EXCEPTION(predrop_ != null && algo !=
"classical",
Exceptions::RuntimeError,
"Dropping function must not be provided for \"" << algo <<
"\" algorithm");
328 TEUCHOS_TEST_FOR_EXCEPTION(algo !=
"classical" && algo !=
"distance laplacian" && algo !=
"signed classical",
Exceptions::RuntimeError,
"\"algorithm\" must be one of (classical|distance laplacian|signed classical)");
330 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
331 std::string distanceLaplacianAlgoStr = pL.get<std::string>(
"aggregation: distance laplacian algo");
332 std::string classicalAlgoStr = pL.get<std::string>(
"aggregation: classical algo");
333 real_type realThreshold = STS::magnitude(threshold);
337#ifdef HAVE_MUELU_DEBUG
338 int distanceLaplacianCutVerbose = 0;
340#ifdef DJS_READ_ENV_VARIABLES
341 if (getenv(
"MUELU_DROP_TOLERANCE_MODE")) {
342 distanceLaplacianAlgoStr = std::string(getenv(
"MUELU_DROP_TOLERANCE_MODE"));
345 if (getenv(
"MUELU_DROP_TOLERANCE_THRESHOLD")) {
346 auto tmp = atoi(getenv(
"MUELU_DROP_TOLERANCE_THRESHOLD"));
347 realThreshold = 1e-4*tmp;
350# ifdef HAVE_MUELU_DEBUG
351 if (getenv(
"MUELU_DROP_TOLERANCE_VERBOSE")) {
352 distanceLaplacianCutVerbose = atoi(getenv(
"MUELU_DROP_TOLERANCE_VERBOSE"));
358 enum decisionAlgoType {defaultAlgo, unscaled_cut, scaled_cut, scaled_cut_symmetric};
360 decisionAlgoType distanceLaplacianAlgo = defaultAlgo;
361 decisionAlgoType classicalAlgo = defaultAlgo;
362 if (algo ==
"distance laplacian") {
363 if (distanceLaplacianAlgoStr ==
"default")
364 distanceLaplacianAlgo = defaultAlgo;
365 else if (distanceLaplacianAlgoStr ==
"unscaled cut")
366 distanceLaplacianAlgo = unscaled_cut;
367 else if (distanceLaplacianAlgoStr ==
"scaled cut")
368 distanceLaplacianAlgo = scaled_cut;
369 else if (distanceLaplacianAlgoStr ==
"scaled cut symmetric")
370 distanceLaplacianAlgo = scaled_cut_symmetric;
372 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"\"aggregation: distance laplacian algo\" must be one of (default|unscaled cut|scaled cut), not \"" << distanceLaplacianAlgoStr <<
"\"");
373 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\" distance laplacian algorithm = \"" << distanceLaplacianAlgoStr <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize()<< std::endl;
374 }
else if (algo ==
"classical") {
375 if (classicalAlgoStr ==
"default")
376 classicalAlgo = defaultAlgo;
377 else if (classicalAlgoStr ==
"unscaled cut")
378 classicalAlgo = unscaled_cut;
379 else if (classicalAlgoStr ==
"scaled cut")
380 classicalAlgo = scaled_cut;
382 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"\"aggregation: classical algo\" must be one of (default|unscaled cut|scaled cut), not \"" << classicalAlgoStr <<
"\"");
383 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\" classical algorithm = \"" << classicalAlgoStr <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
386 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
387 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
389 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<
double>(
"aggregation: Dirichlet threshold")));
393 TEUCHOS_TEST_FOR_EXCEPTION(useSignedClassicalRS && classicalAlgo != defaultAlgo,
Exceptions::RuntimeError,
"\"aggregation: classical algo\" != default is not supported for scalled classical aggregation");
394 TEUCHOS_TEST_FOR_EXCEPTION(useSignedClassicalSA && classicalAlgo != defaultAlgo,
Exceptions::RuntimeError,
"\"aggregation: classical algo\" != default is not supported for scalled classical sa aggregation");
396 GO numDropped = 0, numTotal = 0;
397 std::string graphType =
"unamalgamated";
416 TEUCHOS_TEST_FOR_EXCEPTION(A->GetFixedBlockSize() % A->GetStorageBlockSize() != 0,
Exceptions::RuntimeError,
"A->GetFixedBlockSize() needs to be a multiple of A->GetStorageBlockSize()");
417 const LO BlockSize = A->GetFixedBlockSize() / A->GetStorageBlockSize();
421 if (algo ==
"classical") {
422 if (predrop_ == null) {
427 if (predrop_ != null) {
428 RCP<PreDropFunctionConstVal> predropConstVal = rcp_dynamic_cast<PreDropFunctionConstVal>(predrop_);
430 "MueLu::CoalesceFactory::Build: cast to PreDropFunctionConstVal failed.");
432 SC newt = predropConstVal->GetThreshold();
433 if (newt != threshold) {
434 GetOStream(
Warnings0) <<
"switching threshold parameter from " << threshold <<
" (list) to " << newt <<
" (user function" << std::endl;
441 if ( BlockSize==1 && threshold == STS::zero() && !useSignedClassicalRS && !useSignedClassicalSA && A->hasCrsGraph()) {
443 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
449 graph->SetBoundaryNodeMap(boundaryNodes);
450 numTotal = A->getLocalNumEntries();
453 GO numLocalBoundaryNodes = 0;
454 GO numGlobalBoundaryNodes = 0;
455 for (LO i = 0; i < boundaryNodes.size(); ++i)
456 if (boundaryNodes[i])
457 numLocalBoundaryNodes++;
458 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
459 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
460 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
463 Set(currentLevel,
"DofsPerNode", 1);
464 Set(currentLevel,
"Graph", graph);
466 }
else if ( (BlockSize == 1 && threshold != STS::zero()) ||
467 (BlockSize == 1 && threshold == STS::zero() && !A->hasCrsGraph()) ||
468 (BlockSize == 1 && useSignedClassicalRS) ||
469 (BlockSize == 1 && useSignedClassicalSA) ) {
475 ArrayRCP<LO>
rows (A->getLocalNumRows()+1);
476 ArrayRCP<LO> columns(A->getLocalNumEntries());
478 using MT =
typename STS::magnitudeType;
479 RCP<Vector> ghostedDiag;
480 ArrayRCP<const SC> ghostedDiagVals;
481 ArrayRCP<const MT> negMaxOffDiagonal;
483 if(useSignedClassicalRS) {
484 if(ghostedBlockNumber.is_null()) {
487 GetOStream(
Statistics1) <<
"Calculated max point off-diagonal" << std::endl;
492 GetOStream(
Statistics1) <<
"Calculating max block off-diagonal" << std::endl;
497 ghostedDiagVals = ghostedDiag->getData(0);
500 if (rowSumTol > 0.) {
501 if(ghostedBlockNumber.is_null()) {
503 GetOStream(
Statistics1) <<
"Applying point row sum criterion." << std::endl;
507 GetOStream(
Statistics1) <<
"Applying block row sum criterion." << std::endl;
514 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); ++row) {
515 size_t nnz = A->getNumEntriesInLocalRow(row);
516 bool rowIsDirichlet = boundaryNodes[row];
517 ArrayView<const LO> indices;
518 ArrayView<const SC> vals;
519 A->getLocalRowView(row, indices, vals);
521 if(classicalAlgo == defaultAlgo) {
528 if(useSignedClassicalRS) {
530 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
531 LO col = indices[colID];
532 MT max_neg_aik = realThreshold * STS::real(negMaxOffDiagonal[row]);
533 MT neg_aij = - STS::real(vals[colID]);
538 if ((!rowIsDirichlet && (g_block_id.is_null() || g_block_id[row] == g_block_id[col]) && neg_aij > max_neg_aik) || row == col) {
539 columns[realnnz++] = col;
544 rows[row+1] = realnnz;
546 else if(useSignedClassicalSA) {
548 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
549 LO col = indices[colID];
551 bool is_nonpositive = STS::real(vals[colID]) <= 0;
552 MT aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
553 MT aij = is_nonpositive ? STS::magnitude(vals[colID]*vals[colID]) : (-STS::magnitude(vals[colID]*vals[colID]));
559 if ((!rowIsDirichlet && aij > aiiajj) || row == col) {
560 columns[realnnz++] = col;
565 rows[row+1] = realnnz;
569 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
570 LO col = indices[colID];
571 MT aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
572 MT aij = STS::magnitude(vals[colID]*vals[colID]);
574 if ((!rowIsDirichlet && aij > aiiajj) || row == col) {
575 columns[realnnz++] = col;
580 rows[row+1] = realnnz;
587 std::vector<DropTol> drop_vec;
588 drop_vec.reserve(nnz);
589 const real_type zero = Teuchos::ScalarTraits<real_type>::zero();
590 const real_type one = Teuchos::ScalarTraits<real_type>::one();
595 for (LO colID = 0; colID < (LO)nnz; colID++) {
596 LO col = indices[colID];
598 drop_vec.emplace_back( zero, one, colID,
false);
603 if(boundaryNodes[colID])
continue;
604 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
605 typename STS::magnitudeType aij = STS::magnitude(vals[colID]*vals[colID]);
606 drop_vec.emplace_back(aij, aiiajj, colID,
false);
609 const size_t n = drop_vec.size();
611 if (classicalAlgo == unscaled_cut) {
612 std::sort( drop_vec.begin(), drop_vec.end()
613 , [](DropTol
const& a, DropTol
const& b) {
614 return a.val > b.val;
618 for (
size_t i=1; i<n; ++i) {
620 auto const& x = drop_vec[i-1];
621 auto const& y = drop_vec[i];
624 if (a > realThreshold*b) {
626#ifdef HAVE_MUELU_DEBUG
627 if (distanceLaplacianCutVerbose) {
628 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
633 drop_vec[i].drop = drop;
635 }
else if (classicalAlgo == scaled_cut) {
636 std::sort( drop_vec.begin(), drop_vec.end()
637 , [](DropTol
const& a, DropTol
const& b) {
638 return a.val/a.diag > b.val/b.diag;
643 for (
size_t i=1; i<n; ++i) {
645 auto const& x = drop_vec[i-1];
646 auto const& y = drop_vec[i];
647 auto a = x.val/x.diag;
648 auto b = y.val/y.diag;
649 if (a > realThreshold*b) {
652#ifdef HAVE_MUELU_DEBUG
653 if (distanceLaplacianCutVerbose) {
654 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
661 drop_vec[i].drop = drop;
665 std::sort( drop_vec.begin(), drop_vec.end()
666 , [](DropTol
const& a, DropTol
const& b) {
667 return a.col < b.col;
671 for (LO idxID =0; idxID<(LO)drop_vec.size(); idxID++) {
672 LO col = indices[drop_vec[idxID].col];
675 columns[realnnz++] = col;
680 if (!drop_vec[idxID].drop) {
681 columns[realnnz++] = col;
688 rows[row+1] = realnnz;
693 columns.resize(realnnz);
694 numTotal = A->getLocalNumEntries();
696 if (aggregationMayCreateDirichlet) {
698 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); ++row) {
700 boundaryNodes[row] =
true;
704 RCP<GraphBase> graph = rcp(
new LWGraph(
rows, columns, A->getRowMap(), A->getColMap(),
"thresholded graph of A"));
705 graph->SetBoundaryNodeMap(boundaryNodes);
707 GO numLocalBoundaryNodes = 0;
708 GO numGlobalBoundaryNodes = 0;
709 for (LO i = 0; i < boundaryNodes.size(); ++i)
710 if (boundaryNodes[i])
711 numLocalBoundaryNodes++;
712 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
713 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
714 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
716 Set(currentLevel,
"Graph", graph);
717 Set(currentLevel,
"DofsPerNode", 1);
720 if(generateColoringGraph) {
721 RCP<GraphBase> colorGraph;
722 RCP<const Import> importer = A->getCrsGraph()->getImporter();
723 BlockDiagonalizeGraph(graph,ghostedBlockNumber,colorGraph,importer);
724 Set(currentLevel,
"Coloring Graph",colorGraph);
728 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"m_regular_graph."+std::to_string(currentLevel.
GetLevelID()), *rcp_dynamic_cast<LWGraph>(graph)->GetCrsGraph());
729 Xpetra::IO<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Write(
"m_color_graph."+std::to_string(currentLevel.
GetLevelID()), *rcp_dynamic_cast<LWGraph>(colorGraph)->GetCrsGraph());
745 }
else if (BlockSize > 1 && threshold == STS::zero()) {
747 const RCP<const Map> rowMap = A->getRowMap();
748 const RCP<const Map> colMap = A->getColMap();
750 graphType =
"amalgamated";
756 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
757 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
758 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
759 Array<LO> colTranslation = *(amalInfo->getColTranslation());
762 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
765 ArrayRCP<LO>
rows = ArrayRCP<LO>(numRows+1);
766 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getLocalNumEntries());
768 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
774 ArrayRCP<bool > pointBoundaryNodes;
781 LO blkSize = A->GetFixedBlockSize();
783 LO blkPartSize = A->GetFixedBlockSize();
784 if (A->IsView(
"stridedMaps") ==
true) {
785 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
786 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
788 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
789 blkId = strMap->getStridedBlockId();
791 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
797 Array<LO> indicesExtra;
798 for (LO row = 0; row < numRows; row++) {
799 ArrayView<const LO> indices;
800 indicesExtra.resize(0);
808 bool isBoundary =
false;
809 if (pL.get<
bool>(
"aggregation: greedy Dirichlet") ==
true) {
810 for (LO j = 0; j < blkPartSize; j++) {
811 if (pointBoundaryNodes[row*blkPartSize+j]) {
818 for (LO j = 0; j < blkPartSize; j++) {
819 if (!pointBoundaryNodes[row*blkPartSize+j]) {
829 MergeRows(*A, row, indicesExtra, colTranslation);
831 indicesExtra.push_back(row);
832 indices = indicesExtra;
833 numTotal += indices.size();
837 LO nnz = indices.size(), rownnz = 0;
838 for (LO colID = 0; colID < nnz; colID++) {
839 LO col = indices[colID];
840 columns[realnnz++] = col;
851 amalgBoundaryNodes[row] =
true;
853 rows[row+1] = realnnz;
855 columns.resize(realnnz);
857 RCP<GraphBase> graph = rcp(
new LWGraph(
rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
858 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
861 GO numLocalBoundaryNodes = 0;
862 GO numGlobalBoundaryNodes = 0;
864 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
865 if (amalgBoundaryNodes[i])
866 numLocalBoundaryNodes++;
868 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
869 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
870 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
871 <<
" agglomerated Dirichlet nodes" << std::endl;
874 Set(currentLevel,
"Graph", graph);
875 Set(currentLevel,
"DofsPerNode", blkSize);
877 }
else if (BlockSize > 1 && threshold != STS::zero()) {
879 const RCP<const Map> rowMap = A->getRowMap();
880 const RCP<const Map> colMap = A->getColMap();
881 graphType =
"amalgamated";
887 RCP<const Map> uniqueMap = amalInfo->getNodeRowMap();
888 RCP<const Map> nonUniqueMap = amalInfo->getNodeColMap();
889 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
890 Array<LO> colTranslation = *(amalInfo->getColTranslation());
893 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
896 ArrayRCP<LO>
rows = ArrayRCP<LO>(numRows+1);
897 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getLocalNumEntries());
899 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
905 ArrayRCP<bool > pointBoundaryNodes;
912 LO blkSize = A->GetFixedBlockSize();
914 LO blkPartSize = A->GetFixedBlockSize();
915 if (A->IsView(
"stridedMaps") ==
true) {
916 Teuchos::RCP<const Map> myMap = A->getRowMap(
"stridedMaps");
917 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
919 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
920 blkId = strMap->getStridedBlockId();
922 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
927 const ArrayRCP<const SC> ghostedDiagVals = ghostedDiag->getData(0);
932 Array<LO> indicesExtra;
933 for (LO row = 0; row < numRows; row++) {
934 ArrayView<const LO> indices;
935 indicesExtra.resize(0);
943 bool isBoundary =
false;
944 if (pL.get<
bool>(
"aggregation: greedy Dirichlet") ==
true) {
945 for (LO j = 0; j < blkPartSize; j++) {
946 if (pointBoundaryNodes[row*blkPartSize+j]) {
953 for (LO j = 0; j < blkPartSize; j++) {
954 if (!pointBoundaryNodes[row*blkPartSize+j]) {
964 MergeRowsWithDropping(*A, row, ghostedDiagVals, threshold, indicesExtra, colTranslation);
966 indicesExtra.push_back(row);
967 indices = indicesExtra;
968 numTotal += indices.size();
972 LO nnz = indices.size(), rownnz = 0;
973 for (LO colID = 0; colID < nnz; colID++) {
974 LO col = indices[colID];
975 columns[realnnz++] = col;
986 amalgBoundaryNodes[row] =
true;
988 rows[row+1] = realnnz;
990 columns.resize(realnnz);
992 RCP<GraphBase> graph = rcp(
new LWGraph(
rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
993 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
996 GO numLocalBoundaryNodes = 0;
997 GO numGlobalBoundaryNodes = 0;
999 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
1000 if (amalgBoundaryNodes[i])
1001 numLocalBoundaryNodes++;
1003 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1004 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1005 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
1006 <<
" agglomerated Dirichlet nodes" << std::endl;
1009 Set(currentLevel,
"Graph", graph);
1010 Set(currentLevel,
"DofsPerNode", blkSize);
1013 }
else if (algo ==
"distance laplacian") {
1014 LO blkSize = A->GetFixedBlockSize();
1015 GO indexBase = A->getRowMap()->getIndexBase();
1026 ArrayRCP<bool > pointBoundaryNodes;
1031 if ( (blkSize == 1) && (threshold == STS::zero()) ) {
1033 RCP<GraphBase> graph = rcp(
new Graph(A->getCrsGraph(),
"graph of A"));
1034 graph->SetBoundaryNodeMap(pointBoundaryNodes);
1035 graphType=
"unamalgamated";
1036 numTotal = A->getLocalNumEntries();
1039 GO numLocalBoundaryNodes = 0;
1040 GO numGlobalBoundaryNodes = 0;
1041 for (LO i = 0; i < pointBoundaryNodes.size(); ++i)
1042 if (pointBoundaryNodes[i])
1043 numLocalBoundaryNodes++;
1044 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1045 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1046 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
1049 Set(currentLevel,
"DofsPerNode", blkSize);
1050 Set(currentLevel,
"Graph", graph);
1065 TEUCHOS_TEST_FOR_EXCEPTION(A->getRowMap()->getLocalNumElements()/blkSize != Coords->getLocalLength(),
Exceptions::Incompatible,
1066 "Coordinate vector length (" << Coords->getLocalLength() <<
") is incompatible with number of rows in A (" << A->getRowMap()->getLocalNumElements() <<
") by modulo block size ("<< blkSize <<
").");
1068 const RCP<const Map> colMap = A->getColMap();
1069 RCP<const Map> uniqueMap, nonUniqueMap;
1070 Array<LO> colTranslation;
1072 uniqueMap = A->getRowMap();
1073 nonUniqueMap = A->getColMap();
1074 graphType=
"unamalgamated";
1077 uniqueMap = Coords->getMap();
1079 "Different index bases for matrix and coordinates");
1083 graphType =
"amalgamated";
1085 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getLocalNumElements());
1087 RCP<RealValuedMultiVector> ghostedCoords;
1088 RCP<Vector> ghostedLaplDiag;
1089 Teuchos::ArrayRCP<SC> ghostedLaplDiagData;
1090 if (threshold != STS::zero()) {
1092 RCP<const Import> importer;
1095 if (blkSize == 1 && realA->getCrsGraph()->getImporter() != Teuchos::null) {
1096 GetOStream(
Warnings1) <<
"Using existing importer from matrix graph" << std::endl;
1097 importer = realA->getCrsGraph()->getImporter();
1099 GetOStream(
Warnings0) <<
"Constructing new importer instance" << std::endl;
1100 importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
1103 ghostedCoords = Xpetra::MultiVectorFactory<real_type,LO,GO,NO>::Build(nonUniqueMap, Coords->getNumVectors());
1106 ghostedCoords->doImport(*Coords, *importer, Xpetra::INSERT);
1110 RCP<Vector> localLaplDiag = VectorFactory::Build(uniqueMap);
1111 Array<LO> indicesExtra;
1112 Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
1113 if (threshold != STS::zero()) {
1114 const size_t numVectors = ghostedCoords->getNumVectors();
1115 coordData.reserve(numVectors);
1116 for (
size_t j = 0; j < numVectors; j++) {
1117 Teuchos::ArrayRCP<const real_type> tmpData=ghostedCoords->getData(j);
1118 coordData.push_back(tmpData);
1123 ArrayRCP<SC> localLaplDiagData = localLaplDiag->getDataNonConst(0);
1124 for (LO row = 0; row < numRows; row++) {
1125 ArrayView<const LO> indices;
1128 ArrayView<const SC> vals;
1129 A->getLocalRowView(row, indices, vals);
1133 indicesExtra.resize(0);
1134 MergeRows(*A, row, indicesExtra, colTranslation);
1135 indices = indicesExtra;
1138 LO nnz = indices.size();
1139 bool haveAddedToDiag =
false;
1140 for (LO colID = 0; colID < nnz; colID++) {
1141 const LO col = indices[colID];
1144 if(use_dlap_weights == SINGLE_WEIGHTS) {
1150 else if(use_dlap_weights == BLOCK_WEIGHTS) {
1151 int block_id = row % interleaved_blocksize;
1152 int block_start = block_id * interleaved_blocksize;
1159 haveAddedToDiag =
true;
1164 if (!haveAddedToDiag)
1165 localLaplDiagData[row] = STS::rmax();
1170 ghostedLaplDiag = VectorFactory::Build(nonUniqueMap);
1171 ghostedLaplDiag->doImport(*localLaplDiag, *importer, Xpetra::INSERT);
1172 ghostedLaplDiagData = ghostedLaplDiag->getDataNonConst(0);
1176 GetOStream(
Runtime0) <<
"Skipping distance laplacian construction due to 0 threshold" << std::endl;
1182 ArrayRCP<LO>
rows = ArrayRCP<LO>(numRows+1);
1183 ArrayRCP<LO> columns = ArrayRCP<LO>(A->getLocalNumEntries());
1185#ifdef HAVE_MUELU_DEBUG
1187 for(LO i=0; i<(LO)columns.size(); i++) columns[i]=-666;
1191 ArrayRCP<LO> rows_stop;
1192 bool use_stop_array = threshold != STS::zero() && distanceLaplacianAlgo == scaled_cut_symmetric;
1194 rows_stop.resize(numRows);
1196 const ArrayRCP<bool> amalgBoundaryNodes(numRows,
false);
1201 Array<LO> indicesExtra;
1204 Teuchos::Array<Teuchos::ArrayRCP<const real_type>> coordData;
1205 if (threshold != STS::zero()) {
1206 const size_t numVectors = ghostedCoords->getNumVectors();
1207 coordData.reserve(numVectors);
1208 for (
size_t j = 0; j < numVectors; j++) {
1209 Teuchos::ArrayRCP<const real_type> tmpData=ghostedCoords->getData(j);
1210 coordData.push_back(tmpData);
1214 ArrayView<const SC> vals;
1215 for (LO row = 0; row < numRows; row++) {
1216 ArrayView<const LO> indices;
1217 indicesExtra.resize(0);
1218 bool isBoundary =
false;
1222 A->getLocalRowView(row, indices, vals);
1223 isBoundary = pointBoundaryNodes[row];
1226 for (LO j = 0; j < blkSize; j++) {
1227 if (!pointBoundaryNodes[row*blkSize+j]) {
1235 MergeRows(*A, row, indicesExtra, colTranslation);
1237 indicesExtra.push_back(row);
1238 indices = indicesExtra;
1240 numTotal += indices.size();
1242 LO nnz = indices.size(), rownnz = 0;
1244 if(use_stop_array) {
1246 realnnz =
rows[row];
1249 if (threshold != STS::zero()) {
1251 if (distanceLaplacianAlgo == defaultAlgo) {
1253 for (LO colID = 0; colID < nnz; colID++) {
1255 LO col = indices[colID];
1258 columns[realnnz++] = col;
1264 if(isBoundary)
continue;
1267 if(use_dlap_weights == SINGLE_WEIGHTS) {
1270 else if(use_dlap_weights == BLOCK_WEIGHTS) {
1271 int block_id = row % interleaved_blocksize;
1272 int block_start = block_id * interleaved_blocksize;
1278 real_type aiiajj = STS::magnitude(realThreshold*realThreshold * ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
1279 real_type aij = STS::magnitude(laplVal*laplVal);
1282 columns[realnnz++] = col;
1291 std::vector<DropTol> drop_vec;
1292 drop_vec.reserve(nnz);
1293 const real_type zero = Teuchos::ScalarTraits<real_type>::zero();
1294 const real_type one = Teuchos::ScalarTraits<real_type>::one();
1297 for (LO colID = 0; colID < nnz; colID++) {
1299 LO col = indices[colID];
1302 drop_vec.emplace_back( zero, one, colID,
false);
1306 if(isBoundary)
continue;
1309 if(use_dlap_weights == SINGLE_WEIGHTS) {
1312 else if(use_dlap_weights == BLOCK_WEIGHTS) {
1313 int block_id = row % interleaved_blocksize;
1314 int block_start = block_id * interleaved_blocksize;
1321 real_type aiiajj = STS::magnitude(ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
1322 real_type aij = STS::magnitude(laplVal*laplVal);
1324 drop_vec.emplace_back(aij, aiiajj, colID,
false);
1327 const size_t n = drop_vec.size();
1329 if (distanceLaplacianAlgo == unscaled_cut) {
1331 std::sort( drop_vec.begin(), drop_vec.end()
1332 , [](DropTol
const& a, DropTol
const& b) {
1333 return a.val > b.val;
1338 for (
size_t i=1; i<n; ++i) {
1340 auto const& x = drop_vec[i-1];
1341 auto const& y = drop_vec[i];
1344 if (a > realThreshold*b) {
1346#ifdef HAVE_MUELU_DEBUG
1347 if (distanceLaplacianCutVerbose) {
1348 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
1353 drop_vec[i].drop = drop;
1356 else if (distanceLaplacianAlgo == scaled_cut || distanceLaplacianAlgo == scaled_cut_symmetric) {
1358 std::sort( drop_vec.begin(), drop_vec.end()
1359 , [](DropTol
const& a, DropTol
const& b) {
1360 return a.val/a.diag > b.val/b.diag;
1365 for (
size_t i=1; i<n; ++i) {
1367 auto const& x = drop_vec[i-1];
1368 auto const& y = drop_vec[i];
1369 auto a = x.val/x.diag;
1370 auto b = y.val/y.diag;
1371 if (a > realThreshold*b) {
1373#ifdef HAVE_MUELU_DEBUG
1374 if (distanceLaplacianCutVerbose) {
1375 std::cout <<
"DJS: KEEP, N, ROW: " << i+1 <<
", " << n <<
", " << row << std::endl;
1380 drop_vec[i].drop = drop;
1384 std::sort( drop_vec.begin(), drop_vec.end()
1385 , [](DropTol
const& a, DropTol
const& b) {
1386 return a.col < b.col;
1390 for (LO idxID =0; idxID<(LO)drop_vec.size(); idxID++) {
1391 LO col = indices[drop_vec[idxID].col];
1396 columns[realnnz++] = col;
1402 if (!drop_vec[idxID].drop) {
1403 columns[realnnz++] = col;
1414 for (LO colID = 0; colID < nnz; colID++) {
1415 LO col = indices[colID];
1416 columns[realnnz++] = col;
1428 amalgBoundaryNodes[row] =
true;
1432 rows_stop[row] = rownnz +
rows[row];
1434 rows[row+1] = realnnz;
1439 if (use_stop_array) {
1442 for (LO row = 0; row < numRows; row++) {
1443 for (LO colidx =
rows[row]; colidx < rows_stop[row]; colidx++) {
1444 LO col = columns[colidx];
1445 if(col >= numRows)
continue;
1448 for(LO t_col =
rows[col] ; !found && t_col < rows_stop[col]; t_col++) {
1449 if (columns[t_col] == row)
1454 if(!found && !pointBoundaryNodes[col] && rows_stop[col] <
rows[col+1]) {
1455 LO new_idx = rows_stop[col];
1457 columns[new_idx] = row;
1466 for (LO row = 0; row < numRows; row++) {
1467 LO old_start = current_start;
1468 for (LO col =
rows[row]; col < rows_stop[row]; col++) {
1469 if(current_start != col) {
1470 columns[current_start] = columns[col];
1474 rows[row] = old_start;
1476 rows[numRows] = realnnz = current_start;
1480 columns.resize(realnnz);
1482 RCP<GraphBase> graph;
1485 graph = rcp(
new LWGraph(
rows, columns, uniqueMap, nonUniqueMap,
"amalgamated graph of A"));
1486 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
1490 GO numLocalBoundaryNodes = 0;
1491 GO numGlobalBoundaryNodes = 0;
1493 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
1494 if (amalgBoundaryNodes[i])
1495 numLocalBoundaryNodes++;
1497 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1498 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1499 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" agglomerated Dirichlet nodes"
1500 <<
" using threshold " << dirichletThreshold << std::endl;
1503 Set(currentLevel,
"Graph", graph);
1504 Set(currentLevel,
"DofsPerNode", blkSize);
1508 if ((GetVerbLevel() &
Statistics1) && !(A->GetFixedBlockSize() > 1 && threshold != STS::zero())) {
1509 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1510 GO numGlobalTotal, numGlobalDropped;
1513 GetOStream(
Statistics1) <<
"Number of dropped entries in " << graphType <<
" matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1514 if (numGlobalTotal != 0)
1515 GetOStream(
Statistics1) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
1522 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
1524 GetOStream(
Runtime0) <<
"algorithm = \"" <<
"failsafe" <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
1525 Set<bool>(currentLevel,
"Filtering", (threshold != STS::zero()));
1527 RCP<const Map> rowMap = A->getRowMap();
1528 RCP<const Map> colMap = A->getColMap();
1531 GO indexBase = rowMap->getIndexBase();
1535 if(A->IsView(
"stridedMaps") &&
1536 Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
1537 Xpetra::viewLabel_t oldView = A->SwitchToView(
"stridedMaps");
1538 RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(A->getRowMap());
1539 TEUCHOS_TEST_FOR_EXCEPTION(strMap == Teuchos::null,
Exceptions::BadCast,
"MueLu::CoalesceFactory::Build: cast to strided row map failed.");
1540 blockdim = strMap->getFixedBlockSize();
1541 offset = strMap->getOffset();
1542 oldView = A->SwitchToView(oldView);
1543 GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build():" <<
" found blockdim=" << blockdim <<
" from strided maps. offset=" << offset << std::endl;
1544 }
else GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
1548 RCP<const Map> nodeMap = amalInfo->getNodeRowMap();
1549 GetOStream(
Statistics1) <<
"CoalesceDropFactory: nodeMap " << nodeMap->getLocalNumElements() <<
"/" << nodeMap->getGlobalNumElements() <<
" elements" << std::endl;
1552 RCP<CrsGraph> crsGraph = CrsGraphFactory::Build(nodeMap, A->getLocalMaxNumRowEntries()*blockdim);
1554 LO numRows = A->getRowMap()->getLocalNumElements();
1555 LO numNodes = nodeMap->getLocalNumElements();
1556 const ArrayRCP<bool> amalgBoundaryNodes(numNodes,
false);
1557 const ArrayRCP<int> numberDirichletRowsPerNode(numNodes, 0);
1558 bool bIsDiagonalEntry =
false;
1563 for(LO row=0; row<numRows; row++) {
1565 GO grid = rowMap->getGlobalElement(row);
1568 bIsDiagonalEntry =
false;
1573 size_t nnz = A->getNumEntriesInLocalRow(row);
1574 Teuchos::ArrayView<const LO> indices;
1575 Teuchos::ArrayView<const SC> vals;
1576 A->getLocalRowView(row, indices, vals);
1578 RCP<std::vector<GO> > cnodeIds = Teuchos::rcp(
new std::vector<GO>);
1580 for(LO col=0; col<Teuchos::as<LO>(nnz); col++) {
1581 GO gcid = colMap->getGlobalElement(indices[col]);
1583 if(vals[col]!=STS::zero()) {
1585 cnodeIds->push_back(cnodeId);
1587 if (grid == gcid) bIsDiagonalEntry =
true;
1591 if(realnnz == 1 && bIsDiagonalEntry ==
true) {
1592 LO lNodeId = nodeMap->getLocalElement(nodeId);
1593 numberDirichletRowsPerNode[lNodeId] += 1;
1594 if (numberDirichletRowsPerNode[lNodeId] == blockdim)
1595 amalgBoundaryNodes[lNodeId] =
true;
1598 Teuchos::ArrayRCP<GO> arr_cnodeIds = Teuchos::arcp( cnodeIds );
1600 if(arr_cnodeIds.size() > 0 )
1601 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
1604 crsGraph->fillComplete(nodeMap,nodeMap);
1607 RCP<GraphBase> graph = rcp(
new Graph(crsGraph,
"amalgamated graph of A"));
1610 graph->SetBoundaryNodeMap(amalgBoundaryNodes);
1613 GO numLocalBoundaryNodes = 0;
1614 GO numGlobalBoundaryNodes = 0;
1615 for (LO i = 0; i < amalgBoundaryNodes.size(); ++i)
1616 if (amalgBoundaryNodes[i])
1617 numLocalBoundaryNodes++;
1618 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1619 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1620 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
1625 Set(currentLevel,
"DofsPerNode", blockdim);
1626 Set(currentLevel,
"Graph", graph);
1633 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1635 typedef typename ArrayView<const LO>::size_type size_type;
1638 LO blkSize = A.GetFixedBlockSize();
1639 if (A.IsView(
"stridedMaps") ==
true) {
1640 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
1641 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
1643 if (strMap->getStridedBlockId() > -1)
1644 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1648 size_t nnz = 0, pos = 0;
1649 for (LO j = 0; j < blkSize; j++)
1650 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
1660 ArrayView<const LO> inds;
1661 ArrayView<const SC> vals;
1662 for (LO j = 0; j < blkSize; j++) {
1663 A.getLocalRowView(row*blkSize+j, inds, vals);
1664 size_type numIndices = inds.size();
1666 if (numIndices == 0)
1670 cols[pos++] = translation[inds[0]];
1671 for (size_type k = 1; k < numIndices; k++) {
1672 LO nodeID = translation[inds[k]];
1676 if (nodeID != cols[pos-1])
1677 cols[pos++] = nodeID;
1684 std::sort(cols.begin(), cols.end());
1686 for (
size_t j = 1; j < nnz; j++)
1687 if (cols[j] != cols[pos])
1688 cols[++pos] = cols[j];
1692 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1694 typedef typename ArrayView<const LO>::size_type size_type;
1695 typedef Teuchos::ScalarTraits<SC> STS;
1698 LO blkSize = A.GetFixedBlockSize();
1699 if (A.IsView(
"stridedMaps") ==
true) {
1700 Teuchos::RCP<const Map> myMap = A.getRowMap(
"stridedMaps");
1701 Teuchos::RCP<const StridedMap> strMap = Teuchos::rcp_dynamic_cast<const StridedMap>(myMap);
1703 if (strMap->getStridedBlockId() > -1)
1704 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
1708 size_t nnz = 0, pos = 0;
1709 for (LO j = 0; j < blkSize; j++)
1710 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
1720 ArrayView<const LO> inds;
1721 ArrayView<const SC> vals;
1722 for (LO j = 0; j < blkSize; j++) {
1723 A.getLocalRowView(row*blkSize+j, inds, vals);
1724 size_type numIndices = inds.size();
1726 if (numIndices == 0)
1731 for (size_type k = 0; k < numIndices; k++) {
1733 LO nodeID = translation[inds[k]];
1736 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold*ghostedDiagVals[dofID]*ghostedDiagVals[row*blkSize+j]);
1737 typename STS::magnitudeType aij = STS::magnitude(vals[k]*vals[k]);
1740 if (aij > aiiajj || (row*blkSize+j == dofID)) {
1746 if (nodeID != prevNodeID) {
1747 cols[pos++] = nodeID;
1748 prevNodeID = nodeID;
1757 std::sort(cols.begin(), cols.end());
1759 for (
size_t j = 1; j < nnz; j++)
1760 if (cols[j] != cols[pos])
1761 cols[++pos] = cols[j];
1769 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1771 typedef Teuchos::ScalarTraits<SC> STS;
1773 const ParameterList & pL = GetParameterList();
1774 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<
double>(
"aggregation: Dirichlet threshold")));
1775 const typename STS::magnitudeType rowSumTol = as<typename STS::magnitudeType>(pL.get<
double>(
"aggregation: row sum drop tol"));
1777 RCP<LocalOrdinalVector> BlockNumber = Get<RCP<LocalOrdinalVector> >(currentLevel,
"BlockNumber");
1778 RCP<LocalOrdinalVector> ghostedBlockNumber;
1779 GetOStream(
Statistics1) <<
"Using BlockDiagonal Graph before dropping (with provided blocking)"<<std::endl;
1782 RCP<const Import> importer = A->getCrsGraph()->getImporter();
1783 if(!importer.is_null()) {
1785 ghostedBlockNumber= Xpetra::VectorFactory<LO,LO,GO,NO>::Build(importer->getTargetMap());
1786 ghostedBlockNumber->doImport(*BlockNumber, *importer, Xpetra::INSERT);
1789 ghostedBlockNumber = BlockNumber;
1793 Teuchos::ArrayRCP<const LO> row_block_number = BlockNumber->getData(0);
1794 Teuchos::ArrayRCP<const LO> col_block_number = ghostedBlockNumber->getData(0);
1797 ArrayRCP<size_t> rows_mat;
1798 ArrayRCP<LO> rows_graph,columns;
1799 ArrayRCP<SC> values;
1800 RCP<CrsMatrixWrap> crs_matrix_wrap;
1802 if(generate_matrix) {
1803 crs_matrix_wrap = rcp(
new CrsMatrixWrap(A->getRowMap(), A->getColMap(), 0));
1804 crs_matrix_wrap->getCrsMatrix()->allocateAllValues(A->getLocalNumEntries(), rows_mat, columns, values);
1807 rows_graph.resize(A->getLocalNumRows()+1);
1808 columns.resize(A->getLocalNumEntries());
1809 values.resize(A->getLocalNumEntries());
1813 GO numDropped = 0, numTotal = 0;
1814 for (LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getLocalNumElements()); ++row) {
1815 LO row_block = row_block_number[row];
1816 size_t nnz = A->getNumEntriesInLocalRow(row);
1817 ArrayView<const LO> indices;
1818 ArrayView<const SC> vals;
1819 A->getLocalRowView(row, indices, vals);
1822 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
1823 LO col = indices[colID];
1824 LO col_block = col_block_number[col];
1826 if(row_block == col_block) {
1827 if(generate_matrix) values[realnnz] = vals[colID];
1828 columns[realnnz++] = col;
1833 if(generate_matrix) rows_mat[row+1] = realnnz;
1834 else rows_graph[row+1] = realnnz;
1842 if(!generate_matrix) {
1844 values.resize(realnnz);
1845 columns.resize(realnnz);
1847 numTotal = A->getLocalNumEntries();
1850 GO numLocalBoundaryNodes = 0;
1851 GO numGlobalBoundaryNodes = 0;
1852 for (LO i = 0; i < boundaryNodes.size(); ++i)
1853 if (boundaryNodes[i])
1854 numLocalBoundaryNodes++;
1855 RCP<const Teuchos::Comm<int> > comm = A->getRowMap()->getComm();
1856 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
1857 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
1859 GO numGlobalTotal, numGlobalDropped;
1862 GetOStream(
Statistics1) <<
"Number of dropped entries in block-diagonalized matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1863 if (numGlobalTotal != 0)
1864 GetOStream(
Statistics1) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
1868 Set(currentLevel,
"Filtering",
true);
1870 if(generate_matrix) {
1874 crs_matrix_wrap->getCrsMatrix()->setAllValues(rows_mat,columns,values);
1875 crs_matrix_wrap->getCrsMatrix()->expertStaticFillComplete(A->getColMap(), A->getRowMap());
1878 RCP<GraphBase> graph = rcp(
new LWGraph(rows_graph, columns, A->getRowMap(), A->getColMap(),
"block-diagonalized graph of A"));
1879 graph->SetBoundaryNodeMap(boundaryNodes);
1880 Set(currentLevel,
"Graph", graph);
1884 Set(currentLevel,
"DofsPerNode", 1);
1885 return crs_matrix_wrap;
1889 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1892 TEUCHOS_TEST_FOR_EXCEPTION(ghostedBlockNumber.is_null(),
Exceptions::RuntimeError,
"BlockDiagonalizeGraph(): ghostedBlockNumber is null.");
1893 const ParameterList & pL = GetParameterList();
1895 const bool localizeColoringGraph = pL.get<
bool>(
"aggregation: coloring: localize color graph");
1897 GetOStream(
Statistics1) <<
"Using BlockDiagonal Graph after Dropping (with provided blocking)";
1898 if (localizeColoringGraph)
1899 GetOStream(
Statistics1) <<
", with localization" <<std::endl;
1901 GetOStream(
Statistics1) <<
", without localization" <<std::endl;
1904 Teuchos::ArrayRCP<const LO> row_block_number = ghostedBlockNumber->getData(0);
1905 Teuchos::ArrayRCP<const LO> col_block_number = ghostedBlockNumber->getData(0);
1908 ArrayRCP<size_t> rows_mat;
1909 ArrayRCP<LO> rows_graph,columns;
1911 rows_graph.resize(inputGraph->GetNodeNumVertices()+1);
1912 columns.resize(inputGraph->GetNodeNumEdges());
1915 GO numDropped = 0, numTotal = 0;
1916 const LO numRows = Teuchos::as<LO>(inputGraph->GetDomainMap()->getLocalNumElements());
1917 if (localizeColoringGraph) {
1919 for (LO row = 0; row < numRows; ++row) {
1920 LO row_block = row_block_number[row];
1921 ArrayView<const LO> indices = inputGraph->getNeighborVertices(row);
1924 for (LO colID = 0; colID < Teuchos::as<LO>(indices.size()); colID++) {
1925 LO col = indices[colID];
1926 LO col_block = col_block_number[col];
1928 if((row_block == col_block) && (col < numRows)) {
1929 columns[realnnz++] = col;
1934 rows_graph[row+1] = realnnz;
1938 Teuchos::ArrayRCP<const bool> boundaryNodes = inputGraph->GetBoundaryNodeMap();
1939 auto boundaryNodesVector = Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(inputGraph->GetDomainMap());
1940 for (
size_t i=0; i<inputGraph->GetNodeNumVertices(); i++)
1941 boundaryNodesVector->getDataNonConst(0)[i] = boundaryNodes[i];
1943 auto boundaryColumnVector = Xpetra::VectorFactory<LocalOrdinal,LocalOrdinal,GlobalOrdinal,Node>::Build(inputGraph->GetImportMap());
1944 boundaryColumnVector->doImport(*boundaryNodesVector,*importer, Xpetra::INSERT);
1945 auto boundaryColumn = boundaryColumnVector->getData(0);
1947 for (LO row = 0; row < numRows; ++row) {
1948 LO row_block = row_block_number[row];
1949 ArrayView<const LO> indices = inputGraph->getNeighborVertices(row);
1952 for (LO colID = 0; colID < Teuchos::as<LO>(indices.size()); colID++) {
1953 LO col = indices[colID];
1954 LO col_block = col_block_number[col];
1956 if((row_block == col_block) && ((row == col) || (boundaryColumn[col] == 0))) {
1957 columns[realnnz++] = col;
1962 rows_graph[row+1] = realnnz;
1966 columns.resize(realnnz);
1967 numTotal = inputGraph->GetNodeNumEdges();
1970 RCP<const Teuchos::Comm<int> > comm = inputGraph->GetDomainMap()->getComm();
1971 GO numGlobalTotal, numGlobalDropped;
1974 GetOStream(
Statistics1) <<
"Number of dropped entries in block-diagonalized matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
1975 if (numGlobalTotal != 0)
1976 GetOStream(
Statistics1) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
1980 if (localizeColoringGraph) {
1981 outputGraph = rcp(
new LWGraph(rows_graph, columns, inputGraph->GetDomainMap(), inputGraph->GetImportMap(),
"block-diagonalized graph of A"));
1982 outputGraph->SetBoundaryNodeMap(inputGraph->GetBoundaryNodeMap());
1984 TEUCHOS_ASSERT(inputGraph->GetDomainMap()->lib() == Xpetra::UseTpetra);
1985#ifdef HAVE_XPETRA_TPETRA
1986 auto outputGraph2 = rcp(
new LWGraph(rows_graph, columns, inputGraph->GetDomainMap(), inputGraph->GetImportMap(),
"block-diagonalized graph of A"));
1988 auto tpGraph = Xpetra::toTpetra(rcp_const_cast<const CrsGraph>(outputGraph2->GetCrsGraph()));
1989 auto sym = rcp(
new Tpetra::CrsGraphTransposer<LocalOrdinal,GlobalOrdinal,Node>(tpGraph));
1990 auto tpGraphSym = sym->symmetrize();
1993 Kokkos::Compat::persistingView(tpGraphSym->getLocalIndicesHost());
1995 auto rowsSym = tpGraphSym->getLocalRowPtrsHost();
1996 ArrayRCP<LO> rows_graphSym;
1997 rows_graphSym.resize(rowsSym.size());
1998 for (
size_t row = 0; row < rowsSym.size(); row++)
1999 rows_graphSym[row] = rowsSym[row];
2000 outputGraph = rcp(
new LWGraph(rows_graphSym, colIndsSym, inputGraph->GetDomainMap(), Xpetra::toXpetra(tpGraphSym->getColMap()),
"block-diagonalized graph of A"));
2001 outputGraph->SetBoundaryNodeMap(inputGraph->GetBoundaryNodeMap());
#define SET_VALID_ENTRY(name)
#define MueLu_sumAll(rcpComm, in, out)
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 const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
Translate global (row/column) id to global amalgamation block id.
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
void DeclareInput(Level ¤tLevel) const
Input.
void MergeRowsWithDropping(const Matrix &A, const LO row, const ArrayRCP< const SC > &ghostedDiagVals, SC threshold, Array< LO > &cols, const Array< LO > &translation) const
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
CoalesceDropFactory()
Constructor.
void Build(Level ¤tLevel) const
Build an object with this factory.
void BlockDiagonalizeGraph(const RCP< GraphBase > &inputGraph, const RCP< LocalOrdinalVector > &ghostedBlockNumber, RCP< GraphBase > &outputGraph, RCP< const Import > &importer) const
Exception indicating invalid cast attempted.
Exception throws to report incompatible objects (like maps).
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
MueLu representation of a compressed row storage graph.
Lightweight MueLu representation of a compressed row storage graph.
Class that holds all level-specific information.
int GetLevelID() const
Return level number.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
static Teuchos::ArrayRCP< Magnitude > GetMatrixMaxMinusOffDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Teuchos::Array< Teuchos::ArrayRCP< const Scalar > > &v, LocalOrdinal i0, LocalOrdinal i1)
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Statistics1
Print more statistics.
@ Runtime0
One-liner description of what is happening.
@ Warnings1
Additional warnings.
@ Parameters0
Print class parameters.
DropTol & operator=(DropTol const &)=default
DropTol(DropTol &&)=default
DropTol(real_type val_, real_type diag_, LO col_, bool drop_)
DropTol & operator=(DropTol &&)=default
DropTol(DropTol const &)=default