MueLu Version of the Day
Loading...
Searching...
No Matches
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
59
60#include "MueLu_Level.hpp"
61#include "MueLu_Utilities.hpp"
62#include "MueLu_Monitor.hpp"
63
64namespace 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>
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
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Amesos_BaseSolver * Create(const char *ClassType, const Epetra_LinearProblem &LinearProblem)
Class that encapsulates Amesos direct solvers.
AmesosSmoother(const std::string &type="", const Teuchos::ParameterList &paramList=Teuchos::ParameterList())
Constructor.
std::string description() const
Return a simple one-line description of this object.
std::string type_
amesos-specific key phrase that denote smoother type
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void Setup(Level &currentLevel)
Set up the direct solver. This creates the underlying Amesos solver object according to the parameter...
RCP< SmootherPrototype > Copy() const
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
void DeclareInput(Level &currentLevel) const
Input.
void Apply(MultiVector &X, const MultiVector &B, bool=false) const
Apply the direct solver.
virtual std::string description() const
Return a simple one-line description of this object.
Exception throws to report errors in the internal logical of the program.
Timer to be used in factories. Similar to Monitor but with additional timers.
Class that holds all level-specific information.
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
void declareConstructionOutcome(bool fail, std::string msg)
bool IsSetup() const
Get the state of a smoother prototype.
static RCP< Epetra_MultiVector > MV2NonConstEpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
static RCP< Epetra_CrsMatrix > Op2NonConstEpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
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.
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
Namespace for MueLu classes and methods.
@ Warnings0
Important warning messages (one line)
@ Debug
Print additional debugging information.
@ External
Print external lib objects.
@ Runtime1
Description of what is happening (more verbose)
@ Parameters0
Print class parameters.
@ Parameters1
Print class parameters (more parameters, more verbose)