MueLu  Version of the Day
MueLu_BlockedGaussSeidelSmoother_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_BlockedGaussSeidelSmoother_def.hpp
48  *
49  * Created on: 30.01.2012
50  * Author: tobias
51  */
52 
53 #ifndef MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_
54 #define MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_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 
66 #include "MueLu_Level.hpp"
67 #include "MueLu_Utilities.hpp"
68 #include "MueLu_Monitor.hpp"
69 #include "MueLu_HierarchyUtils.hpp"
70 #include "MueLu_SmootherBase.hpp"
71 
72 namespace MueLu {
73 
74  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
76  : type_("blocked GaussSeidel"), A_(Teuchos::null)
77  {
78  FactManager_.reserve(10); // TODO fix me!
79  }
80 
81  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
83 
84  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
86  RCP<ParameterList> validParamList = rcp(new ParameterList());
87 
88  validParamList->set< RCP<const FactoryBase> >("A", Teuchos::null, "Generating factory of the matrix A");
89  validParamList->set< Scalar > ("Damping factor", 1.0, "Damping/Scaling factor in BGS");
90  validParamList->set< LocalOrdinal > ("Sweeps", 1, "Number of BGS sweeps (default = 1)");
91 
92  return validParamList;
93  }
94 
95  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
97  TEUCHOS_TEST_FOR_EXCEPTION(pos < 0, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::AddFactoryManager: parameter \'pos\' must not be negative! error.");
98 
99  size_t myPos = Teuchos::as<size_t>(pos);
100 
101  if (myPos < FactManager_.size()) {
102  // replace existing entris in FactManager_ vector
103  FactManager_.at(myPos) = FactManager;
104  } else if( myPos == FactManager_.size()) {
105  // add new Factory manager in the end of the vector
106  FactManager_.push_back(FactManager);
107  } else { // if(myPos > FactManager_.size())
108  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
109  *out << "Warning: cannot add new FactoryManager at proper position " << pos << ". The FactoryManager is just appended to the end. Check this!" << std::endl;
110 
111  // add new Factory manager in the end of the vector
112  FactManager_.push_back(FactManager);
113  }
114  }
115 
116  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
118  //this->Input(currentLevel, "A");
119  // TODO: check me: why is this->Input not freeing properly A in release mode?
120  currentLevel.DeclareInput("A",this->GetFactory("A").get());
121 
122  // loop over all factory managers for the subblocks of blocked operator A
123  std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
124  for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
125  SetFactoryManager currentSFM (rcpFromRef(currentLevel), *it);
126 
127  // request "Smoother" for current subblock row.
128  currentLevel.DeclareInput("PreSmoother",(*it)->GetFactory("Smoother").get());
129 
130  // request "A" for current subblock row (only needed for Thyra mode)
131  currentLevel.DeclareInput("A",(*it)->GetFactory("A").get());
132  }
133 
134  //RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
135  }
136 
137  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
139  //typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> BlockedCrsOMatrix;
140 
141  RCP<Teuchos::FancyOStream> out = Teuchos::fancyOStream(Teuchos::rcpFromRef(std::cout));
142 
143  FactoryMonitor m(*this, "Setup blocked Gauss-Seidel Smoother", currentLevel);
144  if (SmootherPrototype::IsSetup() == true) this->GetOStream(Warnings0) << "MueLu::BlockedGaussSeidelSmoother::Setup(): Setup() has already been called";
145 
146  // extract blocked operator A from current level
147  A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A"); // A needed for extracting map extractors
148  RCP<BlockedCrsMatrix> bA = Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(A_);
149  TEUCHOS_TEST_FOR_EXCEPTION(bA==Teuchos::null, Exceptions::BadCast, "MueLu::BlockedPFactory::Build: input matrix A is not of type BlockedCrsMatrix! error.");
150 
151  // plausibility check
152  TEUCHOS_TEST_FOR_EXCEPTION(bA->Rows() != FactManager_.size(), Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Setup: number of block rows of A is " << bA->Rows() << " and does not match number of SubFactoryManagers " << FactManager_.size() << ". error.");
153  TEUCHOS_TEST_FOR_EXCEPTION(bA->Cols() != FactManager_.size(), Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Setup: number of block cols of A is " << bA->Cols() << " and does not match number of SubFactoryManagers " << FactManager_.size() << ". error.");
154 
155  // store map extractors
156  rangeMapExtractor_ = bA->getRangeMapExtractor();
157  domainMapExtractor_ = bA->getDomainMapExtractor();
158 
159  Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
160 
161  // loop over all factory managers for the subblocks of blocked operator A
162  std::vector<Teuchos::RCP<const FactoryManagerBase> >::const_iterator it;
163  for(it = FactManager_.begin(); it!=FactManager_.end(); ++it) {
164  SetFactoryManager currentSFM (rcpFromRef(currentLevel), *it);
165 
166  // extract Smoother for current block row (BGS ordering)
167  RCP<const SmootherBase> Smoo = currentLevel.Get< RCP<SmootherBase> >("PreSmoother",(*it)->GetFactory("Smoother").get());
168  Inverse_.push_back(Smoo);
169 
170  // store whether subblock matrix is blocked or not!
171  RCP<Matrix> Aii = currentLevel.Get< RCP<Matrix> >("A",(*it)->GetFactory("A").get());
172  bIsBlockedOperator_.push_back(Teuchos::rcp_dynamic_cast<BlockedCrsMatrix>(Aii)!=Teuchos::null);
173  }
174 
176  }
177 
178  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
179  void BlockedGaussSeidelSmoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector &X, const MultiVector& B, bool InitialGuessIsZero) const
180  {
181  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::BlockedGaussSeidelSmoother::Apply(): Setup() has not been called");
182 
183 #ifdef HAVE_MUELU_DEBUG
184  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.");
185  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.");
186 #endif
187 
188 #if 0
189  RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > residual = MultiVectorFactory::Build(B.getMap(), B.getNumVectors());
190  RCP<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > tempres = MultiVectorFactory::Build(B.getMap(), B.getNumVectors());
191  RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
192 
193  //Teuchos::RCP<Teuchos::FancyOStream> fos = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cout));
194 
195  // extract parameters from internal parameter list
197  LocalOrdinal nSweeps = pL.get<LocalOrdinal>("Sweeps");
198  Scalar omega = pL.get<Scalar>("Damping factor");
199 
200  // outer Richardson loop
201  for (LocalOrdinal run = 0; run < nSweeps; ++run) {
202  // one BGS sweep
203  // loop over all block rows
204  for(size_t i = 0; i<Inverse_.size(); i++) {
205 
206  // calculate block residual r = B-A*X
207  // note: A_ is the full blocked operator
208  residual->update(1.0,B,0.0); // r = B
209 
210  A_->apply(*rcpX, *residual, Teuchos::NO_TRANS, -1.0, 1.0);
211 
212  // extract corresponding subvectors from X and residual
213  bool bRangeThyraMode = rangeMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] == false);
214  bool bDomainThyraMode = domainMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] == false);
215  Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(rcpX, i, bDomainThyraMode);
216  Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
217  Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode);
218 
219  // apply solver/smoother
220  Inverse_.at(i)->Apply(*tXi, *ri, false);
221 
222  // update vector
223  Xi->update(omega,*tXi,1.0); // X_{i+1} = X_i + omega \Delta X_i
224 
225  // update corresponding part of rhs and lhs
226  domainMapExtractor_->InsertVector(Xi, i, rcpX, bDomainThyraMode); // TODO wrong! fix me
227  }
228  }
229 #else
230  // new implementation uses BlockedMultiVectors
231 
232  // create a new vector for storing the current residual in a blocked multi vector
233  RCP<MultiVector> res = MultiVectorFactory::Build(B.getMap(), B.getNumVectors(), true);
234  RCP<BlockedMultiVector> residual = Teuchos::rcp(new BlockedMultiVector(rangeMapExtractor_,res));
235 
236  // create a new solution vector as a blocked multi vector
237  RCP<MultiVector> rcpX = Teuchos::rcpFromRef(X);
238  RCP<BlockedMultiVector> bX = Teuchos::rcp(new BlockedMultiVector(domainMapExtractor_,rcpX));
239 
240  // create a blocked rhs vector
241  RCP<const MultiVector> rcpB = Teuchos::rcpFromRef(B);
242  RCP<const BlockedMultiVector> bB = Teuchos::rcp(new const BlockedMultiVector(rangeMapExtractor_,rcpB));
243 
244  // extract parameters from internal parameter list
246  LocalOrdinal nSweeps = pL.get<LocalOrdinal>("Sweeps");
247  Scalar omega = pL.get<Scalar>("Damping factor");
248 
249  // outer Richardson loop
250  for (LocalOrdinal run = 0; run < nSweeps; ++run) {
251  // one BGS sweep
252  // loop over all block rows
253  for(size_t i = 0; i<Inverse_.size(); i++) {
254 
255  // calculate block residual r = B-A*X
256  // note: A_ is the full blocked operator
257  residual->update(1.0,*bB,0.0); // r = B
258 
259  A_->apply(*bX, *residual, Teuchos::NO_TRANS, -1.0, 1.0);
260 
261  // extract corresponding subvectors from X and residual
262  bool bRangeThyraMode = rangeMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] == false);
263  bool bDomainThyraMode = domainMapExtractor_->getThyraMode() && (bIsBlockedOperator_[i] == false);
264 
265  Teuchos::RCP<MultiVector> Xi = domainMapExtractor_->ExtractVector(bX, i, bDomainThyraMode);
266  Teuchos::RCP<MultiVector> ri = rangeMapExtractor_->ExtractVector(residual, i, bRangeThyraMode);
267  Teuchos::RCP<MultiVector> tXi = domainMapExtractor_->getVector(i, X.getNumVectors(), bDomainThyraMode, true);
268 
269  // apply solver/smoother
270  Inverse_.at(i)->Apply(*tXi, *ri, false);
271 
272  // update vector
273  Xi->update(omega,*tXi,1.0); // X_{i+1} = X_i + omega \Delta X_i
274 
275  // update corresponding part of rhs and lhs
276 
277  domainMapExtractor_->InsertVector(Xi, i, bX, bDomainThyraMode);
278  }
279  }
280 
281  // write back solution
282  for(size_t r = 0; r < domainMapExtractor_->NumMaps(); ++r) {
283  bool bDomainThyraMode = domainMapExtractor_->getThyraMode() && (bIsBlockedOperator_[r] == false);
284  domainMapExtractor_->InsertVector(bX->getMultiVector(r,bDomainThyraMode),r,rcpX,bDomainThyraMode);
285  }
286 #endif
287  }
288 
289  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
291  return rcp( new BlockedGaussSeidelSmoother(*this) );
292  }
293 
294  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
296  std::ostringstream out;
298  out << "{type = " << type_ << "}";
299  return out.str();
300  }
301 
302  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
305 
306  // extract parameters from internal parameter list
308  LocalOrdinal nSweeps = pL.get<LocalOrdinal>("Sweeps");
309  Scalar omega = pL.get<Scalar>("Damping factor");
310 
311  if (verbLevel & Parameters0) {
312  out0 << "Prec. type: " << type_ << " Sweeps: " << nSweeps << " damping: " << omega << std::endl;
313  }
314 
315  if (verbLevel & Debug) {
316  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
317  }
318  }
319 
320 } // namespace MueLu
321 
322 
323 
324 #endif /* MUELU_BLOCKEDGAUSSSEIDELSMOOTHER_DEF_HPP_ */
Important warning messages (one line)
Exception indicating invalid cast attempted.
bool IsSetup() const
Get the state of a smoother prototype.
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)
block Gauss-Seidel method for blocked matrices
Print additional debugging information.
Namespace for MueLu classes and methods.
std::string description() const
Return a simple one-line description of this object.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
void DeclareInput(Level &currentLevel) const
Input.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
void Setup(Level &currentLevel)
Setup routine In the Setup method the Inverse_ vector is filled with the corresponding SmootherBase o...
RCP< const ParameterList > GetValidParameterList() const
Input.
std::vector< Teuchos::RCP< const FactoryManagerBase > > FactManager_
vector of factory managers
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
virtual const Teuchos::ParameterList & GetParameterList() const
Print class parameters.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the direct solver. Solves the linear system AX=B using the constructed solver.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
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 AddFactoryManager(RCP< const FactoryManagerBase > FactManager, int pos)
Add a factory manager.
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.
std::string toString(const T &t)