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);