46 #ifndef MUELU_AMESOS2SMOOTHER_DEF_HPP
47 #define MUELU_AMESOS2SMOOTHER_DEF_HPP
52 #if defined(HAVE_MUELU_AMESOS2)
56 #include <Amesos2_config.h>
57 #include <Amesos2.hpp>
61 #include "MueLu_Utilities.hpp"
66 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
70 Nullspace->getNumVectors(),
71 Nullspace->getMap()->getIndexBase(),
72 Nullspace->getMap()->getComm(),
80 Kokkos::View<Scalar**, Kokkos::LayoutLeft, Kokkos::HostSpace> Q(
"Q", Nullspace->getNumVectors(), Nullspace->getNumVectors());
83 auto dots = tempMV->getHostLocalView(Xpetra::Access::ReadOnly);
84 Kokkos::deep_copy(Q, dots);
92 lapack.
POTRF(
'L', Nullspace->getNumVectors(), Q.data(), LDQ, &info);
94 lapack.
TRTRI(
'L',
'N', Nullspace->getNumVectors(), Q.data(), LDQ, &info);
99 for (
size_t i = 0; i < Nullspace->getNumVectors(); i++) {
100 for (
size_t j = 0; j <= i; j++) {
101 Nullspace_->getVectorNonConst(i)->update(Q(i, j), *Nullspace->getVector(j), ONE);
106 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
116 auto dots = tempMV->getHostLocalView(Xpetra::Access::ReadOnly);
117 bool doProject =
true;
119 for (
size_t j = 0; j < Nullspace_->getNumVectors(); j++) {
125 for (
size_t j = 0; j < Nullspace_->getNumVectors(); j++) {
132 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
135 , useTransformation_(false) {
138 if (!
type_.empty()) {
141 std::transform(
type_.begin(), ++
type_.begin(),
type_.begin(), ::toupper);
143 if (
type_ ==
"Superlu_dist")
144 type_ =
"Superludist";
149 if (
type_ ==
"" || Amesos2::query(
type_) ==
false) {
150 std::string oldtype =
type_;
151 #if defined(HAVE_AMESOS2_SUPERLU)
153 #elif defined(HAVE_AMESOS2_KLU2)
155 #elif defined(HAVE_AMESOS2_SUPERLUDIST)
156 type_ =
"Superludist";
157 #elif defined(HAVE_AMESOS2_BASKER)
160 this->
declareConstructionOutcome(
true, std::string(
"Amesos2 has been compiled without SuperLU_DIST, SuperLU, Klu, or Basker. By default, MueLu tries") +
161 "to use one of these libraries. Amesos2 must be compiled with one of these solvers, " +
162 "or a valid Amesos2 solver has to be specified explicitly.");
166 this->
GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother: \"" << oldtype <<
"\" is not available. Using \"" <<
type_ <<
"\" instead" << std::endl;
173 "Amesos2 has been compiled without the support of this solver, or the solver name is misspelled.");
176 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
179 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
184 validParamList->
set<
bool>(
"fix nullspace",
false,
"Remove zero eigenvalue by adding rank one correction.");
187 validParamList->
set<
ParameterList>(
"Amesos2", norecurse,
"Parameters that are passed to Amesos2");
188 return validParamList;
191 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
195 this->Input(currentLevel,
"A");
196 if (pL.
get<
bool>(
"fix nullspace"))
197 this->Input(currentLevel,
"Nullspace");
200 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
205 this->GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother::Setup(): Setup() has already been called" << std::endl;
207 RCP<Matrix> A = Factory::Get<RCP<Matrix> >(currentLevel,
"A");
214 if (pL.
get<
bool>(
"fix nullspace")) {
215 this->GetOStream(
Runtime1) <<
"MueLu::Amesos2Smoother::Setup(): fixing nullspace" << std::endl;
217 rowMap = A->getRowMap();
218 size_t gblNumCols = rowMap->getGlobalNumElements();
220 RCP<MultiVector> NullspaceOrig = Factory::Get<RCP<MultiVector> >(currentLevel,
"Nullspace");
228 if (rowMap->getComm()->getSize() > 1) {
229 this->GetOStream(
Warnings0) <<
"MueLu::Amesos2Smoother::Setup(): Applying nullspace fix on distributed matrix. Try rebalancing to single rank!" << std::endl;
231 elements_RCP.
resize(gblNumCols);
233 for (
size_t k = 0; k < gblNumCols; k++)
234 elements[k] = Teuchos::as<GO>(k);
236 importer = ImportFactory::Build(rowMap, colMap);
237 ghostedNullspace = MultiVectorFactory::Build(colMap, Nullspace->getNumVectors());
238 ghostedNullspace->doImport(*Nullspace, *importer,
Xpetra::INSERT);
240 ghostedNullspace = Nullspace;
244 using ATS = Kokkos::ArithTraits<SC>;
245 using impl_Scalar =
typename ATS::val_type;
246 using impl_ATS = Kokkos::ArithTraits<impl_Scalar>;
247 using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
249 typedef typename Matrix::local_matrix_type KCRS;
250 typedef typename KCRS::StaticCrsGraphType graph_t;
251 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
252 typedef typename graph_t::entries_type::non_const_type lno_nnz_view_t;
253 typedef typename KCRS::values_type::non_const_type scalar_view_t;
255 const impl_Scalar impl_SC_ZERO = impl_ATS::zero();
257 size_t lclNumRows = rowMap->getLocalNumElements();
258 LocalOrdinal lclNumCols = Teuchos::as<LocalOrdinal>(gblNumCols);
259 lno_view_t newRowPointers(
"newRowPointers", lclNumRows + 1);
260 lno_nnz_view_t newColIndices(
"newColIndices", lclNumRows * gblNumCols);
261 scalar_view_t newValues(
"newValues", lclNumRows * gblNumCols);
265 RCP<Vector> diag = VectorFactory::Build(A->getRowMap());
266 A->getLocalDiagCopy(*diag);
267 shift = diag->normInf();
272 auto lclNullspace = Nullspace->getDeviceLocalView(Xpetra::Access::ReadOnly);
273 auto lclGhostedNullspace = ghostedNullspace->getDeviceLocalView(Xpetra::Access::ReadOnly);
274 Kokkos::parallel_for(
275 "MueLu:Amesos2Smoother::fixNullspace_1", range_type(0, lclNumRows + 1),
276 KOKKOS_LAMBDA(
const size_t i) {
277 if (i < lclNumRows) {
278 newRowPointers(i) = i * gblNumCols;
280 newColIndices(i * gblNumCols + j) = j;
281 newValues(i * gblNumCols + j) = impl_SC_ZERO;
282 for (
size_t I = 0;
I < lclNullspace.extent(1);
I++)
283 for (
size_t J = 0; J < lclGhostedNullspace.extent(1); J++)
284 newValues(i * gblNumCols + j) += shift * lclNullspace(i,
I) * impl_ATS::conjugate(lclGhostedNullspace(j, J));
287 newRowPointers(lclNumRows) = lclNumRows * gblNumCols;
293 auto lclA = A->getLocalMatrixDevice();
294 auto lclColMapA = A->getColMap()->getLocalMap();
295 auto lclColMapANew = colMap->getLocalMap();
296 Kokkos::parallel_for(
297 "MueLu:Amesos2Smoother::fixNullspace_2", range_type(0, lclNumRows),
298 KOKKOS_LAMBDA(
const size_t i) {
299 for (
size_t jj = lclA.graph.row_map(i); jj < lclA.graph.row_map(i + 1); jj++) {
300 LO j = lclColMapANew.getLocalElement(lclColMapA.getGlobalElement(lclA.graph.entries(jj)));
301 impl_Scalar v = lclA.values(jj);
302 newValues(i * gblNumCols + j) += v;
306 auto lclA = A->getLocalMatrixHost();
307 for (
size_t i = 0; i < lclNumRows; i++) {
308 for (
size_t jj = lclA.graph.row_map(i); jj < lclA.graph.row_map(i + 1); jj++) {
309 LO j = colMap->getLocalElement(A->getColMap()->getGlobalElement(lclA.graph.entries(jj)));
310 SC v = lclA.values(jj);
311 newValues(i * gblNumCols + j) += v;
317 RCP<CrsMatrix> newAcrs = rcp_dynamic_cast<CrsMatrixWrap>(newA)->getCrsMatrix();
318 newAcrs->setAllValues(newRowPointers, newColIndices, newValues);
319 newAcrs->expertStaticFillComplete(A->getDomainMap(), A->getRangeMap(),
320 importer, A->getCrsGraph()->getExporter());
323 rowMap = factorA->getRowMap();
330 prec_ = Amesos2::create<Tpetra_CrsMatrix, Tpetra_MultiVector>(type_, tA);
333 amesos2_params->
setName(
"Amesos2");
334 if ((rowMap->getGlobalNumElements() != as<size_t>((rowMap->getMaxAllGlobalIndex() - rowMap->getMinAllGlobalIndex()) + 1)) ||
335 (!rowMap->isContiguous() && (rowMap->getComm()->getSize() == 1))) {
336 if ((type_ !=
"Cusolver") && !(amesos2_params->
sublist(prec_->name()).
template isType<bool>(
"IsContiguous")))
337 amesos2_params->
sublist(prec_->name()).set(
"IsContiguous",
false,
"Are GIDs Contiguous");
341 prec_->numericFactorization();
346 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
351 if (!useTransformation_) {
356 size_t numVectors = X.getNumVectors();
357 size_t length = X.getLocalLength();
360 "MueLu::Amesos2Smoother::Apply: Fixing coarse matrix for Amesos2 for multivectors has not been implemented yet.");
362 ArrayRCP<SC> X_data = X_->getDataNonConst(0), B_data = B_->getDataNonConst(0);
364 for (
size_t i = 0; i < length; i++) {
365 X_data[i] = Xdata[i];
366 B_data[i] = Bdata[i];
378 prec_->setX(Teuchos::null);
379 prec_->setB(Teuchos::null);
381 if (useTransformation_) {
383 size_t length = X.getLocalLength();
388 for (
size_t i = 0; i < length; i++)
389 Xdata[i] = X_data[i];
394 if (pL.
get<
bool>(
"fix nullspace")) {
395 projection_->projectOut(X);
400 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
407 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
409 std::ostringstream out;
412 out << prec_->description();
416 out <<
"{type = " << type_ <<
"}";
421 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
426 out0 <<
"Prec. type: " << type_ << std::endl;
429 out0 <<
"Parameter list: " << std::endl;
431 out << this->GetParameterList();
434 if ((verbLevel &
External) && prec_ != Teuchos::null) {
436 out << *prec_ << std::endl;
439 if (verbLevel &
Debug)
442 <<
"RCP<prec_>: " << prec_ << std::endl;
445 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
447 if (!prec_.is_null())
448 return prec_->getStatus().getNnzLU();
454 #endif // HAVE_MUELU_AMESOS2
455 #endif // MUELU_AMESOS2SMOOTHER_DEF_HPP
Important warning messages (one line)
MueLu::DefaultLocalOrdinal LocalOrdinal
Teuchos::FancyOStream & GetOStream(MsgType type, int thisProcRankOnly=0) const
Get an output stream for outputting the input message type.
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static Teuchos::RCP< MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Build(const Teuchos::RCP< const Map< LocalOrdinal, GlobalOrdinal, Node >> &map, size_t NumVectors, bool zeroOut=true)
ParameterList & disableRecursiveValidation()
static magnitudeType eps()
Print external lib objects.
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
std::string type_
amesos2-specific key phrase that denote smoother type
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.
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Print additional debugging information.
std::string tolower(const std::string &str)
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.
std::string description() const
Return a simple one-line description of this object.
void resize(const size_type n, const T &val=T())
RCP< SmootherPrototype > Copy() const
virtual void SetParameterList(const Teuchos::ParameterList ¶mList)
Set parameters from a parameter list and return with default values.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
void declareConstructionOutcome(bool fail, std::string msg)
Class that holds all level-specific information.
static Teuchos::RCP< Map< LocalOrdinal, GlobalOrdinal, Node > > Build(UnderlyingLib lib, global_size_t numGlobalElements, GlobalOrdinal indexBase, const Teuchos::RCP< const Teuchos::Comm< int >> &comm, LocalGlobal lg=Xpetra::GloballyDistributed)
Class that encapsulates Amesos2 direct solvers.
void POTRF(const char &UPLO, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV2(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &vec)
bool IsSetup() const
Get the state of a smoother prototype.
void DeclareInput(Level ¤tLevel) const
Input.
#define MUELU_DESCRIBE
Helper macro for implementing Describable::describe() for BaseClass objects.
ParameterList & setParameters(const ParameterList &source)
static magnitudeType magnitude(T a)
virtual ~Amesos2Smoother()
Destructor.
void Setup(Level ¤tLevel)
Set up the direct solver. This creates the underlying Amesos2 solver object according to the paramete...
ParameterList & sublist(const std::string &name, bool mustAlreadyExist=false, const std::string &docString="")
ParameterList & setName(const std::string &name)
Print class parameters (more parameters, more verbose)
virtual Teuchos::RCP< Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > getVectorNonConst(size_t j)=0
Exception throws to report errors in the internal logical of the program.
#define TEUCHOS_ASSERT(assertion_test)
Description of what is happening (more verbose)
virtual size_t getNumVectors() const =0
void TRTRI(const char &UPLO, const char &DIAG, const OrdinalType &n, ScalarType *A, const OrdinalType &lda, OrdinalType *info) const
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
void print(Teuchos::FancyOStream &out, const VerbLevel verbLevel=Default) const
Print the object with some verbosity level to an FancyOStream object.
static RCP< Tpetra::CrsMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraCrs(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node >> Op)
virtual std::string description() const
Return a simple one-line description of this object.
Projection(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node >> &Nullspace)
std::string toString(const T &t)
void projectOut(Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X)