MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_InverseApproximationFactory_def.hpp
Go to the documentation of this file.
1// @HEADER
2//
3// ***********************************************************************
4//
5// MueLu: A package for multigrid based preconditioning
6// Copyright 2012 Sandia Corporation
7//
8// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
9// the U.S. Government retains certain rights in this software.
10//
11// Redistribution and use in source and binary forms, with or without
12// modification, are permitted provided that the following conditions are
13// met:
14//
15// 1. Redistributions of source code must retain the above copyright
16// notice, this list of conditions and the following disclaimer.
17//
18// 2. Redistributions in binary form must reproduce the above copyright
19// notice, this list of conditions and the following disclaimer in the
20// documentation and/or other materials provided with the distribution.
21//
22// 3. Neither the name of the Corporation nor the names of the
23// contributors may be used to endorse or promote products derived from
24// this software without specific prior written permission.
25//
26// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37//
38// Questions? Contact
39// Jonathan Hu (jhu@sandia.gov)
40// Andrey Prokopenko (aprokop@sandia.gov)
41// Ray Tuminaro (rstumin@sandia.gov)
42//
43// ***********************************************************************
44//
45// @HEADER
46#ifndef MUELU_INVERSEAPPROXIMATIONFACTORY_DEF_HPP_
47#define MUELU_INVERSEAPPROXIMATIONFACTORY_DEF_HPP_
48
49#include <Xpetra_BlockedCrsMatrix.hpp>
50#include <Xpetra_CrsGraph.hpp>
51#include <Xpetra_CrsGraphFactory.hpp>
52#include <Xpetra_CrsMatrixWrap.hpp>
53#include <Xpetra_CrsMatrix.hpp>
54#include <Xpetra_MultiVectorFactory.hpp>
55#include <Xpetra_VectorFactory.hpp>
56#include <Xpetra_MatrixFactory.hpp>
57#include <Xpetra_Matrix.hpp>
58#include <Xpetra_MatrixMatrix.hpp>
59#include <Xpetra_TripleMatrixMultiply.hpp>
60
61#include <Teuchos_SerialDenseVector.hpp>
62#include <Teuchos_SerialDenseMatrix.hpp>
63#include <Teuchos_SerialQRDenseSolver.hpp>
64
65#include "MueLu_Level.hpp"
66#include "MueLu_Monitor.hpp"
67#include "MueLu_Utilities.hpp"
69
70namespace MueLu {
71
72 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
74 RCP<ParameterList> validParamList = rcp(new ParameterList());
75 using Magnitude = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
76
77 validParamList->set<RCP<const FactoryBase> >("A", NoFactory::getRCP(), "Matrix to build the approximate inverse on.\n");
78
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");
82
83 return validParamList;
84 }
85
86 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
88 Input(currentLevel, "A");
89 }
90
91 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
93 FactoryMonitor m(*this, "Build", currentLevel);
94
95 using STS = Teuchos::ScalarTraits<SC>;
96 const SC one = STS::one();
97 using Magnitude = typename Teuchos::ScalarTraits<Scalar>::magnitudeType;
98
99 const ParameterList& pL = GetParameterList();
100 const bool fixing = pL.get<bool>("inverse: fixing");
101
102 // check which approximation type to use
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'.");
107
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);
111
112 // if blocked operator is used, defaults to A(0,0)
113 if(isBlocked) A = bA->getMatrix(0,0);
114
115 const Magnitude tol = pL.get<Magnitude>("inverse: drop tolerance");
116 RCP<Matrix> Ainv = Teuchos::null;
117
118 if(method=="diagonal")
119 {
120 const auto diag = VectorFactory::Build(A->getRangeMap(), true);
121 A->getLocalDiagCopy(*diag);
122 const RCP<const Vector> D = (!fixing ? Utilities::GetInverse(diag) : Utilities::GetInverse(diag, tol, one));
123 Ainv = MatrixFactory::Build(D);
124 }
125 else if(method=="lumping")
126 {
127 const auto diag = Utilities::GetLumpedMatrixDiagonal(*A);
128 const RCP<const Vector> D = (!fixing ? Utilities::GetInverse(diag) : Utilities::GetInverse(diag, tol, one));
129 Ainv = MatrixFactory::Build(D);
130 }
131 else if(method=="sparseapproxinverse")
132 {
133 RCP<CrsGraph> sparsityPattern = Utilities::GetThresholdedGraph(A, tol, A->getGlobalMaxNumRowEntries());
134 GetOStream(Statistics1) << "NNZ Graph(A): " << A->getCrsGraph()->getGlobalNumEntries() << " , NNZ Tresholded Graph(A): " << sparsityPattern->getGlobalNumEntries() << std::endl;
135 RCP<Matrix> pAinv = GetSparseInverse(A, sparsityPattern);
136 Ainv = Utilities::GetThresholdedMatrix(pAinv, tol, fixing, pAinv->getGlobalMaxNumRowEntries());
137 GetOStream(Statistics1) << "NNZ Ainv: " << pAinv->getGlobalNumEntries() << ", NNZ Tresholded Ainv (parameter: " << tol << "): " << Ainv->getGlobalNumEntries() << std::endl;
138 }
139
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;
142
143 Set(currentLevel, "Ainv", Ainv);
144 }
145
146 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
147 RCP<Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>>
148 InverseApproximationFactory<Scalar, LocalOrdinal, GlobalOrdinal, Node>::GetSparseInverse(const RCP<Matrix>& Aorg, const RCP<const CrsGraph>& sparsityPattern) const {
149
150 // construct the inverse matrix with the given sparsity pattern
151 RCP<Matrix> Ainv = MatrixFactory::Build(sparsityPattern);
152 Ainv->resumeFill();
153
154 // gather missing rows from other procs to generate an overlapping map
155 RCP<Import> rowImport = ImportFactory::Build(sparsityPattern->getRowMap(), sparsityPattern->getColMap());
156 RCP<Matrix> A = MatrixFactory::Build(Aorg, *rowImport);
157
158 // loop over all rows of the inverse sparsity pattern (this can be done in parallel)
159 for(size_t k=0; k<sparsityPattern->getLocalNumRows(); k++) {
160
161 // 1. get column indices Ik of local row k
162 ArrayView<const LO> Ik;
163 sparsityPattern->getLocalRowView(k, Ik);
164
165 // 2. get all local A(Ik,:) rows
166 Array<ArrayView<const LO>> J(Ik.size());
167 Array<ArrayView<const SC>> Ak(Ik.size());
168 Array<LO> Jk;
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++)
172 Jk.append(J[i][j]);
173 }
174 // set of unique column indices Jk
175 std::sort(Jk.begin(), Jk.end());
176 Jk.erase(std::unique(Jk.begin(), Jk.end()), Jk.end());
177 // create map
178 std::map<LO, LO> G;
179 for (LO i = 0; i < Jk.size(); i++) G.insert(std::pair<LO, LO>(Jk[i], i));
180
181 // 3. merge rows together
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];
186 }
187 }
188
189 // 4. get direction-vector
190 // diagonal needs an entry!
191 Teuchos::SerialDenseVector<LO, SC> ek(Jk.size(), true);
192 ek[std::find(Jk.begin(), Jk.end(), k) - Jk.begin()] = Teuchos::ScalarTraits<Scalar>::one();;
193
194 // 5. solve linear system for x
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();
200 TEUCHOS_TEST_FOR_EXCEPTION(err != 0, Exceptions::RuntimeError,
201 "MueLu::InverseApproximationFactory::GetSparseInverse: Error in serial QR solve.");
202
203 // 6. set calculated row into Ainv
204 ArrayView<const SC> Mk(localX.values(), localX.length());
205 Ainv->replaceLocalValues(k, Ik, Mk);
206
207 }
208 Ainv->fillComplete();
209
210 return Ainv;
211 }
212
213} // namespace MueLu
214
215#endif /* MUELU_INVERSEAPPROXIMATIONFACTORY_DEF_HPP_ */
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
void Build(Level &currentLevel) const
Build an object with this factory.
void DeclareInput(Level &currentLevel) const
Input.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
RCP< Matrix > GetSparseInverse(const RCP< Matrix > &A, const RCP< const CrsGraph > &sparsityPattern) const
Sparse inverse calculation method.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
static const RCP< const NoFactory > getRCP()
Static Get() functions.
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())
Namespace for MueLu classes and methods.
@ Statistics1
Print more statistics.