MueLu  Version of the Day
MueLu_AmesosSmoother.cpp
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 #include <algorithm>
47 
48 #include "MueLu_ConfigDefs.hpp"
49 
50 #if defined(HAVE_MUELU_EPETRA) && defined(HAVE_MUELU_AMESOS)
51 
52 #include <Epetra_LinearProblem.h>
53 
54 #include <Amesos_config.h>
55 #include <Amesos.h>
56 #include <Amesos_BaseSolver.h>
57 
58 #include "MueLu_AmesosSmoother.hpp"
59 
60 #include "MueLu_Level.hpp"
61 #include "MueLu_Utilities.hpp"
62 #include "MueLu_Monitor.hpp"
63 
64 namespace MueLu {
65 
66  template <class Node>
67  AmesosSmoother<Node>::AmesosSmoother(const std::string& type, const Teuchos::ParameterList& paramList)
68  : type_(type) {
69  this->SetParameterList(paramList);
70 
71  if (!type_.empty()) {
72  // Transform string to "Abcde" notation
73  std::transform(type_.begin(), type_.end(), type_.begin(), ::tolower);
74  std::transform(type_.begin(), ++type_.begin(), type_.begin(), ::toupper);
75  }
76  if (type_ == "Amesos_klu") type_ = "Klu";
77  if (type_ == "Klu2") type_ = "Klu";
78  if (type_ == "Amesos_umfpack") type_ = "Umfpack";
79  if (type_ == "Superlu_dist") type_ = "Superludist";
80  if (type_ == "Amesos_mumps") type_ = "Mumps";
81 
82  // Try to come up with something availble
83  // Order corresponds to our preference
84  // TODO: It would be great is Amesos provides directly this kind of logic for us
85  std::string oldtype = type_;
86  if (type_ == "" || Amesos().Query(type_) == false) {
87 #if defined(HAVE_AMESOS_SUPERLU)
88  type_ = "Superlu";
89 #elif defined(HAVE_AMESOS_KLU)
90  type_ = "Klu";
91 #elif defined(HAVE_AMESOS_SUPERLUDIST)
92  type_ = "Superludist";
93 #elif defined(HAVE_AMESOS_UMFPACK)
94  type_ = "Umfpack";
95 #else
96  throw Exceptions::RuntimeError("Amesos has been compiled without SuperLU_DIST, SuperLU, Umfpack or Klu. By default, MueLu tries"
97  "to use one of these libraries. Amesos must be compiled with one of these solvers, "
98  "or a valid Amesos solver has to be specified explicitly.");
99 #endif
100  if (oldtype != "")
101  this->GetOStream(Warnings0) << "MueLu::AmesosSmoother: \"" << oldtype << "\" is not available. Using \"" << type_ << "\" instead" << std::endl;
102  else
103  this->GetOStream(Runtime1) << "MueLu::AmesosSmoother: using \"" << type_ << "\"" << std::endl;
104  }
105  }
106 
107  template <class Node>
108  void AmesosSmoother<Node>::DeclareInput(Level &currentLevel) const {
109  this->Input(currentLevel, "A");
110  }
111 
112  template <class Node>
113  void AmesosSmoother<Node>::Setup(Level& currentLevel) {
114  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
115 
116  if (SmootherPrototype::IsSetup() == true)
117  this->GetOStream(Warnings0) << "MueLu::AmesosSmoother::Setup(): Setup() has already been called" << std::endl;
118 
119  A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");
120 
122  linearProblem_ = rcp( new Epetra_LinearProblem() );
123  linearProblem_->SetOperator(epA.get());
124 
125  Amesos factory;
126  prec_ = rcp(factory.Create(type_, *linearProblem_));
127  TEUCHOS_TEST_FOR_EXCEPTION(prec_ == Teuchos::null, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Setup(): Solver '" + type_ + "' not supported by Amesos");
128 
129  // set Reindex flag, if A is distributed with non-contiguous maps
130  // unfortunately there is no reindex for Amesos2, yet. So, this only works for Epetra based problems
131  if (A_->getRowMap()->isDistributed() == true && A_->getRowMap()->isContiguous() == false)
132  const_cast<ParameterList&>(this->GetParameterList()).set("Reindex", true);
133 
134  const ParameterList& paramList = this->GetParameterList();
135  RCP<ParameterList> precList = this->RemoveFactoriesFromList(paramList);
136 
137  prec_->SetParameters(*precList);
138 
139  const_cast<ParameterList&>(paramList).setParameters(*precList);
140 
141  int r = prec_->NumericFactorization();
142  TEUCHOS_TEST_FOR_EXCEPTION(r != 0, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Setup(): Amesos solver returns value of " +
143  Teuchos::toString(r) + " during NumericFactorization()");
144 
146  }
147 
148  template <class Node>
149  void AmesosSmoother<Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
150  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Apply(): Setup() has not been called");
151 
152  Epetra_MultiVector &epX = Utilities::MV2NonConstEpetraMV(X);
153  Epetra_MultiVector const &epB = Utilities::MV2EpetraMV(B);
154  //Epetra_LinearProblem takes the right-hand side as a non-const pointer.
155  //I think this const_cast is safe because Amesos won't modify the rhs.
156  Epetra_MultiVector &nonconstB = const_cast<Epetra_MultiVector&>(epB);
157 
158  linearProblem_->SetLHS(&epX);
159  linearProblem_->SetRHS(&nonconstB);
160 
161  prec_->Solve();
162 
163  // Don't keep pointers to our vectors in the Epetra_LinearProblem.
164  linearProblem_->SetLHS(0);
165  linearProblem_->SetRHS(0);
166  }
167 
168  template <class Node>
170  return rcp( new AmesosSmoother<Node>(*this) );
171  }
172 
173  template <class Node>
174  std::string AmesosSmoother<Node>::description() const {
175  std::ostringstream out;
177  out << "{type = " << type_ << "}";
178  return out.str();
179  }
180 
181  //using MueLu::Describable::describe; // overloading, not hiding
182  template <class Node>
183  void AmesosSmoother<Node>::print(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
185 
186  if (verbLevel & Parameters0)
187  out0 << "Prec. type: " << type_ << std::endl;
188 
189  if (verbLevel & Parameters1) {
190  out0 << "Parameter list: " << std::endl;
191  Teuchos::OSTab tab2(out);
192  out << this->GetParameterList();
193  }
194 
195  if (verbLevel & External)
196  if (prec_ != Teuchos::null) {
197  prec_->PrintStatus();
198  prec_->PrintTiming();
199  }
200 
201  if (verbLevel & Debug) {
202  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
203  << "-" << std::endl
204  << "RCP<A_>: " << A_ << std::endl
205  << "RCP<linearProblem__>: " << linearProblem_ << std::endl
206  << "RCP<prec_>: " << prec_ << std::endl;
207  }
208  }
209 
210 } // namespace MueLu
211 
212 // The AmesosSmoother is only templated on the Node, since it is an Epetra only object
213 // Therefore we do not need the full ETI instantiations as we do for the other MueLu
214 // objects which are instantiated on all template parameters.
215 #if defined(HAVE_MUELU_EPETRA)
217 #endif
218 
219 
220 #endif // HAVE_MUELU_EPETRA && HAVE_MUELU_AMESOS
Important warning messages (one line)
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the direct solver.
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Print external lib objects.
bool IsSetup() const
Get the state of a smoother prototype.
static RCP< const Epetra_MultiVector > MV2EpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Epetra objects from an Xpetra object.
Timer to be used in factories. Similar to Monitor but with additional timers.
AmesosSmoother(const std::string &type="", const Teuchos::ParameterList &paramList=Teuchos::ParameterList())
Constructor.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Print additional debugging information.
Namespace for MueLu classes and methods.
void Setup(Level &currentLevel)
Set up the direct solver. This creates the underlying Amesos solver object according to the parameter...
std::string type_
amesos-specific key phrase that denote smoother type
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
virtual void SetParameterList(const ParameterList &paramList)
Set parameters from a parameter list and return with default values.
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
RCP< SmootherPrototype > Copy() const
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
void DeclareInput(Level &currentLevel) const
Input.
Print class parameters.
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
std::string description() const
Return a simple one-line description of this object.
Print class parameters (more parameters, more verbose)
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
virtual std::string description() const
Return a simple one-line description of this object.
std::string toString(const T &t)
bool Query(const char *ClassType)
Class that encapsulates Amesos direct solvers.