46 #ifndef MUELU_THRESHOLDAFILTERFACTORY_DEF_HPP 47 #define MUELU_THRESHOLDAFILTERFACTORY_DEF_HPP 49 #include <Xpetra_Matrix.hpp> 50 #include <Xpetra_CrsMatrixWrap.hpp> 59 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
61 : varName_(ename), threshold_(threshold), keepDiagonal_(keepDiagonal), expectedNNZperRow_(expectedNNZperRow)
64 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
69 Input(currentLevel, varName_);
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
78 typedef Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> OMatrix;
79 typedef Xpetra::CrsMatrixWrap<Scalar, LocalOrdinal, GlobalOrdinal, Node> CrsOMatrix;
81 RCP<OMatrix> Ain = Get< RCP<OMatrix> >(currentLevel, varName_);
84 RCP<const Map> rowmap = Ain->getRowMap();
85 RCP<const Map> colmap = Ain->getColMap();
86 RCP<CrsOMatrix> Aout = rcp(
new CrsOMatrix(rowmap, expectedNNZperRow_ <= 0 ? Ain->getGlobalMaxNumRowEntries() : expectedNNZperRow_));
88 for(
size_t row=0; row<Ain->getNodeNumRows(); row++)
90 size_t nnz = Ain->getNumEntriesInLocalRow(row);
92 Teuchos::ArrayView<const LocalOrdinal> indices;
93 Teuchos::ArrayView<const Scalar> vals;
94 Ain->getLocalRowView(row, indices, vals);
96 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<size_t>(indices.size()) != nnz,
Exceptions::RuntimeError,
"MueLu::ThresholdAFilterFactory::Build: number of nonzeros not equal to number of indices? Error.");
98 Teuchos::ArrayRCP<GlobalOrdinal> indout(indices.size(),Teuchos::ScalarTraits<GlobalOrdinal>::zero());
99 Teuchos::ArrayRCP<Scalar> valout(indices.size(),Teuchos::ScalarTraits<Scalar>::zero());
100 size_t nNonzeros = 0;
103 LocalOrdinal lclColIdx = colmap->getLocalElement(glbRow);
104 for(
size_t i=0; i<(size_t)indices.size(); i++) {
105 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(threshold_) || indices[i]==lclColIdx) {
106 indout[nNonzeros] = colmap->getGlobalElement(indices[i]);
107 valout[nNonzeros] = vals[i];
112 for(
size_t i=0; i<(size_t)indices.size(); i++) {
113 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(threshold_)) {
114 indout[nNonzeros] = colmap->getGlobalElement(indices[i]);
115 valout[nNonzeros] = vals[i];
120 indout.resize(nNonzeros);
121 valout.resize(nNonzeros);
123 Aout->insertGlobalValues(Ain->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
126 Aout->fillComplete(Ain->getDomainMap(), Ain->getRangeMap());
128 GetOStream(
Statistics0) <<
"Nonzeros in " << varName_ <<
"(input): " << Ain->getGlobalNumEntries() <<
", Nonzeros after filtering " << varName_ <<
" (parameter: " << threshold_ <<
"): " << Aout->getGlobalNumEntries() << std::endl;
130 currentLevel.Set(varName_, Teuchos::rcp_dynamic_cast<OMatrix>(Aout),
this);
135 #endif // MUELU_THRESHOLDAFILTERFACTORY_DEF_HPP MueLu::DefaultLocalOrdinal LocalOrdinal
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
Print statistics that do not involve significant additional computation.
ThresholdAFilterFactory(const std::string &ename, const Scalar threshold, const bool keepDiagonal=true, const GlobalOrdinal expectedNNZperRow=-1)
Constructor.
void Build(Level ¤tLevel) const
Build an object with this factory.
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
virtual ~ThresholdAFilterFactory()
Destructor.
void DeclareInput(Level ¤tLevel) const
Input.
Exception throws to report errors in the internal logical of the program.