MueLu Version of the Day
Loading...
Searching...
No Matches
MueLu_BraessSarazinSmoother_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
47#ifndef MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
48#define MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
49
50#include <Teuchos_ArrayViewDecl.hpp>
51#include <Teuchos_ScalarTraits.hpp>
52
53#include "MueLu_ConfigDefs.hpp"
54
55#include <Xpetra_Matrix.hpp>
56#include <Xpetra_BlockedCrsMatrix.hpp>
57#include <Xpetra_MultiVectorFactory.hpp>
58#include <Xpetra_VectorFactory.hpp>
59#include <Xpetra_ReorderedBlockedCrsMatrix.hpp>
60
62#include "MueLu_Level.hpp"
63#include "MueLu_Utilities.hpp"
64#include "MueLu_Monitor.hpp"
65#include "MueLu_HierarchyUtils.hpp"
67
68#include "MueLu_FactoryManager.hpp"
69
70namespace MueLu {
71
72 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
73 void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddFactoryManager(RCP<const FactoryManagerBase> FactManager, int pos) {
74 TEUCHOS_TEST_FOR_EXCEPTION(pos != 0, Exceptions::RuntimeError, "MueLu::BraessSarazinSmoother::AddFactoryManager: parameter \'pos\' must be zero! error.");
75 FactManager_ = FactManager;
76 }
77
78 template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
80 RCP<ParameterList> validParamList = rcp(new ParameterList());
81
82 SC one = Teuchos::ScalarTraits<SC>::one();
83
84 validParamList->set<RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
85
86 validParamList->set<bool>("lumping", false, "Use lumping to construct diag(A(0,0)), "
87 "i.e. use row sum of the abs values on the diagonal as approximation of A00 (and A00^{-1})");
88 validParamList->set<SC>("Damping factor", one, "Damping/Scaling factor in BraessSarazin (usually has to be chosen > 1, default = 1.0)");
89 validParamList->set<LO>("Sweeps", 1, "Number of BraessSarazin sweeps (default = 1)");
90 validParamList->set<bool>("q2q1 mode", false, "Run in the mode matching Q2Q1 matlab code");
91
92 return validParamList;
93 }
94
95 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
97 this->Input(currentLevel, "A");
98
99 TEUCHOS_TEST_FOR_EXCEPTION(FactManager_.is_null(), Exceptions::RuntimeError,
100 "MueLu::BraessSarazinSmoother::DeclareInput: FactManager_ must not be null! "
101 "Introduce a FactoryManager for the SchurComplement equation.");
102
103 // carefully call DeclareInput after switching to internal FactoryManager
104 {
105 SetFactoryManager currentSFM(rcpFromRef(currentLevel), FactManager_);
106
107 // request "Smoother" for current subblock row.
108 currentLevel.DeclareInput("PreSmoother", FactManager_->GetFactory("Smoother").get());
109
110 // request Schur matrix just in case
111 currentLevel.DeclareInput("A", FactManager_->GetFactory("A").get());
112 }
113 }
114
115 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
117 FactoryMonitor m(*this, "Setup Smoother", currentLevel);
118
119 if (SmootherPrototype::IsSetup() == true)
120 this->GetOStream(Warnings0) << "MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";
121
122 // Extract blocked operator A from current level
123 A_ = Factory::Get<RCP<Matrix> > (currentLevel, "A");
124 RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
125 TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
126 "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
127
128 // Store map extractors
129 rangeMapExtractor_ = bA->getRangeMapExtractor();
130 domainMapExtractor_ = bA->getDomainMapExtractor();
131
132 // Store the blocks in local member variables
133 A00_ = bA->getMatrix(0,0);
134 A01_ = bA->getMatrix(0,1);
135 A10_ = bA->getMatrix(1,0);
136 A11_ = bA->getMatrix(1,1);
137
138 const ParameterList& pL = Factory::GetParameterList();
139 SC omega = pL.get<SC>("Damping factor");
140
141 // Create the inverse of the diagonal of F
142 // TODO add safety check for zeros on diagonal of F!
143 RCP<Vector> diagFVector = VectorFactory::Build(A00_->getRowMap());
144 if (pL.get<bool>("lumping") == false) {
145 A00_->getLocalDiagCopy(*diagFVector); // extract diagonal of F
146 } else {
147 diagFVector = Utilities::GetLumpedMatrixDiagonal(*A00_);
148 }
149 diagFVector->scale(omega);
150 D_ = Utilities::GetInverse(diagFVector);
151
152 // check whether D_ is a blocked vector with only 1 block
153 RCP<BlockedVector> bD = Teuchos::rcp_dynamic_cast<BlockedVector>(D_);
154 if(bD.is_null() == false && bD->getBlockedMap()->getNumMaps() == 1) {
155 RCP<Vector> nestedVec = bD->getMultiVector(0,bD->getBlockedMap()->getThyraMode())->getVectorNonConst(0);
156 D_.swap(nestedVec);
157 }
158
159 // Set the Smoother
160 // carefully switch to the SubFactoryManagers (defined by the users)
161 {
162 SetFactoryManager currentSFM(rcpFromRef(currentLevel), FactManager_);
163 smoo_ = currentLevel.Get<RCP<SmootherBase> >("PreSmoother", FactManager_->GetFactory("Smoother").get());
164 S_ = currentLevel.Get<RCP<Matrix> > ("A", FactManager_->GetFactory("A").get());
165 }
166
168 }
169
170 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
171 void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
172 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError,
173 "MueLu::BraessSarazinSmoother::Apply(): Setup() has not been called");
174
175 RCP<MultiVector> rcpX = rcpFromRef(X);
176 RCP<const MultiVector> rcpB = rcpFromRef(B);
177
178 // make sure that both rcpX and rcpB are BlockedMultiVector objects
179 bool bCopyResultX = false;
180 RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
181 MUELU_TEST_FOR_EXCEPTION(bA.is_null() == true, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): A_ must be a BlockedCrsMatrix");
182 RCP<BlockedMultiVector> bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
183 RCP<const BlockedMultiVector> bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
184
185 if(bX.is_null() == true) {
186 RCP<MultiVector> test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedDomainMap(),rcpX));
187 rcpX.swap(test);
188 bCopyResultX = true;
189 }
190
191 if(bB.is_null() == true) {
192 RCP<const MultiVector> test = Teuchos::rcp(new BlockedMultiVector(bA->getBlockedRangeMap(),rcpB));
193 rcpB.swap(test);
194 }
195
196 // we now can guarantee that X and B are blocked multi vectors
197 bX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(rcpX);
198 bB = Teuchos::rcp_dynamic_cast<const BlockedMultiVector>(rcpB);
199
200 // check the type of operator
201 RCP<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> > rbA = Teuchos::rcp_dynamic_cast<Xpetra::ReorderedBlockedCrsMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node> >(bA);
202 if(rbA.is_null() == false) {
203 // A is a ReorderedBlockedCrsMatrix
204 Teuchos::RCP<const Xpetra::BlockReorderManager > brm = rbA->getBlockReorderManager();
205
206 // check type of X vector
207 if(bX->getBlockedMap()->getNumMaps() != bA->getDomainMapExtractor()->NumMaps()) {
208 // X is a blocked multi vector but incompatible to the reordered blocked operator A
209 Teuchos::RCP<MultiVector> test =
210 buildReorderedBlockedMultiVector(brm, bX);
211 rcpX.swap(test);
212 }
213 if(bB->getBlockedMap()->getNumMaps() != bA->getRangeMapExtractor()->NumMaps()) {
214 // B is a blocked multi vector but incompatible to the reordered blocked operator A
215 Teuchos::RCP<const MultiVector> test =
216 buildReorderedBlockedMultiVector(brm, bB);
217 rcpB.swap(test);
218 }
219 }
220
221 // use the GIDs of the sub blocks
222 // This is valid as the subblocks actually represent the GIDs (either Thyra, Xpetra or pseudo Xpetra)
223
224 bool bRangeThyraMode = rangeMapExtractor_->getThyraMode();
225 bool bDomainThyraMode = domainMapExtractor_->getThyraMode();
226
227 RCP<MultiVector> deltaX = MultiVectorFactory::Build(rcpX->getMap(), rcpX->getNumVectors());
228 RCP<BlockedMultiVector> bdeltaX = Teuchos::rcp_dynamic_cast<BlockedMultiVector>(deltaX);
229 RCP<MultiVector> deltaX0 = bdeltaX->getMultiVector(0,bDomainThyraMode);
230 RCP<MultiVector> deltaX1 = bdeltaX->getMultiVector(1,bDomainThyraMode);
231
232 RCP<MultiVector> Rtmp = rangeMapExtractor_->getVector(1, rcpB->getNumVectors(), bRangeThyraMode);
233
234
235 typedef Teuchos::ScalarTraits<SC> STS;
236 SC one = STS::one(), zero = STS::zero();
237
238 // extract parameters from internal parameter list
239 const ParameterList& pL = Factory::GetParameterList();
240 LO nSweeps = pL.get<LO>("Sweeps");
241
242 RCP<MultiVector> R;// = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());;
243 if (InitialGuessIsZero) {
244 R = MultiVectorFactory::Build(rcpB->getMap(), rcpB->getNumVectors());
245 R->update(one, *rcpB, zero);
246 } else {
247 R = Utilities::Residual(*A_, *rcpX, *rcpB);
248 }
249
250 // extract diagonal of Schur complement operator
251 RCP<Vector> diagSVector = VectorFactory::Build(S_->getRowMap());
252 S_->getLocalDiagCopy(*diagSVector);
253
254 for (LO run = 0; run < nSweeps; ++run) {
255 // Extract corresponding subvectors from X and R
256 // Automatically detect whether we use Thyra or Xpetra GIDs
257 // The GIDs should be always compatible with the GIDs of A00, A01, etc...
258 RCP<MultiVector> R0 = rangeMapExtractor_ ->ExtractVector(R, 0, bRangeThyraMode);
259 RCP<MultiVector> R1 = rangeMapExtractor_ ->ExtractVector(R, 1, bRangeThyraMode);
260
261 // Calculate Rtmp = R1 - D * deltaX0 (equation 8.14)
262 deltaX0->putScalar(zero);
263 deltaX0->elementWiseMultiply(one, *D_, *R0, zero); // deltaX0 = D * R0 (equation 8.13)
264 A10_->apply(*deltaX0, *Rtmp); // Rtmp = A10*D*deltaX0 (intermediate step)
265 Rtmp->update(one, *R1, -one); // Rtmp = R1 - A10*D*deltaX0
266
267 if (!pL.get<bool>("q2q1 mode")) {
268 deltaX1->putScalar(zero);
269 } else {
270 // special code for q2q1
271 if(Teuchos::rcp_dynamic_cast<BlockedVector>(diagSVector) == Teuchos::null) {
272 ArrayRCP<SC> Sdiag = diagSVector->getDataNonConst(0);
273 ArrayRCP<SC> deltaX1data = deltaX1->getDataNonConst(0);
274 ArrayRCP<SC> Rtmpdata = Rtmp->getDataNonConst(0);
275 for (GO row = 0; row < deltaX1data.size(); row++)
276 deltaX1data[row] = Teuchos::as<SC>(1.1)*Rtmpdata[row] / Sdiag[row];
277 } else {
278 TEUCHOS_TEST_FOR_EXCEPTION(true,MueLu::Exceptions::RuntimeError,"MueLu::BraessSarazinSmoother: q2q1 mode only supported for non-blocked operators.")
279 }
280 }
281
282 smoo_->Apply(*deltaX1,*Rtmp);
283
284 // Compute deltaX0
285 deltaX0->putScalar(zero); // just for safety
286 A01_->apply(*deltaX1, *deltaX0); // deltaX0 = A01*deltaX1
287 deltaX0->update(one, *R0, -one); // deltaX0 = R0 - A01*deltaX1
288 R0.swap(deltaX0);
289 deltaX0->elementWiseMultiply(one, *D_, *R0, zero); // deltaX0 = D*(R0 - A01*deltaX1)
290
291 RCP<MultiVector> X0 = domainMapExtractor_->ExtractVector(rcpX, 0, bDomainThyraMode);
292 RCP<MultiVector> X1 = domainMapExtractor_->ExtractVector(rcpX, 1, bDomainThyraMode);
293
294 // Update solution
295 X0->update(one, *deltaX0, one);
296 X1->update(one, *deltaX1, one);
297
298 domainMapExtractor_->InsertVector(X0, 0, rcpX, bDomainThyraMode);
299 domainMapExtractor_->InsertVector(X1, 1, rcpX, bDomainThyraMode);
300
301 if (run < nSweeps-1) {
302 R = Utilities::Residual(*A_, *rcpX, *rcpB);
303 }
304
305 }
306
307 if (bCopyResultX == true) {
308 RCP<MultiVector> Xmerged = bX->Merge();
309 X.update(one, *Xmerged, zero);
310 }
311
312 }
313
314 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
315 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
319
320 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
322 description () const {
323 std::ostringstream out;
325 out << "{type = " << type_ << "}";
326 return out.str();
327 }
328
329 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
330 void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
332
333 if (verbLevel & Parameters0) {
334 out0 << "Prec. type: " << type_ << /*" Sweeps: " << nSweeps_ << " damping: " << omega_ <<*/ std::endl;
335 }
336
337 if (verbLevel & Debug) {
338 out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
339 }
340 }
341
342 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
344 // FIXME: This is a placeholder
345 return Teuchos::OrdinalTraits<size_t>::invalid();
346 }
347
348} // namespace MueLu
349
350#endif /* MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_ */
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
#define MUELU_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
BraessSarazin smoother for 2x2 block matrices.
RCP< const ParameterList > GetValidParameterList() const
Input.
std::string description() const
Return a simple one-line description of this object.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
void DeclareInput(Level &currentLevel) const
Input.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void Setup(Level &currentLevel)
Setup routine.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos=0)
Add a factory manager for BraessSarazin internal SchurComplement handling.
virtual std::string description() const
Return a simple one-line description of this object.
Exception indicating invalid cast attempted.
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.
Class that holds all level-specific information.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access)....
virtual const Teuchos::ParameterList & GetParameterList() const =0
An exception safe way to call the method 'Level::SetFactoryManager()'.
bool IsSetup() const
Get the state of a smoother prototype.
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::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
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.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ Parameters0
Print class parameters.