95 using STS = Teuchos::ScalarTraits<SC>;
96 const SC one = STS::one();
97 using Magnitude =
typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
99 const ParameterList& pL = GetParameterList();
100 const bool fixing = pL.get<
bool>(
"inverse: fixing");
103 const std::string method = pL.get<std::string>(
"inverse: approximation type");
104 TEUCHOS_TEST_FOR_EXCEPTION(method !=
"diagonal" && method !=
"lumping" && method !=
"sparseapproxinverse",
Exceptions::RuntimeError,
105 "MueLu::InverseApproximationFactory::Build: Approximation type can be 'diagonal' or 'lumping' or "
106 "'sparseapproxinverse'.");
108 RCP<Matrix> A = Get<RCP<Matrix> >(currentLevel,
"A");
109 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(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");
116 RCP<Matrix> Ainv = Teuchos::null;
118 if(method==
"diagonal")
120 const auto diag = VectorFactory::Build(A->getRangeMap(),
true);
121 A->getLocalDiagCopy(*diag);
123 Ainv = MatrixFactory::Build(D);
125 else if(method==
"lumping")
129 Ainv = MatrixFactory::Build(D);
131 else if(method==
"sparseapproxinverse")
134 GetOStream(
Statistics1) <<
"NNZ Graph(A): " << A->getCrsGraph()->getGlobalNumEntries() <<
" , NNZ Tresholded Graph(A): " << sparsityPattern->getGlobalNumEntries() << std::endl;
135 RCP<Matrix> pAinv = GetSparseInverse(A, sparsityPattern);
137 GetOStream(
Statistics1) <<
"NNZ Ainv: " << pAinv->getGlobalNumEntries() <<
", NNZ Tresholded Ainv (parameter: " << tol <<
"): " << Ainv->getGlobalNumEntries() << std::endl;
140 GetOStream(
Statistics1) <<
"Approximate inverse calculated by: " << method <<
"." << std::endl;
141 GetOStream(
Statistics1) <<
"Ainv has " << Ainv->getGlobalNumRows() <<
"x" << Ainv->getGlobalNumCols() <<
" rows and columns." << std::endl;
143 Set(currentLevel,
"Ainv", Ainv);
147 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
151 RCP<Matrix> Ainv = MatrixFactory::Build(sparsityPattern);
155 RCP<Import> rowImport = ImportFactory::Build(sparsityPattern->getRowMap(), sparsityPattern->getColMap());
156 RCP<Matrix> A = MatrixFactory::Build(Aorg, *rowImport);
159 for(
size_t k=0; k<sparsityPattern->getLocalNumRows(); k++) {
162 ArrayView<const LO> Ik;
163 sparsityPattern->getLocalRowView(k, Ik);
166 Array<ArrayView<const LO>> J(Ik.size());
167 Array<ArrayView<const SC>> Ak(Ik.size());
169 for (LO i = 0; i < Ik.size(); i++) {
170 A->getLocalRowView(Ik[i], J[i], Ak[i]);
171 for (LO j = 0; j < J[i].size(); j++)
175 std::sort(Jk.begin(), Jk.end());
176 Jk.erase(std::unique(Jk.begin(), Jk.end()), Jk.end());
179 for (LO i = 0; i < Jk.size(); i++) G.insert(std::pair<LO, LO>(Jk[i], i));
182 Teuchos::SerialDenseMatrix<LO, SC> localA(Jk.size(), Ik.size(),
true);
183 for (LO i = 0; i < Ik.size(); i++) {
184 for (LO j = 0; j < J[i].size(); j++) {
185 localA(G.at(J[i][j]), i) = Ak[i][j];
191 Teuchos::SerialDenseVector<LO, SC> ek(Jk.size(),
true);
192 ek[std::find(Jk.begin(), Jk.end(), k) - Jk.begin()] = Teuchos::ScalarTraits<Scalar>::one();;
195 Teuchos::SerialDenseVector<LO, SC> localX(Ik.size());
196 Teuchos::SerialQRDenseSolver<LO, SC> qrSolver;
197 qrSolver.setMatrix(Teuchos::rcp(&localA,
false));
198 qrSolver.setVectors(Teuchos::rcp(&localX,
false), Teuchos::rcp(&ek,
false));
199 const int err = qrSolver.solve();
201 "MueLu::InverseApproximationFactory::GetSparseInverse: Error in serial QR solve.");
204 ArrayView<const SC> Mk(localX.values(), localX.length());
205 Ainv->replaceLocalValues(k, Ik, Mk);
208 Ainv->fillComplete();
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetLumpedMatrixDiagonal(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > const &A, const bool doReciprocal=false, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero(), const bool replaceSingleEntryRowWithZero=false, const bool useAverageAbsDiagVal=false)
static RCP< Xpetra::CrsGraph< LocalOrdinal, GlobalOrdinal, Node > > GetThresholdedGraph(const RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Magnitude threshold, const GlobalOrdinal expectedNNZperRow)
static RCP< Xpetra::CrsMatrixWrap< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetThresholdedMatrix(const RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Scalar threshold, const bool keepDiagonal, const GlobalOrdinal expectedNNZperRow)
static Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > GetInverse(Teuchos::RCP< const Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > v, Magnitude tol=Teuchos::ScalarTraits< Scalar >::eps() *100, Scalar tolReplacement=Teuchos::ScalarTraits< Scalar >::zero())