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  this->declareConstructionOutcome(true, "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  return;
100 #endif
101  if (oldtype != "")
102  this->GetOStream(Warnings0) << "MueLu::AmesosSmoother: \"" << oldtype << "\" is not available. Using \"" << type_ << "\" instead" << std::endl;
103  else
104  this->GetOStream(Runtime1) << "MueLu::AmesosSmoother: using \"" << type_ << "\"" << std::endl;
105  }
106  this->declareConstructionOutcome(false, "");
107  }
108 
109  template <class Node>
110  void AmesosSmoother<Node>::DeclareInput(Level &currentLevel) const {
111  this->Input(currentLevel, "A");
112  }
113 
114  template <class Node>
115  void AmesosSmoother<Node>::Setup(Level& currentLevel) {
116  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
117 
118  if (SmootherPrototype::IsSetup() == true)
119  this->GetOStream(Warnings0) << "MueLu::AmesosSmoother::Setup(): Setup() has already been called" << std::endl;
120 
121  A_ = Factory::Get< RCP<Matrix> >(currentLevel, "A");
122 
123  RCP<Epetra_CrsMatrix> epA = Utilities::Op2NonConstEpetraCrs(A_);
124  linearProblem_ = rcp( new Epetra_LinearProblem() );
125  linearProblem_->SetOperator(epA.get());
126 
127  Amesos factory;
128  prec_ = rcp(factory.Create(type_, *linearProblem_));
129  TEUCHOS_TEST_FOR_EXCEPTION(prec_ == Teuchos::null, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Setup(): Solver '" + type_ + "' not supported by Amesos");
130 
131  // set Reindex flag, if A is distributed with non-contiguous maps
132  // unfortunately there is no reindex for Amesos2, yet. So, this only works for Epetra based problems
133  if (A_->getRowMap()->isDistributed() == true && A_->getRowMap()->isContiguous() == false)
134  const_cast<ParameterList&>(this->GetParameterList()).set("Reindex", true);
135 
136  const ParameterList& paramList = this->GetParameterList();
137  RCP<ParameterList> precList = this->RemoveFactoriesFromList(paramList);
138 
139  prec_->SetParameters(*precList);
140 
141  const_cast<ParameterList&>(paramList).setParameters(*precList);
142 
143  int r = prec_->NumericFactorization();
144  TEUCHOS_TEST_FOR_EXCEPTION(r != 0, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Setup(): Amesos solver returns value of " +
145  Teuchos::toString(r) + " during NumericFactorization()");
146 
148  }
149 
150  template <class Node>
151  void AmesosSmoother<Node>::Apply(MultiVector& X, const MultiVector& B, bool InitialGuessIsZero) const {
152  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::AmesosSmoother::Apply(): Setup() has not been called");
153 
156  //Epetra_LinearProblem takes the right-hand side as a non-const pointer.
157  //I think this const_cast is safe because Amesos won't modify the rhs.
158  Epetra_MultiVector &nonconstB = const_cast<Epetra_MultiVector&>(epB);
159 
160  linearProblem_->SetLHS(&epX);
161  linearProblem_->SetRHS(&nonconstB);
162 
163  prec_->Solve();
164 
165  // Don't keep pointers to our vectors in the Epetra_LinearProblem.
166  linearProblem_->SetLHS(0);
167  linearProblem_->SetRHS(0);
168  }
169 
170  template <class Node>
171  RCP<MueLu::SmootherPrototype<double,int,int,Node> > AmesosSmoother<Node>::Copy() const {
172  return rcp( new AmesosSmoother<Node>(*this) );
173  }
174 
175  template <class Node>
176  std::string AmesosSmoother<Node>::description() const {
177  std::ostringstream out;
179  out << "{type = " << type_ << "}";
180  return out.str();
181  }
182 
183  //using MueLu::Describable::describe; // overloading, not hiding
184  template <class Node>
185  void AmesosSmoother<Node>::print(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
187 
188  if (verbLevel & Parameters0)
189  out0 << "Prec. type: " << type_ << std::endl;
190 
191  if (verbLevel & Parameters1) {
192  out0 << "Parameter list: " << std::endl;
193  Teuchos::OSTab tab2(out);
194  out << this->GetParameterList();
195  }
196 
197  if (verbLevel & External)
198  if (prec_ != Teuchos::null) {
199  prec_->PrintStatus();
200  prec_->PrintTiming();
201  }
202 
203  if (verbLevel & Debug) {
204  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
205  << "-" << std::endl
206  << "RCP<A_>: " << A_ << std::endl
207  << "RCP<linearProblem__>: " << linearProblem_ << std::endl
208  << "RCP<prec_>: " << prec_ << std::endl;
209  }
210  }
211 
212  template <class Node>
214  // FIXME: This is a placeholder
215  return Teuchos::OrdinalTraits<size_t>::invalid();
216  }
217 
218 
219 
220 } // namespace MueLu
221 
222 // The AmesosSmoother is only templated on the Node, since it is an Epetra only object
223 // Therefore we do not need the full ETI instantiations as we do for the other MueLu
224 // objects which are instantiated on all template parameters.
225 #if defined(HAVE_MUELU_EPETRA)
227 #endif
228 
229 
230 #endif // HAVE_MUELU_EPETRA && HAVE_MUELU_AMESOS
Important warning messages (one line)
std::string toString(const T &what)
Little helper function to convert non-string types to strings.
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
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.
Print additional debugging information.
std::string tolower(const std::string &str)
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...
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
std::string type_
amesos-specific key phrase that denote smoother type
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.
void declareConstructionOutcome(bool fail, std::string msg)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
RCP< SmootherPrototype > Copy() const
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
void DeclareInput(Level &currentLevel) const
Input.
void Apply(MultiVector &X, const MultiVector &B, bool=false) const
Apply the direct solver.
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.
bool Query(const char *ClassType)
Class that encapsulates Amesos direct solvers.