MueLu Version of the Day
Loading...
Searching...
No Matches
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
63namespace 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 ParameterList norecurse;
118 norecurse.disableRecursiveValidation();
119 validParamList->set<ParameterList> ("Amesos2", norecurse, "Parameters that are passed to Amesos2");
120 return validParamList;
121 }
122
123 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
125 ParameterList pL = this->GetParameterList();
126
127 this->Input(currentLevel, "A");
128 if (pL.get<bool>("fix nullspace"))
129 this->Input(currentLevel, "Nullspace");
130 }
131
132 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
134 FactoryMonitor m(*this, "Setup Smoother", currentLevel);
135
136 if (SmootherPrototype::IsSetup() == true)
137 this->GetOStream(Warnings0) << "MueLu::Amesos2Smoother::Setup(): Setup() has already been called" << std::endl;
138
139 RCP<Matrix> A = Factory::Get< RCP<Matrix> >(currentLevel, "A");
140
141 // Do a quick check if we need to modify the matrix
142 RCP<const Map> rowMap = A->getRowMap();
143 RCP<Matrix> factorA;
144 Teuchos::ParameterList pL = this->GetParameterList();
145 if (pL.get<bool>("fix nullspace")) {
146 this->GetOStream(Runtime1) << "MueLu::Amesos2Smoother::Setup(): fixing nullspace" << std::endl;
147
148 rowMap = A->getRowMap();
149 size_t M = rowMap->getGlobalNumElements();
150
151 RCP<MultiVector> Nullspace = Factory::Get< RCP<MultiVector> >(currentLevel, "Nullspace");
152
153 TEUCHOS_TEST_FOR_EXCEPTION(Nullspace->getNumVectors() > 1, Exceptions::RuntimeError,
154 "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 for nullspace of dim > 1 has not been implemented yet.");
155
156 RCP<MultiVector> NullspaceImp;
157 RCP<const Map> colMap;
158 RCP<const Import> importer;
159 if (rowMap->getComm()->getSize() > 1) {
160 this->GetOStream(Warnings0) << "MueLu::Amesos2Smoother::Setup(): Applying nullspace fix on distributed matrix. Try rebalancing to single rank!" << std::endl;
161 ArrayRCP<GO> elements_RCP;
162 elements_RCP.resize(M);
163 ArrayView<GO> elements = elements_RCP();
164 for (size_t k = 0; k<M; k++)
165 elements[k] = Teuchos::as<GO>(k);
166 colMap = MapFactory::Build(rowMap->lib(),M*rowMap->getComm()->getSize(),elements,Teuchos::ScalarTraits<GO>::zero(),rowMap->getComm());
167 importer = ImportFactory::Build(rowMap,colMap);
168 NullspaceImp = MultiVectorFactory::Build(colMap, Nullspace->getNumVectors());
169 NullspaceImp->doImport(*Nullspace,*importer,Xpetra::INSERT);
170 } else {
171 NullspaceImp = Nullspace;
172 colMap = rowMap;
173 }
174
175 ArrayRCP<const SC> nullspaceRCP, nullspaceImpRCP;
176 RCP<CrsMatrixWrap> Acrs = rcp_dynamic_cast<CrsMatrixWrap>(A);
177
178 TEUCHOS_TEST_FOR_EXCEPTION(Acrs.is_null(), Exceptions::RuntimeError,
179 "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 when matrix is not a Crs matrix has not been implemented yet.");
180
181 ArrayRCP<const size_t> rowPointers;
182 ArrayRCP<const LO> colIndices;
183 ArrayRCP<const SC> values;
184 Acrs->getCrsMatrix()->getAllValues(rowPointers, colIndices, values);
185
186 ArrayRCP<size_t> newRowPointers_RCP;
187 ArrayRCP<LO> newColIndices_RCP;
188 ArrayRCP<SC> newValues_RCP;
189
190 size_t N = rowMap->getLocalNumElements();
191 newRowPointers_RCP.resize(N+1);
192 newColIndices_RCP.resize(N*M);
193 newValues_RCP.resize(N*M);
194
195 ArrayView<size_t> newRowPointers = newRowPointers_RCP();
196 ArrayView<LO> newColIndices = newColIndices_RCP();
197 ArrayView<SC> newValues = newValues_RCP();
198
199 SC normalization = Nullspace->getVector(0)->norm2();
200 normalization = Teuchos::ScalarTraits<Scalar>::one()/(normalization*normalization);
201
202 ArrayView<const SC> nullspace, nullspaceImp;
203 nullspaceRCP = Nullspace->getData(0);
204 nullspace = nullspaceRCP();
205 nullspaceImpRCP = NullspaceImp->getData(0);
206 nullspaceImp = nullspaceImpRCP();
207
208 // form nullspace * nullspace^T
209 for (size_t i = 0; i < N; i++) {
210 newRowPointers[i] = i*M;
211 for (size_t j = 0; j < M; j++) {
212 newColIndices[i*M+j] = Teuchos::as<LO>(j);
213 newValues[i*M+j] = normalization * nullspace[i]*Teuchos::ScalarTraits<Scalar>::conjugate(nullspaceImp[j]);
214 }
215 }
216 newRowPointers[N] = N*M;
217
218 // add A
219 for (size_t i = 0; i < N; i++) {
220 for (size_t jj = rowPointers[i]; jj < rowPointers[i+1]; jj++) {
221 LO j = colMap->getLocalElement(A->getColMap()->getGlobalElement(colIndices[jj]));
222 SC v = values[jj];
223 newValues[i*M+j] += v;
224 }
225 }
226
227 RCP<Matrix> newA = rcp(new CrsMatrixWrap(rowMap, colMap, 0));
228 RCP<CrsMatrix> newAcrs = rcp_dynamic_cast<CrsMatrixWrap>(newA)->getCrsMatrix();
229
230 using Teuchos::arcp_const_cast;
231 newAcrs->setAllValues(arcp_const_cast<size_t>(newRowPointers_RCP), arcp_const_cast<LO>(newColIndices_RCP), arcp_const_cast<SC>(newValues_RCP));
232 newAcrs->expertStaticFillComplete(A->getDomainMap(), A->getRangeMap(),
233 importer, A->getCrsGraph()->getExporter());
234
235 factorA = newA;
236 } else {
237 factorA = A;
238 }
239
240 RCP<Tpetra_CrsMatrix> tA = Utilities::Op2NonConstTpetraCrs(factorA);
241
242 prec_ = Amesos2::create<Tpetra_CrsMatrix,Tpetra_MultiVector>(type_, tA);
243 TEUCHOS_TEST_FOR_EXCEPTION(prec_ == Teuchos::null, Exceptions::RuntimeError, "Amesos2::create returns Teuchos::null");
244 RCP<Teuchos::ParameterList> amesos2_params = Teuchos::rcpFromRef(pL.sublist("Amesos2"));
245 amesos2_params->setName("Amesos2");
246 if (rowMap->getGlobalNumElements() != as<size_t>((rowMap->getMaxAllGlobalIndex() - rowMap->getMinAllGlobalIndex())+1)) {
247 if (!(amesos2_params->sublist(prec_->name()).template isType<bool>("IsContiguous")))
248 amesos2_params->sublist(prec_->name()).set("IsContiguous", false, "Are GIDs Contiguous");
249 }
250 prec_->setParameters(amesos2_params);
251
253 }
254
255 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
256 void Amesos2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Apply(MultiVector& X, const MultiVector& B, bool /* InitialGuessIsZero */) const {
257 TEUCHOS_TEST_FOR_EXCEPTION(SmootherPrototype::IsSetup() == false, Exceptions::RuntimeError, "MueLu::Amesos2Smoother::Apply(): Setup() has not been called");
258
259 RCP<Tpetra_MultiVector> tX, tB;
260 if (!useTransformation_) {
262 tB = Utilities::MV2NonConstTpetraMV2(const_cast<MultiVector&>(B));
263 } else {
264 // Copy data of the original vectors into the transformed ones
265 size_t numVectors = X.getNumVectors();
266 size_t length = X.getLocalLength();
267
268 TEUCHOS_TEST_FOR_EXCEPTION(numVectors > 1, Exceptions::RuntimeError,
269 "MueLu::Amesos2Smoother::Apply: Fixing coarse matrix for Amesos2 for multivectors has not been implemented yet.");
270 ArrayRCP<const SC> Xdata = X. getData(0), Bdata = B. getData(0);
271 ArrayRCP<SC> X_data = X_->getDataNonConst(0), B_data = B_->getDataNonConst(0);
272
273 for (size_t i = 0; i < length; i++) {
274 X_data[i] = Xdata[i];
275 B_data[i] = Bdata[i];
276 }
277
280 }
281
282 prec_->setX(tX);
283 prec_->setB(tB);
284
285 prec_->solve();
286
287 prec_->setX(Teuchos::null);
288 prec_->setB(Teuchos::null);
289
290 if (useTransformation_) {
291 // Copy data from the transformed vectors into the original ones
292 size_t length = X.getLocalLength();
293
294 ArrayRCP<SC> Xdata = X. getDataNonConst(0);
295 ArrayRCP<const SC> X_data = X_->getData(0);
296
297 for (size_t i = 0; i < length; i++)
298 Xdata[i] = X_data[i];
299 }
300 }
301
302 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
303 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
305 Copy() const
306 {
307 return rcp (new Amesos2Smoother (*this));
308 }
309
310 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
312 std::ostringstream out;
313
314 if (SmootherPrototype::IsSetup() == true) {
315 out << prec_->description();
316
317 } else {
319 out << "{type = " << type_ << "}";
320 }
321 return out.str();
322 }
323
324 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
325 void Amesos2Smoother<Scalar, LocalOrdinal, GlobalOrdinal, Node>::print(Teuchos::FancyOStream& out, const VerbLevel verbLevel) const {
327
328 if (verbLevel & Parameters0)
329 out0 << "Prec. type: " << type_ << std::endl;
330
331 if (verbLevel & Parameters1) {
332 out0 << "Parameter list: " << std::endl;
333 Teuchos::OSTab tab2(out);
334 out << this->GetParameterList();
335 }
336
337 if ((verbLevel & External) && prec_ != Teuchos::null) {
338 Teuchos::OSTab tab2(out);
339 out << *prec_ << std::endl;
340 }
341
342 if (verbLevel & Debug)
343 out0 << "IsSetup: " << Teuchos::toString(SmootherPrototype::IsSetup()) << std::endl
344 << "-" << std::endl
345 << "RCP<prec_>: " << prec_ << std::endl;
346 }
347
348 template <class Scalar,class LocalOrdinal, class GlobalOrdinal, class Node>
350 if(!prec_.is_null())
351 return prec_->getStatus().getNnzLU();
352 else
353 return 0.0;
354
355 }
356} // namespace MueLu
357
358#endif // HAVE_MUELU_TPETRA && HAVE_MUELU_AMESOS2
359#endif // MUELU_AMESOS2SMOOTHER_DEF_HPP
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Class that encapsulates Amesos2 direct solvers.
void DeclareInput(Level &currentLevel) const
Input.
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
std::string type_
amesos2-specific key phrase that denote smoother type
RCP< SmootherPrototype > Copy() const
std::string description() const
Return a simple one-line description of this object.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
virtual ~Amesos2Smoother()
Destructor.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
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.
Amesos2Smoother(const std::string &type="", const Teuchos::ParameterList &paramList=Teuchos::ParameterList())
Constructor Creates a MueLu interface to the direct solvers in the Amesos2 package....
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.
Definition: MueLu_Level.hpp:99
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< 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)
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)