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)
64 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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<CrsOMatrix> Aout = rcp(
new CrsOMatrix(Ain->getRowMap(),Ain->getGlobalMaxNumRowEntries(),Xpetra::StaticProfile));
87 for(
size_t row=0; row<Ain->getNodeNumRows(); row++)
89 size_t nnz = Ain->getNumEntriesInLocalRow(row);
91 Teuchos::ArrayView<const LocalOrdinal> indices;
92 Teuchos::ArrayView<const Scalar> vals;
93 Ain->getLocalRowView(row, indices, vals);
95 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.");
97 Teuchos::ArrayRCP<GlobalOrdinal> indout(indices.size(),Teuchos::ScalarTraits<GlobalOrdinal>::zero());
98 Teuchos::ArrayRCP<Scalar> valout(indices.size(),Teuchos::ScalarTraits<Scalar>::zero());
100 for(
size_t i=0; i<(size_t)indices.size(); i++) {
101 if(Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > Teuchos::ScalarTraits<Scalar>::magnitude(
threshold_) || indices[i]==(LocalOrdinal)row) {
102 indout[nNonzeros] = Ain->getColMap()->getGlobalElement(indices[i]);
103 valout[nNonzeros] = vals[i];
108 indout.resize(nNonzeros);
109 valout.resize(nNonzeros);
111 Aout->insertGlobalValues(Ain->getRowMap()->getGlobalElement(row), indout.view(0,indout.size()), valout.view(0,valout.size()));
114 Aout->fillComplete(Ain->getDomainMap(), Ain->getRangeMap());
118 currentLevel.Set(
varName_, Teuchos::rcp_dynamic_cast<OMatrix>(Aout),
this);
123 #endif // MUELU_THRESHOLDAFILTERFACTORY_DEF_HPP void DeclareInput(Level ¤tLevel) const
Input.
Timer to be used in factories. Similar to Monitor but with additional timers.
Namespace for MueLu classes and methods.
const Scalar threshold_
threshold parameter
Print statistics that do not involve significant additional computation.
Class that holds all level-specific information.
std::string varName_
name of input and output variable
virtual ~ThresholdAFilterFactory()
Destructor.
ThresholdAFilterFactory(const std::string &ename, const Scalar threshold)
Constructor.
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
void Build(Level ¤tLevel) const
Build an object with this factory.
Exception throws to report errors in the internal logical of the program.
void Input(Level &level, const std::string &varName) const