MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 
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  template <class Node>
212  // FIXME: This is a placeholder
214  }
215 
216 
217 
218 } // namespace MueLu
219 
220 // The AmesosSmoother is only templated on the Node, since it is an Epetra only object
221 // Therefore we do not need the full ETI instantiations as we do for the other MueLu
222 // objects which are instantiated on all template parameters.
223 #if defined(HAVE_MUELU_EPETRA)
225 #endif
226 
227 
228 #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=false) const
Apply the direct solver.
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Print external lib objects.
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)
std::string description() const
Return a simple one-line description of this object.
Print additional debugging information.
std::string tolower(const std::string &str)
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
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
RCP< SmootherPrototype > Copy() const
bool IsSetup() const
Get the state of a smoother prototype.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Print class parameters.
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Print class parameters (more parameters, more verbose)
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Exception throws to report errors in the internal logical of the program.
Description of what is happening (more verbose)
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
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)
void DeclareInput(Level &currentLevel) const
Input.
Class that encapsulates Amesos direct solvers.