46 #ifndef MUELU_COALESCEDROPFACTORY_DEF_HPP 47 #define MUELU_COALESCEDROPFACTORY_DEF_HPP 63 #include "MueLu_AmalgamationFactory.hpp" 64 #include "MueLu_AmalgamationInfo.hpp" 67 #include "MueLu_Graph.hpp" 69 #include "MueLu_LWGraph.hpp" 72 #include "MueLu_PreDropFunctionBaseClass.hpp" 73 #include "MueLu_PreDropFunctionConstVal.hpp" 74 #include "MueLu_Utilities.hpp" 78 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
82 #define SET_VALID_ENTRY(name) validParamList->setEntry(name, MasterList::getEntry(name)) 89 rcp(
new validatorType(Teuchos::tuple<std::string>(
"classical",
"distance laplacian"),
"aggregation: drop scheme")));
91 #undef SET_VALID_ENTRY 92 validParamList->
set<
bool > (
"lightweight wrap",
true,
"Experimental option for lightweight graph access");
95 validParamList->
set<
RCP<const FactoryBase> >(
"UnAmalgamationInfo", Teuchos::null,
"Generating factory for UnAmalgamationInfo");
98 return validParamList;
101 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
104 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
106 Input(currentLevel,
"A");
107 Input(currentLevel,
"UnAmalgamationInfo");
110 if (pL.
get<
bool>(
"lightweight wrap") ==
true) {
111 if (pL.
get<std::string>(
"aggregation: drop scheme") ==
"distance laplacian")
112 Input(currentLevel,
"Coordinates");
117 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
125 if (predrop_ != Teuchos::null)
126 GetOStream(
Parameters0) << predrop_->description();
128 RCP<Matrix> A = Get< RCP<Matrix> >(currentLevel,
"A");
132 bool doExperimentalWrap = pL.get<
bool>(
"lightweight wrap");
134 GetOStream(
Parameters0) <<
"lightweight wrap = " << doExperimentalWrap << std::endl;
147 if (doExperimentalWrap) {
148 std::string algo = pL.get<std::string>(
"aggregation: drop scheme");
153 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
154 GetOStream(
Runtime0) <<
"algorithm = \"" << algo <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
155 Set(currentLevel,
"Filtering", (threshold != STS::zero()));
157 const typename STS::magnitudeType dirichletThreshold = STS::magnitude(as<SC>(pL.get<
double>(
"aggregation: Dirichlet threshold")));
159 GO numDropped = 0, numTotal = 0;
160 std::string graphType =
"unamalgamated";
161 if (algo ==
"classical") {
162 if (predrop_ == null) {
167 if (predrop_ != null) {
170 "MueLu::CoalesceFactory::Build: cast to PreDropFunctionConstVal failed.");
173 if (newt != threshold) {
174 GetOStream(
Warnings0) <<
"switching threshold parameter from " << threshold <<
" (list) to " << newt <<
" (user function" << std::endl;
183 if (A->GetFixedBlockSize() == 1 && threshold == STS::zero()) {
191 numTotal = A->getNodeNumEntries();
194 GO numLocalBoundaryNodes = 0;
195 GO numGlobalBoundaryNodes = 0;
196 for (
LO i = 0; i < boundaryNodes.
size(); ++i)
197 if (boundaryNodes[i])
198 numLocalBoundaryNodes++;
200 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
201 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
204 Set(currentLevel,
"DofsPerNode", 1);
205 Set(currentLevel,
"Graph", graph);
207 }
else if (A->GetFixedBlockSize() == 1 && threshold != STS::zero()) {
222 for (
LO row = 0; row < Teuchos::as<LO>(A->getRowMap()->getNodeNumElements()); ++row) {
223 size_t nnz = A->getNumEntriesInLocalRow(row);
226 A->getLocalRowView(row, indices, vals);
234 for (
LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
235 LO col = indices[colID];
238 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedDiagVals[col]*ghostedDiagVals[row]);
239 typename STS::magnitudeType aij = STS::magnitude(vals[colID]*vals[colID]);
241 if (aij > aiiajj || row == col) {
242 columns[realnnz++] = col;
254 boundaryNodes[row] =
true;
256 rows[row+1] = realnnz;
260 numTotal = A->getNodeNumEntries();
265 GO numLocalBoundaryNodes = 0;
266 GO numGlobalBoundaryNodes = 0;
267 for (
LO i = 0; i < boundaryNodes.size(); ++i)
268 if (boundaryNodes[i])
269 numLocalBoundaryNodes++;
271 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
272 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
274 Set(currentLevel,
"Graph", graph);
275 Set(currentLevel,
"DofsPerNode", 1);
277 }
else if (A->GetFixedBlockSize() > 1 && threshold == STS::zero()) {
283 graphType =
"amalgamated";
291 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
292 Array<LO> colTranslation = *(amalInfo->getColTranslation());
295 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
311 LO blkSize = A->GetFixedBlockSize();
313 LO blkPartSize = A->GetFixedBlockSize();
314 if (A->IsView(
"stridedMaps") ==
true) {
318 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
319 blkId = strMap->getStridedBlockId();
321 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
328 for (
LO row = 0; row < numRows; row++) {
338 bool isBoundary =
false;
340 for (
LO j = 0; j < blkPartSize; j++) {
341 if (!pointBoundaryNodes[row*blkPartSize+j]) {
350 MergeRows(*A, row, indicesExtra, colTranslation);
353 indices = indicesExtra;
354 numTotal += indices.
size();
358 LO nnz = indices.
size(), rownnz = 0;
359 for (
LO colID = 0; colID < nnz; colID++) {
360 LO col = indices[colID];
361 columns[realnnz++] = col;
372 amalgBoundaryNodes[row] =
true;
374 rows[row+1] = realnnz;
382 GO numLocalBoundaryNodes = 0;
383 GO numGlobalBoundaryNodes = 0;
385 for (
LO i = 0; i < amalgBoundaryNodes.
size(); ++i)
386 if (amalgBoundaryNodes[i])
387 numLocalBoundaryNodes++;
390 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
391 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
392 <<
" agglomerated Dirichlet nodes" << std::endl;
395 Set(currentLevel,
"Graph", graph);
396 Set(currentLevel,
"DofsPerNode", blkSize);
398 }
else if (A->GetFixedBlockSize() > 1 && threshold != STS::zero()) {
404 graphType =
"amalgamated";
412 Array<LO> rowTranslation = *(amalInfo->getRowTranslation());
413 Array<LO> colTranslation = *(amalInfo->getColTranslation());
416 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
432 LO blkSize = A->GetFixedBlockSize();
434 LO blkPartSize = A->GetFixedBlockSize();
435 if (A->IsView(
"stridedMaps") ==
true) {
439 blkSize = Teuchos::as<const LO>(strMap->getFixedBlockSize());
440 blkId = strMap->getStridedBlockId();
442 blkPartSize = Teuchos::as<LO>(strMap->getStridingData()[blkId]);
453 for (
LO row = 0; row < numRows; row++) {
463 bool isBoundary =
false;
465 for (
LO j = 0; j < blkPartSize; j++) {
466 if (!pointBoundaryNodes[row*blkPartSize+j]) {
475 MergeRowsWithDropping(*A, row, ghostedDiagVals, threshold, indicesExtra, colTranslation);
478 indices = indicesExtra;
479 numTotal += indices.
size();
483 LO nnz = indices.
size(), rownnz = 0;
484 for (
LO colID = 0; colID < nnz; colID++) {
485 LO col = indices[colID];
486 columns[realnnz++] = col;
497 amalgBoundaryNodes[row] =
true;
499 rows[row+1] = realnnz;
507 GO numLocalBoundaryNodes = 0;
508 GO numGlobalBoundaryNodes = 0;
510 for (
LO i = 0; i < amalgBoundaryNodes.
size(); ++i)
511 if (amalgBoundaryNodes[i])
512 numLocalBoundaryNodes++;
515 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
516 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes
517 <<
" agglomerated Dirichlet nodes" << std::endl;
520 Set(currentLevel,
"Graph", graph);
521 Set(currentLevel,
"DofsPerNode", blkSize);
524 }
else if (algo ==
"distance laplacian") {
525 LO blkSize = A->GetFixedBlockSize();
526 GO indexBase = A->getRowMap()->getIndexBase();
532 RCP<dxMV> Coords = Get< RCP<Xpetra::MultiVector<double,LO,GO,NO> > >(currentLevel,
"Coordinates");
541 if ( (blkSize == 1) && (threshold == STS::zero()) ) {
545 graphType=
"unamalgamated";
546 numTotal = A->getNodeNumEntries();
549 GO numLocalBoundaryNodes = 0;
550 GO numGlobalBoundaryNodes = 0;
551 for (
LO i = 0; i < pointBoundaryNodes.size(); ++i)
552 if (pointBoundaryNodes[i])
553 numLocalBoundaryNodes++;
555 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
556 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
559 Set(currentLevel,
"DofsPerNode", blkSize);
560 Set(currentLevel,
"Graph", graph);
576 "Coordinate vector length (" << Coords->getLocalLength() <<
") is incompatible with number of rows in A (" << A->getRowMap()->getNodeNumElements() <<
") by modulo block size ("<< blkSize <<
").");
582 uniqueMap = A->getRowMap();
583 nonUniqueMap = A->getColMap();
584 graphType=
"unamalgamated";
587 uniqueMap = Coords->getMap();
589 "Different index bases for matrix and coordinates");
593 graphType =
"amalgamated";
595 LO numRows = Teuchos::as<LocalOrdinal>(uniqueMap->getNodeNumElements());
600 if (threshold != STS::zero()) {
605 importer = ImportFactory::Build(uniqueMap, nonUniqueMap);
611 RCP<Vector> localLaplDiag = VectorFactory::Build(uniqueMap);
612 ArrayRCP<SC> localLaplDiagData = localLaplDiag->getDataNonConst(0);
614 for (
LO row = 0; row < numRows; row++) {
620 A->getLocalRowView(row, indices, vals);
624 MergeRows(*A, row, indicesExtra, colTranslation);
625 indices = indicesExtra;
629 for (
LO colID = 0; colID < nnz; colID++) {
630 LO col = indices[colID];
636 ghostedLaplDiag = VectorFactory::Build(nonUniqueMap);
637 ghostedLaplDiag->doImport(*localLaplDiag, *importer,
Xpetra::INSERT);
638 ghostedLaplDiagData = ghostedLaplDiag->getDataNonConst(0);
641 GetOStream(
Runtime0) <<
"Skipping distance laplacian construction due to 0 threshold" << std::endl;
655 for (
LO row = 0; row < numRows; row++) {
661 A->getLocalRowView(row, indices, vals);
665 bool isBoundary =
false;
667 for (
LO j = 0; j < blkSize; j++) {
668 if (!pointBoundaryNodes[row*blkSize+j]) {
676 MergeRows(*A, row, indicesExtra, colTranslation);
679 indices = indicesExtra;
681 numTotal += indices.
size();
683 LO nnz = indices.
size(), rownnz = 0;
684 if (threshold != STS::zero()) {
685 for (
LO colID = 0; colID < nnz; colID++) {
686 LO col = indices[colID];
689 columns[realnnz++] = col;
695 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold * ghostedLaplDiagData[row]*ghostedLaplDiagData[col]);
696 typename STS::magnitudeType aij = STS::magnitude(laplVal*laplVal);
699 columns[realnnz++] = col;
708 for (
LO colID = 0; colID < nnz; colID++) {
709 LO col = indices[colID];
710 columns[realnnz++] = col;
722 amalgBoundaryNodes[row] =
true;
724 rows[row+1] = realnnz;
732 GO numLocalBoundaryNodes = 0;
733 GO numGlobalBoundaryNodes = 0;
735 for (
LO i = 0; i < amalgBoundaryNodes.
size(); ++i)
736 if (amalgBoundaryNodes[i])
737 numLocalBoundaryNodes++;
740 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
741 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" agglomerated Dirichlet nodes" 742 <<
" using threshold " << dirichletThreshold << std::endl;
745 Set(currentLevel,
"Graph", graph);
746 Set(currentLevel,
"DofsPerNode", blkSize);
750 if ((GetVerbLevel() &
Statistics1) && !(A->GetFixedBlockSize() > 1 && threshold != STS::zero())) {
752 GO numGlobalTotal, numGlobalDropped;
755 GetOStream(
Statistics1) <<
"Number of dropped entries in " << graphType <<
" matrix graph: " << numGlobalDropped <<
"/" << numGlobalTotal;
756 if (numGlobalTotal != 0)
757 GetOStream(
Statistics1) <<
" (" << 100*Teuchos::as<double>(numGlobalDropped)/Teuchos::as<double>(numGlobalTotal) <<
"%)";
765 SC threshold = as<SC>(pL.get<
double>(
"aggregation: drop tol"));
767 GetOStream(
Runtime0) <<
"algorithm = \"" <<
"failsafe" <<
"\": threshold = " << threshold <<
", blocksize = " << A->GetFixedBlockSize() << std::endl;
768 Set(currentLevel,
"Filtering", (threshold != STS::zero()));
774 GO indexBase = rowMap->getIndexBase();
778 if(A->IsView(
"stridedMaps") &&
779 Teuchos::rcp_dynamic_cast<
const StridedMap>(A->getRowMap(
"stridedMaps")) != Teuchos::null) {
783 blockdim = strMap->getFixedBlockSize();
784 offset = strMap->getOffset();
785 oldView = A->SwitchToView(oldView);
786 GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build():" <<
" found blockdim=" << blockdim <<
" from strided maps. offset=" << offset << std::endl;
787 }
else GetOStream(
Statistics1) <<
"CoalesceDropFactory::Build(): no striding information available. Use blockdim=1 with offset=0" << std::endl;
792 GetOStream(
Statistics1) <<
"CoalesceDropFactory: nodeMap " << nodeMap->getNodeNumElements() <<
"/" << nodeMap->getGlobalNumElements() <<
" elements" << std::endl;
797 LO numRows = A->getRowMap()->getNodeNumElements();
798 LO numNodes = nodeMap->getNodeNumElements();
801 bool bIsDiagonalEntry =
false;
806 for(
LO row=0; row<numRows; row++) {
808 GO grid = rowMap->getGlobalElement(row);
811 bIsDiagonalEntry =
false;
816 size_t nnz = A->getNumEntriesInLocalRow(row);
819 A->getLocalRowView(row, indices, vals);
823 for(
LO col=0; col<Teuchos::as<LO>(nnz); col++) {
824 GO gcid = colMap->getGlobalElement(indices[col]);
828 cnodeIds->push_back(cnodeId);
830 if (grid == gcid) bIsDiagonalEntry =
true;
834 if(realnnz == 1 && bIsDiagonalEntry ==
true) {
835 LO lNodeId = nodeMap->getLocalElement(nodeId);
836 numberDirichletRowsPerNode[lNodeId] += 1;
837 if (numberDirichletRowsPerNode[lNodeId] == blockdim)
838 amalgBoundaryNodes[lNodeId] =
true;
843 if(arr_cnodeIds.
size() > 0 )
844 crsGraph->insertGlobalIndices(nodeId, arr_cnodeIds());
847 crsGraph->fillComplete(nodeMap,nodeMap);
856 GO numLocalBoundaryNodes = 0;
857 GO numGlobalBoundaryNodes = 0;
858 for (
LO i = 0; i < amalgBoundaryNodes.
size(); ++i)
859 if (amalgBoundaryNodes[i])
860 numLocalBoundaryNodes++;
862 MueLu_sumAll(comm, numLocalBoundaryNodes, numGlobalBoundaryNodes);
863 GetOStream(
Statistics1) <<
"Detected " << numGlobalBoundaryNodes <<
" Dirichlet nodes" << std::endl;
868 Set(currentLevel,
"DofsPerNode", blockdim);
869 Set(currentLevel,
"Graph", graph);
875 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
880 LO blkSize = A.GetFixedBlockSize();
881 if (A.IsView(
"stridedMaps") ==
true) {
885 if (strMap->getStridedBlockId() > -1)
886 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
890 size_t nnz = 0, pos = 0;
891 for (
LO j = 0; j < blkSize; j++)
892 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
904 for (
LO j = 0; j < blkSize; j++) {
905 A.getLocalRowView(row*blkSize+j, inds, vals);
906 size_type numIndices = inds.
size();
912 cols[pos++] = translation[inds[0]];
913 for (size_type k = 1; k < numIndices; k++) {
914 LO nodeID = translation[inds[k]];
918 if (nodeID != cols[pos-1])
919 cols[pos++] = nodeID;
926 std::sort(cols.
begin(), cols.
end());
928 for (
size_t j = 1; j < nnz; j++)
929 if (cols[j] != cols[pos])
930 cols[++pos] = cols[j];
934 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
940 LO blkSize = A.GetFixedBlockSize();
941 if (A.IsView(
"stridedMaps") ==
true) {
945 if (strMap->getStridedBlockId() > -1)
946 blkSize = Teuchos::as<LO>(strMap->getStridingData()[strMap->getStridedBlockId()]);
950 size_t nnz = 0, pos = 0;
951 for (
LO j = 0; j < blkSize; j++)
952 nnz += A.getNumEntriesInLocalRow(row*blkSize+j);
964 for (
LO j = 0; j < blkSize; j++) {
965 A.getLocalRowView(row*blkSize+j, inds, vals);
966 size_type numIndices = inds.
size();
973 for (size_type k = 0; k < numIndices; k++) {
975 LO nodeID = translation[inds[k]];
978 typename STS::magnitudeType aiiajj = STS::magnitude(threshold*threshold*ghostedDiagVals[dofID]*ghostedDiagVals[row*blkSize+j]);
979 typename STS::magnitudeType aij = STS::magnitude(vals[k]*vals[k]);
982 if (aij > aiiajj || (row*blkSize+j == dofID)) {
988 if (nodeID != prevNodeID) {
989 cols[pos++] = nodeID;
999 std::sort(cols.
begin(), cols.
end());
1001 for (
size_t j = 1; j < nnz; j++)
1002 if (cols[j] != cols[pos])
1003 cols[++pos] = cols[j];
1010 #endif // MUELU_COALESCEDROPFACTORY_DEF_HPP Important warning messages (one line)
Exception indicating invalid cast attempted.
#define MueLu_sumAll(rcpComm, in, out)
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)
MueLu representation of a compressed row storage graph.
static Teuchos::ScalarTraits< Scalar >::magnitudeType Distance2(const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &v, LocalOrdinal i0, LocalOrdinal i1)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
One-liner description of what is happening.
Namespace for MueLu classes and methods.
Exception throws to report incompatible objects (like maps).
static const GlobalOrdinal DOFGid2NodeId(GlobalOrdinal gid, LocalOrdinal blockSize, const GlobalOrdinal offset, const GlobalOrdinal indexBase)
translate global (row/column) id to global amalgamation block id
virtual void SetBoundaryNodeMap(const ArrayRCP< const bool > &boundaryArray)=0
void resize(const size_type n, const T &val=T())
Scalar GetThreshold() const
Return threshold value.
CoalesceDropFactory()
Constructor.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void DeclareInput(Level ¤tLevel) const
Input.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
#define SET_VALID_ENTRY(name)
static void AmalgamateMap(const Map &sourceMap, const Matrix &A, RCP< const Map > &amalgamatedMap, Array< LO > &translation)
Method to create merged map for systems of PDEs.
void resize(size_type new_size, const value_type &x=value_type())
void MergeRowsWithDropping(const Matrix &A, const LO row, const ArrayRCP< const SC > &ghostedDiagVals, SC threshold, Array< LO > &cols, const Array< LO > &translation) const
void Build(Level ¤tLevel) const
Build an object with this factory.
void 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.
static RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetMatrixOverlappedDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Lightweight MueLu representation of a compressed row storage graph.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node > > &map, size_t NumVectors, bool zeroOut=true)
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.))
Exception throws to report errors in the internal logical of the program.
ParameterEntry & getEntry(const std::string &name)