MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_DirectSolver_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_DIRECTSOLVER_DEF_HPP
47 #define MUELU_DIRECTSOLVER_DEF_HPP
48 
49 #include <Xpetra_Utils.hpp>
50 #include <Xpetra_Matrix.hpp>
51 
53 
54 #include "MueLu_Exceptions.hpp"
55 #include "MueLu_Level.hpp"
56 
57 #include "MueLu_Amesos2Smoother.hpp"
58 #include "MueLu_AmesosSmoother.hpp"
59 
60 namespace MueLu {
61 
62  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
64  : type_(type) {
65  // The original idea behind all smoothers was to use prototype pattern. However, it does not fully work of the dependencies
66  // calculation. Particularly, we need to propagate DeclareInput to proper prototypes. Therefore, both TrilinosSmoother and
67  // DirectSolver do not follow the pattern exactly.
68  // The difference is that in order to propagate the calculation of dependencies, we need to call a DeclareInput on a
69  // constructed object (or objects, as we have two different code branches for Epetra and Tpetra). The only place where we
70  // could construct these objects is the constructor. Thus, we need to store RCPs, and both TrilinosSmoother and DirectSolver
71  // obtain a state: they contain RCP to smoother prototypes.
72  sEpetra_ = Teuchos::null;
73  sTpetra_ = Teuchos::null;
74 
75  ParameterList paramList = paramListIn;
76 
77  // We want DirectSolver to be able to work with both Epetra and Tpetra objects, therefore we try to construct both
78  // Amesos and Amesos2 solver prototypes. The construction really depends on configuration options.
79  triedEpetra_ = triedTpetra_ = false;
80 #if defined(HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_AMESOS2)
81  try {
82  sTpetra_ = rcp(new Amesos2Smoother(type_, paramList));
83  TEUCHOS_TEST_FOR_EXCEPTION(sTpetra_.is_null(), Exceptions::RuntimeError, "Unable to construct Amesos2 direct solver");
84  } catch (Exceptions::RuntimeError& e) {
85  errorTpetra_ = e.what();
86  } catch (Exceptions::BadCast& e) {
87  errorTpetra_ = e.what();
88  }
89  triedTpetra_ = true;
90 #endif
91 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_AMESOS)
92  try {
93  // GetAmesosSmoother masks the template argument matching, and simply throws if template arguments are incompatible with Epetra
94  sEpetra_ = GetAmesosSmoother<SC,LO,GO,NO>(type_, paramList);
95  TEUCHOS_TEST_FOR_EXCEPTION(sEpetra_.is_null(), Exceptions::RuntimeError, "Unable to construct Amesos direct solver");
96  } catch (Exceptions::RuntimeError& e) {
97  // AmesosSmoother throws if Scalar != double, LocalOrdinal != int, GlobalOrdinal != int
98  errorEpetra_ = e.what();
99  }
100  triedEpetra_ = true;
101 #endif
102 
103  // Check if we were able to construct at least one solver. In many cases that's all we need, for instance if a user
104  // simply wants to use Tpetra only stack, never enables Amesos, and always runs Tpetra objects.
105  TEUCHOS_TEST_FOR_EXCEPTION(!triedEpetra_ && !triedTpetra_, Exceptions::RuntimeError, "Unable to construct any direct solver."
106  "Plase enable (TPETRA and AMESOS2) or (EPETRA and AMESOS)");
107 
109  "Could not enable any direct solver:\n"
110  << (triedEpetra_ ? "Epetra mode was disabled due to an error:\n" : "")
111  << (triedEpetra_ ? errorEpetra_ : "")
112  << (triedTpetra_ ? "Tpetra mode was disabled due to an error:\n" : "")
113  << (triedTpetra_ ? errorTpetra_ : ""));
114 
115  this->SetParameterList(paramList);
116  }
117 
118  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
120  // We need to propagate SetFactory to proper place
121  if (!sEpetra_.is_null()) sEpetra_->SetFactory(varName, factory);
122  if (!sTpetra_.is_null()) sTpetra_->SetFactory(varName, factory);
123  }
124 
125  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
127  // Decide whether we are running in Epetra or Tpetra mode
128  //
129  // Theoretically, we could make this decision in the constructor, and create only
130  // one of the smoothers. But we want to be able to reuse, so one can imagine a scenario
131  // where one first runs hierarchy with tpetra matrix, and then with epetra.
132  bool useTpetra = (currentLevel.lib() == Xpetra::UseTpetra);
133  s_ = (useTpetra ? sTpetra_ : sEpetra_);
134  if (s_.is_null()) {
135  if (useTpetra) {
136 #if not defined(HAVE_MUELU_AMESOS2)
138  "Error: running in Tpetra mode, but MueLu with Amesos2 was disabled during the configure stage.\n"
139  "Please make sure that:\n"
140  " - Amesos2 is enabled (Trilinos_ENABLE_Amesos2=ON),\n"
141  " - Amesos2 is available for MueLu to use (MueLu_ENABLE_Amesos2=ON)\n");
142 #else
143  if (triedTpetra_)
144  this->GetOStream(Errors) << "Tpetra mode was disabled due to an error:\n" << errorTpetra_ << std::endl;
145 #endif
146  }
147  if (!useTpetra) {
148 #if not defined(HAVE_MUELU_AMESOS)
150  "Error: running in Epetra mode, but MueLu with Amesos was disabled during the configure stage.\n"
151  "Please make sure that:\n"
152  " - Amesos is enabled (you can do that with Trilinos_ENABLE_Amesos=ON),\n"
153  " - Amesos is available for MueLu to use (MueLu_ENABLE_Amesos=ON)\n");
154 #else
155  if (triedEpetra_)
156  this->GetOStream(Errors) << "Epetra mode was disabled due to an error:\n" << errorEpetra_ << std::endl;
157 #endif
158  }
160  "Direct solver for " << (useTpetra ? "Tpetra" : "Epetra") << " was not constructed");
161  }
162 
163  s_->DeclareInput(currentLevel);
164  }
165 
166  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
168  if (SmootherPrototype::IsSetup() == true)
169  this->GetOStream(Warnings0) << "MueLu::DirectSolver::Setup(): Setup() has already been called" << std::endl;
170 
171  int oldRank = s_->SetProcRankVerbose(this->GetProcRankVerbose());
172 
173  s_->Setup(currentLevel);
174 
175  s_->SetProcRankVerbose(oldRank);
176 
178 
179  this->SetParameterList(s_->GetParameterList());
180  }
181 
182  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
183  void DirectSolver<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
184  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Apply(): Setup() has not been called");
185 
186  s_->Apply(X, B, InitialGuessIsZero);
187  }
188 
189  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
191  RCP<DirectSolver> newSmoo = rcp(new DirectSolver(*this));
192 
193  // We need to be quite careful with Copy
194  // We still want DirectSolver to follow Prototype Pattern, so we need to hide the fact that we do have some state
195  if (!sEpetra_.is_null())
196  newSmoo->sEpetra_ = sEpetra_->Copy();
197  if (!sTpetra_.is_null())
198  newSmoo->sTpetra_ = sTpetra_->Copy();
199 
200  // Copy the default mode
201  newSmoo->s_ = (s_.get() == sTpetra_.get() ? newSmoo->sTpetra_ : newSmoo->sEpetra_);
202  newSmoo->SetParameterList(this->GetParameterList());
203 
204  return newSmoo;
205  }
206 
207  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
209  std::ostringstream out;
210  if (s_ != Teuchos::null) {
211  out << s_->description();
212  } else {
214  out << "{type = " << type_ << "}";
215  }
216  return out.str();
217  }
218 
219  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
222 
223  if (verbLevel & Parameters0)
224  out0 << "Prec. type: " << type_ << std::endl;
225 
226  if (verbLevel & Parameters1) {
227  out0 << "Parameter list: " << std::endl;
228  Teuchos::OSTab tab3(out);
229  out << this->GetParameterList();
230  }
231 
232  if (verbLevel & Debug)
233  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl;
234  }
235 
236 } // namespace MueLu
237 
238 #endif // MUELU_DIRECTSOLVER_DEF_HPP
Important warning messages (one line)
Exception indicating invalid cast attempted.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
DirectSolver cannot be applied. Apply() always returns a RuntimeError exception.
void Setup(Level &currentLevel)
DirectSolver cannot be turned into a smoother using Setup(). Setup() always returns a RuntimeError ex...
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void DeclareInput(Level &currentLevel) const
Input.
Print additional debugging information.
std::string type_
amesos1/2-specific key phrase that denote smoother type
DirectSolver(const std::string &type="", const Teuchos::ParameterList &paramList=Teuchos::ParameterList())
Constructor Note: only parameters shared by Amesos and Amesos2 should be used for type and paramList ...
Class that encapsulates direct solvers. Autoselection of AmesosSmoother or Amesos2Smoother according ...
RCP< SmootherPrototype > Copy() const
When this prototype is cloned using Copy(), the clone is an Amesos or an Amesos2 smoother.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
std::string description() const
Return a simple one-line description of this object.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Class that encapsulates Amesos2 direct solvers.
bool IsSetup() const
Get the state of a smoother prototype.
Xpetra::UnderlyingLib lib()
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
RCP< SmootherPrototype > sTpetra_
Print class parameters.
void SetFactory(const std::string &varName, const RCP< const FactoryBase > &factory)
Custom SetFactory.
RCP< SmootherPrototype > sEpetra_
Smoother.
RCP< SmootherPrototype > s_
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
virtual std::string description() const
Return a simple one-line description of this object.
std::string toString(const T &t)