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> 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 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 return validParamList;
120 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
122 ParameterList pL = this->GetParameterList();
124 this->Input(currentLevel,
"A");
125 if (pL.get<
bool>(
"fix nullspace"))
126 this->Input(currentLevel,
"Nullspace");
129 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
134 this->GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother::Setup(): Setup() has already been called" << std::endl;
136 RCP<Matrix> A = Factory::Get< RCP<Matrix> >(currentLevel,
"A");
139 RCP<const Map> rowMap = A->getRowMap();
141 Teuchos::ParameterList pL = this->GetParameterList();
142 if (pL.get<
bool>(
"fix nullspace")) {
143 this->GetOStream(
Runtime1) <<
"MueLu::Amesos2Smoother::Setup(): fixing nullspace" << std::endl;
145 rowMap = A->getRowMap();
146 size_t M = rowMap->getGlobalNumElements();
148 RCP<MultiVector> Nullspace = Factory::Get< RCP<MultiVector> >(currentLevel,
"Nullspace");
151 "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 for nullspace of dim > 1 has not been implemented yet.");
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);
168 NullspaceImp = Nullspace;
172 ArrayRCP<const SC> nullspaceRCP, nullspaceImpRCP;
173 RCP<CrsMatrixWrap> Acrs = rcp_dynamic_cast<CrsMatrixWrap>(A);
176 "MueLu::Amesos2Smoother::Setup Fixing nullspace for coarse matrix for Amesos2 when matrix is not a Crs matrix has not been implemented yet.");
178 ArrayRCP<const size_t> rowPointers;
179 ArrayRCP<const LO> colIndices;
180 ArrayRCP<const SC> values;
181 Acrs->getCrsMatrix()->getAllValues(rowPointers, colIndices, values);
183 ArrayRCP<size_t> newRowPointers_RCP;
184 ArrayRCP<LO> newColIndices_RCP;
185 ArrayRCP<SC> newValues_RCP;
187 size_t N = rowMap->getNodeNumElements();
188 newRowPointers_RCP.resize(N+1);
189 newColIndices_RCP.resize(N*M);
190 newValues_RCP.resize(N*M);
192 ArrayView<size_t> newRowPointers = newRowPointers_RCP();
193 ArrayView<LO> newColIndices = newColIndices_RCP();
194 ArrayView<SC> newValues = newValues_RCP();
196 SC normalization = Nullspace->getVector(0)->norm2();
197 normalization = Teuchos::ScalarTraits<Scalar>::one()/(normalization*normalization);
199 ArrayView<const SC> nullspace, nullspaceImp;
200 nullspaceRCP = Nullspace->getData(0);
201 nullspace = nullspaceRCP();
202 nullspaceImpRCP = NullspaceImp->getData(0);
203 nullspaceImp = nullspaceImpRCP();
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]);
213 newRowPointers[N] = N*M;
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]));
220 newValues[i*M+j] += v;
224 RCP<Matrix> newA = rcp(
new CrsMatrixWrap(rowMap, colMap, 0));
225 RCP<CrsMatrix> newAcrs = rcp_dynamic_cast<CrsMatrixWrap>(newA)->getCrsMatrix();
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());
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);
250 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
254 RCP<Tpetra_MultiVector> tX, tB;
255 if (!useTransformation_) {
260 size_t numVectors = X.getNumVectors();
261 size_t length = X.getLocalLength();
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);
268 for (
size_t i = 0; i < length; i++) {
269 X_data[i] = Xdata[i];
270 B_data[i] = Bdata[i];
282 prec_->setX(Teuchos::null);
283 prec_->setB(Teuchos::null);
285 if (useTransformation_) {
287 size_t length = X.getLocalLength();
289 ArrayRCP<SC> Xdata = X. getDataNonConst(0);
290 ArrayRCP<const SC> X_data = X_->getData(0);
292 for (
size_t i = 0; i < length; i++)
293 Xdata[i] = X_data[i];
297 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
298 RCP<MueLu::SmootherPrototype<Scalar, LocalOrdinal, GlobalOrdinal, Node> >
305 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
307 std::ostringstream out;
310 out << prec_->description();
314 out <<
"{type = " << type_ <<
"}";
319 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
324 out0 <<
"Prec. type: " << type_ << std::endl;
327 out0 <<
"Parameter list: " << std::endl;
328 Teuchos::OSTab tab2(out);
329 out << this->GetParameterList();
332 if ((verbLevel &
External) && prec_ != Teuchos::null) {
333 Teuchos::OSTab tab2(out);
334 out << *prec_ << std::endl;
337 if (verbLevel &
Debug)
340 <<
"RCP<prec_>: " << prec_ << std::endl;
343 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
346 return prec_->getStatus().getNnzLU();
353 #endif // HAVE_MUELU_TPETRA && HAVE_MUELU_AMESOS2 354 #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.
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 ¶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 Teuchos::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.
void declareConstructionOutcome(bool fail, std::string msg)
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)
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.
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.