46 #ifndef MUELU_INVERSEAPPROXIMATIONFACTORY_DEF_HPP_
47 #define MUELU_INVERSEAPPROXIMATIONFACTORY_DEF_HPP_
52 #include <Xpetra_CrsMatrixWrap.hpp>
54 #include <Xpetra_MultiVectorFactory.hpp>
67 #include "MueLu_Utilities.hpp"
72 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
79 validParamList->
set<std::string>(
"inverse: approximation type",
"diagonal",
"Method used to approximate the inverse.");
80 validParamList->
set<Magnitude>(
"inverse: drop tolerance", 0.0,
"Values below this threshold are dropped from the matrix (or fixed if diagonal fixing is active).");
81 validParamList->
set<
bool>(
"inverse: fixing",
false,
"Keep diagonal and fix small entries with 1.0");
83 return validParamList;
86 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
88 Input(currentLevel,
"A");
91 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 const SC one = STS::one();
100 const bool fixing = pL.
get<
bool>(
"inverse: fixing");
103 const std::string method = pL.
get<std::string>(
"inverse: approximation type");
105 "MueLu::InverseApproximationFactory::Build: Approximation type can be 'diagonal' or 'lumping' or "
106 "'sparseapproxinverse'.");
108 RCP<Matrix> A = Get<RCP<Matrix>>(currentLevel,
"A");
110 const bool isBlocked = (bA == Teuchos::null ?
false :
true);
113 if (isBlocked) A = bA->getMatrix(0, 0);
115 const Magnitude tol = pL.
get<Magnitude>(
"inverse: drop tolerance");
118 if (method ==
"diagonal") {
119 const auto diag = VectorFactory::Build(A->getRangeMap(),
true);
120 A->getLocalDiagCopy(*diag);
122 Ainv = MatrixFactory::Build(D);
123 }
else if (method ==
"lumping") {
126 Ainv = MatrixFactory::Build(D);
127 }
else if (method ==
"sparseapproxinverse") {
129 GetOStream(
Statistics1) <<
"NNZ Graph(A): " << A->getCrsGraph()->getGlobalNumEntries() <<
" , NNZ Tresholded Graph(A): " << sparsityPattern->getGlobalNumEntries() << std::endl;
130 RCP<Matrix> pAinv = GetSparseInverse(A, sparsityPattern);
132 GetOStream(
Statistics1) <<
"NNZ Ainv: " << pAinv->getGlobalNumEntries() <<
", NNZ Tresholded Ainv (parameter: " << tol <<
"): " << Ainv->getGlobalNumEntries() << std::endl;
135 GetOStream(
Statistics1) <<
"Approximate inverse calculated by: " << method <<
"." << std::endl;
136 GetOStream(
Statistics1) <<
"Ainv has " << Ainv->getGlobalNumRows() <<
"x" << Ainv->getGlobalNumCols() <<
" rows and columns." << std::endl;
138 Set(currentLevel,
"Ainv", Ainv);
141 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
145 RCP<Matrix> Ainv = MatrixFactory::Build(sparsityPattern);
149 RCP<Import> rowImport = ImportFactory::Build(sparsityPattern->getRowMap(), sparsityPattern->getColMap());
150 RCP<Matrix> A = MatrixFactory::Build(Aorg, *rowImport);
153 for (
size_t k = 0; k < sparsityPattern->getLocalNumRows(); k++) {
156 sparsityPattern->getLocalRowView(k, Ik);
162 for (
LO i = 0; i < Ik.
size(); i++) {
163 A->getLocalRowView(Ik[i], J[i], Ak[i]);
164 for (
LO j = 0; j < J[i].size(); j++)
172 for (
LO i = 0; i < Jk.
size(); i++) G.insert(std::pair<LO, LO>(Jk[i], i));
176 for (
LO i = 0; i < Ik.
size(); i++) {
177 for (
LO j = 0; j < J[i].size(); j++) {
178 localA(G.at(J[i][j]), i) = Ak[i][j];
193 const int err = qrSolver.solve();
195 "MueLu::InverseApproximationFactory::GetSparseInverse: Error in serial QR solve.");
199 Ainv->replaceLocalValues(k, Ik, Mk);
201 Ainv->fillComplete();
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static RCP< CrsMatrixWrap > GetThresholdedMatrix(const RCP< Matrix > &Ain, const Scalar threshold, const bool keepDiagonal=true, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a matrix.
Array< T > & append(const T &x)
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)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
static RCP< Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > GetThresholdedGraph(const RCP< Matrix > &A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow=-1)
Threshold a graph.
iterator erase(iterator position)
void Build(Level ¤tLevel) const
Build an object with this factory.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Class that holds all level-specific information.
RCP< Matrix > GetSparseInverse(const RCP< Matrix > &A, const RCP< const CrsGraph > &sparsityPattern) const
Sparse inverse calculation method.
void DeclareInput(Level ¤tLevel) const
Input.
Exception throws to report errors in the internal logical of the program.
static Teuchos::RCP< Vector > GetInverse(Teuchos::RCP< const Vector > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps()*100, Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero())
Return vector containing inverse of input vector.
static Teuchos::RCP< Vector > GetLumpedMatrixDiagonal(Matrix const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::magnitude(Teuchos::ScalarTraits< Scalar >::zero()), Scalar valReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
Extract Matrix Diagonal of lumped matrix.
static const RCP< const NoFactory > getRCP()
Static Get() functions.