46 #ifndef MUELU_REFMAXWELL_DEF_HPP 47 #define MUELU_REFMAXWELL_DEF_HPP 53 #include "Xpetra_Map.hpp" 54 #include "Xpetra_MatrixMatrix.hpp" 55 #include "Xpetra_TripleMatrixMultiply.hpp" 56 #include "Xpetra_CrsMatrixUtils.hpp" 57 #include "Xpetra_MatrixUtils.hpp" 61 #include "MueLu_AmalgamationFactory.hpp" 62 #include "MueLu_RAPFactory.hpp" 63 #include "MueLu_ThresholdAFilterFactory.hpp" 64 #include "MueLu_TransPFactory.hpp" 65 #include "MueLu_SmootherFactory.hpp" 67 #include "MueLu_CoalesceDropFactory.hpp" 68 #include "MueLu_CoarseMapFactory.hpp" 69 #include "MueLu_CoordinatesTransferFactory.hpp" 70 #include "MueLu_UncoupledAggregationFactory.hpp" 71 #include "MueLu_TentativePFactory.hpp" 72 #include "MueLu_SaPFactory.hpp" 73 #include "MueLu_AggregationExportFactory.hpp" 74 #include "MueLu_Utilities.hpp" 76 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 77 #include "MueLu_AmalgamationFactory_kokkos.hpp" 78 #include "MueLu_CoalesceDropFactory_kokkos.hpp" 79 #include "MueLu_CoarseMapFactory_kokkos.hpp" 80 #include "MueLu_CoordinatesTransferFactory_kokkos.hpp" 81 #include "MueLu_UncoupledAggregationFactory_kokkos.hpp" 82 #include "MueLu_TentativePFactory_kokkos.hpp" 83 #include "MueLu_SaPFactory_kokkos.hpp" 84 #include "MueLu_Utilities_kokkos.hpp" 85 #include <Kokkos_Core.hpp> 86 #include <KokkosSparse_CrsMatrix.hpp> 89 #include "MueLu_ZoltanInterface.hpp" 90 #include "MueLu_Zoltan2Interface.hpp" 91 #include "MueLu_RepartitionHeuristicFactory.hpp" 92 #include "MueLu_RepartitionFactory.hpp" 93 #include "MueLu_RebalanceAcFactory.hpp" 94 #include "MueLu_RebalanceTransferFactory.hpp" 101 #ifdef HAVE_MUELU_CUDA 102 #include "cuda_profiler_api.h" 108 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
110 Teuchos::ArrayRCP<bool> nonzeros) {
111 TEUCHOS_ASSERT(vals.size() == nonzeros.size());
112 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType magnitudeType;
113 const magnitudeType eps = 2.0*Teuchos::ScalarTraits<magnitudeType>::eps();
114 for(
size_t i=0; i<static_cast<size_t>(vals.size()); i++) {
115 nonzeros[i] = (Teuchos::ScalarTraits<Scalar>::magnitude(vals[i]) > eps);
120 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
122 const Teuchos::ArrayRCP<bool>& dirichletRows,
123 Teuchos::ArrayRCP<bool> dirichletCols,
124 Teuchos::ArrayRCP<bool> dirichletDomain) {
125 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
126 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
127 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
128 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
129 TEUCHOS_ASSERT(static_cast<size_t>(dirichletRows.size()) == rowMap->getNodeNumElements());
130 TEUCHOS_ASSERT(static_cast<size_t>(dirichletCols.size()) == colMap->getNodeNumElements());
131 TEUCHOS_ASSERT(static_cast<size_t>(dirichletDomain.size()) == domMap->getNodeNumElements());
132 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap, 1,
true);
134 for(
size_t i=0; i<(size_t) dirichletRows.size(); i++) {
135 if (dirichletRows[i]) {
136 ArrayView<const LocalOrdinal> indices;
137 ArrayView<const Scalar> values;
138 A.getLocalRowView(i,indices,values);
139 for(
size_t j=0; j<static_cast<size_t>(indices.size()); j++)
140 myColsToZero->replaceLocalValue(indices[j],0,one);
144 RCP<Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
145 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
146 if (!importer.is_null()) {
147 globalColsToZero = Xpetra::MultiVectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap, 1,
true);
149 globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
151 myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
154 globalColsToZero = myColsToZero;
156 FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(globalColsToZero->getData(0),dirichletDomain);
157 FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(myColsToZero->getData(0),dirichletCols);
161 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 163 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
164 void FindNonZeros(
const typename Xpetra::MultiVector<Scalar,LocalOrdinal,GlobalOrdinal,Node>::dual_view_type::t_dev_um vals,
166 using ATS = Kokkos::ArithTraits<Scalar>;
167 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
169 TEUCHOS_ASSERT(vals.extent(0) == nonzeros.
extent(0));
170 const typename ATS::magnitudeType eps = 2.0*impl_ATS::eps();
173 KOKKOS_LAMBDA (
const size_t i) {
174 nonzeros(i) = (impl_ATS::magnitude(vals(i,0)) > eps);
178 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
184 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
185 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > domMap = A.getDomainMap();
186 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > rowMap = A.getRowMap();
187 RCP<const Xpetra::Map<LocalOrdinal,GlobalOrdinal,Node> > colMap = A.getColMap();
188 TEUCHOS_ASSERT(dirichletRows.
extent(0) == rowMap->getNodeNumElements());
189 TEUCHOS_ASSERT(dirichletCols.
extent(0) == colMap->getNodeNumElements());
190 TEUCHOS_ASSERT(dirichletDomain.
extent(0) == domMap->getNodeNumElements());
191 RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > myColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(colMap,
true);
193 auto myColsToZeroView = myColsToZero->template getLocalView<typename Node::device_type>();
194 auto localMatrix = A.getLocalMatrix();
195 Kokkos::parallel_for(
"MueLu:RefMaxwell::DetectDirichletCols", range_type(0,rowMap->getNodeNumElements()),
197 if (dirichletRows(row)) {
198 auto rowView = localMatrix.row(row);
199 auto length = rowView.length;
201 for (decltype(length) colID = 0; colID < length; colID++)
202 myColsToZeroView(rowView.colidx(colID),0) = one;
206 RCP<Xpetra::Vector<Scalar,LocalOrdinal,GlobalOrdinal,Node> > globalColsToZero;
207 RCP<const Xpetra::Import<LocalOrdinal,GlobalOrdinal,Node> > importer = A.getCrsGraph()->getImporter();
208 if (!importer.is_null()) {
209 globalColsToZero = Xpetra::VectorFactory<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Build(domMap,
true);
211 globalColsToZero->doExport(*myColsToZero,*importer,Xpetra::ADD);
213 myColsToZero->doImport(*globalColsToZero,*importer,Xpetra::INSERT);
216 globalColsToZero = myColsToZero;
217 FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(globalColsToZero->template getLocalView<typename Node::device_type>(),dirichletDomain);
218 FindNonZeros<Scalar,LocalOrdinal,GlobalOrdinal,Node>(myColsToZero->template getLocalView<typename Node::device_type>(),dirichletCols);
223 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
225 return SM_Matrix_->getDomainMap();
229 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
231 return SM_Matrix_->getRangeMap();
235 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
238 if (list.isType<std::string>(
"parameterlist: syntax") && list.get<std::string>(
"parameterlist: syntax") ==
"ml") {
240 if(list.isSublist(
"refmaxwell: 11list") && list.sublist(
"refmaxwell: 11list").isSublist(
"edge matrix free: coarse"))
242 if(list.isSublist(
"refmaxwell: 22list"))
247 parameterList_ = list;
248 std::string verbosityLevel = parameterList_.get<std::string>(
"verbosity", MasterList::getDefault<std::string>(
"verbosity"));
250 std::string outputFilename = parameterList_.get<std::string>(
"output filename", MasterList::getDefault<std::string>(
"output filename"));
251 if (outputFilename !=
"")
253 if (parameterList_.isType<Teuchos::RCP<Teuchos::FancyOStream> >(
"output stream"))
256 if (parameterList_.get(
"print initial parameters",MasterList::getDefault<bool>(
"print initial parameters")))
257 GetOStream(static_cast<MsgType>(
Runtime1), 0) << parameterList_ << std::endl;
258 disable_addon_ = list.get(
"refmaxwell: disable addon", MasterList::getDefault<bool>(
"refmaxwell: disable addon"));
259 mode_ = list.get(
"refmaxwell: mode", MasterList::getDefault<std::string>(
"refmaxwell: mode"));
260 use_as_preconditioner_ = list.get(
"refmaxwell: use as preconditioner", MasterList::getDefault<bool>(
"refmaxwell: use as preconditioner"));
261 dump_matrices_ = list.get(
"refmaxwell: dump matrices", MasterList::getDefault<bool>(
"refmaxwell: dump matrices"));
262 implicitTranspose_ = list.get(
"transpose: use implicit", MasterList::getDefault<bool>(
"transpose: use implicit"));
263 fuseProlongationAndUpdate_ = list.get(
"fuse prolongation and update", MasterList::getDefault<bool>(
"fuse prolongation and update"));
264 syncTimers_ = list.get(
"sync timers",
false);
265 numItersH_ = list.get(
"refmaxwell: num iters H", 1);
266 numIters22_ = list.get(
"refmaxwell: num iters 22", 1);
268 if(list.isSublist(
"refmaxwell: 11list"))
269 precList11_ = list.sublist(
"refmaxwell: 11list");
270 if(!precList11_.isType<std::string>(
"smoother: type") && !precList11_.isType<std::string>(
"smoother: pre type") && !precList11_.isType<std::string>(
"smoother: post type")) {
271 precList11_.set(
"smoother: type",
"CHEBYSHEV");
272 precList11_.sublist(
"smoother: params").set(
"chebyshev: degree",2);
273 precList11_.sublist(
"smoother: params").set(
"chebyshev: ratio eigenvalue",5.4);
274 precList11_.sublist(
"smoother: params").set(
"chebyshev: eigenvalue max iterations",30);
277 if(list.isSublist(
"refmaxwell: 22list"))
278 precList22_ = list.sublist(
"refmaxwell: 22list");
279 if(!precList22_.isType<std::string>(
"smoother: type") && !precList22_.isType<std::string>(
"smoother: pre type") && !precList22_.isType<std::string>(
"smoother: post type")) {
280 precList22_.set(
"smoother: type",
"CHEBYSHEV");
281 precList22_.sublist(
"smoother: params").set(
"chebyshev: degree",2);
282 precList22_.sublist(
"smoother: params").set(
"chebyshev: ratio eigenvalue",7.0);
283 precList22_.sublist(
"smoother: params").set(
"chebyshev: eigenvalue max iterations",30);
286 if(!list.isType<std::string>(
"smoother: type") && !list.isType<std::string>(
"smoother: pre type") && !list.isType<std::string>(
"smoother: post type")) {
287 list.set(
"smoother: type",
"CHEBYSHEV");
288 list.sublist(
"smoother: params").set(
"chebyshev: degree",2);
289 list.sublist(
"smoother: params").set(
"chebyshev: ratio eigenvalue",20.0);
290 list.sublist(
"smoother: params").set(
"chebyshev: eigenvalue max iterations",30);
293 if(list.isSublist(
"smoother: params")) {
294 smootherList_ = list.sublist(
"smoother: params");
297 #if !defined(HAVE_MUELU_KOKKOS_REFACTOR) 300 # ifdef HAVE_MUELU_SERIAL 301 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosSerialWrapperNode).name())
304 # ifdef HAVE_MUELU_OPENMP 305 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosOpenMPWrapperNode).name())
308 # ifdef HAVE_MUELU_CUDA 309 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosCudaWrapperNode).name())
312 useKokkos_ = list.get(
"use kokkos refactor",useKokkos_);
317 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
320 #ifdef HAVE_MUELU_CUDA 321 if (parameterList_.get<
bool>(
"refmaxwell: cuda profile setup",
false)) cudaProfilerStart();
324 std::string timerLabel;
326 timerLabel =
"MueLu RefMaxwell: compute (reuse)";
328 timerLabel =
"MueLu RefMaxwell: compute";
329 RCP<Teuchos::TimeMonitor> tmCompute = getTimer(timerLabel);
334 bool defaultFilter =
false;
339 if (parameterList_.get<
bool>(
"refmaxwell: filter D0",
true) && D0_Matrix_->getNodeMaxNumRowEntries()>2) {
343 fineLevel.
Set(
"A",D0_Matrix_);
344 fineLevel.
setlib(D0_Matrix_->getDomainMap()->lib());
347 fineLevel.
Request(
"A",ThreshFact.get());
348 ThreshFact->Build(fineLevel);
349 D0_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
352 defaultFilter =
true;
355 if (parameterList_.get<
bool>(
"refmaxwell: filter SM", defaultFilter)) {
356 RCP<Vector> diag = VectorFactory::Build(SM_Matrix_->getRowMap());
358 SM_Matrix_->getLocalDiagCopy(*diag);
364 fineLevel.
Set(
"A",SM_Matrix_);
365 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
367 fineLevel.
Request(
"A",ThreshFact.get());
368 ThreshFact->Build(fineLevel);
369 SM_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
372 if (parameterList_.get<
bool>(
"refmaxwell: filter M1", defaultFilter)) {
373 RCP<Vector> diag = VectorFactory::Build(M1_Matrix_->getRowMap());
375 M1_Matrix_->getLocalDiagCopy(*diag);
381 fineLevel.
Set(
"A",M1_Matrix_);
382 fineLevel.
setlib(M1_Matrix_->getDomainMap()->lib());
384 fineLevel.
Request(
"A",ThreshFact.get());
385 ThreshFact->Build(fineLevel);
386 M1_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
389 if (parameterList_.get<
bool>(
"refmaxwell: filter Ms", defaultFilter)) {
390 RCP<Vector> diag = VectorFactory::Build(Ms_Matrix_->getRowMap());
392 Ms_Matrix_->getLocalDiagCopy(*diag);
398 fineLevel.
Set(
"A",Ms_Matrix_);
399 fineLevel.
setlib(Ms_Matrix_->getDomainMap()->lib());
401 fineLevel.
Request(
"A",ThreshFact.get());
402 ThreshFact->Build(fineLevel);
403 Ms_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
407 RCP<ParameterList> params = rcp(
new ParameterList());;
408 params->set(
"printLoadBalancingInfo",
true);
409 params->set(
"printCommInfo",
true);
421 int BCedgesLocal = 0;
422 int BCnodesLocal = 0;
423 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 427 double rowsumTol = parameterList_.get(
"refmaxwell: row sum drop tol",-1.0);
428 if (rowsumTol > 0.) {
429 typedef Teuchos::ScalarTraits<Scalar> STS;
430 RCP<const Map> rowmap = SM_Matrix_->getRowMap();
431 for (LO row = 0; row < Teuchos::as<LO>(SM_Matrix_->getRowMap()->getNodeNumElements()); ++row) {
432 size_t nnz = SM_Matrix_->getNumEntriesInLocalRow(row);
433 ArrayView<const LO> indices;
434 ArrayView<const SC> vals;
435 SM_Matrix_->getLocalRowView(row, indices, vals);
437 SC rowsum = STS::zero();
438 SC diagval = STS::zero();
439 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
440 LO col = indices[colID];
442 diagval = vals[colID];
443 rowsum += vals[colID];
445 if (STS::real(rowsum) > STS::magnitude(diagval) * rowsumTol)
446 BCrowsKokkos_(row) =
true;
452 DetectDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(*D0_Matrix_,BCrowsKokkos_,BCcolsKokkos_,BCdomainKokkos_);
454 dump(BCrowsKokkos_,
"BCrows.m");
455 dump(BCcolsKokkos_,
"BCcols.m");
456 dump(BCdomainKokkos_,
"BCdomain.m");
458 for (
size_t i = 0; i<BCrowsKokkos_.size(); i++)
459 if (BCrowsKokkos_(i))
461 for (
size_t i = 0; i<BCdomainKokkos_.size(); i++)
462 if (BCdomainKokkos_(i))
465 #endif // HAVE_MUELU_KOKKOS_REFACTOR 469 double rowsumTol = parameterList_.get(
"refmaxwell: row sum drop tol",-1.0);
470 if (rowsumTol > 0.) {
471 typedef Teuchos::ScalarTraits<Scalar> STS;
472 RCP<const Map> rowmap = SM_Matrix_->getRowMap();
473 for (LO row = 0; row < Teuchos::as<LO>(SM_Matrix_->getRowMap()->getNodeNumElements()); ++row) {
474 size_t nnz = SM_Matrix_->getNumEntriesInLocalRow(row);
475 ArrayView<const LO> indices;
476 ArrayView<const SC> vals;
477 SM_Matrix_->getLocalRowView(row, indices, vals);
479 SC rowsum = STS::zero();
480 SC diagval = STS::zero();
481 for (LO colID = 0; colID < Teuchos::as<LO>(nnz); colID++) {
482 LO col = indices[colID];
484 diagval = vals[colID];
485 rowsum += vals[colID];
487 if (STS::real(rowsum) > STS::magnitude(diagval) * rowsumTol)
492 BCcols_.resize(D0_Matrix_->getColMap()->getNodeNumElements());
493 BCdomain_.resize(D0_Matrix_->getDomainMap()->getNodeNumElements());
494 DetectDirichletCols<Scalar,LocalOrdinal,GlobalOrdinal,Node>(*D0_Matrix_,BCrows_,BCcols_,BCdomain_);
496 dump(BCrows_,
"BCrows.m");
497 dump(BCcols_,
"BCcols.m");
498 dump(BCdomain_,
"BCdomain.m");
500 for (
auto it = BCrows_.begin(); it != BCrows_.end(); ++it)
503 for (
auto it = BCdomain_.begin(); it != BCdomain_.end(); ++it)
509 MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCedgesLocal, BCedges_);
510 MueLu_sumAll(SM_Matrix_->getRowMap()->getComm(), BCnodesLocal, BCnodes_);
512 BCedges_ = BCedgesLocal;
513 BCnodes_ = BCnodesLocal;
516 GetOStream(
Statistics2) <<
"MueLu::RefMaxwell::compute(): Detected " << BCedges_ <<
" BC rows and " << BCnodes_ <<
" BC columns." << std::endl;
519 TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::as<Xpetra::global_size_t>(BCedges_) >= D0_Matrix_->getRangeMap()->getGlobalNumElements(),
Exceptions::RuntimeError,
520 "All edges are detected as boundary edges!");
527 if(Nullspace_ != null) {
529 TEUCHOS_ASSERT(Nullspace_->getMap()->isCompatible(*(SM_Matrix_->getRowMap())));
531 else if(Nullspace_ == null && Coords_ != null) {
533 Array<coordinateType> norms(Coords_->getNumVectors());
534 Coords_->norm2(norms);
535 for (
size_t i=0;i<Coords_->getNumVectors();i++)
537 Nullspace_ = MultiVectorFactory::Build(SM_Matrix_->getRowMap(),Coords_->getNumVectors());
540 Array<Scalar> normsSC(Coords_->getNumVectors());
541 for (
size_t i=0;i<Coords_->getNumVectors();i++)
542 normsSC[i] = (SC) norms[i];
543 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 544 RCP<MultiVector> CoordsSC;
546 CoordsSC = Utilities_kokkos::RealValuedToScalarMultiVector(Coords_);
552 D0_Matrix_->apply(*CoordsSC,*Nullspace_);
555 ArrayRCP<ArrayRCP<const Scalar> > localNullspace(Nullspace_->getNumVectors());
556 for (
size_t i = 0; i < Nullspace_->getNumVectors(); i++)
557 localNullspace[i] = Nullspace_->getData(i);
558 coordinateType localMinLen = Teuchos::ScalarTraits<coordinateType>::rmax();
559 coordinateType localMeanLen = Teuchos::ScalarTraits<coordinateType>::zero();
560 coordinateType localMaxLen = Teuchos::ScalarTraits<coordinateType>::zero();
561 for (
size_t j=0; j < Nullspace_->getMap()->getNodeNumElements(); j++) {
562 Scalar lenSC = Teuchos::ScalarTraits<Scalar>::zero();
563 for (
size_t i=0; i < Nullspace_->getNumVectors(); i++)
564 lenSC += localNullspace[i][j]*localNullspace[i][j];
565 coordinateType len = sqrt(Teuchos::ScalarTraits<Scalar>::real(lenSC));
566 localMinLen = std::min(localMinLen, len);
567 localMaxLen = std::max(localMaxLen, len);
572 RCP<const Teuchos::Comm<int> > comm = Nullspace_->getMap()->getComm();
577 minLen = localMinLen;
578 meanLen = localMeanLen;
579 maxLen = localMaxLen;
581 meanLen /= Nullspace_->getMap()->getGlobalNumElements();
582 GetOStream(
Statistics0) <<
"Edge length (min/mean/max): " << minLen <<
" / " << meanLen <<
" / " << maxLen << std::endl;
584 Nullspace_->scale(normsSC());
587 GetOStream(
Errors) <<
"MueLu::RefMaxwell::compute(): either the nullspace or the nodal coordinates must be provided." << std::endl;
592 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 600 dump(*Nullspace_,
"nullspace.m");
608 Level fineLevel, coarseLevel;
614 fineLevel.
Set(
"A",Ms_Matrix_);
615 coarseLevel.
Set(
"P",D0_Matrix_);
616 coarseLevel.
setlib(Ms_Matrix_->getDomainMap()->lib());
617 fineLevel.
setlib(Ms_Matrix_->getDomainMap()->lib());
618 coarseLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
619 fineLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
621 RCP<RAPFactory> rapFact = rcp(
new RAPFactory());
622 ParameterList rapList = *(rapFact->GetValidParameterList());
623 rapList.set(
"transpose: use implicit",
true);
624 rapList.set(
"rap: fix zero diagonals", parameterList_.get<
bool>(
"rap: fix zero diagonals",
true));
625 rapList.set(
"rap: triple product", parameterList_.get<
bool>(
"rap: triple product",
false));
626 rapFact->SetParameterList(rapList);
629 coarseLevel.
Request(
"A", rapFact.get());
631 A_nodal_Matrix_ = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
634 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 640 dump(*A_nodal_Matrix_,
"A_nodal.m");
643 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): building special prolongator" << std::endl;
648 bool doRebalancing =
false;
649 doRebalancing = parameterList_.get<
bool>(
"refmaxwell: subsolves on subcommunicators", MasterList::getDefault<bool>(
"refmaxwell: subsolves on subcommunicators"));
650 int rebalanceStriding = parameterList_.get<
int>(
"refmaxwell: subsolves striding", -1);
651 int numProcsAH, numProcsA22;
658 int numProcs = SM_Matrix_->getDomainMap()->getComm()->getSize();
659 if (doRebalancing && numProcs > 1) {
668 ParameterList repartheurParams;
669 repartheurParams.set(
"repartition: start level", 0);
671 int defaultTargetRows = 10000;
672 repartheurParams.set(
"repartition: min rows per proc", precList11_.get<
int>(
"repartition: target rows per proc", defaultTargetRows));
673 repartheurParams.set(
"repartition: target rows per proc", precList11_.get<
int>(
"repartition: target rows per proc", defaultTargetRows));
674 repartheurParams.set(
"repartition: min rows per thread", precList11_.get<
int>(
"repartition: target rows per thread", defaultTargetRows));
675 repartheurParams.set(
"repartition: target rows per thread", precList11_.get<
int>(
"repartition: target rows per thread", defaultTargetRows));
676 repartheurParams.set(
"repartition: max imbalance", precList11_.get<
double>(
"repartition: max imbalance", 1.1));
677 repartheurFactory->SetParameterList(repartheurParams);
679 level.
Request(
"number of partitions", repartheurFactory.get());
680 repartheurFactory->Build(level);
681 numProcsAH = level.
Get<
int>(
"number of partitions", repartheurFactory.get());
682 numProcsAH = std::min(numProcsAH,numProcs);
690 level.
Set(
"Map",D0_Matrix_->getDomainMap());
693 ParameterList repartheurParams;
694 repartheurParams.set(
"repartition: start level", 0);
695 repartheurParams.set(
"repartition: use map",
true);
697 int defaultTargetRows = 10000;
698 repartheurParams.set(
"repartition: min rows per proc", precList22_.get<
int>(
"repartition: target rows per proc", defaultTargetRows));
699 repartheurParams.set(
"repartition: target rows per proc", precList22_.get<
int>(
"repartition: target rows per proc", defaultTargetRows));
700 repartheurParams.set(
"repartition: min rows per thread", precList22_.get<
int>(
"repartition: target rows per thread", defaultTargetRows));
701 repartheurParams.set(
"repartition: target rows per thread", precList22_.get<
int>(
"repartition: target rows per thread", defaultTargetRows));
703 repartheurFactory->SetParameterList(repartheurParams);
705 level.
Request(
"number of partitions", repartheurFactory.get());
706 repartheurFactory->Build(level);
707 numProcsA22 = level.
Get<
int>(
"number of partitions", repartheurFactory.get());
708 numProcsA22 = std::min(numProcsA22,numProcs);
711 if (rebalanceStriding >= 1) {
712 TEUCHOS_ASSERT(rebalanceStriding*numProcsAH<=numProcs);
713 TEUCHOS_ASSERT(rebalanceStriding*numProcsA22<=numProcs);
714 if (rebalanceStriding*(numProcsAH+numProcsA22)>numProcs) {
715 GetOStream(
Warnings0) <<
"RefMaxwell::compute(): Disabling striding = " << rebalanceStriding <<
", since AH needs " << numProcsAH
716 <<
" procs and A22 needs " << numProcsA22 <<
" procs."<< std::endl;
717 rebalanceStriding = -1;
722 doRebalancing =
false;
726 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu RefMaxwell: Rebalance AH");
728 Level fineLevel, coarseLevel;
734 coarseLevel.
Set(
"A",AH_);
735 coarseLevel.
Set(
"P",P11_);
736 coarseLevel.
Set(
"Coordinates",CoordsH_);
737 coarseLevel.
Set(
"number of partitions", numProcsAH);
738 coarseLevel.
Set(
"repartition: heuristic target rows per process", 1000);
740 coarseLevel.
setlib(AH_->getDomainMap()->lib());
741 fineLevel.
setlib(AH_->getDomainMap()->lib());
742 coarseLevel.setObjectLabel(
"RefMaxwell coarse (1,1)");
743 fineLevel.setObjectLabel(
"RefMaxwell coarse (1,1)");
745 std::string partName = precList11_.get<std::string>(
"repartition: partitioner",
"zoltan2");
746 RCP<Factory> partitioner;
747 if (partName ==
"zoltan") {
748 #ifdef HAVE_MUELU_ZOLTAN 755 }
else if (partName ==
"zoltan2") {
756 #ifdef HAVE_MUELU_ZOLTAN2 758 ParameterList partParams;
759 RCP<const ParameterList> partpartParams = rcp(
new ParameterList(precList11_.sublist(
"repartition: params",
false)));
760 partParams.set(
"ParameterList", partpartParams);
761 partitioner->SetParameterList(partParams);
769 ParameterList repartParams;
770 repartParams.set(
"repartition: print partition distribution", precList11_.get<
bool>(
"repartition: print partition distribution",
false));
771 repartParams.set(
"repartition: remap parts", precList11_.get<
bool>(
"repartition: remap parts",
true));
772 if (rebalanceStriding >= 1) {
773 bool acceptPart = (SM_Matrix_->getDomainMap()->getComm()->getRank() % rebalanceStriding) == 0;
774 if (SM_Matrix_->getDomainMap()->getComm()->getRank() >= numProcsAH*rebalanceStriding)
776 repartParams.set(
"repartition: remap accept partition", acceptPart);
778 repartFactory->SetParameterList(repartParams);
780 repartFactory->SetFactory(
"Partition", partitioner);
783 ParameterList newPparams;
784 newPparams.set(
"type",
"Interpolation");
785 newPparams.set(
"repartition: rebalance P and R", precList11_.get<
bool>(
"repartition: rebalance P and R",
false));
786 newPparams.set(
"repartition: use subcommunicators",
true);
787 newPparams.set(
"repartition: rebalance Nullspace",
false);
789 newP->SetParameterList(newPparams);
790 newP->SetFactory(
"Importer", repartFactory);
793 ParameterList rebAcParams;
794 rebAcParams.set(
"repartition: use subcommunicators",
true);
795 newA->SetParameterList(rebAcParams);
796 newA->SetFactory(
"Importer", repartFactory);
798 coarseLevel.
Request(
"P", newP.get());
799 coarseLevel.
Request(
"Importer", repartFactory.get());
800 coarseLevel.
Request(
"A", newA.get());
801 coarseLevel.
Request(
"Coordinates", newP.get());
802 repartFactory->Build(coarseLevel);
804 if (!precList11_.get<
bool>(
"repartition: rebalance P and R",
false))
805 ImporterH_ = coarseLevel.
Get< RCP<const Import> >(
"Importer", repartFactory.get());
806 P11_ = coarseLevel.
Get< RCP<Matrix> >(
"P", newP.get());
807 AH_ = coarseLevel.
Get< RCP<Matrix> >(
"A", newA.get());
808 CoordsH_ = coarseLevel.
Get< RCP<RealValuedMultiVector> >(
"Coordinates", newP.get());
811 ParameterList XpetraList;
812 XpetraList.set(
"Restrict Communicator",
true);
813 XpetraList.set(
"Timer Label",
"MueLu RefMaxwell::RebalanceAH");
814 RCP<const Map> targetMap = ImporterH_->getTargetMap();
815 AH_ = MatrixFactory::Build(AH_, *ImporterH_, *ImporterH_, targetMap, targetMap, rcp(&XpetraList,
false));
818 AH_->setObjectLabel(
"RefMaxwell coarse (1,1)");
822 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 826 if (useKokkos_ && precList11_.isParameter(
"aggregation: drop tol") && precList11_.get<
double>(
"aggregation: drop tol") != 0.0) {
827 GetOStream(
Warnings0) <<
"RefMaxwell::compute(): Setting \"aggregation: drop tol\". to 0.0, since CoalesceDropFactory_kokkos does not " 828 <<
"support BlockSize > 1 and drop tol != 0.0" << std::endl;
829 precList11_.set(
"aggregation: drop tol", 0.0);
832 dump(*P11_,
"P11.m");
834 if (!implicitTranspose_) {
835 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 837 R11_ = Utilities_kokkos::Transpose(*P11_);
843 dump(*R11_,
"R11.m");
847 if (!AH_.is_null()) {
849 dumpCoords(*CoordsH_,
"coordsH.m");
850 int oldRank = SetProcRankVerbose(AH_->getDomainMap()->getComm()->getRank());
852 RCP<ParameterList> params = rcp(
new ParameterList());;
853 params->set(
"printLoadBalancingInfo",
true);
854 params->set(
"printCommInfo",
true);
858 ParameterList& userParamList = precList11_.sublist(
"user data");
859 userParamList.set<RCP<RealValuedMultiVector> >(
"Coordinates", CoordsH_);
862 RCP<MueLu::Level> level0 = HierarchyH_->GetLevel(0);
863 level0->Set(
"A", AH_);
864 HierarchyH_->SetupRe();
866 SetProcRankVerbose(oldRank);
873 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): nuking BC nodes of D0" << std::endl;
875 D0_Matrix_->resumeFill();
877 if (D0_Matrix_->getRowMap()->lib() == Xpetra::UseEpetra)
878 replaceWith= Teuchos::ScalarTraits<SC>::eps();
880 replaceWith = Teuchos::ScalarTraits<SC>::zero();
881 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 890 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
896 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): building MG for (2,2)-block" << std::endl;
899 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu RefMaxwell: Build A22");
901 Level fineLevel, coarseLevel;
907 fineLevel.
Set(
"A",SM_Matrix_);
908 coarseLevel.
Set(
"P",D0_Matrix_);
909 coarseLevel.
Set(
"Coordinates",Coords_);
911 coarseLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
912 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
913 coarseLevel.setObjectLabel(
"RefMaxwell (2,2)");
914 fineLevel.setObjectLabel(
"RefMaxwell (2,2)");
916 RCP<RAPFactory> rapFact = rcp(
new RAPFactory());
917 ParameterList rapList = *(rapFact->GetValidParameterList());
918 rapList.set(
"transpose: use implicit",
true);
919 rapList.set(
"rap: fix zero diagonals", parameterList_.get<
bool>(
"rap: fix zero diagonals",
true));
920 rapList.set(
"rap: triple product", parameterList_.get<
bool>(
"rap: triple product",
false));
921 rapFact->SetParameterList(rapList);
923 if (!A22_AP_reuse_data_.is_null()) {
924 coarseLevel.
AddKeepFlag(
"AP reuse data", rapFact.get());
925 coarseLevel.
Set<Teuchos::RCP<Teuchos::ParameterList> >(
"AP reuse data", A22_AP_reuse_data_, rapFact.get());
927 if (!A22_RAP_reuse_data_.is_null()) {
928 coarseLevel.
AddKeepFlag(
"RAP reuse data", rapFact.get());
929 coarseLevel.
Set<Teuchos::RCP<Teuchos::ParameterList> >(
"RAP reuse data", A22_RAP_reuse_data_, rapFact.get());
937 coarseLevel.
Set(
"number of partitions", numProcsA22);
938 coarseLevel.
Set(
"repartition: heuristic target rows per process", 1000);
940 std::string partName = precList22_.get<std::string>(
"repartition: partitioner",
"zoltan2");
941 RCP<Factory> partitioner;
942 if (partName ==
"zoltan") {
943 #ifdef HAVE_MUELU_ZOLTAN 945 partitioner->SetFactory(
"A", rapFact);
951 }
else if (partName ==
"zoltan2") {
952 #ifdef HAVE_MUELU_ZOLTAN2 954 ParameterList partParams;
955 RCP<const ParameterList> partpartParams = rcp(
new ParameterList(precList22_.sublist(
"repartition: params",
false)));
956 partParams.set(
"ParameterList", partpartParams);
957 partitioner->SetParameterList(partParams);
958 partitioner->SetFactory(
"A", rapFact);
966 ParameterList repartParams;
967 repartParams.set(
"repartition: print partition distribution", precList22_.get<
bool>(
"repartition: print partition distribution",
false));
968 repartParams.set(
"repartition: remap parts", precList22_.get<
bool>(
"repartition: remap parts",
true));
969 if (rebalanceStriding >= 1) {
970 bool acceptPart = ((SM_Matrix_->getDomainMap()->getComm()->getSize()-1-SM_Matrix_->getDomainMap()->getComm()->getRank()) % rebalanceStriding) == 0;
971 if (SM_Matrix_->getDomainMap()->getComm()->getSize()-1-SM_Matrix_->getDomainMap()->getComm()->getRank() >= numProcsA22*rebalanceStriding)
974 TEUCHOS_ASSERT(AH_.is_null());
975 repartParams.set(
"repartition: remap accept partition", acceptPart);
977 repartParams.set(
"repartition: remap accept partition", AH_.is_null());
978 repartFactory->SetParameterList(repartParams);
979 repartFactory->SetFactory(
"A", rapFact);
981 repartFactory->SetFactory(
"Partition", partitioner);
984 ParameterList newPparams;
985 newPparams.set(
"type",
"Interpolation");
986 newPparams.set(
"repartition: rebalance P and R", precList22_.get<
bool>(
"repartition: rebalance P and R",
false));
987 newPparams.set(
"repartition: use subcommunicators",
true);
988 newPparams.set(
"repartition: rebalance Nullspace",
false);
990 newP->SetParameterList(newPparams);
991 newP->SetFactory(
"Importer", repartFactory);
994 ParameterList rebAcParams;
995 rebAcParams.set(
"repartition: use subcommunicators",
true);
996 newA->SetParameterList(rebAcParams);
997 newA->SetFactory(
"A", rapFact);
998 newA->SetFactory(
"Importer", repartFactory);
1000 coarseLevel.
Request(
"P", newP.get());
1001 coarseLevel.
Request(
"Importer", repartFactory.get());
1002 coarseLevel.
Request(
"A", newA.get());
1003 if (precList22_.isType<std::string>(
"reuse: type") && precList22_.get<std::string>(
"reuse: type") !=
"none") {
1004 if (!parameterList_.get<
bool>(
"rap: triple product",
false))
1005 coarseLevel.
Request(
"AP reuse data", rapFact.get());
1006 coarseLevel.
Request(
"RAP reuse data", rapFact.get());
1008 coarseLevel.
Request(
"Coordinates", newP.get());
1009 rapFact->Build(fineLevel,coarseLevel);
1010 repartFactory->Build(coarseLevel);
1012 if (!precList22_.get<
bool>(
"repartition: rebalance P and R",
false))
1013 Importer22_ = coarseLevel.
Get< RCP<const Import> >(
"Importer", repartFactory.get());
1014 D0_Matrix_ = coarseLevel.
Get< RCP<Matrix> >(
"P", newP.get());
1015 A22_ = coarseLevel.
Get< RCP<Matrix> >(
"A", newA.get());
1016 if (precList22_.isType<std::string>(
"reuse: type") && precList22_.get<std::string>(
"reuse: type") !=
"none") {
1017 if (!parameterList_.get<
bool>(
"rap: triple product",
false))
1018 A22_AP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data", rapFact.get());
1019 A22_RAP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"RAP reuse data", rapFact.get());
1021 Coords_ = coarseLevel.
Get< RCP<RealValuedMultiVector> >(
"Coordinates", newP.get());
1023 coarseLevel.
Request(
"A", rapFact.get());
1024 if (precList22_.isType<std::string>(
"reuse: type") && precList22_.get<std::string>(
"reuse: type") !=
"none") {
1025 coarseLevel.
Request(
"AP reuse data", rapFact.get());
1026 coarseLevel.
Request(
"RAP reuse data", rapFact.get());
1029 A22_ = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
1030 if (precList22_.isType<std::string>(
"reuse: type") && precList22_.get<std::string>(
"reuse: type") !=
"none") {
1031 if (!parameterList_.get<
bool>(
"rap: triple product",
false))
1032 A22_AP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data", rapFact.get());
1033 A22_RAP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"RAP reuse data", rapFact.get());
1036 ParameterList XpetraList;
1037 XpetraList.set(
"Restrict Communicator",
true);
1038 XpetraList.set(
"Timer Label",
"MueLu RefMaxwell::RebalanceA22");
1039 RCP<const Map> targetMap = Importer22_->getTargetMap();
1040 A22_ = MatrixFactory::Build(A22_, *Importer22_, *Importer22_, targetMap, targetMap, rcp(&XpetraList,
false));
1045 coarseLevel.
Request(
"A", rapFact.get());
1046 if (precList22_.isType<std::string>(
"reuse: type") && precList22_.get<std::string>(
"reuse: type") !=
"none") {
1047 coarseLevel.
Request(
"AP reuse data", rapFact.get());
1048 coarseLevel.
Request(
"RAP reuse data", rapFact.get());
1051 A22_ = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
1053 if (precList22_.isType<std::string>(
"reuse: type") && precList22_.get<std::string>(
"reuse: type") !=
"none") {
1054 if (!parameterList_.get<
bool>(
"rap: triple product",
false))
1055 A22_AP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data", rapFact.get());
1056 A22_RAP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"RAP reuse data", rapFact.get());
1061 if (!implicitTranspose_) {
1062 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1064 D0_T_Matrix_ = Utilities_kokkos::Transpose(*D0_Matrix_);
1073 if (!A22_.is_null()) {
1074 dump(*A22_,
"A22.m");
1075 A22_->setObjectLabel(
"RefMaxwell (2,2)");
1076 int oldRank = SetProcRankVerbose(A22_->getDomainMap()->getComm()->getRank());
1078 RCP<ParameterList> params = rcp(
new ParameterList());;
1079 params->set(
"printLoadBalancingInfo",
true);
1080 params->set(
"printCommInfo",
true);
1084 ParameterList& userParamList = precList22_.sublist(
"user data");
1085 userParamList.set<RCP<RealValuedMultiVector> >(
"Coordinates", Coords_);
1088 std::string coarseType =
"";
1089 if (precList22_.isParameter(
"coarse: type")) {
1090 coarseType = precList22_.get<std::string>(
"coarse: type");
1092 std::transform(coarseType.begin(), coarseType.end(), coarseType.begin(),
::tolower);
1093 std::transform(coarseType.begin(), ++coarseType.begin(), coarseType.begin(), ::toupper);
1095 if (BCedges_ == 0 &&
1096 (coarseType ==
"" ||
1097 coarseType ==
"Klu" ||
1098 coarseType ==
"Klu2") &&
1099 (!precList22_.isSublist(
"coarse: params") ||
1100 !precList22_.sublist(
"coarse: params").isParameter(
"fix nullspace")))
1101 precList22_.sublist(
"coarse: params").set(
"fix nullspace",
true);
1104 RCP<MueLu::Level> level0 = Hierarchy22_->GetLevel(0);
1105 level0->Set(
"A", A22_);
1106 Hierarchy22_->SetupRe();
1108 SetProcRankVerbose(oldRank);
1115 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): nuking BC edges of D0" << std::endl;
1117 D0_Matrix_->resumeFill();
1119 if (D0_Matrix_->getRowMap()->lib() == Xpetra::UseEpetra)
1120 replaceWith= Teuchos::ScalarTraits<SC>::eps();
1122 replaceWith = Teuchos::ScalarTraits<SC>::zero();
1123 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1132 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
1133 dump(*D0_Matrix_,
"D0_nuked.m");
1137 if (parameterList_.isType<std::string>(
"smoother: type") &&
1138 parameterList_.get<std::string>(
"smoother: type") ==
"hiptmair" &&
1139 SM_Matrix_->getDomainMap()->lib() == Xpetra::UseTpetra &&
1140 A22_->getDomainMap()->lib() == Xpetra::UseTpetra &&
1141 D0_Matrix_->getDomainMap()->lib() == Xpetra::UseTpetra) {
1142 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 1143 ParameterList hiptmairPreList, hiptmairPostList, smootherPreList, smootherPostList;
1145 if (smootherList_.isSublist(
"smoother: pre params"))
1146 smootherPreList = smootherList_.sublist(
"smoother: pre params");
1147 else if (smootherList_.isSublist(
"smoother: params"))
1148 smootherPreList = smootherList_.sublist(
"smoother: params");
1149 hiptmairPreList.set(
"hiptmair: smoother type 1",
1150 smootherPreList.get<std::string>(
"hiptmair: smoother type 1",
"CHEBYSHEV"));
1151 hiptmairPreList.set(
"hiptmair: smoother type 2",
1152 smootherPreList.get<std::string>(
"hiptmair: smoother type 2",
"CHEBYSHEV"));
1153 if(smootherPreList.isSublist(
"hiptmair: smoother list 1"))
1154 hiptmairPreList.set(
"hiptmair: smoother list 1", smootherPreList.sublist(
"hiptmair: smoother list 1"));
1155 if(smootherPreList.isSublist(
"hiptmair: smoother list 2"))
1156 hiptmairPreList.set(
"hiptmair: smoother list 2", smootherPreList.sublist(
"hiptmair: smoother list 2"));
1157 hiptmairPreList.set(
"hiptmair: pre or post",
1158 smootherPreList.get<std::string>(
"hiptmair: pre or post",
"pre"));
1159 hiptmairPreList.set(
"hiptmair: zero starting solution",
1160 smootherPreList.get<
bool>(
"hiptmair: zero starting solution",
true));
1162 if (smootherList_.isSublist(
"smoother: post params"))
1163 smootherPostList = smootherList_.sublist(
"smoother: post params");
1164 else if (smootherList_.isSublist(
"smoother: params"))
1165 smootherPostList = smootherList_.sublist(
"smoother: params");
1166 hiptmairPostList.set(
"hiptmair: smoother type 1",
1167 smootherPostList.get<std::string>(
"hiptmair: smoother type 1",
"CHEBYSHEV"));
1168 hiptmairPostList.set(
"hiptmair: smoother type 2",
1169 smootherPostList.get<std::string>(
"hiptmair: smoother type 2",
"CHEBYSHEV"));
1170 if(smootherPostList.isSublist(
"hiptmair: smoother list 1"))
1171 hiptmairPostList.set(
"hiptmair: smoother list 1", smootherPostList.sublist(
"hiptmair: smoother list 1"));
1172 if(smootherPostList.isSublist(
"hiptmair: smoother list 2"))
1173 hiptmairPostList.set(
"hiptmair: smoother list 2", smootherPostList.sublist(
"hiptmair: smoother list 2"));
1174 hiptmairPostList.set(
"hiptmair: pre or post",
1175 smootherPostList.get<std::string>(
"hiptmair: pre or post",
"post"));
1176 hiptmairPostList.set(
"hiptmair: zero starting solution",
1177 smootherPostList.get<
bool>(
"hiptmair: zero starting solution",
false));
1179 typedef Tpetra::RowMatrix<SC, LO, GO, NO> TROW;
1185 hiptmairPreSmoother_ -> setParameters(hiptmairPreList);
1186 hiptmairPreSmoother_ -> initialize();
1187 hiptmairPreSmoother_ -> compute();
1189 hiptmairPostSmoother_ -> setParameters(hiptmairPostList);
1190 hiptmairPostSmoother_ -> initialize();
1191 hiptmairPostSmoother_ -> compute();
1192 useHiptmairSmoothing_ =
true;
1194 throw(Xpetra::Exceptions::RuntimeError(
"MueLu must be compiled with Ifpack2 for Hiptmair smoothing."));
1195 #endif // defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 1197 if (parameterList_.isType<std::string>(
"smoother: pre type") && parameterList_.isType<std::string>(
"smoother: post type")) {
1198 std::string preSmootherType = parameterList_.get<std::string>(
"smoother: pre type");
1199 std::string postSmootherType = parameterList_.get<std::string>(
"smoother: post type");
1201 ParameterList preSmootherList, postSmootherList;
1202 if (parameterList_.isSublist(
"smoother: pre params"))
1203 preSmootherList = parameterList_.sublist(
"smoother: pre params");
1204 if (parameterList_.isSublist(
"smoother: post params"))
1205 postSmootherList = parameterList_.sublist(
"smoother: post params");
1208 RCP<MueLu::FactoryManagerBase> factoryHandler = rcp(
new FactoryManager());
1211 level.setObjectLabel(
"RefMaxwell (1,1)");
1212 level.
Set(
"A",SM_Matrix_);
1213 level.
setlib(SM_Matrix_->getDomainMap()->lib());
1215 RCP<SmootherPrototype> preSmootherPrototype = rcp(
new TrilinosSmoother(preSmootherType, preSmootherList));
1216 RCP<SmootherFactory> preSmootherFact = rcp(
new SmootherFactory(preSmootherPrototype));
1218 RCP<SmootherPrototype> postSmootherPrototype = rcp(
new TrilinosSmoother(postSmootherType, postSmootherList));
1219 RCP<SmootherFactory> postSmootherFact = rcp(
new SmootherFactory(postSmootherPrototype));
1221 level.
Request(
"PreSmoother",preSmootherFact.get());
1222 preSmootherFact->Build(level);
1223 PreSmoother_ = level.
Get<RCP<SmootherBase> >(
"PreSmoother",preSmootherFact.get());
1225 level.
Request(
"PostSmoother",postSmootherFact.get());
1226 postSmootherFact->Build(level);
1227 PostSmoother_ = level.
Get<RCP<SmootherBase> >(
"PostSmoother",postSmootherFact.get());
1229 std::string smootherType = parameterList_.get<std::string>(
"smoother: type",
"CHEBYSHEV");
1231 RCP<MueLu::FactoryManagerBase> factoryHandler = rcp(
new FactoryManager());
1232 level.SetFactoryManager(factoryHandler);
1233 level.SetLevelID(0);
1234 level.setObjectLabel(
"RefMaxwell (1,1)");
1235 level.Set(
"A",SM_Matrix_);
1236 level.setlib(SM_Matrix_->getDomainMap()->lib());
1237 RCP<SmootherPrototype> smootherPrototype = rcp(
new TrilinosSmoother(smootherType, smootherList_));
1238 RCP<SmootherFactory> SmootherFact = rcp(
new SmootherFactory(smootherPrototype));
1239 level.Request(
"PreSmoother",SmootherFact.get());
1240 SmootherFact->Build(level);
1241 PreSmoother_ = level.Get<RCP<SmootherBase> >(
"PreSmoother",SmootherFact.get());
1242 PostSmoother_ = PreSmoother_;
1244 useHiptmairSmoothing_ =
false;
1248 if (!ImporterH_.is_null()) {
1249 RCP<const Import> ImporterP11 = ImportFactory::Build(ImporterH_->getTargetMap(),P11_->getColMap());
1250 rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix()->replaceDomainMapAndImporter(ImporterH_->getTargetMap(), ImporterP11);
1253 if (!Importer22_.is_null()) {
1254 RCP<const Import> ImporterD0 = ImportFactory::Build(Importer22_->getTargetMap(),D0_Matrix_->getColMap());
1255 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->replaceDomainMapAndImporter(Importer22_->getTargetMap(), ImporterD0);
1258 if ((!D0_T_Matrix_.is_null()) &&
1259 (!R11_.is_null()) &&
1260 (!rcp_dynamic_cast<CrsMatrixWrap>(D0_T_Matrix_)->getCrsMatrix()->getCrsGraph()->getImporter().is_null()) &&
1261 (!rcp_dynamic_cast<CrsMatrixWrap>(R11_)->getCrsMatrix()->getCrsGraph()->getImporter().is_null()) &&
1262 (D0_T_Matrix_->getColMap()->lib() == Xpetra::UseTpetra) &&
1263 (R11_->getColMap()->lib() == Xpetra::UseTpetra))
1264 D0_T_R11_colMapsMatch_ = D0_T_Matrix_->getColMap()->isSameAs(*R11_->getColMap());
1266 D0_T_R11_colMapsMatch_ =
false;
1267 if (D0_T_R11_colMapsMatch_)
1268 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): D0_T and R11 have matching colMaps" << std::endl;
1273 if (parameterList_.isSublist(
"matvec params"))
1275 RCP<ParameterList> matvecParams = rcpFromRef(parameterList_.sublist(
"matvec params"));
1276 setMatvecParams(*SM_Matrix_, matvecParams);
1277 setMatvecParams(*D0_Matrix_, matvecParams);
1278 setMatvecParams(*P11_, matvecParams);
1279 if (!D0_T_Matrix_.is_null()) setMatvecParams(*D0_T_Matrix_, matvecParams);
1280 if (!R11_.is_null()) setMatvecParams(*R11_, matvecParams);
1281 if (!ImporterH_.is_null()) ImporterH_->setDistributorParameters(matvecParams);
1282 if (!Importer22_.is_null()) Importer22_->setDistributorParameters(matvecParams);
1284 if (!ImporterH_.is_null() && parameterList_.isSublist(
"refmaxwell: ImporterH params")){
1285 RCP<ParameterList> importerParams = rcpFromRef(parameterList_.sublist(
"refmaxwell: ImporterH params"));
1286 ImporterH_->setDistributorParameters(importerParams);
1288 if (!Importer22_.is_null() && parameterList_.isSublist(
"refmaxwell: Importer22 params")){
1289 RCP<ParameterList> importerParams = rcpFromRef(parameterList_.sublist(
"refmaxwell: Importer22 params"));
1290 Importer22_->setDistributorParameters(importerParams);
1295 #ifdef HAVE_MUELU_CUDA 1296 if (parameterList_.get<
bool>(
"refmaxwell: cuda profile setup",
false)) cudaProfilerStop();
1301 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1303 if (!R11_.is_null())
1304 P11res_ = MultiVectorFactory::Build(R11_->getRangeMap(), numVectors);
1306 P11res_ = MultiVectorFactory::Build(P11_->getDomainMap(), numVectors);
1307 if (D0_T_R11_colMapsMatch_)
1308 D0TR11Tmp_ = MultiVectorFactory::Build(R11_->getColMap(), numVectors);
1309 if (!ImporterH_.is_null()) {
1310 P11resTmp_ = MultiVectorFactory::Build(ImporterH_->getTargetMap(), numVectors);
1311 P11xTmp_ = MultiVectorFactory::Build(ImporterH_->getSourceMap(), numVectors);
1312 P11x_ = MultiVectorFactory::Build(ImporterH_->getTargetMap(), numVectors);
1314 P11x_ = MultiVectorFactory::Build(P11_->getDomainMap(), numVectors);
1315 if (!D0_T_Matrix_.is_null())
1316 D0res_ = MultiVectorFactory::Build(D0_T_Matrix_->getRangeMap(), numVectors);
1318 D0res_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
1319 if (!Importer22_.is_null()) {
1320 D0resTmp_ = MultiVectorFactory::Build(Importer22_->getTargetMap(), numVectors);
1321 D0xTmp_ = MultiVectorFactory::Build(Importer22_->getSourceMap(), numVectors);
1322 D0x_ = MultiVectorFactory::Build(Importer22_->getTargetMap(), numVectors);
1324 D0x_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
1325 residual_ = MultiVectorFactory::Build(SM_Matrix_->getDomainMap(), numVectors);
1329 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1331 if (dump_matrices_) {
1332 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1333 Xpetra::IO<SC, LO, GO, NO>::Write(name, A);
1338 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1340 if (dump_matrices_) {
1341 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1342 Xpetra::IO<SC, LO, GO, NO>::Write(name, X);
1347 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1349 if (dump_matrices_) {
1350 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1351 Xpetra::IO<coordinateType, LO, GO, NO>::Write(name, X);
1356 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1358 if (dump_matrices_) {
1359 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1360 std::ofstream out(name);
1361 for (
size_t i = 0; i < Teuchos::as<size_t>(v.size()); i++)
1362 out << v[i] <<
"\n";
1366 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1367 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1369 if (dump_matrices_) {
1370 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1371 std::ofstream out(name);
1372 auto vH = Kokkos::create_mirror_view (v);
1374 for (
size_t i = 0; i < vH.size(); i++)
1375 out << vH[i] <<
"\n";
1380 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1384 return Teuchos::rcp(
new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name)));
1387 return Teuchos::rcp(
new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), SM_Matrix_->getRowMap()->getComm().ptr()));
1389 return Teuchos::rcp(
new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), comm.ptr()));
1392 return Teuchos::null;
1396 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1409 const SC SC_ZERO = Teuchos::ScalarTraits<SC>::zero();
1410 const SC SC_ONE = Teuchos::ScalarTraits<SC>::one();
1411 const Scalar half = SC_ONE / (SC_ONE + SC_ONE);
1412 size_t dim = Nullspace_->getNumVectors();
1413 size_t numLocalRows = SM_Matrix_->getNodeNumRows();
1417 RCP<Matrix> P_nodal;
1418 bool read_P_from_file = parameterList_.get(
"refmaxwell: read_P_from_file",
false);
1419 if (read_P_from_file) {
1422 std::string P_filename = parameterList_.get(
"refmaxwell: P_filename",std::string(
"P.m"));
1423 std::string domainmap_filename = parameterList_.get(
"refmaxwell: P_domainmap_filename",std::string(
"domainmap_P.m"));
1424 std::string colmap_filename = parameterList_.get(
"refmaxwell: P_colmap_filename",std::string(
"colmap_P.m"));
1425 std::string coords_filename = parameterList_.get(
"refmaxwell: CoordsH",std::string(
"coordsH.m"));
1426 RCP<const Map> colmap = Xpetra::IO<SC, LO, GO, NO>::ReadMap(colmap_filename, A_nodal_Matrix_->getDomainMap()->lib(),A_nodal_Matrix_->getDomainMap()->getComm());
1427 RCP<const Map> domainmap = Xpetra::IO<SC, LO, GO, NO>::ReadMap(domainmap_filename, A_nodal_Matrix_->getDomainMap()->lib(),A_nodal_Matrix_->getDomainMap()->getComm());
1428 P_nodal = Xpetra::IO<SC, LO, GO, NO>::Read(P_filename, A_nodal_Matrix_->getDomainMap(), colmap, domainmap, A_nodal_Matrix_->getDomainMap());
1429 CoordsH_ = Xpetra::IO<typename RealValuedMultiVector::scalar_type, LO, GO, NO>::ReadMultiVector(coords_filename, domainmap);
1431 Level fineLevel, coarseLevel;
1437 fineLevel.
Set(
"A",A_nodal_Matrix_);
1438 fineLevel.
Set(
"Coordinates",Coords_);
1439 fineLevel.
Set(
"DofsPerNode",1);
1440 coarseLevel.
setlib(A_nodal_Matrix_->getDomainMap()->lib());
1441 fineLevel.
setlib(A_nodal_Matrix_->getDomainMap()->lib());
1442 coarseLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
1443 fineLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
1446 RCP<MultiVector> nullSpace = MultiVectorFactory::Build(A_nodal_Matrix_->getRowMap(),NSdim);
1447 nullSpace->putScalar(SC_ONE);
1448 fineLevel.
Set(
"Nullspace",nullSpace);
1450 RCP<Factory> amalgFact, dropFact, UncoupledAggFact, coarseMapFact, TentativePFact, Tfact, SaPFact;
1451 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1453 amalgFact = rcp(
new AmalgamationFactory_kokkos());
1454 dropFact = rcp(
new CoalesceDropFactory_kokkos());
1455 UncoupledAggFact = rcp(
new UncoupledAggregationFactory_kokkos());
1456 coarseMapFact = rcp(
new CoarseMapFactory_kokkos());
1457 TentativePFact = rcp(
new TentativePFactory_kokkos());
1458 if (parameterList_.get(
"multigrid algorithm",
"unsmoothed") ==
"sa")
1459 SaPFact = rcp(
new SaPFactory_kokkos());
1460 Tfact = rcp(
new CoordinatesTransferFactory_kokkos());
1469 if (parameterList_.get(
"multigrid algorithm",
"unsmoothed") ==
"sa")
1473 dropFact->SetFactory(
"UnAmalgamationInfo", amalgFact);
1474 double dropTol = parameterList_.get(
"aggregation: drop tol",0.0);
1475 std::string dropScheme = parameterList_.get(
"aggregation: drop scheme",
"classical");
1476 std::string distLaplAlgo = parameterList_.get(
"aggregation: distance laplacian algo",
"default");
1477 dropFact->SetParameter(
"aggregation: drop tol",Teuchos::ParameterEntry(dropTol));
1478 dropFact->SetParameter(
"aggregation: drop scheme",Teuchos::ParameterEntry(dropScheme));
1480 dropFact->SetParameter(
"aggregation: distance laplacian algo",Teuchos::ParameterEntry(distLaplAlgo));
1482 UncoupledAggFact->SetFactory(
"Graph", dropFact);
1483 int minAggSize = parameterList_.get(
"aggregation: min agg size",2);
1484 UncoupledAggFact->SetParameter(
"aggregation: min agg size",Teuchos::ParameterEntry(minAggSize));
1485 int maxAggSize = parameterList_.get(
"aggregation: max agg size",-1);
1486 UncoupledAggFact->SetParameter(
"aggregation: max agg size",Teuchos::ParameterEntry(maxAggSize));
1488 coarseMapFact->SetFactory(
"Aggregates", UncoupledAggFact);
1490 TentativePFact->SetFactory(
"Aggregates", UncoupledAggFact);
1491 TentativePFact->SetFactory(
"UnAmalgamationInfo", amalgFact);
1492 TentativePFact->SetFactory(
"CoarseMap", coarseMapFact);
1494 Tfact->SetFactory(
"Aggregates", UncoupledAggFact);
1495 Tfact->SetFactory(
"CoarseMap", coarseMapFact);
1497 if (parameterList_.get(
"multigrid algorithm",
"unsmoothed") ==
"sa") {
1498 SaPFact->SetFactory(
"P", TentativePFact);
1499 coarseLevel.
Request(
"P", SaPFact.get());
1501 coarseLevel.
Request(
"P",TentativePFact.get());
1502 coarseLevel.
Request(
"Coordinates",Tfact.get());
1504 RCP<AggregationExportFactory> aggExport;
1505 if (parameterList_.get(
"aggregation: export visualization data",
false)) {
1507 ParameterList aggExportParams;
1508 aggExportParams.set(
"aggregation: output filename",
"aggs.vtk");
1509 aggExportParams.set(
"aggregation: output file: agg style",
"Jacks");
1510 aggExport->SetParameterList(aggExportParams);
1512 aggExport->SetFactory(
"Aggregates", UncoupledAggFact);
1513 aggExport->SetFactory(
"UnAmalgamationInfo", amalgFact);
1514 fineLevel.
Request(
"Aggregates",UncoupledAggFact.get());
1515 fineLevel.
Request(
"UnAmalgamationInfo",amalgFact.get());
1518 if (parameterList_.get(
"multigrid algorithm",
"unsmoothed") ==
"sa")
1519 coarseLevel.
Get(
"P",P_nodal,SaPFact.get());
1521 coarseLevel.
Get(
"P",P_nodal,TentativePFact.get());
1522 coarseLevel.
Get(
"Coordinates",CoordsH_,Tfact.get());
1524 if (parameterList_.get(
"aggregation: export visualization data",
false))
1525 aggExport->Build(fineLevel, coarseLevel);
1527 dump(*P_nodal,
"P_nodal.m");
1529 RCP<CrsMatrix> D0Crs = rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix();
1532 RCP<CrsMatrix> P_nodal_imported;
1533 int numProcs = P_nodal->getDomainMap()->getComm()->getSize();
1535 RCP<CrsMatrixWrap> P_nodal_temp;
1536 RCP<const Map> targetMap = D0Crs->getColMap();
1537 P_nodal_temp = rcp(
new CrsMatrixWrap(targetMap));
1538 RCP<const Import> importer = D0Crs->getCrsGraph()->getImporter();
1539 P_nodal_temp->doImport(*P_nodal, *importer, Xpetra::INSERT);
1540 P_nodal_temp->fillComplete(rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix()->getDomainMap(),
1541 rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix()->getRangeMap());
1542 P_nodal_imported = P_nodal_temp->getCrsMatrix();
1543 dump(*P_nodal_temp,
"P_nodal_imported.m");
1545 P_nodal_imported = rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix();
1547 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1550 using ATS = Kokkos::ArithTraits<SC>;
1551 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1554 typedef typename Matrix::local_matrix_type KCRS;
1555 typedef typename KCRS::device_type device_t;
1556 typedef typename KCRS::StaticCrsGraphType graph_t;
1557 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
1558 typedef typename graph_t::entries_type::non_const_type lno_nnz_view_t;
1559 typedef typename KCRS::values_type::non_const_type scalar_view_t;
1562 auto localP = P_nodal_imported->getLocalMatrix();
1563 auto localD0 = D0_Matrix_->getLocalMatrix();
1568 std::string defaultAlgo =
"mat-mat";
1569 std::string algo = parameterList_.get(
"refmaxwell: prolongator compute algorithm",defaultAlgo);
1571 if (algo ==
"mat-mat") {
1572 RCP<Matrix> D0_P_nodal = MatrixFactory::Build(SM_Matrix_->getRowMap());
1573 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
false,*P_nodal,
false,*D0_P_nodal,
true,
true);
1575 #ifdef HAVE_MUELU_DEBUG 1576 TEUCHOS_ASSERT(D0_P_nodal->getColMap()->isSameAs(*P_nodal_imported->getColMap()));
1580 auto localD0P = D0_P_nodal->getLocalMatrix();
1583 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
1584 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
1586 size_t nnzEstimate = dim*localD0P.graph.entries.size();
1587 lno_view_t P11rowptr(
"P11_rowptr", numLocalRows+1);
1588 lno_nnz_view_t P11colind(
"P11_colind",nnzEstimate);
1589 scalar_view_t P11vals(
"P11_vals",nnzEstimate);
1592 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_adjustRowptr", range_type(0,numLocalRows+1),
1593 KOKKOS_LAMBDA(
const size_t i) {
1594 P11rowptr(i) = dim*localD0P.graph.row_map(i);
1598 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_adjustColind", range_type(0,localD0P.graph.entries.size()),
1599 KOKKOS_LAMBDA(
const size_t jj) {
1600 for (
size_t k = 0; k < dim; k++) {
1601 P11colind(dim*jj+k) = dim*localD0P.graph.entries(jj)+k;
1602 P11vals(dim*jj+k) = SC_ZERO;
1606 auto localNullspace = Nullspace_->template getLocalView<device_t>();
1609 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
1613 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
1615 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
1617 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_enterValues_D0wZeros", range_type(0,numLocalRows),
1618 KOKKOS_LAMBDA(
const size_t i) {
1619 for (
size_t ll = localD0.graph.row_map(i); ll < localD0.graph.row_map(i+1); ll++) {
1620 LO l = localD0.graph.entries(ll);
1621 SC p = localD0.values(ll);
1622 if (impl_ATS::magnitude(p) < tol)
1624 for (
size_t jj = localP.graph.row_map(l); jj < localP.graph.row_map(l+1); jj++) {
1625 LO j = localP.graph.entries(jj);
1626 SC v = localP.values(jj);
1627 for (
size_t k = 0; k < dim; k++) {
1629 SC n = localNullspace(i,k);
1631 for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
1632 if (P11colind(m) == jNew)
1634 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA) 1635 TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
1637 P11vals(m) += half * v * n;
1644 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_enterValues", range_type(0,numLocalRows),
1645 KOKKOS_LAMBDA(
const size_t i) {
1646 for (
size_t ll = localD0.graph.row_map(i); ll < localD0.graph.row_map(i+1); ll++) {
1647 LO l = localD0.graph.entries(ll);
1648 for (
size_t jj = localP.graph.row_map(l); jj < localP.graph.row_map(l+1); jj++) {
1649 LO j = localP.graph.entries(jj);
1650 SC v = localP.values(jj);
1651 for (
size_t k = 0; k < dim; k++) {
1653 SC n = localNullspace(i,k);
1655 for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
1656 if (P11colind(m) == jNew)
1658 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA) 1659 TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
1661 P11vals(m) += half * v * n;
1668 P11_ = rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0));
1669 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
1670 P11Crs->setAllValues(P11rowptr, P11colind, P11vals);
1671 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
1674 TEUCHOS_TEST_FOR_EXCEPTION(
false,std::invalid_argument,algo <<
" is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
1676 #endif // ifdef(HAVE_MUELU_KOKKOS_REFACTOR) 1679 ArrayRCP<ArrayRCP<const SC> > nullspaceRCP(dim);
1680 ArrayRCP<ArrayView<const SC> > nullspace(dim);
1681 for(
size_t i=0; i<dim; i++) {
1682 nullspaceRCP[i] = Nullspace_->getData(i);
1683 nullspace[i] = nullspaceRCP[i]();
1687 ArrayRCP<const size_t> Prowptr_RCP, D0rowptr_RCP;
1688 ArrayRCP<const LO> Pcolind_RCP, D0colind_RCP;
1689 ArrayRCP<const SC> Pvals_RCP, D0vals_RCP;
1690 ArrayRCP<size_t> P11rowptr_RCP;
1691 ArrayRCP<LO> P11colind_RCP;
1692 ArrayRCP<SC> P11vals_RCP;
1694 P_nodal_imported->getAllValues(Prowptr_RCP, Pcolind_RCP, Pvals_RCP);
1695 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->getAllValues(D0rowptr_RCP, D0colind_RCP, D0vals_RCP);
1700 ArrayView<const size_t> Prowptr, D0rowptr;
1701 ArrayView<const LO> Pcolind, D0colind;
1702 ArrayView<const SC> Pvals, D0vals;
1703 Prowptr = Prowptr_RCP(); Pcolind = Pcolind_RCP(); Pvals = Pvals_RCP();
1704 D0rowptr = D0rowptr_RCP(); D0colind = D0colind_RCP(); D0vals = D0vals_RCP();
1712 std::string defaultAlgo =
"mat-mat";
1713 std::string algo = parameterList_.get(
"refmaxwell: prolongator compute algorithm",defaultAlgo);
1715 if (algo ==
"mat-mat") {
1716 RCP<Matrix> D0_P_nodal = MatrixFactory::Build(SM_Matrix_->getRowMap());
1717 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
false,*P_nodal,
false,*D0_P_nodal,
true,
true);
1720 ArrayRCP<const size_t> D0Prowptr_RCP;
1721 ArrayRCP<const LO> D0Pcolind_RCP;
1722 ArrayRCP<const SC> D0Pvals_RCP;
1723 rcp_dynamic_cast<CrsMatrixWrap>(D0_P_nodal)->getCrsMatrix()->getAllValues(D0Prowptr_RCP, D0Pcolind_RCP, D0Pvals_RCP);
1728 ArrayView<const size_t> D0Prowptr;
1729 ArrayView<const LO> D0Pcolind;
1730 D0Prowptr = D0Prowptr_RCP(); D0Pcolind = D0Pcolind_RCP();
1733 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
1734 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
1735 P11_ = rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0));
1736 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
1737 size_t nnzEstimate = dim*D0Prowptr[numLocalRows];
1738 P11Crs->allocateAllValues(nnzEstimate, P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1740 ArrayView<size_t> P11rowptr = P11rowptr_RCP();
1741 ArrayView<LO> P11colind = P11colind_RCP();
1742 ArrayView<SC> P11vals = P11vals_RCP();
1745 for (
size_t i = 0; i < numLocalRows+1; i++) {
1746 P11rowptr[i] = dim*D0Prowptr[i];
1750 for (
size_t jj = 0; jj < (size_t) D0Prowptr[numLocalRows]; jj++)
1751 for (
size_t k = 0; k < dim; k++) {
1752 P11colind[dim*jj+k] = dim*D0Pcolind[jj]+k;
1753 P11vals[dim*jj+k] = SC_ZERO;
1756 RCP<const Map> P_nodal_imported_colmap = P_nodal_imported->getColMap();
1757 RCP<const Map> D0_P_nodal_colmap = D0_P_nodal->getColMap();
1759 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
1763 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
1765 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
1766 for (
size_t i = 0; i < numLocalRows; i++) {
1767 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1768 LO l = D0colind[ll];
1770 if (Teuchos::ScalarTraits<Scalar>::magnitude(p) < tol)
1772 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1774 j = D0_P_nodal_colmap->getLocalElement(P_nodal_imported_colmap->getGlobalElement(j));
1776 for (
size_t k = 0; k < dim; k++) {
1778 SC n = nullspace[k][i];
1780 for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
1781 if (P11colind[m] == jNew)
1783 #ifdef HAVE_MUELU_DEBUG 1784 TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
1786 P11vals[m] += half * v * n;
1793 for (
size_t i = 0; i < numLocalRows; i++) {
1794 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1795 LO l = D0colind[ll];
1796 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1798 j = D0_P_nodal_colmap->getLocalElement(P_nodal_imported_colmap->getGlobalElement(j));
1800 for (
size_t k = 0; k < dim; k++) {
1802 SC n = nullspace[k][i];
1804 for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
1805 if (P11colind[m] == jNew)
1807 #ifdef HAVE_MUELU_DEBUG 1808 TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
1810 P11vals[m] += half * v * n;
1817 P11Crs->setAllValues(P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1818 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
1820 }
else if (algo ==
"gustavson") {
1822 LO maxP11col = dim * P_nodal_imported->getColMap()->getMaxLocalIndex();
1823 const size_t ST_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1824 Array<size_t> P11_status(dim*maxP11col, ST_INVALID);
1826 size_t nnz_alloc = dim*D0vals_RCP.size();
1829 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
1830 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
1831 P11_ = rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0));
1832 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
1833 P11Crs->allocateAllValues(nnz_alloc, P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1835 ArrayView<size_t> P11rowptr = P11rowptr_RCP();
1836 ArrayView<LO> P11colind = P11colind_RCP();
1837 ArrayView<SC> P11vals = P11vals_RCP();
1840 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
1844 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
1846 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
1849 for (
size_t i = 0; i < numLocalRows; i++) {
1851 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1852 LO l = D0colind[ll];
1854 if (Teuchos::ScalarTraits<Scalar>::magnitude(p) < tol)
1856 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1859 for (
size_t k = 0; k < dim; k++) {
1861 SC n = nullspace[k][i];
1863 if (P11_status[jNew] == ST_INVALID || P11_status[jNew] < nnz_old) {
1864 P11_status[jNew] = nnz;
1865 P11colind[nnz] = jNew;
1866 P11vals[nnz] = half * v * n;
1871 P11vals[P11_status[jNew]] += half * v * n;
1880 P11rowptr[numLocalRows] = nnz;
1884 for (
size_t i = 0; i < numLocalRows; i++) {
1886 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1887 LO l = D0colind[ll];
1888 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1891 for (
size_t k = 0; k < dim; k++) {
1893 SC n = nullspace[k][i];
1895 if (P11_status[jNew] == ST_INVALID || P11_status[jNew] < nnz_old) {
1896 P11_status[jNew] = nnz;
1897 P11colind[nnz] = jNew;
1898 P11vals[nnz] = half * v * n;
1903 P11vals[P11_status[jNew]] += half * v * n;
1912 P11rowptr[numLocalRows] = nnz;
1915 if (blockDomainMap->lib() == Xpetra::UseTpetra) {
1919 P11vals_RCP.resize(nnz);
1920 P11colind_RCP.resize(nnz);
1923 P11Crs->setAllValues(P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1924 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
1926 TEUCHOS_TEST_FOR_EXCEPTION(
false,std::invalid_argument,algo <<
" is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
1927 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1932 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1934 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu RefMaxwell: Build coarse (1,1) matrix");
1936 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
1939 RCP<Matrix> Matrix1;
1941 Level fineLevel, coarseLevel;
1947 fineLevel.
Set(
"A",SM_Matrix_);
1948 coarseLevel.
Set(
"P",P11_);
1950 coarseLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
1951 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
1952 coarseLevel.setObjectLabel(
"RefMaxwell (1,1)");
1953 fineLevel.setObjectLabel(
"RefMaxwell (1,1)");
1955 RCP<RAPFactory> rapFact = rcp(
new RAPFactory());
1956 ParameterList rapList = *(rapFact->GetValidParameterList());
1957 rapList.set(
"transpose: use implicit",
true);
1958 rapList.set(
"rap: fix zero diagonals", parameterList_.get<
bool>(
"rap: fix zero diagonals",
true));
1959 rapList.set(
"rap: triple product", parameterList_.get<
bool>(
"rap: triple product",
false));
1960 rapFact->SetParameterList(rapList);
1962 if (precList11_.isType<std::string>(
"reuse: type") && precList11_.get<std::string>(
"reuse: type") !=
"none") {
1963 if (!parameterList_.get<
bool>(
"rap: triple product",
false))
1964 coarseLevel.
Request(
"AP reuse data", rapFact.get());
1965 coarseLevel.
Request(
"RAP reuse data", rapFact.get());
1968 if (!AH_AP_reuse_data_.is_null()) {
1969 coarseLevel.
AddKeepFlag(
"AP reuse data", rapFact.get());
1970 coarseLevel.
Set<Teuchos::RCP<Teuchos::ParameterList> >(
"AP reuse data", AH_AP_reuse_data_, rapFact.get());
1972 if (!AH_RAP_reuse_data_.is_null()) {
1973 coarseLevel.
AddKeepFlag(
"RAP reuse data", rapFact.get());
1974 coarseLevel.
Set<Teuchos::RCP<Teuchos::ParameterList> >(
"RAP reuse data", AH_RAP_reuse_data_, rapFact.get());
1977 coarseLevel.
Request(
"A", rapFact.get());
1979 Matrix1 = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
1980 if (precList11_.isType<std::string>(
"reuse: type") && precList11_.get<std::string>(
"reuse: type") !=
"none") {
1981 if (!parameterList_.get<
bool>(
"rap: triple product",
false))
1982 AH_AP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data", rapFact.get());
1983 AH_RAP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"RAP reuse data", rapFact.get());
1988 AH_ = Teuchos::null;
1990 if(disable_addon_==
true) {
1995 if (Addon_Matrix_.is_null()) {
1996 RCP<Teuchos::TimeMonitor> tmAddon = getTimer(
"MueLu RefMaxwell: Build coarse addon matrix");
1998 TEUCHOS_TEST_FOR_EXCEPTION(M0inv_Matrix_==Teuchos::null,std::invalid_argument,
1999 "MueLu::RefMaxwell::formCoarseMatrix(): Inverse of " 2000 "lumped mass matrix required for add-on (i.e. M0inv_Matrix is null)");
2003 RCP<Matrix> Zaux = MatrixFactory::Build(M1_Matrix_->getRowMap());
2004 RCP<Matrix> Z = MatrixFactory::Build(D0_Matrix_->getDomainMap());
2007 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M1_Matrix_,
false,*P11_,
false,*Zaux,
true,
true);
2009 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
true,*Zaux,
false,*Z,
true,
true);
2012 RCP<Matrix> Matrix2 = MatrixFactory::Build(Z->getDomainMap());
2013 if (M0inv_Matrix_->getGlobalMaxNumRowEntries()<=1) {
2016 RCP<Vector> diag = VectorFactory::Build(M0inv_Matrix_->getRowMap());
2017 M0inv_Matrix_->getLocalDiagCopy(*diag);
2018 ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
2019 for (
size_t j=0; j < diag->getMap()->getNodeNumElements(); j++) {
2020 diagVals[j] = Teuchos::ScalarTraits<Scalar>::squareroot(diagVals[j]);
2022 if (Z->getRowMap()->isSameAs(*(diag->getMap())))
2023 Z->leftScale(*diag);
2025 RCP<Import> importer = ImportFactory::Build(diag->getMap(),Z->getRowMap());
2026 RCP<Vector> diag2 = VectorFactory::Build(Z->getRowMap());
2027 diag2->doImport(*diag,*importer,Xpetra::INSERT);
2028 Z->leftScale(*diag2);
2030 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*Z,
true,*Z,
false,*Matrix2,
true,
true);
2031 }
else if (parameterList_.get<
bool>(
"rap: triple product",
false) ==
false) {
2032 RCP<Matrix> C2 = MatrixFactory::Build(M0inv_Matrix_->getRowMap());
2034 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M0inv_Matrix_,
false,*Z,
false,*C2,
true,
true);
2036 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*Z,
true,*C2,
false,*Matrix2,
true,
true);
2039 Xpetra::TripleMatrixMultiply<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
2040 MultiplyRAP(*Z,
true, *M0inv_Matrix_,
false, *Z,
false, *Matrix2,
true,
true);
2043 if (precList11_.isType<std::string>(
"reuse: type") && precList11_.get<std::string>(
"reuse: type") !=
"none")
2044 Addon_Matrix_ = Matrix2;
2047 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*Matrix1,
false,one,*Matrix2,
false,one,AH_,GetOStream(
Runtime0));
2048 AH_->fillComplete();
2051 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*Matrix1,
false,one,*Addon_Matrix_,
false,one,AH_,GetOStream(
Runtime0));
2052 AH_->fillComplete();
2057 size_t dim = Nullspace_->getNumVectors();
2058 AH_->SetFixedBlockSize(dim);
2059 AH_->setObjectLabel(
"RefMaxwell coarse (1,1)");
2064 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2066 bool reuse = !SM_Matrix_.is_null();
2067 SM_Matrix_ = SM_Matrix_new;
2068 dump(*SM_Matrix_,
"SM.m");
2074 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2077 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
2081 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu RefMaxwell: residual calculation");
2087 if (implicitTranspose_) {
2089 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu RefMaxwell: restriction coarse (1,1) (implicit)");
2090 P11_->apply(*residual_,*P11res_,Teuchos::TRANS);
2093 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: restriction (2,2) (implicit)");
2094 D0_Matrix_->apply(*residual_,*D0res_,Teuchos::TRANS);
2097 if (D0_T_R11_colMapsMatch_) {
2100 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: restrictions import");
2101 D0TR11Tmp_->doImport(*residual_, *rcp_dynamic_cast<CrsMatrixWrap>(R11_)->getCrsMatrix()->getCrsGraph()->getImporter(), Xpetra::INSERT);
2104 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: restriction (2,2) (explicit)");
2105 rcp_dynamic_cast<TpetraCrsMatrix>(rcp_dynamic_cast<CrsMatrixWrap>(D0_T_Matrix_)->getCrsMatrix())->getTpetra_CrsMatrix()->localApply(toTpetra(*D0TR11Tmp_),toTpetra(*D0res_),Teuchos::NO_TRANS);
2108 RCP<Teuchos::TimeMonitor> tmP11 = getTimer(
"MueLu RefMaxwell: restriction coarse (1,1) (explicit)");
2109 rcp_dynamic_cast<TpetraCrsMatrix>(rcp_dynamic_cast<CrsMatrixWrap>(R11_)->getCrsMatrix())->getTpetra_CrsMatrix()->localApply(toTpetra(*D0TR11Tmp_),toTpetra(*P11res_),Teuchos::NO_TRANS);
2113 RCP<Teuchos::TimeMonitor> tmP11 = getTimer(
"MueLu RefMaxwell: restriction coarse (1,1) (explicit)");
2114 R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
2117 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: restriction (2,2) (explicit)");
2118 D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
2125 RCP<Teuchos::TimeMonitor> tmSubSolves = getTimer(
"MueLu RefMaxwell: subsolves");
2129 if (!ImporterH_.is_null() && !implicitTranspose_) {
2130 RCP<Teuchos::TimeMonitor> tmH = getTimer(
"MueLu RefMaxwell: import coarse (1,1)");
2131 P11resTmp_->doImport(*P11res_, *ImporterH_, Xpetra::INSERT);
2132 P11res_.swap(P11resTmp_);
2134 if (!Importer22_.is_null() && !implicitTranspose_) {
2135 RCP<Teuchos::TimeMonitor> tm22 = getTimer(
"MueLu RefMaxwell: import (2,2)");
2136 D0resTmp_->doImport(*D0res_, *Importer22_, Xpetra::INSERT);
2137 D0res_.swap(D0resTmp_);
2141 if (!AH_.is_null()) {
2142 RCP<Teuchos::TimeMonitor> tmH = getTimer(
"MueLu RefMaxwell: solve coarse (1,1)", AH_->getRowMap()->getComm());
2144 RCP<const Map> origXMap = P11x_->getMap();
2145 RCP<const Map> origRhsMap = P11res_->getMap();
2148 P11res_->replaceMap(AH_->getRangeMap());
2149 P11x_ ->replaceMap(AH_->getDomainMap());
2150 HierarchyH_->Iterate(*P11res_, *P11x_, numItersH_,
true);
2151 P11x_ ->replaceMap(origXMap);
2152 P11res_->replaceMap(origRhsMap);
2156 if (!A22_.is_null()) {
2157 RCP<Teuchos::TimeMonitor> tm22 = getTimer(
"MueLu RefMaxwell: solve (2,2)", A22_->getRowMap()->getComm());
2159 RCP<const Map> origXMap = D0x_->getMap();
2160 RCP<const Map> origRhsMap = D0res_->getMap();
2163 D0res_->replaceMap(A22_->getRangeMap());
2164 D0x_ ->replaceMap(A22_->getDomainMap());
2165 Hierarchy22_->Iterate(*D0res_, *D0x_, numIters22_,
true);
2166 D0x_ ->replaceMap(origXMap);
2167 D0res_->replaceMap(origRhsMap);
2172 if (fuseProlongationAndUpdate_) {
2174 RCP<Teuchos::TimeMonitor> tmP11 = getTimer(
"MueLu RefMaxwell: prolongation coarse (1,1) (fused)");
2175 P11_->apply(*P11x_,X,Teuchos::NO_TRANS,one,one);
2179 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: prolongation (2,2) (fused)");
2180 D0_Matrix_->apply(*D0x_,X,Teuchos::NO_TRANS,one,one);
2184 RCP<Teuchos::TimeMonitor> tmP11 = getTimer(
"MueLu RefMaxwell: prolongation coarse (1,1) (unfused)");
2185 P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
2189 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: prolongation (2,2) (unfused)");
2190 D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS,one,one);
2194 RCP<Teuchos::TimeMonitor> tmUpdate = getTimer(
"MueLu RefMaxwell: update");
2195 X.update(one, *residual_, one);
2199 if (!ImporterH_.is_null()) {
2200 P11res_.swap(P11resTmp_);
2202 if (!Importer22_.is_null()) {
2203 D0res_.swap(D0resTmp_);
2209 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2222 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2234 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2237 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
2240 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu RefMaxwell: residual calculation");
2242 if (implicitTranspose_)
2243 P11_->apply(*residual_,*P11res_,Teuchos::TRANS);
2245 R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
2249 if (!ImporterH_.is_null() && !implicitTranspose_) {
2250 RCP<Teuchos::TimeMonitor> tmH = getTimer(
"MueLu RefMaxwell: import coarse (1,1)");
2251 P11resTmp_->doImport(*P11res_, *ImporterH_, Xpetra::INSERT);
2252 P11res_.swap(P11resTmp_);
2254 if (!AH_.is_null()) {
2255 RCP<Teuchos::TimeMonitor> tmH = getTimer(
"MueLu RefMaxwell: solve coarse (1,1)", AH_->getRowMap()->getComm());
2257 RCP<const Map> origXMap = P11x_->getMap();
2258 RCP<const Map> origRhsMap = P11res_->getMap();
2261 P11res_->replaceMap(AH_->getRangeMap());
2262 P11x_ ->replaceMap(AH_->getDomainMap());
2263 HierarchyH_->Iterate(*P11res_, *P11x_, numItersH_,
true);
2264 P11x_ ->replaceMap(origXMap);
2265 P11res_->replaceMap(origRhsMap);
2270 RCP<Teuchos::TimeMonitor> tmUp = getTimer(
"MueLu RefMaxwell: update");
2271 P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
2272 X.update(one, *residual_, one);
2274 if (!ImporterH_.is_null()) {
2275 P11res_.swap(P11resTmp_);
2281 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2284 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
2287 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu RefMaxwell: residual calculation");
2289 if (implicitTranspose_)
2290 D0_Matrix_->apply(*residual_,*D0res_,Teuchos::TRANS);
2292 D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
2296 if (!Importer22_.is_null() && !implicitTranspose_) {
2297 RCP<Teuchos::TimeMonitor> tm22 = getTimer(
"MueLu RefMaxwell: import (2,2)");
2298 D0resTmp_->doImport(*D0res_, *Importer22_, Xpetra::INSERT);
2299 D0res_.swap(D0resTmp_);
2301 if (!A22_.is_null()) {
2302 RCP<Teuchos::TimeMonitor> tm22 = getTimer(
"MueLu RefMaxwell: solve (2,2)", A22_->getRowMap()->getComm());
2304 RCP<const Map> origXMap = D0x_->getMap();
2305 RCP<const Map> origRhsMap = D0res_->getMap();
2308 D0res_->replaceMap(A22_->getRangeMap());
2309 D0x_ ->replaceMap(A22_->getDomainMap());
2310 Hierarchy22_->Iterate(*D0res_, *D0x_, numIters22_,
true);
2311 D0x_ ->replaceMap(origXMap);
2312 D0res_->replaceMap(origRhsMap);
2317 RCP<Teuchos::TimeMonitor> tmUp = getTimer(
"MueLu RefMaxwell: update");
2318 D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS);
2319 X.update(one, *residual_, one);
2321 if (!Importer22_.is_null()) {
2322 D0res_.swap(D0resTmp_);
2328 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2334 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu RefMaxwell: solve");
2337 if (X.getNumVectors() != P11res_->getNumVectors())
2338 allocateMemory(X.getNumVectors());
2342 RCP<Teuchos::TimeMonitor> tmSm = getTimer(
"MueLu RefMaxwell: smoothing");
2344 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 2345 if (useHiptmairSmoothing_) {
2348 hiptmairPreSmoother_->apply(tRHS, tX);
2352 PreSmoother_->Apply(X, RHS, use_as_preconditioner_);
2356 if(mode_==
"additive")
2357 applyInverseAdditive(RHS,X);
2358 else if(mode_==
"121")
2359 applyInverse121(RHS,X);
2360 else if(mode_==
"212")
2361 applyInverse212(RHS,X);
2366 else if(mode_==
"none") {
2370 applyInverseAdditive(RHS,X);
2374 RCP<Teuchos::TimeMonitor> tmSm = getTimer(
"MueLu RefMaxwell: smoothing");
2376 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 2377 if (useHiptmairSmoothing_)
2381 hiptmairPostSmoother_->apply(tRHS, tX);
2385 PostSmoother_->Apply(X, RHS,
false);
2391 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2397 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2400 const Teuchos::RCP<Matrix> & Ms_Matrix,
2401 const Teuchos::RCP<Matrix> & M0inv_Matrix,
2402 const Teuchos::RCP<Matrix> & M1_Matrix,
2403 const Teuchos::RCP<MultiVector> & Nullspace,
2404 const Teuchos::RCP<RealValuedMultiVector> & Coords,
2405 Teuchos::ParameterList& List)
2408 TEUCHOS_ASSERT(D0_Matrix!=Teuchos::null);
2409 TEUCHOS_ASSERT(Ms_Matrix!=Teuchos::null);
2410 TEUCHOS_ASSERT(M1_Matrix!=Teuchos::null);
2412 #ifdef HAVE_MUELU_DEBUG 2413 TEUCHOS_ASSERT(M1_Matrix->getDomainMap()->isSameAs(*M1_Matrix->getRangeMap()));
2414 TEUCHOS_ASSERT(M1_Matrix->getDomainMap()->isSameAs(*M1_Matrix->getRowMap()));
2415 TEUCHOS_ASSERT(M1_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getRangeMap()));
2417 TEUCHOS_ASSERT(Ms_Matrix->getDomainMap()->isSameAs(*Ms_Matrix->getRangeMap()));
2418 TEUCHOS_ASSERT(Ms_Matrix->getDomainMap()->isSameAs(*Ms_Matrix->getRowMap()));
2419 TEUCHOS_ASSERT(Ms_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getRangeMap()));
2421 TEUCHOS_ASSERT(D0_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getRowMap()));
2423 if (!M0inv_Matrix) {
2424 TEUCHOS_ASSERT(M0inv_Matrix->getDomainMap()->isSameAs(*M0inv_Matrix->getRowMap()));
2425 TEUCHOS_ASSERT(M0inv_Matrix->getDomainMap()->isSameAs(*M0inv_Matrix->getRangeMap()));
2426 TEUCHOS_ASSERT(M0inv_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getDomainMap()));
2430 HierarchyH_ = Teuchos::null;
2431 Hierarchy22_ = Teuchos::null;
2432 PreSmoother_ = Teuchos::null;
2433 PostSmoother_ = Teuchos::null;
2434 disable_addon_ =
false;
2438 setParameters(List);
2440 if (D0_Matrix->getRowMap()->lib() == Xpetra::UseTpetra) {
2445 RCP<Matrix> D0copy = MatrixFactory::Build(D0_Matrix->getRowMap(), D0_Matrix->getColMap(), 0);
2446 RCP<CrsMatrix> D0copyCrs = rcp_dynamic_cast<CrsMatrixWrap>(D0copy,
true)->getCrsMatrix();
2447 ArrayRCP<const size_t> D0rowptr_RCP;
2448 ArrayRCP<const LO> D0colind_RCP;
2449 ArrayRCP<const SC> D0vals_RCP;
2450 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getAllValues(D0rowptr_RCP,
2454 ArrayRCP<size_t> D0copyrowptr_RCP;
2455 ArrayRCP<LO> D0copycolind_RCP;
2456 ArrayRCP<SC> D0copyvals_RCP;
2457 D0copyCrs->allocateAllValues(D0vals_RCP.size(),D0copyrowptr_RCP,D0copycolind_RCP,D0copyvals_RCP);
2458 D0copyrowptr_RCP.deepCopy(D0rowptr_RCP());
2459 D0copycolind_RCP.deepCopy(D0colind_RCP());
2460 D0copyvals_RCP.deepCopy(D0vals_RCP());
2461 D0copyCrs->setAllValues(D0copyrowptr_RCP,
2464 D0copyCrs->expertStaticFillComplete(D0_Matrix->getDomainMap(), D0_Matrix->getRangeMap(),
2465 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getImporter(),
2466 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getExporter());
2467 D0_Matrix_ = D0copy;
2469 D0_Matrix_ = MatrixFactory::BuildCopy(D0_Matrix);
2472 M0inv_Matrix_ = M0inv_Matrix;
2473 Ms_Matrix_ = Ms_Matrix;
2474 M1_Matrix_ = M1_Matrix;
2476 Nullspace_ = Nullspace;
2478 dump(*D0_Matrix_,
"D0_clean.m");
2479 dump(*Ms_Matrix_,
"Ms.m");
2480 dump(*M1_Matrix_,
"M1.m");
2481 if (!M0inv_Matrix_.is_null()) dump(*M0inv_Matrix_,
"M0inv.m");
2482 if (!Coords_.is_null()) dumpCoords(*Coords_,
"coords.m");
2487 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2489 describe(Teuchos::FancyOStream& out,
const Teuchos::EVerbosityLevel )
const {
2491 std::ostringstream oss;
2493 RCP<const Teuchos::Comm<int> > comm = SM_Matrix_->getDomainMap()->getComm();
2497 if (!A22_.is_null())
2498 root = comm->getRank();
2503 reduceAll(*comm, Teuchos::REDUCE_MAX, root, Teuchos::ptr(&actualRoot));
2508 oss <<
"\n--------------------------------------------------------------------------------\n" <<
2509 "--- RefMaxwell Summary ---\n" 2510 "--------------------------------------------------------------------------------" << std::endl;
2516 SM_Matrix_->getRowMap()->getComm()->barrier();
2518 numRows = SM_Matrix_->getGlobalNumRows();
2519 nnz = SM_Matrix_->getGlobalNumEntries();
2521 Xpetra::global_size_t tt = numRows;
2522 int rowspacer = 3;
while (tt != 0) { tt /= 10; rowspacer++; }
2524 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
2526 oss <<
"block " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " << std::setw(9) <<
" nnz/row" << std::endl;
2527 oss <<
"(1, 1)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
2529 if (!A22_.is_null()) {
2531 numRows = A22_->getGlobalNumRows();
2532 nnz = A22_->getGlobalNumEntries();
2534 oss <<
"(2, 2)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
2539 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 2540 if (useHiptmairSmoothing_) {
2541 if (hiptmairPreSmoother_ != null && hiptmairPreSmoother_ == hiptmairPostSmoother_)
2542 oss <<
"Smoother both : " << hiptmairPreSmoother_->description() << std::endl;
2544 oss <<
"Smoother pre : " 2545 << (hiptmairPreSmoother_ != null ? hiptmairPreSmoother_->description() :
"no smoother") << std::endl;
2546 oss <<
"Smoother post : " 2547 << (hiptmairPostSmoother_ != null ? hiptmairPostSmoother_->description() :
"no smoother") << std::endl;
2553 if (PreSmoother_ != null && PreSmoother_ == PostSmoother_)
2554 oss <<
"Smoother both : " << PreSmoother_->description() << std::endl;
2556 oss <<
"Smoother pre : " 2557 << (PreSmoother_ != null ? PreSmoother_->description() :
"no smoother") << std::endl;
2558 oss <<
"Smoother post : " 2559 << (PostSmoother_ != null ? PostSmoother_->description() :
"no smoother") << std::endl;
2564 std::string outstr = oss.str();
2567 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
2568 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
2570 int strLength = outstr.size();
2571 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
2572 if (comm->getRank() != root)
2573 outstr.resize(strLength);
2574 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
2579 if (!HierarchyH_.is_null())
2580 HierarchyH_->describe(out, GetVerbLevel());
2582 if (!Hierarchy22_.is_null())
2583 Hierarchy22_->describe(out, GetVerbLevel());
2587 std::ostringstream oss2;
2589 oss2 <<
"Sub-solver distribution over ranks" << std::endl;
2590 oss2 <<
"( (1,1) block only is indicated by '1', (2,2) block only by '2', and both blocks by 'B' and none by '.')" << std::endl;
2592 int numProcs = comm->getSize();
2594 RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
2595 TEUCHOS_TEST_FOR_EXCEPTION(tmpic == Teuchos::null,
Exceptions::RuntimeError,
"Cannot cast base Teuchos::Comm to Teuchos::MpiComm object.");
2596 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
2602 if (!A22_.is_null())
2604 std::vector<char> states(numProcs, 0);
2606 MPI_Gather(&status, 1, MPI_CHAR, &states[0], 1, MPI_CHAR, 0, *rawMpiComm);
2608 states.push_back(status);
2611 int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
2612 for (
int proc = 0; proc < numProcs; proc += rowWidth) {
2613 for (
int j = 0; j < rowWidth; j++)
2614 if (proc + j < numProcs)
2615 if (states[proc+j] == 0)
2617 else if (states[proc+j] == 1)
2619 else if (states[proc+j] == 2)
2626 oss2 <<
" " << proc <<
":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
2637 #define MUELU_REFMAXWELL_SHORT 2638 #endif //ifdef MUELU_REFMAXWELL_DEF_HPP Important warning messages (one line)
Generic Smoother Factory for generating the smoothers of the MG hierarchy.
This class specifies the default factory that should generate some data on a Level if the data does n...
#define MueLu_sumAll(rcpComm, in, out)
static void SetMueLuOFileStream(const std::string &filename)
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
MueLu::DefaultLocalOrdinal LocalOrdinal
Factory for determing the number of partitions for rebalancing.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
Factory for generating coarse level map. Used by TentativePFactory.
static RCP< const Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2TpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > const vec)
Helper utility to pull out the underlying Tpetra objects from an Xpetra object.
void solveH(const MultiVector &RHS, MultiVector &X) const
apply solve to 1-1 block only
static void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows)
#define MueLu_maxAll(rcpComm, in, out)
static void SetDefaultVerbLevel(const VerbLevel defaultVerbLevel)
Set the default (global) verbosity level.
KOKKOS_INLINE_FUNCTION constexpr std::enable_if< std::is_integral< iType >::value, size_t >::type extent(const iType &r) const noexcept
void FindNonZeros(const Teuchos::ArrayRCP< const Scalar > vals, Teuchos::ArrayRCP< bool > nonzeros)
Class that encapsulates external library smoothers.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
Factory for building permutation matrix that can be be used to shuffle data (matrices, vectors) among processes.
void ZeroDirichletCols(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletCols, Scalar replaceWith)
void ZeroDirichletRows(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows, Scalar replaceWith)
One-liner description of what is happening.
std::string tolower(const std::string &str)
Namespace for MueLu classes and methods.
void SetPreviousLevel(const RCP< Level > &previousLevel)
Interface to Zoltan library.This interface provides access to partitioning methods in Zoltan...
void ApplyOAZToMatrixRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletRows)
static VerbLevel GetDefaultVerbLevel()
Get the default (global) verbosity level.
void SetFactoryManager(const RCP< const FactoryManagerBase > &factoryManager)
Set default factories (used internally by Hierarchy::SetLevel()).
void setParameters(Teuchos::ParameterList &list)
Set parameters.
void applyInverseAdditive(const MultiVector &RHS, MultiVector &X) const
apply additive algorithm for 2x2 solve
void buildProlongator()
Setup the prolongator for the (1,1)-block.
static RCP< Tpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > MV2NonConstTpetraMV(RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > vec)
#define MueLu_minAll(rcpComm, in, out)
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
Factory for building tentative prolongator.
MsgType toVerbLevel(const std::string &verbLevelStr)
void allocateMemory(int numVectors) const
allocate multivectors for solve
Print even more statistics.
void setlib(Xpetra::UnderlyingLib lib2)
Print statistics that do not involve significant additional computation.
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Residual(const Xpetra::Operator< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &X, const Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > &RHS)
void applyInverse121(const MultiVector &RHS, MultiVector &X) const
apply 1-2-1 algorithm for 2x2 solve
Kokkos::View< bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol, const bool count_twos_as_dirichlet)
void dumpCoords(const RealValuedMultiVector &X, std::string name) const
dump out real-valued multivector
MueLu::DefaultScalar Scalar
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > RealValuedToScalarMultiVector(RCP< Xpetra::MultiVector< typename Teuchos::ScalarTraits< Scalar >::magnitudeType, LocalOrdinal, GlobalOrdinal, Node > > X)
Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
Factory to export aggregation info or visualize aggregates using VTK.
Interface to Zoltan2 library.This interface provides access to partitioning methods in Zoltan2...
static void SetMueLuOStream(const Teuchos::RCP< Teuchos::FancyOStream > &mueluOStream)
static RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Transpose(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, bool optimizeTranspose=false, const std::string &label=std::string(), const Teuchos::RCP< Teuchos::ParameterList > ¶ms=Teuchos::null)
Xpetra::MultiVector< coordinateType, LO, GO, NO > RealValuedMultiVector
AmalgamationFactory for subblocks of strided map based amalgamation data.
void AddKeepFlag(const std::string &ename, const FactoryBase *factory=NoFactory::get(), KeepType keep=MueLu::Keep)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
void compute(bool reuse=false)
Setup the preconditioner.
Applies permutation to grid transfer operators.
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Print all timing information.
void initialize(const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &Ms_Matrix, const Teuchos::RCP< Matrix > &M0inv_Matrix, const Teuchos::RCP< Matrix > &M1_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List)
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
void applyInverse212(const MultiVector &RHS, MultiVector &X) const
apply 2-1-2 algorithm for 2x2 solve
Factory for creating a graph base on a given matrix.
Teuchos::RCP< Teuchos::TimeMonitor > getTimer(std::string name, RCP< const Teuchos::Comm< int > > comm=Teuchos::null) const
get a (synced) timer
void parallel_for(const ExecPolicy &policy, const FunctorType &functor, const std::string &str="", typename std::enable_if< Kokkos::Impl::is_execution_policy< ExecPolicy >::value >::type *=nullptr)
void SetLevelID(int levelID)
Set level number.
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList)
Helper function to create a MueLu preconditioner that can be used by Xpetra.Given an Xpetra::Matrix...
Class for transferring coordinates from a finer level to a coarser one.
void dump(const Matrix &A, std::string name) const
dump out matrix
Factory for building coarse matrices.
void formCoarseMatrix()
Compute P11^{T}*A*P11 efficiently.
Exception throws to report errors in the internal logical of the program.
void apply(const MultiVector &X, MultiVector &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, Scalar alpha=Teuchos::ScalarTraits< Scalar >::one(), Scalar beta=Teuchos::ScalarTraits< Scalar >::zero()) const
void deep_copy(const View< DT, DP... > &dst, typename ViewTraits< DT, DP... >::const_value_type &value, typename std::enable_if< std::is_same< typename ViewTraits< DT, DP... >::specialize, void >::value >::type *=nullptr)
Description of what is happening (more verbose)
static void ZeroDirichletCols(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Teuchos::ArrayRCP< const bool > &dirichletCols, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
Factory for building coarse matrices.
Teuchos::ScalarTraits< Scalar >::coordinateType coordinateType
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Factory for building Smoothed Aggregation prolongators.
Factory for building uncoupled aggregates.
Kokkos::View< bool *, typename NO::device_type > DetectDirichletCols(const Xpetra::Matrix< SC, LO, GO, NO > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows)
static std::string translate(Teuchos::ParameterList ¶mList, const std::string &defaultVals="")
: Translate ML parameters to MueLu parameter XML string
void solve22(const MultiVector &RHS, MultiVector &X) const
apply solve to 2-2 block only
Factory for building a thresholded operator.
Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
static const RCP< const NoFactory > getRCP()
Static Get() functions.