MueLu  Version of the Day
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  * MueLu_BraessSarazinSmoother_def.hpp
48  *
49  * Created on: Apr 16, 2012
50  * Author: wiesner
51  */
52 
53 #ifndef MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
54 #define MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_
55 
56 #include "Teuchos_ArrayViewDecl.hpp"
57 #include "Teuchos_ScalarTraits.hpp"
58 
59 #include "MueLu_ConfigDefs.hpp"
60 
61 #include <Xpetra_Matrix.hpp>
62 #include <Xpetra_BlockedCrsMatrix.hpp>
63 #include <Xpetra_MultiVectorFactory.hpp>
64 #include <Xpetra_VectorFactory.hpp>
65 
67 #include "MueLu_Level.hpp"
68 #include "MueLu_Utilities.hpp"
69 #include "MueLu_Monitor.hpp"
70 #include "MueLu_HierarchyUtils.hpp"
71 #include "MueLu_SmootherBase.hpp"
72 
73 // include files for default FactoryManager
74 #include "MueLu_SchurComplementFactory.hpp"
75 #include "MueLu_DirectSolver.hpp"
76 #include "MueLu_SmootherFactory.hpp"
77 #include "MueLu_FactoryManager.hpp"
78 
79 namespace MueLu {
80 
81  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
83  : type_("Braess Sarazin"), A_(null)
84  {
85  //Factory::SetParameter("Sweeps", ParameterEntry(sweeps));
86  //Factory::SetParameter("Damping factor",ParameterEntry(omega));
87 
88 #if 0
89  // when declaring default factories without overwriting them leads to a multipleCallCheck exception
90  // TODO: debug into this
91  // workaround: always define your factory managers outside either using the C++ API or the XML files
92  RCP<SchurComplementFactory> SchurFact = rcp(new SchurComplementFactory());
93  SchurFact->SetParameter("omega",ParameterEntry(omega));
94  SchurFact->SetFactory("A", this->GetFactory("A"));
95 
96  // define smoother/solver for BraessSarazin
97  ParameterList SCparams;
98  std::string SCtype;
99  RCP<SmootherPrototype> smoProtoSC = rcp( new DirectSolver(SCtype,SCparams) );
100  smoProtoSC->SetFactory("A", SchurFact);
101 
102  RCP<SmootherFactory> SmooSCFact = rcp( new SmootherFactory(smoProtoSC) );
103 
104  RCP<FactoryManager> FactManager = rcp(new FactoryManager());
105  FactManager->SetFactory("A", SchurFact);
106  FactManager->SetFactory("Smoother", SmooSCFact);
107  FactManager->SetIgnoreUserData(true);
108 
109  AddFactoryManager(FactManager,0);
110 #endif
111  }
112 
113  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
115 
117  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
118  void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::AddFactoryManager(RCP<const FactoryManagerBase> FactManager, int pos) {
119  TEUCHOS_TEST_FOR_EXCEPTION(pos != 0, Exceptions::RuntimeError, "MueLu::BraessSarazinSmoother::AddFactoryManager: parameter \'pos\' must be zero! error.");
120  FactManager_ = FactManager;
121  }
122 
123  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
125  RCP<ParameterList> validParamList = rcp(new ParameterList());
126 
127  SC one = Teuchos::ScalarTraits<SC>::one();
128 
129  validParamList->set<RCP<const FactoryBase> >("A", null, "Generating factory of the matrix A");
130 
131  validParamList->set<bool> ("lumping", false, "Use lumping to construct diag(A(0,0)), i.e. use row sum of the abs values on the diagonal "
132  "as approximation of A00 (and A00^{-1})");
133  validParamList->set<SC> ("Damping factor", one, "Damping/Scaling factor in BraessSarazin (usually has to be chosen > 1, default = 1.0)");
134  validParamList->set<LO> ("Sweeps", 1, "Number of BraessSarazin sweeps (default = 1)");
135  //validParamList->set<ParameterList> ("block1", ParameterList(), "Sublist for parameters for SchurComplement block. At least has to contain some information about a smoother \"Smoother\" for variable \"A\" which is generated by a SchurComplementFactory.");
136  validParamList->set<bool> ("q2q1 mode", false, "Run in the mode matching Q2Q1 matlab code");
137 
138  return validParamList;
139  }
140 
141  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
143  this->Input(currentLevel, "A");
144 
145  TEUCHOS_TEST_FOR_EXCEPTION(FactManager_.is_null(), Exceptions::RuntimeError,
146  "MueLu::BraessSarazinSmoother::DeclareInput: FactManager_ must not be null! "
147  "Introduce a FactoryManager for the SchurComplement equation.");
148 
149  // carefully call DeclareInput after switching to internal FactoryManager
150  {
151  SetFactoryManager currentSFM(rcpFromRef(currentLevel), FactManager_);
152 
153  // request "Smoother" for current subblock row.
154  currentLevel.DeclareInput("PreSmoother", FactManager_->GetFactory("Smoother").get());
155 
156  // request Schur matrix just in case
157  currentLevel.DeclareInput("A", FactManager_->GetFactory("A").get());
158  }
159  }
160 
161  // Setup routine can be summarized in 4 steps:
162  // - set the map extractors
163  // - set the blocks
164  // - create and set the inverse of the diagonal of F
165  // - set the smoother for the Schur Complement
166  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
168  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
169 
170  if (SmootherPrototype::IsSetup() == true)
171  this->GetOStream(Warnings0) << "MueLu::BreaessSarazinSmoother::Setup(): Setup() has already been called";
172 
173  // Extract blocked operator A from current level
174  A_ = Factory::Get<RCP<Matrix> > (currentLevel, "A");
175  RCP<BlockedCrsMatrix> bA = rcp_dynamic_cast<BlockedCrsMatrix>(A_);
176  TEUCHOS_TEST_FOR_EXCEPTION(bA.is_null(), Exceptions::BadCast,
177  "MueLu::BraessSarazinSmoother::Setup: input matrix A is not of type BlockedCrsMatrix! error.");
178 
179  // Store map extractors
180  rangeMapExtractor_ = bA->getRangeMapExtractor();
181  domainMapExtractor_ = bA->getDomainMapExtractor();
182 
183  // Store the blocks in local member variables
184  A00_ = bA->getMatrix(0,0);
185  A01_ = bA->getMatrix(0,1);
186  A10_ = bA->getMatrix(1,0);
187  A11_ = bA->getMatrix(1,1);
188 
189  const ParameterList& pL = Factory::GetParameterList();
190  SC omega = pL.get<SC>("Damping factor");
191 
192 #if 0 // old code
193  // Create the inverse of the diagonal of F
194  D_ = VectorFactory::Build(A00_->getRowMap());
195 
196  ArrayRCP<SC> diag;
197  if (pL.get<bool>("lumping") == false)
198  diag = Utilities::GetMatrixDiagonal (*A00_);
199  else
201 
202  SC one = Teuchos::ScalarTraits<SC>::one();
203 
204  ArrayRCP<SC> Ddata = D_->getDataNonConst(0);
205  for (GO row = 0; row < Ddata.size(); row++)
206  Ddata[row] = one / (diag[row]*omega);*/
207 #else
208  // Create the inverse of the diagonal of F
209  // TODO add safety check for zeros on diagonal of F!
210  RCP<Vector> diagFVector = VectorFactory::Build(A00_->getRowMap());
211  if (pL.get<bool>("lumping") == false) {
212  A00_->getLocalDiagCopy(*diagFVector); // extract diagonal of F
213  } else {
214  diagFVector = Utilities::GetLumpedMatrixDiagonal(A00_);
215  }
216  diagFVector->scale(omega);
217  D_ = Utilities::GetInverse(diagFVector);
218 #endif
219 
220  // Set the Smoother
221  // carefully switch to the SubFactoryManagers (defined by the users)
222  {
223  SetFactoryManager currentSFM(rcpFromRef(currentLevel), FactManager_);
224  smoo_ = currentLevel.Get<RCP<SmootherBase> >("PreSmoother", FactManager_->GetFactory("Smoother").get());
225  S_ = currentLevel.Get<RCP<Matrix> > ("A", FactManager_->GetFactory("A").get());
226  }
227 
229  }
230 
231  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
232  void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
233  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError,
234  "MueLu::BraessSarazinSmoother::Apply(): Setup() has not been called");
235 
236 #ifdef HAVE_MUELU_DEBUG
237  TEUCHOS_TEST_FOR_EXCEPTION(A_->getRangeMap()->isSameAs(*(B.getMap())) == false, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): The map of RHS vector B is not the same as range map of the blocked operator A. Please check the map of B and A.");
238  TEUCHOS_TEST_FOR_EXCEPTION(A_->getDomainMap()->isSameAs(*(X.getMap())) == false, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): The map of the solution vector X is not the same as domain map of the blocked operator A. Please check the map of X and A.");
239 #endif
240 
241 
242  // The following boolean flags catch the case where we need special transformation
243  // for the GIDs when calling the subsmoothers.
244  RCP<BlockedCrsMatrix> bA11 = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A11_);
245  bool bA11ThyraSpecialTreatment = false;
246  if (bA11 != Teuchos::null) {
247  if(bA11->Rows() == 1 && bA11->Cols() == 1 && rangeMapExtractor_->getThyraMode() == true) bA11ThyraSpecialTreatment = true;
248  }
249 
250  RCP<MultiVector> rcpX = rcpFromRef(X);
251 
252  // use the GIDs of the sub blocks
253  // This is valid as the subblocks actually represent the GIDs (either Thyra, Xpetra or pseudo Xpetra)
254  RCP<MultiVector> deltaX0 = MultiVectorFactory::Build(A00_->getRowMap(), X.getNumVectors());
255  RCP<MultiVector> deltaX1 = MultiVectorFactory::Build(A10_->getRowMap(), X.getNumVectors());
256  RCP<MultiVector> Rtmp = MultiVectorFactory::Build(A10_->getRowMap(), B.getNumVectors());
257 
258  typedef Teuchos::ScalarTraits<SC> STS;
259  SC one = STS::one(), zero = STS::zero();
260 
261  // extract parameters from internal parameter list
262  const ParameterList& pL = Factory::GetParameterList();
263  LO nSweeps = pL.get<LO>("Sweeps");
264 
265  RCP<MultiVector> R;
266  if (InitialGuessIsZero) {
267  R = MultiVectorFactory::Build(B.getMap(), B.getNumVectors());
268  R->update(one, B, zero);
269  } else {
270  R = Utilities::Residual(*A_, X, B);
271  }
272 
273  // extract diagonal of Schur complement operator
274  RCP<Vector> diagSVector = VectorFactory::Build(S_->getRowMap());
275  S_->getLocalDiagCopy(*diagSVector);
276  ArrayRCP<SC> Sdiag = diagSVector->getDataNonConst(0);
277 
278  for (LO run = 0; run < nSweeps; ++run) {
279  // Extract corresponding subvectors from X and R
280  // Automatically detect whether we use Thyra or Xpetra GIDs
281  // The GIDs should be always compatible with the GIDs of A00, A01, etc...
282  RCP<MultiVector> R0 = rangeMapExtractor_ ->ExtractVector(R, 0);
283  RCP<MultiVector> R1 = rangeMapExtractor_ ->ExtractVector(R, 1);
284 
285  RCP<MultiVector> X0 = domainMapExtractor_->ExtractVector(rcpX, 0);
286  RCP<MultiVector> X1 = domainMapExtractor_->ExtractVector(rcpX, 1);
287 
288  // Calculate Rtmp = R1 - D * deltaX0 (equation 8.14)
289  deltaX0->putScalar(zero);
290  deltaX0->elementWiseMultiply(one, *D_, *R0, zero); // deltaX0 = D * R0 (equation 8.13)
291  A10_->apply(*deltaX0, *Rtmp); // Rtmp = A10*D*deltaX0 (intermediate step)
292  Rtmp->update(one, *R1, -one); // Rtmp = R1 - A10*D*deltaX0
293 
294  if (!pL.get<bool>("q2q1 mode")) {
295  deltaX1->putScalar(zero);
296  } else {
297  ArrayRCP<SC> deltaX1data = deltaX1->getDataNonConst(0);
298  ArrayRCP<SC> Rtmpdata = Rtmp->getDataNonConst(0);
299  for (GO row = 0; row < deltaX1data.size(); row++)
300  deltaX1data[row] = 1.1*Rtmpdata[row] / Sdiag[row];
301  }
302 
303  // Special handling if SchurComplement operator was a 1x1 blocked operator in Thyra mode
304  // Then, we have to translate the Xpetra offset GIDs to plain Thyra GIDs and vice versa
305  if(bA11ThyraSpecialTreatment == true) {
306  RCP<MultiVector> deltaX1_thyra = domainMapExtractor_->getVector(1, X.getNumVectors(), true);
307  RCP<MultiVector> Rtmp_thyra = rangeMapExtractor_->getVector(1, B.getNumVectors(), true);
308  // transform vector
309  for(size_t k=0; k < Rtmp->getNumVectors(); k++) {
310  Teuchos::ArrayRCP<const Scalar> xpetraVecData = Rtmp->getData(k);
311  Teuchos::ArrayRCP<Scalar> thyraVecData = Rtmp_thyra->getDataNonConst(k);
312  for(size_t i=0; i < Rtmp->getLocalLength(); i++) {
313  thyraVecData[i] = xpetraVecData[i];
314  }
315  }
316 
317  smoo_->Apply(*deltaX1_thyra,*Rtmp_thyra);
318 
319  for(size_t k=0; k < deltaX1_thyra->getNumVectors(); k++) {
320  Teuchos::ArrayRCP<Scalar> xpetraVecData = deltaX1->getDataNonConst(k);
321  Teuchos::ArrayRCP<const Scalar> thyraVecData = deltaX1_thyra->getData(k);
322  for(size_t i=0; i < deltaX1_thyra->getLocalLength(); i++) {
323  xpetraVecData[i] = thyraVecData[i];
324  }
325  }
326  } else {
327  // Compute deltaX1 (pressure correction)
328  // We use user provided preconditioner
329  smoo_->Apply(*deltaX1,*Rtmp);
330  }
331 
332  // Compute deltaX0
333  deltaX0->putScalar(zero); // just for safety
334  A01_->apply(*deltaX1, *deltaX0); // deltaX0 = A01*deltaX1
335  deltaX0->update(one, *R0, -one); // deltaX0 = R0 - A01*deltaX1
336  R0.swap(deltaX0);
337  deltaX0->elementWiseMultiply(one, *D_, *R0, zero); // deltaX0 = D*(R0 - A01*deltaX1)
338 
339  // Update solution
340  X0->update(one, *deltaX0, one);
341  X1->update(one, *deltaX1, one);
342 
343  domainMapExtractor_->InsertVector(X0, 0, rcpX);
344  domainMapExtractor_->InsertVector(X1, 1, rcpX);
345 
346  if (run < nSweeps-1)
347  R = Utilities::Residual(*A_, X, B);
348  }
349  }
350 
351  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
352  RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
354  return rcp (new BraessSarazinSmoother (*this));
355  }
356 
357  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
359  description () const {
360  std::ostringstream out;
362  out << "{type = " << type_ << "}";
363  return out.str();
364  }
365 
366  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
367  void BraessSarazinSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream &out, const VerbLevel verbLevel) const {
369 
370  if (verbLevel & Parameters0) {
371  out0 << "Prec. type: " << type_ << /*" Sweeps: " << nSweeps_ << " damping: " << omega_ <<*/ std::endl;
372  }
373 
374  if (verbLevel & Debug) {
375  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
376  }
377  }
378 
379 } // namespace MueLu
380 
381 #endif /* MUELU_BRAESSSARAZINSMOOTHER_DEF_HPP_ */
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Important warning messages (one line)
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
Exception indicating invalid cast attempted.
This class specifies the default factory that should generate some data on a Level if the data does n...
RCP< SmootherPrototype > Copy() const
BraessSarazin smoother for 2x2 block matrices.
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
bool IsSetup() const
Get the state of a smoother prototype.
Timer to be used in factories. Similar to Monitor but with additional timers.
Print additional debugging information.
static Teuchos::ArrayRCP< Scalar > GetLumpedMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
Namespace for MueLu classes and methods.
std::string description() const
Return a simple one-line description of this object.
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
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)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the Braess Sarazin smoother.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
virtual const Teuchos::ParameterList & GetParameterList() const
Print class parameters.
Factory for building the Schur Complement for a 2x2 block matrix.
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())
void AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos=0)
Add a factory manager for BraessSarazin internal SchurComplement handling.
void DeclareInput(Level &currentLevel) const
Input.
void Setup(Level &currentLevel)
Setup routine.
const RCP< const FactoryBase > GetFactory(const std::string &varName) const
Default implementation of FactoryAcceptor::GetFactory()
Exception throws to report errors in the internal logical of the program.
An exception safe way to call the method &#39;Level::SetFactoryManager()&#39;.
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
virtual std::string description() const
Return a simple one-line description of this object.
RCP< const ParameterList > GetValidParameterList() const
Input.