MueLu  Version of the Day
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>
66  Amesos2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Amesos2Smoother(const std::string& type, const Teuchos::ParameterList& paramList)
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  this->declareConstructionOutcome(true, std::string("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  return;
96 #endif
97  if (oldtype != "")
98  this->GetOStream(Warnings0) << "MueLu::Amesos2Smoother: \"" << oldtype << "\" is not available. Using \"" << type_ << "\" instead" << std::endl;
99  else
100  this->GetOStream(Runtime1) << "MueLu::Amesos2Smoother: using \"" << type_ << "\"" << std::endl;
101  }
102 
103  // Check the validity of the solver type parameter
104  this->declareConstructionOutcome(Amesos2::query(type_) == false, "The Amesos2 library reported that the solver '" + type_ + "' is not available. " +
105  "Amesos2 has been compiled without the support of this solver, or the solver name is misspelled.");
106  }
107 
108  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
110 
111  template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
113  RCP<ParameterList> validParamList = rcp(new ParameterList());
114  validParamList->set< RCP<const FactoryBase> >("A", null, "Factory of the coarse matrix");
115  validParamList->set< RCP<const FactoryBase> >("Nullspace", null, "Factory of the nullspace");
116  validParamList->set<bool>("fix nullspace", false, "Remove zero eigenvalue by adding rank one correction.");
117  return validParamList;
118  }
119 
120  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
122  ParameterList pL = this->GetParameterList();
123 
124  this->Input(currentLevel, "A");
125  if (pL.get<bool>("fix nullspace"))
126  this->Input(currentLevel, "Nullspace");
127  }
128 
129  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
131  FactoryMonitor m(*this, "Setup Smoother", currentLevel);
132 
133  if (SmootherPrototype::IsSetup() == true)
134  this->GetOStream(Warnings0) << "MueLu::Amesos2Smoother::Setup(): Setup() has already been called" << std::endl;
135 
136  RCP<Matrix> A = Factory::Get< RCP<Matrix> >(currentLevel, "A");
137 
138  // Do a quick check if we need to modify the matrix
139  RCP<const Map> rowMap = A->getRowMap();
140  RCP<Matrix> factorA;
141  Teuchos::ParameterList pL = this->GetParameterList();
142  if (pL.get<bool>("fix nullspace")) {
143  this->GetOStream(Runtime1) << "MueLu::Amesos2Smoother::Setup(): fixing nullspace" << std::endl;
144 
145  rowMap = A->getRowMap();
146  size_t M = rowMap->getGlobalNumElements();
147 
148  RCP<MultiVector> Nullspace = Factory::Get< RCP<MultiVector> >(currentLevel, "Nullspace");
149 
150  TEUCHOS_TEST_FOR_EXCEPTION(Nullspace->getNumVectors() > 1, Exceptions::RuntimeError,
151  "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 for nullspace of dim > 1 has not been implemented yet.");
152 
153  RCP<MultiVector> NullspaceImp;
154  RCP<const Map> colMap;
155  RCP<const Import> importer;
156  if (rowMap->getComm()->getSize() > 1) {
157  this->GetOStream(Warnings0) << "MueLu::Amesos2Smoother::Setup(): Applying nullspace fix on distributed matrix. Try rebalancing to single rank!" << std::endl;
158  ArrayRCP<GO> elements_RCP;
159  elements_RCP.resize(M);
160  ArrayView<GO> elements = elements_RCP();
161  for (size_t k = 0; k<M; k++)
162  elements[k] = Teuchos::as<GO>(k);
163  colMap = MapFactory::Build(rowMap->lib(),M*rowMap->getComm()->getSize(),elements,Teuchos::ScalarTraits<GO>::zero(),rowMap->getComm());
164  importer = ImportFactory::Build(rowMap,colMap);
165  NullspaceImp = MultiVectorFactory::Build(colMap, Nullspace->getNumVectors());
166  NullspaceImp->doImport(*Nullspace,*importer,Xpetra::INSERT);
167  } else {
168  NullspaceImp = Nullspace;
169  colMap = rowMap;
170  }
171 
172  ArrayRCP<const SC> nullspaceRCP, nullspaceImpRCP;
173  RCP<CrsMatrixWrap> Acrs = rcp_dynamic_cast<CrsMatrixWrap>(A);
174 
175  TEUCHOS_TEST_FOR_EXCEPTION(Acrs.is_null(), Exceptions::RuntimeError,
176  "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 when matrix is not a Crs matrix has not been implemented yet.");
177 
178  ArrayRCP<const size_t> rowPointers;
179  ArrayRCP<const LO> colIndices;
180  ArrayRCP<const SC> values;
181  Acrs->getCrsMatrix()->getAllValues(rowPointers, colIndices, values);
182 
183  ArrayRCP<size_t> newRowPointers_RCP;
184  ArrayRCP<LO> newColIndices_RCP;
185  ArrayRCP<SC> newValues_RCP;
186 
187  size_t N = rowMap->getNodeNumElements();
188  newRowPointers_RCP.resize(N+1);
189  newColIndices_RCP.resize(N*M);
190  newValues_RCP.resize(N*M);
191 
192  ArrayView<size_t> newRowPointers = newRowPointers_RCP();
193  ArrayView<LO> newColIndices = newColIndices_RCP();
194  ArrayView<SC> newValues = newValues_RCP();
195 
196  SC normalization = Nullspace->getVector(0)->norm2();
197  normalization = Teuchos::ScalarTraits<Scalar>::one()/(normalization*normalization);
198 
199  ArrayView<const SC> nullspace, nullspaceImp;
200  nullspaceRCP = Nullspace->getData(0);
201  nullspace = nullspaceRCP();
202  nullspaceImpRCP = NullspaceImp->getData(0);
203  nullspaceImp = nullspaceImpRCP();
204 
205  // form nullspace * nullspace^T
206  for (size_t i = 0; i < N; i++) {
207  newRowPointers[i] = i*M;
208  for (size_t j = 0; j < M; j++) {
209  newColIndices[i*M+j] = Teuchos::as<LO>(j);
210  newValues[i*M+j] = normalization * nullspace[i]*Teuchos::ScalarTraits<Scalar>::conjugate(nullspaceImp[j]);
211  }
212  }
213  newRowPointers[N] = N*M;
214 
215  // add A
216  for (size_t i = 0; i < N; i++) {
217  for (size_t jj = rowPointers[i]; jj < rowPointers[i+1]; jj++) {
218  LO j = colMap->getLocalElement(A->getColMap()->getGlobalElement(colIndices[jj]));
219  SC v = values[jj];
220  newValues[i*M+j] += v;
221  }
222  }
223 
224  RCP<Matrix> newA = rcp(new CrsMatrixWrap(rowMap, colMap, 0));
225  RCP<CrsMatrix> newAcrs = rcp_dynamic_cast<CrsMatrixWrap>(newA)->getCrsMatrix();
226 
227  using Teuchos::arcp_const_cast;
228  newAcrs->setAllValues(arcp_const_cast<size_t>(newRowPointers_RCP), arcp_const_cast<LO>(newColIndices_RCP), arcp_const_cast<SC>(newValues_RCP));
229  newAcrs->expertStaticFillComplete(A->getDomainMap(), A->getRangeMap(),
230  importer, A->getCrsGraph()->getExporter());
231 
232  factorA = newA;
233  } else {
234  factorA = A;
235  }
236 
237  RCP<Tpetra_CrsMatrix> tA = Utilities::Op2NonConstTpetraCrs(factorA);
238 
239  prec_ = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(type_, tA);
240  TEUCHOS_TEST_FOR_EXCEPTION(prec_ == Teuchos::null, Exceptions::RuntimeError, "Amesos2::create returns Teuchos::null");
241  if (rowMap->getGlobalNumElements() != as<size_t>((rowMap->getMaxAllGlobalIndex() - rowMap->getMinAllGlobalIndex())+1)) {
242  auto amesos2_params = Teuchos::rcp(new Teuchos::ParameterList("Amesos2"));
243  amesos2_params->sublist(prec_->name()).set("IsContiguous", false, "Are GIDs Contiguous");
244  prec_->setParameters(amesos2_params);
245  }
246 
248  }
249 
250  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
251  void Amesos2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool /* InitialGuessIsZero */) const {
252  TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::Amesos2Smoother::Apply(): Setup() has not been called");
253 
254  RCP<Tpetra_MultiVector> tX, tB;
255  if (!useTransformation_) {
257  tB = Utilities::MV2NonConstTpetraMV2(const_cast<MultiVector&>(B));
258  } else {
259  // Copy data of the original vectors into the transformed ones
260  size_t numVectors = X.getNumVectors();
261  size_t length = X.getLocalLength();
262 
263  TEUCHOS_TEST_FOR_EXCEPTION(numVectors > 1, Exceptions::RuntimeError,
264  "MueLu::Amesos2Smoother::Apply: Fixing coarse matrix for Amesos2 for multivectors has not been implemented yet.");
265  ArrayRCP<const SC> Xdata = X. getData(0), Bdata = B. getData(0);
266  ArrayRCP<SC> X_data = X_->getDataNonConst(0), B_data = B_->getDataNonConst(0);
267 
268  for (size_t i = 0; i < length; i++) {
269  X_data[i] = Xdata[i];
270  B_data[i] = Bdata[i];
271  }
272 
275  }
276 
277  prec_->setX(tX);
278  prec_->setB(tB);
279 
280  prec_->solve();
281 
282  prec_->setX(Teuchos::null);
283  prec_->setB(Teuchos::null);
284 
285  if (useTransformation_) {
286  // Copy data from the transformed vectors into the original ones
287  size_t length = X.getLocalLength();
288 
289  ArrayRCP<SC> Xdata = X. getDataNonConst(0);
290  ArrayRCP<const SC> X_data = X_->getData(0);
291 
292  for (size_t i = 0; i < length; i++)
293  Xdata[i] = X_data[i];
294  }
295  }
296 
297  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
298  RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
300  Copy() const
301  {
302  return rcp (new Amesos2Smoother (*this));
303  }
304 
305  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
307  std::ostringstream out;
308 
309  if (SmootherPrototype::IsSetup() == true) {
310  out << prec_->description();
311 
312  } else {
314  out << "{type = " << type_ << "}";
315  }
316  return out.str();
317  }
318 
319  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
320  void Amesos2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
322 
323  if (verbLevel & Parameters0)
324  out0 << "Prec. type: " << type_ << std::endl;
325 
326  if (verbLevel & Parameters1) {
327  out0 << "Parameter list: " << std::endl;
328  Teuchos::OSTab tab2(out);
329  out << this->GetParameterList();
330  }
331 
332  if ((verbLevel & External) && prec_ != Teuchos::null) {
333  Teuchos::OSTab tab2(out);
334  out << *prec_ << std::endl;
335  }
336 
337  if (verbLevel & Debug)
338  out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
339  << "-" << std::endl
340  << "RCP<prec_>: " << prec_ << std::endl;
341  }
342 
343  template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
345  if(!prec_.is_null())
346  return prec_->getStatus().getNnzLU();
347  else
348  return 0.0;
349 
350  }
351 } // namespace MueLu
352 
353 #endif // HAVE_MUELU_TPETRA && HAVE_MUELU_AMESOS2
354 #endif // MUELU_AMESOS2SMOOTHER_DEF_HPP
Important warning messages (one line)
void DeclareInput(Level &currentLevel) const
Input.
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.
Print external lib objects.
bool IsSetup() const
Get the state of a smoother prototype.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
std::string type_
amesos2-specific key phrase that denote smoother type
Timer to be used in factories. Similar to Monitor but with additional timers.
Print additional debugging information.
std::string tolower(const std::string &str)
Namespace for MueLu classes and methods.
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=="", then either SuperLU or KLU2 are used by default.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
virtual void SetParameterList(const Teuchos::ParameterList &paramList)
Set parameters from a parameter list and return with default values.
std::string description() const
Return a simple one-line description of this object.
void declareConstructionOutcome(bool fail, std::string msg)
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)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
#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...
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.
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)
RCP< SmootherPrototype > Copy() const
virtual std::string description() const
Return a simple one-line description of this object.