45 #ifndef MUELU_MAXWELL1_DEF_HPP 46 #define MUELU_MAXWELL1_DEF_HPP 52 #include "Xpetra_Map.hpp" 53 #include "Xpetra_MatrixMatrix.hpp" 54 #include "Xpetra_TripleMatrixMultiply.hpp" 55 #include "Xpetra_CrsMatrixUtils.hpp" 56 #include "Xpetra_MatrixUtils.hpp" 59 #include "MueLu_Maxwell_Utils.hpp" 61 #include "MueLu_ReitzingerPFactory.hpp" 62 #include "MueLu_SaPFactory.hpp" 63 #include "MueLu_AggregationExportFactory.hpp" 64 #include "MueLu_Utilities.hpp" 66 #include "MueLu_Hierarchy.hpp" 67 #include "MueLu_RAPFactory.hpp" 68 #include "MueLu_PerfUtils.hpp" 69 #include "MueLu_ParameterListInterpreter.hpp" 72 #if defined(HAVE_MUELU_KOKKOS_REFACTOR) 73 #include "MueLu_Utilities_kokkos.hpp" 82 #ifdef HAVE_MUELU_CUDA 83 #include "cuda_profiler_api.h" 90 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 return SM_Matrix_->getDomainMap();
96 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
98 return SM_Matrix_->getRangeMap();
102 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
105 std::string mode_string = list.get(
"maxwell1: mode", MasterList::getDefault<std::string>(
"maxwell1: mode"));
106 applyBCsTo22_ = list.get(
"maxwell1: apply BCs to 22",
true);
109 Teuchos::ParameterList defaultSmootherList;
110 defaultSmootherList.set(
"smoother: type",
"CHEBYSHEV");
111 defaultSmootherList.sublist(
"smoother: params").set(
"chebyshev: degree",2);
112 defaultSmootherList.sublist(
"smoother: params").set(
"chebyshev: ratio eigenvalue",7.0);
113 defaultSmootherList.sublist(
"smoother: params").set(
"chebyshev: eigenvalue max iterations",30);
116 std::string verbosity = list.get(
"verbosity",
"high");
119 if(mode_string ==
"standard") mode_ = MODE_STANDARD;
120 else if(mode_string ==
"refmaxwell") mode_ = MODE_REFMAXWELL;
121 else if(mode_string ==
"edge only") mode_ = MODE_EDGE_ONLY;
123 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"Must use mode 'standard', 'refmaxwell' or 'edge only'.");
128 if(list.isSublist(
"maxwell1: 22list"))
129 precList22_ = list.sublist(
"maxwell1: 22list");
130 else if(list.isSublist(
"refmaxwell: 22list"))
131 precList22_ = list.sublist(
"refmaxwell: 22list");
132 if(mode_ == MODE_EDGE_ONLY || mode_ == MODE_STANDARD)
133 precList22_.set(
"smoother: pre or post",
"none");
134 else if(!precList22_.isType<std::string>(
"Preconditioner Type") &&
135 !precList22_.isType<std::string>(
"smoother: type") &&
136 !precList22_.isType<std::string>(
"smoother: pre type") &&
137 !precList22_.isType<std::string>(
"smoother: post type")) {
138 precList22_ = defaultSmootherList;
140 precList22_.set(
"verbosity",precList22_.get(
"verbosity",verbosity));
145 if(list.isSublist(
"maxwell1: 11list"))
146 precList11_ = list.sublist(
"maxwell1: 11list");
147 else if(list.isSublist(
"refmaxwell: 11list"))
148 precList11_ = list.sublist(
"refmaxwell: 11list");
149 if(!precList11_.isType<std::string>(
"Preconditioner Type") &&
150 !precList11_.isType<std::string>(
"smoother: type") &&
151 !precList11_.isType<std::string>(
"smoother: pre type") &&
152 !precList11_.isType<std::string>(
"smoother: post type")) {
153 if(mode_ == MODE_EDGE_ONLY || mode_ == MODE_REFMAXWELL) {
154 precList11_ = defaultSmootherList;
157 precList11_.set(
"smoother: type",
"HIPTMAIR");
158 precList11_.sublist(
"hiptmair: smoother type 1",
"CHEBYSHEV");
159 precList11_.sublist(
"hiptmair: smoother type 2",
"CHEBYSHEV");
160 precList11_.sublist(
"hiptmair: smoother list 1") = defaultSmootherList;
161 precList11_.sublist(
"hiptmair: smoother list 2") = defaultSmootherList;
164 precList11_.set(
"verbosity",precList11_.get(
"verbosity",verbosity));
168 !precList11_.isType<std::string>(
"Preconditioner Type") &&
169 !precList11_.isParameter(
"reuse: type"))
170 precList11_.set(
"reuse: type",
"full");
172 !precList22_.isType<std::string>(
"Preconditioner Type") &&
173 !precList22_.isParameter(
"reuse: type"))
174 precList22_.set(
"reuse: type",
"full");
178 #if !defined(HAVE_MUELU_KOKKOS_REFACTOR) 181 # ifdef HAVE_MUELU_SERIAL 182 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosSerialWrapperNode).name())
185 # ifdef HAVE_MUELU_OPENMP 186 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosOpenMPWrapperNode).name())
189 # ifdef HAVE_MUELU_CUDA 190 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosCudaWrapperNode).name())
193 useKokkos_ = list.get(
"use kokkos refactor",useKokkos_);
199 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
202 #ifdef HAVE_MUELU_CUDA 203 if (parameterList_.get<
bool>(
"maxwell1: cuda profile setup",
false)) cudaProfilerStart();
206 std::string timerLabel;
208 timerLabel =
"MueLu Maxwell1: compute (reuse)";
210 timerLabel =
"MueLu Maxwell1: compute";
211 RCP<Teuchos::TimeMonitor> tmCompute = getTimer(timerLabel);
219 RCP<ParameterList> params = rcp(
new ParameterList());;
220 params->set(
"printLoadBalancingInfo",
true);
221 params->set(
"printCommInfo",
true);
228 magnitudeType rowSumTol = precList11_.get(
"aggregation: row sum drop tol",-1.0);
231 useKokkos_,BCrowsKokkos_,BCcolsKokkos_,BCdomainKokkos_,
233 BCedges_,BCnodes_,BCrows_,BCcols_,BCdomain_,
234 allEdgesBoundary_,allNodesBoundary_);
236 GetOStream(
Statistics2) <<
"MueLu::Maxwell1::compute(): Detected " << BCedges_ <<
" BC rows and " << BCnodes_ <<
" BC columns." << std::endl;
240 if (allEdgesBoundary_) {
243 GetOStream(
Warnings0) <<
"All edges are detected as boundary edges!" << std::endl;
244 mode_ = MODE_EDGE_ONLY;
247 if(precList11_.get<std::string>(
"smoother: type") ==
"HIPTMAIR") {
248 precList11_.set(
"smoother: type",precList11_.get(
"hiptmair: smoother type 1",
"CHEBYSHEV"));
249 precList11_.sublist(
"smoother: sublist") = precList11_.sublist(
"hiptmair: smoother list 1");
253 throw std::runtime_error(
"Maxwell1: Not yet supported");
258 if (allNodesBoundary_) {
261 GetOStream(
Warnings0) <<
"All nodes are detected as boundary edges!" << std::endl;
262 mode_ = MODE_EDGE_ONLY;
264 if(precList11_.get<std::string>(
"smoother: type") ==
"HIPTMAIR") {
265 precList11_.set(
"smoother: type",precList11_.get(
"hiptmair: smoother type 1",
"CHEBYSHEV"));
266 precList11_.sublist(
"smoother: sublist") = precList11_.sublist(
"hiptmair: smoother list 1");
270 throw std::runtime_error(
"Maxwell1: Not yet supported");
278 if(!reuse && applyBCsTo22_) {
279 GetOStream(
Runtime0) <<
"Maxwell1::compute(): nuking BC nodes of D0" << std::endl;
280 D0_Matrix_->resumeFill();
282 if (D0_Matrix_->getRowMap()->lib() == Xpetra::UseEpetra)
283 replaceWith= Teuchos::ScalarTraits<SC>::eps();
285 replaceWith = Teuchos::ScalarTraits<SC>::zero();
286 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 295 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
300 if(Kn_Matrix_.is_null()) {
301 Kn_Matrix_ = generate_kn();
313 for(
int i=0; i<Hierarchy22_->GetNumLevels(); i++) {
314 Hierarchy11_->AddNewLevel();
315 RCP<Level> NodeL = Hierarchy22_->GetLevel(i);
316 RCP<Level> EdgeL = Hierarchy11_->GetLevel(i);
317 EdgeL->Set(
"NodeMatrix",NodeL->Get<RCP<Matrix> >(
"A"));
319 EdgeL->Set(
"A", SM_Matrix_);
320 EdgeL->Set(
"D0", D0_Matrix_);
323 EdgeL->Set(
"Pnodal",NodeL->Get<RCP<Matrix> >(
"P"));
331 SM_Matrix_->setObjectLabel(
"A(1,1)");
332 precList11_.set(
"coarse: max size",1);
333 precList11_.set(
"max levels",Hierarchy22_->GetNumLevels());
336 Hierarchy11_->setlib(SM_Matrix_->getDomainMap()->lib());
337 Hierarchy11_->SetProcRankVerbose(SM_Matrix_->getDomainMap()->getComm()->getRank());
338 mueLuFactory->SetupHierarchy(*Hierarchy11_);
340 if(mode_ == MODE_REFMAXWELL) {
341 if(Hierarchy11_->GetNumLevels() > 1) {
342 RCP<Level> EdgeL = Hierarchy11_->GetLevel(1);
343 P11_ = EdgeL->Get<RCP<Matrix> >(
"P");
355 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
361 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu Maxwell1: Build Kn");
363 Level fineLevel, coarseLevel;
369 fineLevel.
Set(
"A",SM_Matrix_);
370 coarseLevel.
Set(
"P",D0_Matrix_);
373 coarseLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
374 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
375 coarseLevel.setObjectLabel(
"Maxwell1 (2,2)");
376 fineLevel.setObjectLabel(
"Maxwell1 (2,2)");
378 RCP<RAPFactory> rapFact = rcp(
new RAPFactory());
379 ParameterList rapList = *(rapFact->GetValidParameterList());
380 rapList.set(
"transpose: use implicit",
true);
381 rapList.set(
"rap: fix zero diagonals", parameterList_.get<
bool>(
"rap: fix zero diagonals",
true));
382 rapList.set(
"rap: fix zero diagonals threshold", parameterList_.get<
double>(
"rap: fix zero diagonals threshold", Teuchos::ScalarTraits<double>::eps()));
383 rapList.set(
"rap: triple product", parameterList_.get<
bool>(
"rap: triple product",
false));
384 rapFact->SetParameterList(rapList);
385 coarseLevel.
Request(
"A", rapFact.get());
387 coarseLevel.
Request(
"AP reuse data", rapFact.get());
388 coarseLevel.
Request(
"RAP reuse data", rapFact.get());
391 RCP<Matrix> Kn_Matrix = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
392 Kn_Matrix->setObjectLabel(
"A(2,2)");
399 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
401 if(mode_ == MODE_REFMAXWELL) {
402 RCP<Teuchos::TimeMonitor> tmAlloc = getTimer(
"MueLu Maxwell1: Allocate MVs");
404 residualFine_ = MultiVectorFactory::Build(SM_Matrix_->getRangeMap(), numVectors);
406 if(!Hierarchy11_.is_null() && Hierarchy11_->GetNumLevels() > 1) {
407 RCP<Level> EdgeL = Hierarchy11_->GetLevel(1);
408 RCP<Matrix> A = EdgeL->Get<RCP<Matrix> >(
"A");
409 residual11c_ = MultiVectorFactory::Build(A->getRangeMap(), numVectors);
410 update11c_ = MultiVectorFactory::Build(A->getDomainMap(), numVectors);
413 if(!Hierarchy22_.is_null()) {
414 residual22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
415 update22_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
423 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
425 if (dump_matrices_) {
426 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
427 Xpetra::IO<SC, LO, GO, NO>::Write(name, A);
432 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
434 if (dump_matrices_) {
435 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
436 Xpetra::IO<SC, LO, GO, NO>::Write(name, X);
441 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
443 if (dump_matrices_) {
444 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
445 Xpetra::IO<coordinateType, LO, GO, NO>::Write(name, X);
450 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
452 if (dump_matrices_) {
453 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
454 std::ofstream out(name);
455 for (
size_t i = 0; i < Teuchos::as<size_t>(v.size()); i++)
460 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 461 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
463 if (dump_matrices_) {
464 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
465 std::ofstream out(name);
466 auto vH = Kokkos::create_mirror_view (v);
467 Kokkos::deep_copy(vH , v);
468 for (
size_t i = 0; i < vH.size(); i++)
469 out << vH[i] <<
"\n";
474 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
478 return Teuchos::rcp(
new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name)));
481 return Teuchos::rcp(
new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), SM_Matrix_->getRowMap()->getComm().ptr()));
483 return Teuchos::rcp(
new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), comm.ptr()));
486 return Teuchos::null;
492 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
494 bool reuse = !SM_Matrix_.is_null();
495 SM_Matrix_ = SM_Matrix_new;
496 dump(*SM_Matrix_,
"SM.m");
502 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
505 const SC one = Teuchos::ScalarTraits<SC>::one();
506 if (!allEdgesBoundary_ && X.getNumVectors() != residualFine_->getNumVectors())
507 allocateMemory(X.getNumVectors());
509 TEUCHOS_TEST_FOR_EXCEPTION(Hierarchy11_.is_null() || Hierarchy11_->GetNumLevels() == 0,
Exceptions::RuntimeError,
"(1,1) Hiearchy is null.");
512 RCP<Level> Fine = Hierarchy11_->GetLevel(0);
513 if (Fine->IsAvailable(
"PreSmoother")) {
514 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: PreSmoother");
515 RCP<SmootherBase> preSmoo = Fine->Get< RCP<SmootherBase> >(
"PreSmoother");
516 preSmoo->Apply(X, RHS,
true);
521 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: residual calculation");
526 if(!P11_.is_null()) {
527 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: (1,1) correction");
528 P11_->apply(*residualFine_,*residual11c_,Teuchos::TRANS);
529 Hierarchy11_->Iterate(*residual11c_,*update11c_,
true,1);
533 if (!allNodesBoundary_) {
534 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: (2,2) correction");
535 D0_Matrix_->apply(*residualFine_,*residual22_,Teuchos::TRANS);
536 Hierarchy22_->Iterate(*residual22_,*update22_,
true,0);
541 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: Orolongation");
543 P11_->apply(*update11c_,X,Teuchos::NO_TRANS,one,one);
544 if (!allNodesBoundary_)
545 D0_Matrix_->apply(*update22_,X,Teuchos::NO_TRANS,one,one);
549 if (Fine->IsAvailable(
"PostSmoother")) {
550 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu Maxwell1: PostSmoother");
551 RCP<SmootherBase> postSmoo = Fine->Get< RCP<SmootherBase> >(
"PostSmoother");
552 postSmoo->Apply(X, RHS,
false);
558 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
560 Hierarchy11_->Iterate(RHS,X);
563 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
568 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu Maxwell1: solve");
569 if(mode_ == MODE_STANDARD || mode_ == MODE_EDGE_ONLY)
570 applyInverseStandard(RHS,X);
571 else if(mode_ == MODE_REFMAXWELL)
572 applyInverseRefMaxwellAdditive(RHS,X);
574 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"Must use mode 'standard', 'refmaxwell' or 'edge only'.");
579 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
585 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
588 const Teuchos::RCP<Matrix> & Kn_Matrix,
589 const Teuchos::RCP<MultiVector> & Nullspace,
590 const Teuchos::RCP<RealValuedMultiVector> & Coords,
591 Teuchos::ParameterList& List)
594 TEUCHOS_ASSERT(D0_Matrix!=Teuchos::null);
596 #ifdef HAVE_MUELU_DEBUG 597 if(!Kn_Matrix.is_null()) {
598 TEUCHOS_ASSERT(Kn_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getDomainMap()));
599 TEUCHOS_ASSERT(Kn_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getDomainMap()));
603 TEUCHOS_ASSERT(D0_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getRowMap()));
606 Hierarchy11_ = Teuchos::null;
607 Hierarchy22_ = Teuchos::null;
608 mode_ = MODE_STANDARD;
612 allEdgesBoundary_=
false;
613 allNodesBoundary_=
false;
614 dump_matrices_ =
false;
617 applyBCsTo22_ =
true;
624 if (D0_Matrix->getRowMap()->lib() == Xpetra::UseTpetra) {
629 RCP<Matrix> D0copy = MatrixFactory::Build(D0_Matrix->getRowMap(), D0_Matrix->getColMap(), 0);
630 RCP<CrsMatrix> D0copyCrs = rcp_dynamic_cast<CrsMatrixWrap>(D0copy,
true)->getCrsMatrix();
631 ArrayRCP<const size_t> D0rowptr_RCP;
632 ArrayRCP<const LO> D0colind_RCP;
633 ArrayRCP<const SC> D0vals_RCP;
634 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getAllValues(D0rowptr_RCP,
638 ArrayRCP<size_t> D0copyrowptr_RCP;
639 ArrayRCP<LO> D0copycolind_RCP;
640 ArrayRCP<SC> D0copyvals_RCP;
641 D0copyCrs->allocateAllValues(D0vals_RCP.size(),D0copyrowptr_RCP,D0copycolind_RCP,D0copyvals_RCP);
642 D0copyrowptr_RCP.deepCopy(D0rowptr_RCP());
643 D0copycolind_RCP.deepCopy(D0colind_RCP());
644 D0copyvals_RCP.deepCopy(D0vals_RCP());
645 D0copyCrs->setAllValues(D0copyrowptr_RCP,
648 D0copyCrs->expertStaticFillComplete(D0_Matrix->getDomainMap(), D0_Matrix->getRangeMap(),
649 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getImporter(),
650 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getExporter());
653 D0_Matrix_ = MatrixFactory::BuildCopy(D0_Matrix);
656 Kn_Matrix_ = Kn_Matrix;
658 Nullspace_ = Nullspace;
660 if(!Kn_Matrix_.is_null()) dump(*Kn_Matrix_,
"Kn.m");
661 if (!Nullspace_.is_null()) dump(*Nullspace_,
"nullspace.m");
662 if (!Coords_.is_null()) dumpCoords(*Coords_,
"coords.m");
668 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
670 describe(Teuchos::FancyOStream& out,
const Teuchos::EVerbosityLevel )
const {
672 std::ostringstream oss;
674 RCP<const Teuchos::Comm<int> > comm = SM_Matrix_->getDomainMap()->getComm();
678 if (!Kn_Matrix_.is_null())
679 root = comm->getRank();
684 reduceAll(*comm, Teuchos::REDUCE_MAX, root, Teuchos::ptr(&actualRoot));
689 oss <<
"\n--------------------------------------------------------------------------------\n" <<
690 "--- Maxwell1 Summary ---\n" 691 "--------------------------------------------------------------------------------" << std::endl;
697 SM_Matrix_->getRowMap()->getComm()->barrier();
699 numRows = SM_Matrix_->getGlobalNumRows();
700 nnz = SM_Matrix_->getGlobalNumEntries();
702 Xpetra::global_size_t tt = numRows;
703 int rowspacer = 3;
while (tt != 0) { tt /= 10; rowspacer++; }
705 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
707 oss <<
"block " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " << std::setw(9) <<
" nnz/row" << std::endl;
708 oss <<
"(1, 1)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
710 if (!Kn_Matrix_.is_null()) {
712 numRows = Kn_Matrix_->getGlobalNumRows();
713 nnz = Kn_Matrix_->getGlobalNumEntries();
715 oss <<
"(2, 2)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
721 std::string outstr = oss.str();
724 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
725 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
727 int strLength = outstr.size();
728 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
729 if (comm->getRank() != root)
730 outstr.resize(strLength);
731 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
736 if (!Hierarchy11_.is_null())
737 Hierarchy11_->describe(out, GetVerbLevel());
739 if (!Hierarchy22_.is_null())
740 Hierarchy22_->describe(out, GetVerbLevel());
744 std::ostringstream oss2;
746 oss2 <<
"Sub-solver distribution over ranks" << std::endl;
747 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;
749 int numProcs = comm->getSize();
751 RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
753 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
757 if (!Kn_Matrix_.is_null())
759 std::vector<char> states(numProcs, 0);
761 MPI_Gather(&status, 1, MPI_CHAR, &states[0], 1, MPI_CHAR, 0, *rawMpiComm);
763 states.push_back(status);
766 int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
767 for (
int proc = 0; proc < numProcs; proc += rowWidth) {
768 for (
int j = 0; j < rowWidth; j++)
769 if (proc + j < numProcs)
770 if (states[proc+j] == 0)
772 else if (states[proc+j] == 1)
774 else if (states[proc+j] == 2)
781 oss2 <<
" " << proc <<
":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
794 #define MUELU_MAXWELL1_SHORT 795 #endif //ifdef MUELU_MAXWELL1_DEF_HPP Important warning messages (one line)
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
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.
#define HAVE_MUELU_KOKKOS_REFACTOR
static void removeExplicitZeros(Teuchos::ParameterList ¶meterList, RCP< Matrix > &D0_Matrix, RCP< Matrix > &SM_Matrix, RCP< Matrix > &M1_Matrix, RCP< Matrix > &Ms_Matrix)
Remove explicit zeros.
Teuchos::RCP< const Map > getDomainMap() const
Returns the Xpetra::Map object associated with the domain of this operator.
void applyInverseStandard(const MultiVector &RHS, MultiVector &X) const
apply standard Maxwell1 cycle
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new, bool ComputePrec=true)
Reset system matrix.
void applyInverseRefMaxwellAdditive(const MultiVector &RHS, MultiVector &X) const
apply RefMaxwell additive 2x2 style cycle
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
void ZeroDirichletCols(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletCols, Scalar replaceWith)
void initialize(const Teuchos::RCP< Matrix > &D0_Matrix, const Teuchos::RCP< Matrix > &Kn_Matrix, const Teuchos::RCP< MultiVector > &Nullspace, const Teuchos::RCP< RealValuedMultiVector > &Coords, Teuchos::ParameterList &List)
One-liner description of what is happening.
Namespace for MueLu classes and methods.
void SetPreviousLevel(const RCP< Level > &previousLevel)
void SetFactoryManager(const RCP< const FactoryManagerBase > &factoryManager)
Set default factories (used internally by Hierarchy::SetLevel()).
Print even more statistics.
void setlib(Xpetra::UnderlyingLib lib2)
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
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 setParameters(Teuchos::ParameterList &list)
Set parameters.
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
Class that holds all level-specific information.
Teuchos::RCP< Teuchos::TimeMonitor > getTimer(std::string name, RCP< const Teuchos::Comm< int > > comm=Teuchos::null) const
get a (synced) timer
static void detectBoundaryConditionsSM(RCP< Matrix > &SM_Matrix, RCP< Matrix > &D0_Matrix, magnitudeType rowSumTol, int &BCedges, int &BCnodes, Teuchos::ArrayRCP< bool > &BCrows, Teuchos::ArrayRCP< bool > &BCcols, Teuchos::ArrayRCP< bool > &BCdomain, bool &allEdgesBoundary, bool &allNodesBoundary)
Detect Dirichlet boundary conditions.
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
Print all timing information.
void dump(const Matrix &A, std::string name) const
dump out matrix
Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
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...
void compute(bool reuse=false)
Setup the preconditioner.
Teuchos::RCP< Matrix > generate_kn() const
Generates the Kn matrix.
Xpetra::MultiVector< coordinateType, LO, GO, NO > RealValuedMultiVector
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
Exception throws to report errors in the internal logical of the program.
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.
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
void dumpCoords(const RealValuedMultiVector &X, std::string name) const
dump out real-valued multivector
Provides methods to build a multigrid hierarchy and apply multigrid cycles.
void Request(const FactoryBase &factory)
Increment the storage counter for all the inputs of a factory.
void allocateMemory(int numVectors) const
allocate multivectors for solve