MueLu  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MueLu_Amesos2Smoother_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_AMESOS2SMOOTHER_DEF_HPP
47 #define MUELU_AMESOS2SMOOTHER_DEF_HPP
48 
49 #include <algorithm>
50 
51 #include "MueLu_ConfigDefs.hpp"
52 #if defined (HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_AMESOS2)
53 #include <Xpetra_Matrix.hpp>
54 
55 #include <Amesos2_config.h>
56 #include <Amesos2.hpp>
57 
59 #include "MueLu_Level.hpp"
60 #include "MueLu_Utilities.hpp"
61 #include "MueLu_Monitor.hpp"
62 
63 namespace MueLu {
64 
65  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
67  : type_(type), useTransformation_(false) {
68  this->SetParameterList(paramList);
69 
70  if (!type_.empty()) {
71  // Transform string to "Abcde" notation
72  std::transform(type_.begin(), type_.end(), type_.begin(), ::tolower);
73  std::transform(type_.begin(), ++type_.begin(), type_.begin(), ::toupper);
74  }
75  if (type_ == "Superlu_dist")
76  type_ = "Superludist";
77 
78  // Try to come up with something availble
79  // Order corresponds to our preference
80  // TODO: It would be great is Amesos2 provides directly this kind of logic for us
81  if (type_ == "" || Amesos2::query(type_) == false) {
82  std::string oldtype = type_;
83 #if defined(HAVE_AMESOS2_SUPERLU)
84  type_ = "Superlu";
85 #elif defined(HAVE_AMESOS2_KLU2)
86  type_ = "Klu";
87 #elif defined(HAVE_AMESOS2_SUPERLUDIST)
88  type_ = "Superludist";
89 #elif defined(HAVE_AMESOS2_BASKER)
90  type_ = "Basker";
91 #else
92  throw Exceptions::RuntimeError("Amesos2 has been compiled without SuperLU_DIST, SuperLU, Klu, or Basker. By default, MueLu tries"
93  "to use one of these libraries. Amesos2 must be compiled with one of these solvers, "
94  "or a valid Amesos2 solver has to be specified explicitly.");
95 #endif
96  if (oldtype != "")
97  this->GetOStream(Warnings0) << "MueLu::Amesos2Smoother: \"" << oldtype << "\" is not available. Using \"" << type_ << "\" instead" << std::endl;
98  else
99  this->GetOStream(Runtime1) << "MueLu::Amesos2Smoother: using \"" << type_ << "\"" << std::endl;
100  }
101 
102  // Check the validity of the solver type parameter
103  TEUCHOS_TEST_FOR_EXCEPTION(Amesos2::query(type_) == false, Exceptions::RuntimeError, "The Amesos2 library reported that the solver '" << type_ << "' is not available. "
104  "Amesos2 has been compiled without the support of this solver, or the solver name is misspelled.");
105  }
106 
107  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
109 
110  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
112  this->Input(currentLevel, "A");
113  }
114 
115  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
117  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
118 
119  if (SmootherPrototype::IsSetup() == true)
120  this->GetOStream(Warnings0) << "MueLu::Amesos2Smoother::Setup(): Setup() has already been called" << std::endl;
121 
122  RCP<Matrix> A = Factory::Get< RCP<Matrix> >(currentLevel, "A");
123 
124  // Do a quick check if we need to modify the matrix
125  RCP<const Map> rowMap = A->getRowMap();
126  if (rowMap->getGlobalNumElements() != as<size_t>((rowMap->getMaxAllGlobalIndex() - rowMap->getMinAllGlobalIndex())+1)) {
127  // If our system is non-conventional, here is the place where we try to fix it
128  // One example is: if our maps contain a gap in them, for instance GIDs
129  // are [0,..., 100, 10000, ..., 10010], then Amesos2 breaks down.
130  //
131  // The approach we take is to construct a second system with maps
132  // replaced by their continuous versions.
133  //
134  // FIXME: for the moment, this functionality works for selected limited scenarios
135  this->GetOStream(Runtime1) << "MueLu::Amesos2Smoother::Setup(): using system transformation" << std::endl;
136 
137  TEUCHOS_TEST_FOR_EXCEPTION(rowMap->getComm()->getSize() > 1, Exceptions::RuntimeError,
138  "MueLu::Amesos2Smoother::Setup Fixing coarse matrix for Amesos2 for multiple processors has not been implemented yet.");
139  TEUCHOS_TEST_FOR_EXCEPTION(!rowMap->isSameAs(*A->getColMap()), Exceptions::RuntimeError,
140  "MueLu::Amesos2Smoother::Setup Fixing coarse matrix for Amesos2 when row map is different from column map has not been implemented yet.");
141 
142  RCP<CrsMatrixWrap> Acrs = rcp_dynamic_cast<CrsMatrixWrap>(A);
143  TEUCHOS_TEST_FOR_EXCEPTION(Acrs.is_null(), Exceptions::RuntimeError,
144  "MueLu::Amesos2Smoother::Setup Fixing coarse matrix for Amesos2 when matrix is not a Crs matrix has not been implemented yet.");
145 
146  useTransformation_ = true;
147 
148  ArrayRCP<const size_t> rowPointers;
149  ArrayRCP<const LO> colIndices;
150  ArrayRCP<const SC> values;
151  Acrs->getCrsMatrix()->getAllValues(rowPointers, colIndices, values);
152 
153  // Create new map
154  RCP<Map> map = MapFactory::Build(rowMap->lib(), rowMap->getGlobalNumElements(), 0, rowMap->getComm());
155  RCP<Matrix> newA = rcp(new CrsMatrixWrap(map, map, 0, Xpetra::StaticProfile));
156  RCP<CrsMatrix> newAcrs = rcp_dynamic_cast<CrsMatrixWrap>(newA)->getCrsMatrix();
157 
158  using Teuchos::arcp_const_cast;
159  newAcrs->setAllValues(arcp_const_cast<size_t>(rowPointers), arcp_const_cast<LO>(colIndices), arcp_const_cast<SC>(values));
160  newAcrs->expertStaticFillComplete(map, map);
161 
162  A.swap(newA);
163 
164  X_ = MultiVectorFactory::Build(map, 1);
165  B_ = MultiVectorFactory::Build(map, 1);
166  }
167 
169 
170  prec_ = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(type_, tA);
171  TEUCHOS_TEST_FOR_EXCEPTION(prec_ == Teuchos::null, Exceptions::RuntimeError, "Amesos2::create returns Teuchos::null");
172 
174  }
175 
176  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
177  void Amesos2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool /* InitialGuessIsZero */) const {
178  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::Amesos2Smoother::Apply(): Setup() has not been called");
179 
181  if (!useTransformation_) {
183  tB = Utilities::MV2NonConstTpetraMV2(const_cast<MultiVector&>(B));
184  } else {
185  // Copy data of the original vectors into the transformed ones
186  size_t numVectors = X.getNumVectors();
187  size_t length = X.getLocalLength();
188 
190  "MueLu::Amesos2Smoother::Apply: Fixing coarse matrix for Amesos2 for multivectors has not been implemented yet.");
191  ArrayRCP<const SC> Xdata = X. getData(0), Bdata = B. getData(0);
192  ArrayRCP<SC> X_data = X_->getDataNonConst(0), B_data = B_->getDataNonConst(0);
193 
194  for (size_t i = 0; i < length; i++) {
195  X_data[i] = Xdata[i];
196  B_data[i] = Bdata[i];
197  }
198 
201  }
202 
203  prec_->setX(tX);
204  prec_->setB(tB);
205 
206  prec_->solve();
207 
208  prec_->setX(Teuchos::null);
209  prec_->setB(Teuchos::null);
210 
211  if (useTransformation_) {
212  // Copy data from the transformed vectors into the original ones
213  size_t length = X.getLocalLength();
214 
215  ArrayRCP<SC> Xdata = X. getDataNonConst(0);
216  ArrayRCP<const SC> X_data = X_->getData(0);
217 
218  for (size_t i = 0; i < length; i++)
219  Xdata[i] = X_data[i];
220  }
221  }
222 
223  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
226  Copy() const
227  {
228  return rcp (new Amesos2Smoother (*this));
229  }
230 
231  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
233  std::ostringstream out;
234 
235  if (SmootherPrototype::IsSetup() == true) {
236  out << prec_->description();
237 
238  } else {
240  out << "{type = " << type_ << "}";
241  }
242  return out.str();
243  }
244 
245  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
248 
249  if (verbLevel & Parameters0)
250  out0 << "Prec. type: " << type_ << std::endl;
251 
252  if (verbLevel & Parameters1) {
253  out0 << "Parameter list: " << std::endl;
254  Teuchos::OSTab tab2(out);
255  out << this->GetParameterList();
256  }
257 
258  if ((verbLevel & External) && prec_ != Teuchos::null) {
259  Teuchos::OSTab tab2(out);
260  out << *prec_ << std::endl;
261  }
262 
263  if (verbLevel & Debug)
264  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
265  << "-" << std::endl
266  << "RCP<prec_>: " << prec_ << std::endl;
267  }
268 
269  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
271  if(!prec_.is_null())
272  return prec_->getStatus().getNnzLU();
273  else
274  return 0.0;
275 
276  }
277 } // namespace MueLu
278 
279 #endif // HAVE_MUELU_TPETRA && HAVE_MUELU_AMESOS2
280 #endif // MUELU_AMESOS2SMOOTHER_DEF_HPP
Important warning messages (one line)
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
Print external lib objects.
std::string type_
amesos2-specific key phrase that denote smoother type
void Apply(MultiVector &X, const MultiVector &B, bool InitialGuessIsZero=false) const
Apply the direct solver. Solves the linear system AX=B using the constructed solver.
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
void swap(RCP< T > &r_ptr)
Print additional debugging information.
std::string tolower(const std::string &str)
Amesos2Smoother(const std::string &type="", const Teuchos::ParameterList &paramList=Teuchos::ParameterList())
Constructor Creates a MueLu interface to the direct solvers in the Amesos2 package. If you are using type==&quot;&quot;, then either SuperLU or KLU2 are used by default.
std::string description() const
Return a simple one-line description of this object.
RCP< SmootherPrototype > Copy() 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)
Class that holds all level-specific information.
Definition: MueLu_Level.hpp:99
Class that encapsulates Amesos2 direct solvers.
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
bool IsSetup() const
Get the state of a smoother prototype.
void DeclareInput(Level &currentLevel) const
Input.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Print class parameters.
virtual ~Amesos2Smoother()
Destructor.
void Setup(Level &currentLevel)
Set up the direct solver. This creates the underlying Amesos2 solver object according to the paramete...
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)
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
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 is_null() const