46#ifndef MUELU_AMESOS2SMOOTHER_DEF_HPP
47#define MUELU_AMESOS2SMOOTHER_DEF_HPP
52#if defined (HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_AMESOS2)
53#include <Xpetra_Matrix.hpp>
55#include <Amesos2_config.h>
60#include "MueLu_Utilities.hpp"
65 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
67 : type_(type), useTransformation_(false) {
73 std::transform(
type_.begin(), ++
type_.begin(),
type_.begin(), ::toupper);
75 if (
type_ ==
"Superlu_dist")
76 type_ =
"Superludist";
81 if (
type_ ==
"" || Amesos2::query(
type_) ==
false) {
82 std::string oldtype =
type_;
83#if defined(HAVE_AMESOS2_SUPERLU)
85#elif defined(HAVE_AMESOS2_KLU2)
87#elif defined(HAVE_AMESOS2_SUPERLUDIST)
88 type_ =
"Superludist";
89#elif defined(HAVE_AMESOS2_BASKER)
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.");
98 this->
GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother: \"" << oldtype <<
"\" is not available. Using \"" <<
type_ <<
"\" instead" << std::endl;
105 "Amesos2 has been compiled without the support of this solver, or the solver name is misspelled.");
108 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
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;
123 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
125 ParameterList pL = this->GetParameterList();
127 this->Input(currentLevel,
"A");
128 if (pL.get<
bool>(
"fix nullspace"))
129 this->Input(currentLevel,
"Nullspace");
132 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
137 this->GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother::Setup(): Setup() has already been called" << std::endl;
139 RCP<Matrix> A = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
142 RCP<const Map> rowMap = A->getRowMap();
144 Teuchos::ParameterList pL = this->GetParameterList();
145 if (pL.get<
bool>(
"fix nullspace")) {
146 this->GetOStream(
Runtime1) <<
"MueLu::Amesos2Smoother::Setup(): fixing nullspace" << std::endl;
148 rowMap = A->getRowMap();
149 size_t M = rowMap->getGlobalNumElements();
151 RCP<MultiVector> Nullspace = Factory::Get< RCP<MultiVector> >(currentLevel,
"Nullspace");
154 "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 for nullspace of dim > 1 has not been implemented yet.");
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);
171 NullspaceImp = Nullspace;
175 ArrayRCP<const SC> nullspaceRCP, nullspaceImpRCP;
176 RCP<CrsMatrixWrap> Acrs = rcp_dynamic_cast<CrsMatrixWrap>(A);
179 "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 when matrix is not a Crs matrix has not been implemented yet.");
181 ArrayRCP<const size_t> rowPointers;
182 ArrayRCP<const LO> colIndices;
183 ArrayRCP<const SC> values;
184 Acrs->getCrsMatrix()->getAllValues(rowPointers, colIndices, values);
186 ArrayRCP<size_t> newRowPointers_RCP;
187 ArrayRCP<LO> newColIndices_RCP;
188 ArrayRCP<SC> newValues_RCP;
190 size_t N = rowMap->getLocalNumElements();
191 newRowPointers_RCP.resize(N+1);
192 newColIndices_RCP.resize(N*M);
193 newValues_RCP.resize(N*M);
195 ArrayView<size_t> newRowPointers = newRowPointers_RCP();
196 ArrayView<LO> newColIndices = newColIndices_RCP();
197 ArrayView<SC> newValues = newValues_RCP();
199 SC normalization = Nullspace->getVector(0)->norm2();
200 normalization = Teuchos::ScalarTraits<Scalar>::one()/(normalization*normalization);
202 ArrayView<const SC> nullspace, nullspaceImp;
203 nullspaceRCP = Nullspace->getData(0);
204 nullspace = nullspaceRCP();
205 nullspaceImpRCP = NullspaceImp->getData(0);
206 nullspaceImp = nullspaceImpRCP();
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]);
216 newRowPointers[N] = N*M;
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]));
223 newValues[i*M+j] += v;
227 RCP<Matrix> newA = rcp(
new CrsMatrixWrap(rowMap, colMap, 0));
228 RCP<CrsMatrix> newAcrs = rcp_dynamic_cast<CrsMatrixWrap>(newA)->getCrsMatrix();
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());
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");
250 prec_->setParameters(amesos2_params);
255 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
259 RCP<Tpetra_MultiVector> tX, tB;
260 if (!useTransformation_) {
265 size_t numVectors = X.getNumVectors();
266 size_t length = X.getLocalLength();
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);
273 for (
size_t i = 0; i < length; i++) {
274 X_data[i] = Xdata[i];
275 B_data[i] = Bdata[i];
287 prec_->setX(Teuchos::null);
288 prec_->setB(Teuchos::null);
290 if (useTransformation_) {
292 size_t length = X.getLocalLength();
294 ArrayRCP<SC> Xdata = X. getDataNonConst(0);
295 ArrayRCP<const SC> X_data = X_->getData(0);
297 for (
size_t i = 0; i < length; i++)
298 Xdata[i] = X_data[i];
302 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
303 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
310 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
312 std::ostringstream out;
315 out << prec_->description();
319 out <<
"{type = " << type_ <<
"}";
324 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
329 out0 <<
"Prec. type: " << type_ << std::endl;
332 out0 <<
"Parameter list: " << std::endl;
333 Teuchos::OSTab tab2(out);
334 out << this->GetParameterList();
337 if ((verbLevel &
External) && prec_ != Teuchos::null) {
338 Teuchos::OSTab tab2(out);
339 out << *prec_ << std::endl;
342 if (verbLevel &
Debug)
345 <<
"RCP<prec_>: " << prec_ << std::endl;
348 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
351 return prec_->getStatus().getNnzLU();
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
Class that encapsulates Amesos2 direct solvers.
void DeclareInput(Level ¤tLevel) 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 ¤tLevel)
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 ¶mList=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.
virtual void SetParameterList(const Teuchos::ParameterList ¶mList)
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)