46 #ifndef MUELU_AMESOS2SMOOTHER_DEF_HPP 47 #define MUELU_AMESOS2SMOOTHER_DEF_HPP 52 #if defined (HAVE_MUELU_TPETRA) && defined(HAVE_MUELU_AMESOS2) 55 #include <Amesos2_config.h> 56 #include <Amesos2.hpp> 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 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.");
97 this->
GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother: \"" << oldtype <<
"\" is not available. Using \"" <<
type_ <<
"\" instead" << std::endl;
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.");
107 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
110 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
112 this->Input(currentLevel,
"A");
115 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
120 this->GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother::Setup(): Setup() has already been called" << std::endl;
122 RCP<Matrix> A = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
125 RCP<const Map> rowMap = A->getRowMap();
126 if (rowMap->getGlobalNumElements() != as<size_t>((rowMap->getMaxAllGlobalIndex() - rowMap->getMinAllGlobalIndex())+1)) {
135 this->GetOStream(
Runtime1) <<
"MueLu::Amesos2Smoother::Setup(): using system transformation" << std::endl;
138 "MueLu::Amesos2Smoother::Setup Fixing coarse matrix for Amesos2 for multiple processors has not been implemented yet.");
140 "MueLu::Amesos2Smoother::Setup Fixing coarse matrix for Amesos2 when row map is different from column map has not been implemented yet.");
142 RCP<CrsMatrixWrap> Acrs = rcp_dynamic_cast<CrsMatrixWrap>(A);
144 "MueLu::Amesos2Smoother::Setup Fixing coarse matrix for Amesos2 when matrix is not a Crs matrix has not been implemented yet.");
146 useTransformation_ =
true;
148 ArrayRCP<const size_t> rowPointers;
149 ArrayRCP<const LO> colIndices;
150 ArrayRCP<const SC> values;
151 Acrs->getCrsMatrix()->getAllValues(rowPointers, colIndices, values);
154 RCP<Map> map = MapFactory::Build(rowMap->lib(), rowMap->getGlobalNumElements(), 0, rowMap->getComm());
156 RCP<CrsMatrix> newAcrs = rcp_dynamic_cast<CrsMatrixWrap>(newA)->getCrsMatrix();
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);
164 X_ = MultiVectorFactory::Build(map, 1);
165 B_ = MultiVectorFactory::Build(map, 1);
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");
176 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
180 RCP<Tpetra_MultiVector> tX, tB;
181 if (!useTransformation_) {
186 size_t numVectors = X.getNumVectors();
187 size_t length = X.getLocalLength();
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);
194 for (
size_t i = 0; i < length; i++) {
195 X_data[i] = Xdata[i];
196 B_data[i] = Bdata[i];
208 prec_->setX(Teuchos::null);
209 prec_->setB(Teuchos::null);
211 if (useTransformation_) {
213 size_t length = X.getLocalLength();
215 ArrayRCP<SC> Xdata = X. getDataNonConst(0);
216 ArrayRCP<const SC> X_data = X_->getData(0);
218 for (
size_t i = 0; i < length; i++)
219 Xdata[i] = X_data[i];
223 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
224 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
231 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
233 std::ostringstream out;
236 out << prec_->description();
240 out <<
"{type = " << type_ <<
"}";
245 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
250 out0 <<
"Prec. type: " << type_ << std::endl;
253 out0 <<
"Parameter list: " << std::endl;
254 Teuchos::OSTab tab2(out);
255 out << this->GetParameterList();
258 if ((verbLevel &
External) && prec_ != Teuchos::null) {
259 Teuchos::OSTab tab2(out);
260 out << *prec_ << std::endl;
263 if (verbLevel &
Debug)
266 <<
"RCP<prec_>: " << prec_ << std::endl;
271 #endif // HAVE_MUELU_TPETRA && HAVE_MUELU_AMESOS2 272 #endif // MUELU_AMESOS2SMOOTHER_DEF_HPP Important warning messages (one line)
void DeclareInput(Level ¤tLevel) 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.
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.
Namespace for MueLu classes and methods.
Amesos2Smoother(const std::string &type="", const Teuchos::ParameterList ¶mList=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 ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
std::string description() const
Return a simple one-line description of this object.
Class that holds all level-specific information.
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)
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
virtual ~Amesos2Smoother()
Destructor.
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.
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.