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" 75 #include "MueLu_Maxwell_Utils.hpp" 77 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 78 #include "MueLu_AmalgamationFactory_kokkos.hpp" 79 #include "MueLu_CoalesceDropFactory_kokkos.hpp" 80 #include "MueLu_CoarseMapFactory_kokkos.hpp" 81 #include "MueLu_CoordinatesTransferFactory_kokkos.hpp" 82 #include "MueLu_UncoupledAggregationFactory_kokkos.hpp" 83 #include "MueLu_TentativePFactory_kokkos.hpp" 84 #include "MueLu_SaPFactory_kokkos.hpp" 85 #include "MueLu_Utilities_kokkos.hpp" 86 #include <Kokkos_Core.hpp> 87 #include <KokkosSparse_CrsMatrix.hpp> 90 #include "MueLu_ZoltanInterface.hpp" 91 #include "MueLu_Zoltan2Interface.hpp" 92 #include "MueLu_RepartitionHeuristicFactory.hpp" 93 #include "MueLu_RepartitionFactory.hpp" 94 #include "MueLu_RebalanceAcFactory.hpp" 95 #include "MueLu_RebalanceTransferFactory.hpp" 102 #ifdef HAVE_MUELU_CUDA 103 #include "cuda_profiler_api.h" 107 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 109 #include <Thyra_VectorBase.hpp> 110 #include <Thyra_SolveSupportTypes.hpp> 112 #include <Stratimikos_DefaultLinearSolverBuilder.hpp> 114 #include "Teuchos_AbstractFactoryStd.hpp" 116 #ifdef HAVE_MUELU_IFPACK2 117 #include <Thyra_Ifpack2PreconditionerFactory.hpp> 124 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
126 return SM_Matrix_->getDomainMap();
130 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
132 return SM_Matrix_->getRangeMap();
136 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
139 if (list.isType<std::string>(
"parameterlist: syntax") && list.get<std::string>(
"parameterlist: syntax") ==
"ml") {
141 if(list.isSublist(
"refmaxwell: 11list") && list.sublist(
"refmaxwell: 11list").isSublist(
"edge matrix free: coarse"))
143 if(list.isSublist(
"refmaxwell: 22list"))
148 parameterList_ = list;
149 std::string verbosityLevel = parameterList_.get<std::string>(
"verbosity", MasterList::getDefault<std::string>(
"verbosity"));
151 std::string outputFilename = parameterList_.get<std::string>(
"output filename", MasterList::getDefault<std::string>(
"output filename"));
152 if (outputFilename !=
"")
154 if (parameterList_.isType<Teuchos::RCP<Teuchos::FancyOStream> >(
"output stream"))
157 if (parameterList_.get(
"print initial parameters",MasterList::getDefault<bool>(
"print initial parameters")))
158 GetOStream(static_cast<MsgType>(
Runtime1), 0) << parameterList_ << std::endl;
159 disable_addon_ = list.get(
"refmaxwell: disable addon", MasterList::getDefault<bool>(
"refmaxwell: disable addon"));
160 mode_ = list.get(
"refmaxwell: mode", MasterList::getDefault<std::string>(
"refmaxwell: mode"));
161 use_as_preconditioner_ = list.get(
"refmaxwell: use as preconditioner", MasterList::getDefault<bool>(
"refmaxwell: use as preconditioner"));
162 dump_matrices_ = list.get(
"refmaxwell: dump matrices", MasterList::getDefault<bool>(
"refmaxwell: dump matrices"));
163 enable_reuse_ = list.get(
"refmaxwell: enable reuse", MasterList::getDefault<bool>(
"refmaxwell: enable reuse"));
164 implicitTranspose_ = list.get(
"transpose: use implicit", MasterList::getDefault<bool>(
"transpose: use implicit"));
165 fuseProlongationAndUpdate_ = list.get(
"fuse prolongation and update", MasterList::getDefault<bool>(
"fuse prolongation and update"));
166 skipFirstLevel_ = list.get(
"refmaxwell: skip first (1,1) level", MasterList::getDefault<bool>(
"refmaxwell: skip first (1,1) level"));
167 syncTimers_ = list.get(
"sync timers",
false);
168 numItersH_ = list.get(
"refmaxwell: num iters H", 1);
169 numIters22_ = list.get(
"refmaxwell: num iters 22", 1);
170 applyBCsToAnodal_ = list.get(
"refmaxwell: apply BCs to Anodal",
false);
171 applyBCsToH_ = list.get(
"refmaxwell: apply BCs to H",
true);
172 applyBCsTo22_ = list.get(
"refmaxwell: apply BCs to 22",
true);
174 if(list.isSublist(
"refmaxwell: 11list"))
175 precList11_ = list.sublist(
"refmaxwell: 11list");
176 if(!precList11_.isType<std::string>(
"Preconditioner Type") &&
177 !precList11_.isType<std::string>(
"smoother: type") &&
178 !precList11_.isType<std::string>(
"smoother: pre type") &&
179 !precList11_.isType<std::string>(
"smoother: post type")) {
180 precList11_.set(
"smoother: type",
"CHEBYSHEV");
181 precList11_.sublist(
"smoother: params").set(
"chebyshev: degree",2);
182 precList11_.sublist(
"smoother: params").set(
"chebyshev: ratio eigenvalue",5.4);
183 precList11_.sublist(
"smoother: params").set(
"chebyshev: eigenvalue max iterations",30);
186 if(list.isSublist(
"refmaxwell: 22list"))
187 precList22_ = list.sublist(
"refmaxwell: 22list");
188 if(!precList22_.isType<std::string>(
"Preconditioner Type") &&
189 !precList22_.isType<std::string>(
"smoother: type") &&
190 !precList22_.isType<std::string>(
"smoother: pre type") &&
191 !precList22_.isType<std::string>(
"smoother: post type")) {
192 precList22_.set(
"smoother: type",
"CHEBYSHEV");
193 precList22_.sublist(
"smoother: params").set(
"chebyshev: degree",2);
194 precList22_.sublist(
"smoother: params").set(
"chebyshev: ratio eigenvalue",7.0);
195 precList22_.sublist(
"smoother: params").set(
"chebyshev: eigenvalue max iterations",30);
198 if(!list.isType<std::string>(
"smoother: type") && !list.isType<std::string>(
"smoother: pre type") && !list.isType<std::string>(
"smoother: post type")) {
199 list.set(
"smoother: type",
"CHEBYSHEV");
200 list.sublist(
"smoother: params").set(
"chebyshev: degree",2);
201 list.sublist(
"smoother: params").set(
"chebyshev: ratio eigenvalue",20.0);
202 list.sublist(
"smoother: params").set(
"chebyshev: eigenvalue max iterations",30);
205 if(list.isSublist(
"smoother: params")) {
206 smootherList_ = list.sublist(
"smoother: params");
210 !precList11_.isType<std::string>(
"Preconditioner Type") &&
211 !precList11_.isParameter(
"reuse: type"))
212 precList11_.set(
"reuse: type",
"full");
214 !precList22_.isType<std::string>(
"Preconditioner Type") &&
215 !precList22_.isParameter(
"reuse: type"))
216 precList22_.set(
"reuse: type",
"full");
218 #if !defined(HAVE_MUELU_KOKKOS_REFACTOR) 221 # ifdef HAVE_MUELU_SERIAL 222 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosSerialWrapperNode).name())
225 # ifdef HAVE_MUELU_OPENMP 226 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosOpenMPWrapperNode).name())
229 # ifdef HAVE_MUELU_CUDA 230 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosCudaWrapperNode).name())
233 # ifdef HAVE_MUELU_HIP 234 if (
typeid(
Node).name() ==
typeid(Kokkos::Compat::KokkosHIPWrapperNode).name())
237 useKokkos_ = list.get(
"use kokkos refactor",useKokkos_);
243 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
246 #ifdef HAVE_MUELU_CUDA 247 if (parameterList_.get<
bool>(
"refmaxwell: cuda profile setup",
false)) cudaProfilerStart();
250 std::string timerLabel;
252 timerLabel =
"MueLu RefMaxwell: compute (reuse)";
254 timerLabel =
"MueLu RefMaxwell: compute";
255 RCP<Teuchos::TimeMonitor> tmCompute = getTimer(timerLabel);
262 RCP<ParameterList> params = rcp(
new ParameterList());;
263 params->set(
"printLoadBalancingInfo",
true);
264 params->set(
"printCommInfo",
true);
271 magnitudeType rowSumTol = parameterList_.get(
"refmaxwell: row sum drop tol (1,1)",-1.0);
274 useKokkos_,BCrowsKokkos_,BCcolsKokkos_,BCdomainKokkos_,
276 BCedges_,BCnodes_,BCrows_,BCcols_,BCdomain_,
277 allEdgesBoundary_,allNodesBoundary_);
279 GetOStream(
Statistics2) <<
"MueLu::RefMaxwell::compute(): Detected " << BCedges_ <<
" BC rows and " << BCnodes_ <<
" BC columns." << std::endl;
283 if (allEdgesBoundary_) {
286 GetOStream(
Warnings0) <<
"All edges are detected as boundary edges!" << std::endl;
288 setFineLevelSmoother();
295 if(Nullspace_ != null) {
297 TEUCHOS_ASSERT(Nullspace_->getMap()->isCompatible(*(SM_Matrix_->getRowMap())));
299 else if(Nullspace_ == null && Coords_ != null) {
300 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 301 RCP<MultiVector> CoordsSC;
303 CoordsSC = Utilities_kokkos::RealValuedToScalarMultiVector(Coords_);
309 Nullspace_ = MultiVectorFactory::Build(SM_Matrix_->getRowMap(),Coords_->getNumVectors());
310 D0_Matrix_->apply(*CoordsSC,*Nullspace_);
312 bool normalize = parameterList_.get<
bool>(
"refmaxwell: normalize nullspace", MasterList::getDefault<bool>(
"refmaxwell: normalize nullspace"));
317 ArrayRCP<ArrayRCP<const Scalar> > localNullspace(Nullspace_->getNumVectors());
318 for (
size_t i = 0; i < Nullspace_->getNumVectors(); i++)
319 localNullspace[i] = Nullspace_->getData(i);
320 coordinateType localMinLen = Teuchos::ScalarTraits<coordinateType>::rmax();
321 coordinateType localMeanLen = Teuchos::ScalarTraits<coordinateType>::zero();
322 coordinateType localMaxLen = Teuchos::ScalarTraits<coordinateType>::zero();
323 for (
size_t j=0; j < Nullspace_->getMap()->getNodeNumElements(); j++) {
324 Scalar lenSC = Teuchos::ScalarTraits<Scalar>::zero();
325 for (
size_t i=0; i < Nullspace_->getNumVectors(); i++)
326 lenSC += localNullspace[i][j]*localNullspace[i][j];
327 coordinateType len = sqrt(Teuchos::ScalarTraits<Scalar>::real(lenSC));
328 localMinLen = std::min(localMinLen, len);
329 localMaxLen = std::max(localMaxLen, len);
333 RCP<const Teuchos::Comm<int> > comm = Nullspace_->getMap()->getComm();
337 meanLen /= Nullspace_->getMap()->getGlobalNumElements();
341 GetOStream(
Statistics0) <<
"Edge length (min/mean/max): " << minLen <<
" / " << meanLen <<
" / " << maxLen << std::endl;
346 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): normalizing nullspace" << std::endl;
348 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
350 Array<Scalar> normsSC(Coords_->getNumVectors(), one / Teuchos::as<Scalar>(meanLen));
351 Nullspace_->scale(normsSC());
355 GetOStream(
Errors) <<
"MueLu::RefMaxwell::compute(): either the nullspace or the nodal coordinates must be provided." << std::endl;
358 if (!reuse && skipFirstLevel_) {
360 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 368 dump(*Nullspace_,
"nullspace.m");
375 if (skipFirstLevel_) {
377 Level fineLevel, coarseLevel;
383 fineLevel.
Set(
"A",Ms_Matrix_);
384 coarseLevel.
Set(
"P",D0_Matrix_);
385 coarseLevel.
setlib(Ms_Matrix_->getDomainMap()->lib());
386 fineLevel.
setlib(Ms_Matrix_->getDomainMap()->lib());
387 coarseLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
388 fineLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
390 RCP<RAPFactory> rapFact = rcp(
new RAPFactory());
391 ParameterList rapList = *(rapFact->GetValidParameterList());
392 rapList.set(
"transpose: use implicit",
true);
393 rapList.set(
"rap: fix zero diagonals", parameterList_.get<
bool>(
"rap: fix zero diagonals",
true));
394 rapList.set(
"rap: fix zero diagonals threshold", parameterList_.get<
double>(
"rap: fix zero diagonals threshold", Teuchos::ScalarTraits<double>::eps()));
395 rapList.set(
"rap: triple product", parameterList_.get<
bool>(
"rap: triple product",
false));
396 rapFact->SetParameterList(rapList);
399 coarseLevel.
Request(
"A", rapFact.get());
401 A_nodal_Matrix_ = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
403 if (applyBCsToAnodal_) {
405 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 412 dump(*A_nodal_Matrix_,
"A_nodal.m");
416 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): building special prolongator" << std::endl;
421 bool doRebalancing =
false;
422 doRebalancing = parameterList_.get<
bool>(
"refmaxwell: subsolves on subcommunicators", MasterList::getDefault<bool>(
"refmaxwell: subsolves on subcommunicators"));
423 int rebalanceStriding = parameterList_.get<
int>(
"refmaxwell: subsolves striding", -1);
424 int numProcsAH, numProcsA22;
425 int numProcs = SM_Matrix_->getDomainMap()->getComm()->getSize();
427 doRebalancing =
false;
446 ParameterList repartheurParams;
447 repartheurParams.set(
"repartition: start level", 0);
449 int defaultTargetRows = 10000;
450 repartheurParams.set(
"repartition: min rows per proc", precList11_.get<
int>(
"repartition: target rows per proc", defaultTargetRows));
451 repartheurParams.set(
"repartition: target rows per proc", precList11_.get<
int>(
"repartition: target rows per proc", defaultTargetRows));
452 repartheurParams.set(
"repartition: min rows per thread", precList11_.get<
int>(
"repartition: target rows per thread", defaultTargetRows));
453 repartheurParams.set(
"repartition: target rows per thread", precList11_.get<
int>(
"repartition: target rows per thread", defaultTargetRows));
454 repartheurParams.set(
"repartition: max imbalance", precList11_.get<
double>(
"repartition: max imbalance", 1.1));
455 repartheurFactory->SetParameterList(repartheurParams);
457 level.
Request(
"number of partitions", repartheurFactory.get());
458 repartheurFactory->Build(level);
459 numProcsAH = level.
Get<
int>(
"number of partitions", repartheurFactory.get());
460 numProcsAH = std::min(numProcsAH,numProcs);
470 level.
Set(
"Map",D0_Matrix_->getDomainMap());
473 ParameterList repartheurParams;
474 repartheurParams.set(
"repartition: start level", 0);
475 repartheurParams.set(
"repartition: use map",
true);
477 int defaultTargetRows = 10000;
478 repartheurParams.set(
"repartition: min rows per proc", precList22_.get<
int>(
"repartition: target rows per proc", defaultTargetRows));
479 repartheurParams.set(
"repartition: target rows per proc", precList22_.get<
int>(
"repartition: target rows per proc", defaultTargetRows));
480 repartheurParams.set(
"repartition: min rows per thread", precList22_.get<
int>(
"repartition: target rows per thread", defaultTargetRows));
481 repartheurParams.set(
"repartition: target rows per thread", precList22_.get<
int>(
"repartition: target rows per thread", defaultTargetRows));
483 repartheurFactory->SetParameterList(repartheurParams);
485 level.
Request(
"number of partitions", repartheurFactory.get());
486 repartheurFactory->Build(level);
487 numProcsA22 = level.
Get<
int>(
"number of partitions", repartheurFactory.get());
488 numProcsA22 = std::min(numProcsA22,numProcs);
491 if (rebalanceStriding >= 1) {
492 TEUCHOS_ASSERT(rebalanceStriding*numProcsAH<=numProcs);
493 TEUCHOS_ASSERT(rebalanceStriding*numProcsA22<=numProcs);
494 if (rebalanceStriding*(numProcsAH+numProcsA22)>numProcs) {
495 GetOStream(
Warnings0) <<
"RefMaxwell::compute(): Disabling striding = " << rebalanceStriding <<
", since AH needs " << numProcsAH
496 <<
" procs and A22 needs " << numProcsA22 <<
" procs."<< std::endl;
497 rebalanceStriding = -1;
505 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu RefMaxwell: Rebalance AH");
507 Level fineLevel, coarseLevel;
513 coarseLevel.
Set(
"A",AH_);
514 coarseLevel.
Set(
"P",P11_);
515 coarseLevel.
Set(
"Coordinates",CoordsH_);
516 if (!NullspaceH_.is_null())
517 coarseLevel.
Set(
"Nullspace",NullspaceH_);
518 coarseLevel.
Set(
"number of partitions", numProcsAH);
519 coarseLevel.
Set(
"repartition: heuristic target rows per process", 1000);
521 coarseLevel.
setlib(AH_->getDomainMap()->lib());
522 fineLevel.
setlib(AH_->getDomainMap()->lib());
523 coarseLevel.setObjectLabel(
"RefMaxwell coarse (1,1)");
524 fineLevel.setObjectLabel(
"RefMaxwell coarse (1,1)");
526 std::string partName = precList11_.get<std::string>(
"repartition: partitioner",
"zoltan2");
527 RCP<Factory> partitioner;
528 if (partName ==
"zoltan") {
529 #ifdef HAVE_MUELU_ZOLTAN 536 }
else if (partName ==
"zoltan2") {
537 #ifdef HAVE_MUELU_ZOLTAN2 539 ParameterList partParams;
540 RCP<const ParameterList> partpartParams = rcp(
new ParameterList(precList11_.sublist(
"repartition: params",
false)));
541 partParams.set(
"ParameterList", partpartParams);
542 partitioner->SetParameterList(partParams);
550 ParameterList repartParams;
551 repartParams.set(
"repartition: print partition distribution", precList11_.get<
bool>(
"repartition: print partition distribution",
false));
552 repartParams.set(
"repartition: remap parts", precList11_.get<
bool>(
"repartition: remap parts",
true));
553 if (rebalanceStriding >= 1) {
554 bool acceptPart = (SM_Matrix_->getDomainMap()->getComm()->getRank() % rebalanceStriding) == 0;
555 if (SM_Matrix_->getDomainMap()->getComm()->getRank() >= numProcsAH*rebalanceStriding)
557 repartParams.set(
"repartition: remap accept partition", acceptPart);
559 repartFactory->SetParameterList(repartParams);
561 repartFactory->SetFactory(
"Partition", partitioner);
564 ParameterList newPparams;
565 newPparams.set(
"type",
"Interpolation");
566 newPparams.set(
"repartition: rebalance P and R", precList11_.get<
bool>(
"repartition: rebalance P and R",
false));
567 newPparams.set(
"repartition: use subcommunicators",
true);
568 newPparams.set(
"repartition: rebalance Nullspace", !NullspaceH_.is_null());
570 if (!NullspaceH_.is_null())
572 newP->SetParameterList(newPparams);
573 newP->SetFactory(
"Importer", repartFactory);
576 ParameterList rebAcParams;
577 rebAcParams.set(
"repartition: use subcommunicators",
true);
578 newA->SetParameterList(rebAcParams);
579 newA->SetFactory(
"Importer", repartFactory);
581 coarseLevel.
Request(
"P", newP.get());
582 coarseLevel.
Request(
"Importer", repartFactory.get());
583 coarseLevel.
Request(
"A", newA.get());
584 coarseLevel.
Request(
"Coordinates", newP.get());
585 if (!NullspaceH_.is_null())
586 coarseLevel.
Request(
"Nullspace", newP.get());
587 repartFactory->Build(coarseLevel);
589 if (!precList11_.get<
bool>(
"repartition: rebalance P and R",
false))
590 ImporterH_ = coarseLevel.
Get< RCP<const Import> >(
"Importer", repartFactory.get());
591 P11_ = coarseLevel.
Get< RCP<Matrix> >(
"P", newP.get());
592 AH_ = coarseLevel.
Get< RCP<Matrix> >(
"A", newA.get());
593 CoordsH_ = coarseLevel.
Get< RCP<RealValuedMultiVector> >(
"Coordinates", newP.get());
594 if (!NullspaceH_.is_null())
595 NullspaceH_ = coarseLevel.
Get< RCP<MultiVector> >(
"Nullspace", newP.get());
597 AH_AP_reuse_data_ = Teuchos::null;
598 AH_RAP_reuse_data_ = Teuchos::null;
600 if (!disable_addon_ && enable_reuse_) {
602 RCP<const Import> ImporterH = coarseLevel.
Get< RCP<const Import> >(
"Importer", repartFactory.get());
603 RCP<const Map> targetMap = ImporterH->getTargetMap();
604 ParameterList XpetraList;
605 XpetraList.set(
"Restrict Communicator",
true);
606 Addon_Matrix_ = MatrixFactory::Build(Addon_Matrix_, *ImporterH, *ImporterH, targetMap, targetMap, rcp(&XpetraList,
false));
611 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 615 if (useKokkos_ && precList11_.isParameter(
"aggregation: drop tol") && precList11_.get<
double>(
"aggregation: drop tol") != 0.0) {
616 GetOStream(
Warnings0) <<
"RefMaxwell::compute(): Setting \"aggregation: drop tol\". to 0.0, since CoalesceDropFactory_kokkos does not " 617 <<
"support BlockSize > 1 and drop tol != 0.0" << std::endl;
618 precList11_.set(
"aggregation: drop tol", 0.0);
621 dump(*P11_,
"P11.m");
623 if (!implicitTranspose_) {
624 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 626 R11_ = Utilities_kokkos::Transpose(*P11_);
632 dump(*R11_,
"R11.m");
637 if (!AH_.is_null()) {
639 size_t dim = Nullspace_->getNumVectors();
640 AH_->SetFixedBlockSize(dim);
641 AH_->setObjectLabel(
"RefMaxwell coarse (1,1)");
643 int oldRank = SetProcRankVerbose(AH_->getDomainMap()->getComm()->getRank());
645 RCP<ParameterList> params = rcp(
new ParameterList());;
646 params->set(
"printLoadBalancingInfo",
true);
647 params->set(
"printCommInfo",
true);
650 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 651 if (precList11_.isType<std::string>(
"Preconditioner Type")) {
653 if (precList11_.get<std::string>(
"Preconditioner Type") ==
"MueLu") {
654 ParameterList& userParamList = precList11_.sublist(
"Preconditioner Types").sublist(
"MueLu").sublist(
"user data");
655 userParamList.set<RCP<RealValuedMultiVector> >(
"Coordinates", CoordsH_);
664 dumpCoords(*CoordsH_,
"coordsH.m");
665 if (!NullspaceH_.is_null())
666 dump(*NullspaceH_,
"NullspaceH.m");
667 ParameterList& userParamList = precList11_.sublist(
"user data");
668 userParamList.set<RCP<RealValuedMultiVector> >(
"Coordinates", CoordsH_);
669 if (!NullspaceH_.is_null())
670 userParamList.set<RCP<MultiVector> >(
"Nullspace", NullspaceH_);
673 RCP<MueLu::Level> level0 = HierarchyH_->GetLevel(0);
674 level0->Set(
"A", AH_);
675 HierarchyH_->SetupRe();
678 SetProcRankVerbose(oldRank);
684 if(!reuse && applyBCsTo22_) {
685 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): nuking BC nodes of D0" << std::endl;
687 D0_Matrix_->resumeFill();
689 if (D0_Matrix_->getRowMap()->lib() == Xpetra::UseEpetra)
690 replaceWith= Teuchos::ScalarTraits<SC>::eps();
692 replaceWith = Teuchos::ScalarTraits<SC>::zero();
693 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 702 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
707 if (!allNodesBoundary_) {
708 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): building MG for (2,2)-block" << std::endl;
711 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu RefMaxwell: Build A22");
713 Level fineLevel, coarseLevel;
719 fineLevel.
Set(
"A",SM_Matrix_);
720 coarseLevel.
Set(
"P",D0_Matrix_);
721 coarseLevel.
Set(
"Coordinates",Coords_);
723 coarseLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
724 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
725 coarseLevel.setObjectLabel(
"RefMaxwell (2,2)");
726 fineLevel.setObjectLabel(
"RefMaxwell (2,2)");
728 RCP<RAPFactory> rapFact = rcp(
new RAPFactory());
729 ParameterList rapList = *(rapFact->GetValidParameterList());
730 rapList.set(
"transpose: use implicit",
true);
731 rapList.set(
"rap: fix zero diagonals", parameterList_.get<
bool>(
"rap: fix zero diagonals",
true));
732 rapList.set(
"rap: fix zero diagonals threshold", parameterList_.get<
double>(
"rap: fix zero diagonals threshold", Teuchos::ScalarTraits<double>::eps()));
733 rapList.set(
"rap: triple product", parameterList_.get<
bool>(
"rap: triple product",
false));
734 rapFact->SetParameterList(rapList);
736 if (!A22_AP_reuse_data_.is_null()) {
737 coarseLevel.
AddKeepFlag(
"AP reuse data", rapFact.get());
738 coarseLevel.
Set<Teuchos::RCP<Teuchos::ParameterList> >(
"AP reuse data", A22_AP_reuse_data_, rapFact.get());
740 if (!A22_RAP_reuse_data_.is_null()) {
741 coarseLevel.
AddKeepFlag(
"RAP reuse data", rapFact.get());
742 coarseLevel.
Set<Teuchos::RCP<Teuchos::ParameterList> >(
"RAP reuse data", A22_RAP_reuse_data_, rapFact.get());
748 coarseLevel.
Set(
"number of partitions", numProcsA22);
749 coarseLevel.
Set(
"repartition: heuristic target rows per process", 1000);
751 std::string partName = precList22_.get<std::string>(
"repartition: partitioner",
"zoltan2");
752 RCP<Factory> partitioner;
753 if (partName ==
"zoltan") {
754 #ifdef HAVE_MUELU_ZOLTAN 756 partitioner->SetFactory(
"A", rapFact);
762 }
else if (partName ==
"zoltan2") {
763 #ifdef HAVE_MUELU_ZOLTAN2 765 ParameterList partParams;
766 RCP<const ParameterList> partpartParams = rcp(
new ParameterList(precList22_.sublist(
"repartition: params",
false)));
767 partParams.set(
"ParameterList", partpartParams);
768 partitioner->SetParameterList(partParams);
769 partitioner->SetFactory(
"A", rapFact);
777 ParameterList repartParams;
778 repartParams.set(
"repartition: print partition distribution", precList22_.get<
bool>(
"repartition: print partition distribution",
false));
779 repartParams.set(
"repartition: remap parts", precList22_.get<
bool>(
"repartition: remap parts",
true));
780 if (rebalanceStriding >= 1) {
781 bool acceptPart = ((SM_Matrix_->getDomainMap()->getComm()->getSize()-1-SM_Matrix_->getDomainMap()->getComm()->getRank()) % rebalanceStriding) == 0;
782 if (SM_Matrix_->getDomainMap()->getComm()->getSize()-1-SM_Matrix_->getDomainMap()->getComm()->getRank() >= numProcsA22*rebalanceStriding)
785 TEUCHOS_ASSERT(AH_.is_null());
786 repartParams.set(
"repartition: remap accept partition", acceptPart);
788 repartParams.set(
"repartition: remap accept partition", AH_.is_null());
789 repartFactory->SetParameterList(repartParams);
790 repartFactory->SetFactory(
"A", rapFact);
792 repartFactory->SetFactory(
"Partition", partitioner);
795 ParameterList newPparams;
796 newPparams.set(
"type",
"Interpolation");
797 newPparams.set(
"repartition: rebalance P and R", precList22_.get<
bool>(
"repartition: rebalance P and R",
false));
798 newPparams.set(
"repartition: use subcommunicators",
true);
799 newPparams.set(
"repartition: rebalance Nullspace",
false);
801 newP->SetParameterList(newPparams);
802 newP->SetFactory(
"Importer", repartFactory);
805 ParameterList rebAcParams;
806 rebAcParams.set(
"repartition: use subcommunicators",
true);
807 newA->SetParameterList(rebAcParams);
808 newA->SetFactory(
"A", rapFact);
809 newA->SetFactory(
"Importer", repartFactory);
811 coarseLevel.
Request(
"P", newP.get());
812 coarseLevel.
Request(
"Importer", repartFactory.get());
813 coarseLevel.
Request(
"A", newA.get());
814 coarseLevel.
Request(
"Coordinates", newP.get());
815 rapFact->Build(fineLevel,coarseLevel);
816 repartFactory->Build(coarseLevel);
818 if (!precList22_.get<
bool>(
"repartition: rebalance P and R",
false))
819 Importer22_ = coarseLevel.
Get< RCP<const Import> >(
"Importer", repartFactory.get());
820 D0_Matrix_ = coarseLevel.
Get< RCP<Matrix> >(
"P", newP.get());
821 A22_ = coarseLevel.
Get< RCP<Matrix> >(
"A", newA.get());
822 Coords_ = coarseLevel.
Get< RCP<RealValuedMultiVector> >(
"Coordinates", newP.get());
827 coarseLevel.
Request(
"A", rapFact.get());
829 coarseLevel.
Request(
"AP reuse data", rapFact.get());
830 coarseLevel.
Request(
"RAP reuse data", rapFact.get());
833 A22_ = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
836 if (!parameterList_.get<
bool>(
"rap: triple product",
false))
837 A22_AP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"AP reuse data", rapFact.get());
838 A22_RAP_reuse_data_ = coarseLevel.
Get< RCP<ParameterList> >(
"RAP reuse data", rapFact.get());
842 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu RefMaxwell: Build A22");
843 if (Importer22_.is_null()) {
845 temp = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*SM_Matrix_,
false,*D0_Matrix_,
false,temp,GetOStream(
Runtime0),
true,
true);
846 if (!implicitTranspose_)
847 A22_ = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_T_Matrix_,
false,*temp,
false,A22_,GetOStream(
Runtime0),
true,
true);
849 A22_ = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
true,*temp,
false,A22_,GetOStream(
Runtime0),
true,
true);
852 RCP<const Import> D0importer = rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->getCrsGraph()->getImporter();
853 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->replaceDomainMapAndImporter(D0origDomainMap_, D0origImporter_);
855 RCP<Matrix> temp, temp2;
856 temp = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*SM_Matrix_,
false,*D0_Matrix_,
false,temp,GetOStream(
Runtime0),
true,
true);
857 if (!implicitTranspose_)
858 temp2 = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_T_Matrix_,
false,*temp,
false,temp2,GetOStream(
Runtime0),
true,
true);
860 temp2 = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
true,*temp,
false,temp2,GetOStream(
Runtime0),
true,
true);
863 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->replaceDomainMapAndImporter(Importer22_->getTargetMap(), D0importer);
865 ParameterList XpetraList;
866 XpetraList.set(
"Restrict Communicator",
true);
867 XpetraList.set(
"Timer Label",
"MueLu::RebalanceA22");
868 RCP<const Map> targetMap = Importer22_->getTargetMap();
869 A22_ = MatrixFactory::Build(temp2, *Importer22_, *Importer22_, targetMap, targetMap, rcp(&XpetraList,
false));
873 if (!implicitTranspose_ && !reuse) {
874 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 876 D0_T_Matrix_ = Utilities_kokkos::Transpose(*D0_Matrix_);
885 if (!A22_.is_null()) {
886 dump(*A22_,
"A22.m");
887 A22_->setObjectLabel(
"RefMaxwell (2,2)");
888 int oldRank = SetProcRankVerbose(A22_->getDomainMap()->getComm()->getRank());
890 RCP<ParameterList> params = rcp(
new ParameterList());;
891 params->set(
"printLoadBalancingInfo",
true);
892 params->set(
"printCommInfo",
true);
895 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 896 if (precList22_.isType<std::string>(
"Preconditioner Type")) {
898 if (precList22_.get<std::string>(
"Preconditioner Type") ==
"MueLu") {
899 ParameterList& userParamList = precList22_.sublist(
"Preconditioner Types").sublist(
"MueLu").sublist(
"user data");
900 userParamList.set<RCP<RealValuedMultiVector> >(
"Coordinates", Coords_);
908 ParameterList& userParamList = precList22_.sublist(
"user data");
909 userParamList.set<RCP<RealValuedMultiVector> >(
"Coordinates", Coords_);
912 std::string coarseType =
"";
913 if (precList22_.isParameter(
"coarse: type")) {
914 coarseType = precList22_.get<std::string>(
"coarse: type");
916 std::transform(coarseType.begin(), coarseType.end(), coarseType.begin(),
::tolower);
917 std::transform(coarseType.begin(), ++coarseType.begin(), coarseType.begin(), ::toupper);
921 coarseType ==
"Klu" ||
922 coarseType ==
"Klu2") &&
923 (!precList22_.isSublist(
"coarse: params") ||
924 !precList22_.sublist(
"coarse: params").isParameter(
"fix nullspace")))
925 precList22_.sublist(
"coarse: params").set(
"fix nullspace",
true);
928 RCP<MueLu::Level> level0 = Hierarchy22_->GetLevel(0);
929 level0->Set(
"A", A22_);
930 Hierarchy22_->SetupRe();
933 SetProcRankVerbose(oldRank);
939 if(!reuse && !allNodesBoundary_ && applyBCsTo22_) {
940 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): nuking BC edges of D0" << std::endl;
942 D0_Matrix_->resumeFill();
944 if (D0_Matrix_->getRowMap()->lib() == Xpetra::UseEpetra)
945 replaceWith= Teuchos::ScalarTraits<SC>::eps();
947 replaceWith = Teuchos::ScalarTraits<SC>::zero();
948 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 957 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
958 dump(*D0_Matrix_,
"D0_nuked.m");
961 setFineLevelSmoother();
964 if (!ImporterH_.is_null()) {
965 RCP<const Import> ImporterP11 = ImportFactory::Build(ImporterH_->getTargetMap(),P11_->getColMap());
966 rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix()->replaceDomainMapAndImporter(ImporterH_->getTargetMap(), ImporterP11);
969 if (!Importer22_.is_null()) {
971 D0origDomainMap_ = D0_Matrix_->getDomainMap();
972 D0origImporter_ = rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->getCrsGraph()->getImporter();
974 RCP<const Import> ImporterD0 = ImportFactory::Build(Importer22_->getTargetMap(),D0_Matrix_->getColMap());
975 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->replaceDomainMapAndImporter(Importer22_->getTargetMap(), ImporterD0);
978 #ifdef HAVE_MUELU_TPETRA 979 if ((!D0_T_Matrix_.is_null()) &&
981 (!rcp_dynamic_cast<CrsMatrixWrap>(D0_T_Matrix_)->getCrsMatrix()->getCrsGraph()->getImporter().is_null()) &&
982 (!rcp_dynamic_cast<CrsMatrixWrap>(R11_)->getCrsMatrix()->getCrsGraph()->getImporter().is_null()) &&
983 (D0_T_Matrix_->getColMap()->lib() == Xpetra::UseTpetra) &&
984 (R11_->getColMap()->lib() == Xpetra::UseTpetra))
985 D0_T_R11_colMapsMatch_ = D0_T_Matrix_->getColMap()->isSameAs(*R11_->getColMap());
988 D0_T_R11_colMapsMatch_ =
false;
989 if (D0_T_R11_colMapsMatch_)
990 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): D0_T and R11 have matching colMaps" << std::endl;
995 if (parameterList_.isSublist(
"matvec params"))
997 RCP<ParameterList> matvecParams = rcpFromRef(parameterList_.sublist(
"matvec params"));
1002 if (!ImporterH_.is_null()) ImporterH_->setDistributorParameters(matvecParams);
1003 if (!Importer22_.is_null()) Importer22_->setDistributorParameters(matvecParams);
1005 if (!ImporterH_.is_null() && parameterList_.isSublist(
"refmaxwell: ImporterH params")){
1006 RCP<ParameterList> importerParams = rcpFromRef(parameterList_.sublist(
"refmaxwell: ImporterH params"));
1007 ImporterH_->setDistributorParameters(importerParams);
1009 if (!Importer22_.is_null() && parameterList_.isSublist(
"refmaxwell: Importer22 params")){
1010 RCP<ParameterList> importerParams = rcpFromRef(parameterList_.sublist(
"refmaxwell: Importer22 params"));
1011 Importer22_->setDistributorParameters(importerParams);
1017 #ifdef HAVE_MUELU_CUDA 1018 if (parameterList_.get<
bool>(
"refmaxwell: cuda profile setup",
false)) cudaProfilerStop();
1023 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1025 if (parameterList_.isType<std::string>(
"smoother: type") &&
1026 parameterList_.get<std::string>(
"smoother: type") ==
"hiptmair" &&
1027 SM_Matrix_->getDomainMap()->lib() == Xpetra::UseTpetra &&
1028 A22_->getDomainMap()->lib() == Xpetra::UseTpetra &&
1029 D0_Matrix_->getDomainMap()->lib() == Xpetra::UseTpetra) {
1030 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 1031 ParameterList hiptmairPreList, hiptmairPostList, smootherPreList, smootherPostList;
1033 if (smootherList_.isSublist(
"smoother: pre params"))
1034 smootherPreList = smootherList_.sublist(
"smoother: pre params");
1035 else if (smootherList_.isSublist(
"smoother: params"))
1036 smootherPreList = smootherList_.sublist(
"smoother: params");
1037 hiptmairPreList.set(
"hiptmair: smoother type 1",
1038 smootherPreList.get<std::string>(
"hiptmair: smoother type 1",
"CHEBYSHEV"));
1039 hiptmairPreList.set(
"hiptmair: smoother type 2",
1040 smootherPreList.get<std::string>(
"hiptmair: smoother type 2",
"CHEBYSHEV"));
1041 if(smootherPreList.isSublist(
"hiptmair: smoother list 1"))
1042 hiptmairPreList.set(
"hiptmair: smoother list 1", smootherPreList.sublist(
"hiptmair: smoother list 1"));
1043 if(smootherPreList.isSublist(
"hiptmair: smoother list 2"))
1044 hiptmairPreList.set(
"hiptmair: smoother list 2", smootherPreList.sublist(
"hiptmair: smoother list 2"));
1045 hiptmairPreList.set(
"hiptmair: pre or post",
1046 smootherPreList.get<std::string>(
"hiptmair: pre or post",
"pre"));
1047 hiptmairPreList.set(
"hiptmair: zero starting solution",
1048 smootherPreList.get<
bool>(
"hiptmair: zero starting solution",
true));
1050 if (smootherList_.isSublist(
"smoother: post params"))
1051 smootherPostList = smootherList_.sublist(
"smoother: post params");
1052 else if (smootherList_.isSublist(
"smoother: params"))
1053 smootherPostList = smootherList_.sublist(
"smoother: params");
1054 hiptmairPostList.set(
"hiptmair: smoother type 1",
1055 smootherPostList.get<std::string>(
"hiptmair: smoother type 1",
"CHEBYSHEV"));
1056 hiptmairPostList.set(
"hiptmair: smoother type 2",
1057 smootherPostList.get<std::string>(
"hiptmair: smoother type 2",
"CHEBYSHEV"));
1058 if(smootherPostList.isSublist(
"hiptmair: smoother list 1"))
1059 hiptmairPostList.set(
"hiptmair: smoother list 1", smootherPostList.sublist(
"hiptmair: smoother list 1"));
1060 if(smootherPostList.isSublist(
"hiptmair: smoother list 2"))
1061 hiptmairPostList.set(
"hiptmair: smoother list 2", smootherPostList.sublist(
"hiptmair: smoother list 2"));
1062 hiptmairPostList.set(
"hiptmair: pre or post",
1063 smootherPostList.get<std::string>(
"hiptmair: pre or post",
"post"));
1064 hiptmairPostList.set(
"hiptmair: zero starting solution",
1065 smootherPostList.get<
bool>(
"hiptmair: zero starting solution",
false));
1067 typedef Tpetra::RowMatrix<SC, LO, GO, NO> TROW;
1073 hiptmairPreSmoother_ -> setParameters(hiptmairPreList);
1074 hiptmairPreSmoother_ -> initialize();
1075 hiptmairPreSmoother_ -> compute();
1077 hiptmairPostSmoother_ -> setParameters(hiptmairPostList);
1078 hiptmairPostSmoother_ -> initialize();
1079 hiptmairPostSmoother_ -> compute();
1080 useHiptmairSmoothing_ =
true;
1082 throw(Xpetra::Exceptions::RuntimeError(
"MueLu must be compiled with Ifpack2 for Hiptmair smoothing."));
1083 #endif // defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 1087 RCP<MueLu::FactoryManagerBase> factoryHandler = rcp(
new FactoryManager());
1090 level.setObjectLabel(
"RefMaxwell (1,1)");
1091 level.
Set(
"A",SM_Matrix_);
1092 level.
setlib(SM_Matrix_->getDomainMap()->lib());
1094 if (parameterList_.isType<std::string>(
"smoother: pre type") && parameterList_.isType<std::string>(
"smoother: post type")) {
1095 std::string preSmootherType = parameterList_.get<std::string>(
"smoother: pre type");
1096 std::string postSmootherType = parameterList_.get<std::string>(
"smoother: post type");
1098 ParameterList preSmootherList, postSmootherList;
1099 if (parameterList_.isSublist(
"smoother: pre params"))
1100 preSmootherList = parameterList_.sublist(
"smoother: pre params");
1101 if (parameterList_.isSublist(
"smoother: post params"))
1102 postSmootherList = parameterList_.sublist(
"smoother: post params");
1104 RCP<SmootherPrototype> preSmootherPrototype = rcp(
new TrilinosSmoother(preSmootherType, preSmootherList));
1105 RCP<SmootherPrototype> postSmootherPrototype = rcp(
new TrilinosSmoother(postSmootherType, postSmootherList));
1106 RCP<SmootherFactory> smootherFact = rcp(
new SmootherFactory(preSmootherPrototype, postSmootherPrototype));
1108 level.
Request(
"PreSmoother",smootherFact.get());
1109 level.
Request(
"PostSmoother",smootherFact.get());
1110 if (enable_reuse_) {
1111 ParameterList smootherFactoryParams;
1112 smootherFactoryParams.set(
"keep smoother data",
true);
1113 smootherFact->SetParameterList(smootherFactoryParams);
1114 level.
Request(
"PreSmoother data", smootherFact.get());
1115 level.
Request(
"PostSmoother data", smootherFact.get());
1116 if (!PreSmootherData_.is_null())
1117 level.
Set(
"PreSmoother data", PreSmootherData_, smootherFact.get());
1118 if (!PostSmootherData_.is_null())
1119 level.
Set(
"PostSmoother data", PostSmootherData_, smootherFact.get());
1121 smootherFact->Build(level);
1122 PreSmoother_ = level.
Get<RCP<SmootherBase> >(
"PreSmoother",smootherFact.get());
1123 PostSmoother_ = level.
Get<RCP<SmootherBase> >(
"PostSmoother",smootherFact.get());
1124 if (enable_reuse_) {
1125 PreSmootherData_ = level.
Get<RCP<SmootherPrototype> >(
"PreSmoother data",smootherFact.get());
1126 PostSmootherData_ = level.
Get<RCP<SmootherPrototype> >(
"PostSmoother data",smootherFact.get());
1129 std::string smootherType = parameterList_.get<std::string>(
"smoother: type",
"CHEBYSHEV");
1131 RCP<SmootherPrototype> smootherPrototype = rcp(
new TrilinosSmoother(smootherType, smootherList_));
1132 RCP<SmootherFactory> smootherFact = rcp(
new SmootherFactory(smootherPrototype));
1133 level.
Request(
"PreSmoother",smootherFact.get());
1134 if (enable_reuse_) {
1135 ParameterList smootherFactoryParams;
1136 smootherFactoryParams.set(
"keep smoother data",
true);
1137 smootherFact->SetParameterList(smootherFactoryParams);
1138 level.
Request(
"PreSmoother data", smootherFact.get());
1139 if (!PreSmootherData_.is_null())
1140 level.
Set(
"PreSmoother data", PreSmootherData_, smootherFact.get());
1142 smootherFact->Build(level);
1143 PreSmoother_ = level.
Get<RCP<SmootherBase> >(
"PreSmoother",smootherFact.get());
1144 PostSmoother_ = PreSmoother_;
1146 PreSmootherData_ = level.
Get<RCP<SmootherPrototype> >(
"PreSmoother data",smootherFact.get());
1148 useHiptmairSmoothing_ =
false;
1153 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1156 RCP<Teuchos::TimeMonitor> tmAlloc = getTimer(
"MueLu RefMaxwell: Allocate MVs");
1158 if (!R11_.is_null())
1159 P11res_ = MultiVectorFactory::Build(R11_->getRangeMap(), numVectors);
1161 P11res_ = MultiVectorFactory::Build(P11_->getDomainMap(), numVectors);
1162 P11res_->setObjectLabel(
"P11res");
1163 if (D0_T_R11_colMapsMatch_) {
1164 D0TR11Tmp_ = MultiVectorFactory::Build(R11_->getColMap(), numVectors);
1165 D0TR11Tmp_->setObjectLabel(
"D0TR11Tmp");
1167 if (!ImporterH_.is_null()) {
1168 P11resTmp_ = MultiVectorFactory::Build(ImporterH_->getTargetMap(), numVectors);
1169 P11resTmp_->setObjectLabel(
"P11resTmp");
1170 P11x_ = MultiVectorFactory::Build(ImporterH_->getTargetMap(), numVectors);
1172 P11x_ = MultiVectorFactory::Build(P11_->getDomainMap(), numVectors);
1173 P11x_->setObjectLabel(
"P11x");
1174 if (!D0_T_Matrix_.is_null())
1175 D0res_ = MultiVectorFactory::Build(D0_T_Matrix_->getRangeMap(), numVectors);
1177 D0res_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
1178 D0res_->setObjectLabel(
"D0res");
1179 if (!Importer22_.is_null()) {
1180 D0resTmp_ = MultiVectorFactory::Build(Importer22_->getTargetMap(), numVectors);
1181 D0resTmp_->setObjectLabel(
"D0resTmp");
1182 D0x_ = MultiVectorFactory::Build(Importer22_->getTargetMap(), numVectors);
1184 D0x_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(), numVectors);
1185 D0x_->setObjectLabel(
"D0x");
1186 if (!AH_.is_null()) {
1187 if (!ImporterH_.is_null() && !implicitTranspose_)
1188 P11resSubComm_ = MultiVectorFactory::Build(P11resTmp_, Teuchos::View);
1190 P11resSubComm_ = MultiVectorFactory::Build(P11res_, Teuchos::View);
1191 P11resSubComm_->replaceMap(AH_->getRangeMap());
1192 P11resSubComm_->setObjectLabel(
"P11resSubComm");
1194 P11xSubComm_ = MultiVectorFactory::Build(P11x_, Teuchos::View);
1195 P11xSubComm_->replaceMap(AH_->getDomainMap());
1196 P11xSubComm_->setObjectLabel(
"P11xSubComm");
1198 if (!A22_.is_null()) {
1199 if (!Importer22_.is_null() && !implicitTranspose_)
1200 D0resSubComm_ = MultiVectorFactory::Build(D0resTmp_, Teuchos::View);
1202 D0resSubComm_ = MultiVectorFactory::Build(D0res_, Teuchos::View);
1203 D0resSubComm_->replaceMap(A22_->getRangeMap());
1204 D0resSubComm_->setObjectLabel(
"D0resSubComm");
1206 D0xSubComm_ = MultiVectorFactory::Build(D0x_, Teuchos::View);
1207 D0xSubComm_->replaceMap(A22_->getDomainMap());
1208 D0xSubComm_->setObjectLabel(
"D0xSubComm");
1210 residual_ = MultiVectorFactory::Build(SM_Matrix_->getDomainMap(), numVectors);
1211 residual_->setObjectLabel(
"residual");
1215 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1217 if (dump_matrices_) {
1218 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1219 Xpetra::IO<SC, LO, GO, NO>::Write(name, A);
1224 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1226 if (dump_matrices_) {
1227 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1228 Xpetra::IO<SC, LO, GO, NO>::Write(name, X);
1233 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1235 if (dump_matrices_) {
1236 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1237 Xpetra::IO<coordinateType, LO, GO, NO>::Write(name, X);
1242 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1244 if (dump_matrices_) {
1245 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1246 std::ofstream out(name);
1247 for (
size_t i = 0; i < Teuchos::as<size_t>(v.size()); i++)
1248 out << v[i] <<
"\n";
1252 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1253 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1255 if (dump_matrices_) {
1256 GetOStream(
Runtime0) <<
"Dumping to " << name << std::endl;
1257 std::ofstream out(name);
1258 auto vH = Kokkos::create_mirror_view (v);
1259 Kokkos::deep_copy(vH , v);
1260 for (
size_t i = 0; i < vH.size(); i++)
1261 out << vH[i] <<
"\n";
1266 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1270 return Teuchos::rcp(
new Teuchos::TimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name)));
1273 return Teuchos::rcp(
new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), SM_Matrix_->getRowMap()->getComm().ptr()));
1275 return Teuchos::rcp(
new Teuchos::SyncTimeMonitor(*Teuchos::TimeMonitor::getNewTimer(name), comm.ptr()));
1278 return Teuchos::null;
1282 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1295 const SC SC_ZERO = Teuchos::ScalarTraits<SC>::zero();
1296 const SC SC_ONE = Teuchos::ScalarTraits<SC>::one();
1297 const Scalar half = SC_ONE / (SC_ONE + SC_ONE);
1298 size_t dim = Nullspace_->getNumVectors();
1299 size_t numLocalRows = SM_Matrix_->getNodeNumRows();
1301 RCP<Matrix> P_nodal;
1302 RCP<CrsMatrix> P_nodal_imported;
1303 RCP<MultiVector> Nullspace_nodal;
1304 if (skipFirstLevel_) {
1307 bool read_P_from_file = parameterList_.get(
"refmaxwell: read_P_from_file",
false);
1308 if (read_P_from_file) {
1311 std::string P_filename = parameterList_.get(
"refmaxwell: P_filename",std::string(
"P.m"));
1312 std::string domainmap_filename = parameterList_.get(
"refmaxwell: P_domainmap_filename",std::string(
"domainmap_P.m"));
1313 std::string colmap_filename = parameterList_.get(
"refmaxwell: P_colmap_filename",std::string(
"colmap_P.m"));
1314 std::string coords_filename = parameterList_.get(
"refmaxwell: CoordsH",std::string(
"coordsH.m"));
1315 RCP<const Map> colmap = Xpetra::IO<SC, LO, GO, NO>::ReadMap(colmap_filename, A_nodal_Matrix_->getDomainMap()->lib(),A_nodal_Matrix_->getDomainMap()->getComm());
1316 RCP<const Map> domainmap = Xpetra::IO<SC, LO, GO, NO>::ReadMap(domainmap_filename, A_nodal_Matrix_->getDomainMap()->lib(),A_nodal_Matrix_->getDomainMap()->getComm());
1317 P_nodal = Xpetra::IO<SC, LO, GO, NO>::Read(P_filename, A_nodal_Matrix_->getDomainMap(), colmap, domainmap, A_nodal_Matrix_->getDomainMap());
1318 CoordsH_ = Xpetra::IO<typename RealValuedMultiVector::scalar_type, LO, GO, NO>::ReadMultiVector(coords_filename, domainmap);
1320 Level fineLevel, coarseLevel;
1326 fineLevel.
Set(
"A",A_nodal_Matrix_);
1327 fineLevel.
Set(
"Coordinates",Coords_);
1328 fineLevel.
Set(
"DofsPerNode",1);
1329 coarseLevel.
setlib(A_nodal_Matrix_->getDomainMap()->lib());
1330 fineLevel.
setlib(A_nodal_Matrix_->getDomainMap()->lib());
1331 coarseLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
1332 fineLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
1335 RCP<MultiVector> nullSpace = MultiVectorFactory::Build(A_nodal_Matrix_->getRowMap(),NSdim);
1336 nullSpace->putScalar(SC_ONE);
1337 fineLevel.
Set(
"Nullspace",nullSpace);
1339 RCP<Factory> amalgFact, dropFact, UncoupledAggFact, coarseMapFact, TentativePFact, Tfact, SaPFact;
1340 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1342 amalgFact = rcp(
new AmalgamationFactory_kokkos());
1343 dropFact = rcp(
new CoalesceDropFactory_kokkos());
1344 UncoupledAggFact = rcp(
new UncoupledAggregationFactory_kokkos());
1345 coarseMapFact = rcp(
new CoarseMapFactory_kokkos());
1346 TentativePFact = rcp(
new TentativePFactory_kokkos());
1347 if (parameterList_.get(
"multigrid algorithm",
"unsmoothed") ==
"sa")
1348 SaPFact = rcp(
new SaPFactory_kokkos());
1349 Tfact = rcp(
new CoordinatesTransferFactory_kokkos());
1358 if (parameterList_.get(
"multigrid algorithm",
"unsmoothed") ==
"sa")
1362 dropFact->SetFactory(
"UnAmalgamationInfo", amalgFact);
1363 double dropTol = parameterList_.get(
"aggregation: drop tol",0.0);
1364 std::string dropScheme = parameterList_.get(
"aggregation: drop scheme",
"classical");
1365 std::string distLaplAlgo = parameterList_.get(
"aggregation: distance laplacian algo",
"default");
1366 dropFact->SetParameter(
"aggregation: drop tol",Teuchos::ParameterEntry(dropTol));
1367 dropFact->SetParameter(
"aggregation: drop scheme",Teuchos::ParameterEntry(dropScheme));
1369 dropFact->SetParameter(
"aggregation: distance laplacian algo",Teuchos::ParameterEntry(distLaplAlgo));
1371 UncoupledAggFact->SetFactory(
"Graph", dropFact);
1372 int minAggSize = parameterList_.get(
"aggregation: min agg size",2);
1373 UncoupledAggFact->SetParameter(
"aggregation: min agg size",Teuchos::ParameterEntry(minAggSize));
1374 int maxAggSize = parameterList_.get(
"aggregation: max agg size",-1);
1375 UncoupledAggFact->SetParameter(
"aggregation: max agg size",Teuchos::ParameterEntry(maxAggSize));
1377 coarseMapFact->SetFactory(
"Aggregates", UncoupledAggFact);
1379 TentativePFact->SetFactory(
"Aggregates", UncoupledAggFact);
1380 TentativePFact->SetFactory(
"UnAmalgamationInfo", amalgFact);
1381 TentativePFact->SetFactory(
"CoarseMap", coarseMapFact);
1383 Tfact->SetFactory(
"Aggregates", UncoupledAggFact);
1384 Tfact->SetFactory(
"CoarseMap", coarseMapFact);
1386 if (parameterList_.get(
"multigrid algorithm",
"unsmoothed") ==
"sa") {
1387 SaPFact->SetFactory(
"P", TentativePFact);
1388 coarseLevel.
Request(
"P", SaPFact.get());
1390 coarseLevel.
Request(
"P",TentativePFact.get());
1391 coarseLevel.
Request(
"Nullspace",TentativePFact.get());
1392 coarseLevel.
Request(
"Coordinates",Tfact.get());
1394 RCP<AggregationExportFactory> aggExport;
1395 if (parameterList_.get(
"aggregation: export visualization data",
false)) {
1397 ParameterList aggExportParams;
1398 aggExportParams.set(
"aggregation: output filename",
"aggs.vtk");
1399 aggExportParams.set(
"aggregation: output file: agg style",
"Jacks");
1400 aggExport->SetParameterList(aggExportParams);
1402 aggExport->SetFactory(
"Aggregates", UncoupledAggFact);
1403 aggExport->SetFactory(
"UnAmalgamationInfo", amalgFact);
1404 fineLevel.
Request(
"Aggregates",UncoupledAggFact.get());
1405 fineLevel.
Request(
"UnAmalgamationInfo",amalgFact.get());
1408 if (parameterList_.get(
"multigrid algorithm",
"unsmoothed") ==
"sa")
1409 coarseLevel.
Get(
"P",P_nodal,SaPFact.get());
1411 coarseLevel.
Get(
"P",P_nodal,TentativePFact.get());
1412 coarseLevel.
Get(
"Nullspace",Nullspace_nodal,TentativePFact.get());
1413 coarseLevel.
Get(
"Coordinates",CoordsH_,Tfact.get());
1416 if (parameterList_.get(
"aggregation: export visualization data",
false))
1417 aggExport->Build(fineLevel, coarseLevel);
1419 dump(*P_nodal,
"P_nodal.m");
1420 dump(*Nullspace_nodal,
"Nullspace_nodal.m");
1423 RCP<CrsMatrix> D0Crs = rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix();
1426 int numProcs = P_nodal->getDomainMap()->getComm()->getSize();
1428 RCP<CrsMatrixWrap> P_nodal_temp;
1429 RCP<const Map> targetMap = D0Crs->getColMap();
1430 P_nodal_temp = rcp(
new CrsMatrixWrap(targetMap));
1431 RCP<const Import> importer = D0Crs->getCrsGraph()->getImporter();
1432 P_nodal_temp->doImport(*P_nodal, *importer, Xpetra::INSERT);
1433 P_nodal_temp->fillComplete(rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix()->getDomainMap(),
1434 rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix()->getRangeMap());
1435 P_nodal_imported = P_nodal_temp->getCrsMatrix();
1436 dump(*P_nodal_temp,
"P_nodal_imported.m");
1438 P_nodal_imported = rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix();
1442 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1445 using ATS = Kokkos::ArithTraits<SC>;
1446 using impl_ATS = Kokkos::ArithTraits<typename ATS::val_type>;
1449 typedef typename Matrix::local_matrix_type KCRS;
1450 typedef typename KCRS::device_type device_t;
1451 typedef typename KCRS::StaticCrsGraphType graph_t;
1452 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
1453 typedef typename graph_t::entries_type::non_const_type lno_nnz_view_t;
1454 typedef typename KCRS::values_type::non_const_type scalar_view_t;
1460 std::string defaultAlgo =
"mat-mat";
1461 std::string algo = parameterList_.get(
"refmaxwell: prolongator compute algorithm",defaultAlgo);
1463 if (skipFirstLevel_) {
1466 if (algo ==
"mat-mat") {
1467 RCP<Matrix> D0_P_nodal = MatrixFactory::Build(SM_Matrix_->getRowMap());
1468 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
false,*P_nodal,
false,*D0_P_nodal,
true,
true);
1470 #ifdef HAVE_MUELU_DEBUG 1471 TEUCHOS_ASSERT(D0_P_nodal->getColMap()->isSameAs(*P_nodal_imported->getColMap()));
1475 auto localD0P = D0_P_nodal->getLocalMatrixDevice();
1478 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
1479 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
1481 size_t nnzEstimate = dim*localD0P.graph.entries.size();
1482 lno_view_t P11rowptr(
"P11_rowptr", numLocalRows+1);
1483 lno_nnz_view_t P11colind(
"P11_colind",nnzEstimate);
1484 scalar_view_t P11vals(
"P11_vals",nnzEstimate);
1487 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_adjustRowptr", range_type(0,numLocalRows+1),
1488 KOKKOS_LAMBDA(
const size_t i) {
1489 P11rowptr(i) = dim*localD0P.graph.row_map(i);
1493 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_adjustColind", range_type(0,localD0P.graph.entries.size()),
1494 KOKKOS_LAMBDA(
const size_t jj) {
1495 for (
size_t k = 0; k < dim; k++) {
1496 P11colind(dim*jj+k) = dim*localD0P.graph.entries(jj)+k;
1497 P11vals(dim*jj+k) = SC_ZERO;
1501 auto localNullspace = Nullspace_->getDeviceLocalView(Xpetra::Access::ReadOnly);
1504 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
1508 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
1510 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
1512 auto localD0 = D0_Matrix_->getLocalMatrixDevice();
1513 auto localP = P_nodal_imported->getLocalMatrixDevice();
1514 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_enterValues_D0wZeros", range_type(0,numLocalRows),
1515 KOKKOS_LAMBDA(
const size_t i) {
1516 for (
size_t ll = localD0.graph.row_map(i); ll < localD0.graph.row_map(i+1); ll++) {
1517 LO l = localD0.graph.entries(ll);
1518 SC p = localD0.values(ll);
1519 if (impl_ATS::magnitude(p) < tol)
1521 for (
size_t jj = localP.graph.row_map(l); jj < localP.graph.row_map(l+1); jj++) {
1522 LO j = localP.graph.entries(jj);
1523 SC v = localP.values(jj);
1524 for (
size_t k = 0; k < dim; k++) {
1526 SC n = localNullspace(i,k);
1528 for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
1529 if (P11colind(m) == jNew)
1531 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA) && !defined(HAVE_MUELU_HIP) 1532 TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
1534 P11vals(m) += half * v * n;
1541 auto localD0 = D0_Matrix_->getLocalMatrixDevice();
1542 auto localP = P_nodal_imported->getLocalMatrixDevice();
1543 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_enterValues", range_type(0,numLocalRows),
1544 KOKKOS_LAMBDA(
const size_t i) {
1545 for (
size_t ll = localD0.graph.row_map(i); ll < localD0.graph.row_map(i+1); ll++) {
1546 LO l = localD0.graph.entries(ll);
1547 for (
size_t jj = localP.graph.row_map(l); jj < localP.graph.row_map(l+1); jj++) {
1548 LO j = localP.graph.entries(jj);
1549 SC v = localP.values(jj);
1550 for (
size_t k = 0; k < dim; k++) {
1552 SC n = localNullspace(i,k);
1554 for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
1555 if (P11colind(m) == jNew)
1557 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA) && !defined(HAVE_MUELU_HIP) 1558 TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
1560 P11vals(m) += half * v * n;
1567 P11_ = rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0));
1568 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
1569 P11Crs->setAllValues(P11rowptr, P11colind, P11vals);
1570 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
1573 TEUCHOS_TEST_FOR_EXCEPTION(
false,std::invalid_argument,algo <<
" is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
1575 NullspaceH_ = MultiVectorFactory::Build(P11_->getDomainMap(), dim);
1577 auto localNullspace_nodal = Nullspace_nodal->getDeviceLocalView(Xpetra::Access::ReadOnly);
1578 auto localNullspaceH = NullspaceH_->getDeviceLocalView(Xpetra::Access::ReadWrite);
1579 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_nullspace", range_type(0,Nullspace_nodal->getLocalLength()),
1580 KOKKOS_LAMBDA(
const size_t i) {
1581 Scalar val = localNullspace_nodal(i,0);
1582 for (
size_t j = 0; j < dim; j++)
1583 localNullspaceH(dim*i+j, j) = val;
1588 auto localD0 = D0_Matrix_->getLocalMatrixDevice();
1592 if (algo ==
"mat-mat") {
1595 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(D0_Matrix_->getColMap(), dim);
1596 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(D0_Matrix_->getDomainMap(), dim);
1598 size_t nnzEstimate = dim*localD0.graph.entries.size();
1599 lno_view_t P11rowptr(
"P11_rowptr", numLocalRows+1);
1600 lno_nnz_view_t P11colind(
"P11_colind",nnzEstimate);
1601 scalar_view_t P11vals(
"P11_vals",nnzEstimate);
1604 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_adjustRowptr", range_type(0,numLocalRows+1),
1605 KOKKOS_LAMBDA(
const size_t i) {
1606 P11rowptr(i) = dim*localD0.graph.row_map(i);
1610 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_adjustColind", range_type(0,localD0.graph.entries.size()),
1611 KOKKOS_LAMBDA(
const size_t jj) {
1612 for (
size_t k = 0; k < dim; k++) {
1613 P11colind(dim*jj+k) = dim*localD0.graph.entries(jj)+k;
1614 P11vals(dim*jj+k) = SC_ZERO;
1618 auto localNullspace = Nullspace_->getDeviceLocalView(Xpetra::Access::ReadOnly);
1621 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
1625 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
1627 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
1629 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_enterValues_D0wZeros", range_type(0,numLocalRows),
1630 KOKKOS_LAMBDA(
const size_t i) {
1631 for (
size_t jj = localD0.graph.row_map(i); jj < localD0.graph.row_map(i+1); jj++) {
1632 LO j = localD0.graph.entries(jj);
1633 SC p = localD0.values(jj);
1634 if (impl_ATS::magnitude(p) < tol)
1636 for (
size_t k = 0; k < dim; k++) {
1638 SC n = localNullspace(i,k);
1640 for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
1641 if (P11colind(m) == jNew)
1643 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA) && !defined(HAVE_MUELU_HIP) 1644 TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
1646 P11vals(m) += half * n;
1652 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_enterValues", range_type(0,numLocalRows),
1653 KOKKOS_LAMBDA(
const size_t i) {
1654 for (
size_t jj = localD0.graph.row_map(i); jj < localD0.graph.row_map(i+1); jj++) {
1655 LO j = localD0.graph.entries(jj);
1656 for (
size_t k = 0; k < dim; k++) {
1658 SC n = localNullspace(i,k);
1660 for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
1661 if (P11colind(m) == jNew)
1663 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA) && !defined(HAVE_MUELU_HIP) 1664 TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
1666 P11vals(m) += half * n;
1672 P11_ = rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0));
1673 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
1674 P11Crs->setAllValues(P11rowptr, P11colind, P11vals);
1675 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
1677 TEUCHOS_TEST_FOR_EXCEPTION(
false,std::invalid_argument,algo <<
" is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
1681 #endif // ifdef(HAVE_MUELU_KOKKOS_REFACTOR) 1684 ArrayRCP<ArrayRCP<const SC> > nullspaceRCP(dim);
1685 ArrayRCP<ArrayView<const SC> > nullspace(dim);
1686 for(
size_t i=0; i<dim; i++) {
1687 nullspaceRCP[i] = Nullspace_->getData(i);
1688 nullspace[i] = nullspaceRCP[i]();
1692 ArrayRCP<size_t> P11rowptr_RCP;
1693 ArrayRCP<LO> P11colind_RCP;
1694 ArrayRCP<SC> P11vals_RCP;
1703 std::string defaultAlgo =
"mat-mat";
1704 std::string algo = parameterList_.get(
"refmaxwell: prolongator compute algorithm",defaultAlgo);
1706 if (skipFirstLevel_) {
1708 if (algo ==
"mat-mat") {
1709 RCP<Matrix> D0_P_nodal = MatrixFactory::Build(SM_Matrix_->getRowMap());
1710 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
false,*P_nodal,
false,*D0_P_nodal,
true,
true);
1713 ArrayRCP<const size_t> D0rowptr_RCP;
1714 ArrayRCP<const LO> D0colind_RCP;
1715 ArrayRCP<const SC> D0vals_RCP;
1716 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->getAllValues(D0rowptr_RCP, D0colind_RCP, D0vals_RCP);
1720 ArrayView<const size_t> D0rowptr;
1721 ArrayView<const LO> D0colind;
1722 ArrayView<const SC> D0vals;
1723 D0rowptr = D0rowptr_RCP(); D0colind = D0colind_RCP(); D0vals = D0vals_RCP();
1726 ArrayRCP<const size_t> Prowptr_RCP;
1727 ArrayRCP<const LO> Pcolind_RCP;
1728 ArrayRCP<const SC> Pvals_RCP;
1729 P_nodal_imported->getAllValues(Prowptr_RCP, Pcolind_RCP, Pvals_RCP);
1730 ArrayView<const size_t> Prowptr;
1731 ArrayView<const LO> Pcolind;
1732 ArrayView<const SC> Pvals;
1733 Prowptr = Prowptr_RCP(); Pcolind = Pcolind_RCP(); Pvals = Pvals_RCP();
1736 ArrayRCP<const size_t> D0Prowptr_RCP;
1737 ArrayRCP<const LO> D0Pcolind_RCP;
1738 ArrayRCP<const SC> D0Pvals_RCP;
1739 rcp_dynamic_cast<CrsMatrixWrap>(D0_P_nodal)->getCrsMatrix()->getAllValues(D0Prowptr_RCP, D0Pcolind_RCP, D0Pvals_RCP);
1744 ArrayView<const size_t> D0Prowptr;
1745 ArrayView<const LO> D0Pcolind;
1746 D0Prowptr = D0Prowptr_RCP(); D0Pcolind = D0Pcolind_RCP();
1749 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
1750 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
1751 P11_ = rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0));
1752 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
1753 size_t nnzEstimate = dim*D0Prowptr[numLocalRows];
1754 P11Crs->allocateAllValues(nnzEstimate, P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1756 ArrayView<size_t> P11rowptr = P11rowptr_RCP();
1757 ArrayView<LO> P11colind = P11colind_RCP();
1758 ArrayView<SC> P11vals = P11vals_RCP();
1761 for (
size_t i = 0; i < numLocalRows+1; i++) {
1762 P11rowptr[i] = dim*D0Prowptr[i];
1766 for (
size_t jj = 0; jj < (size_t) D0Prowptr[numLocalRows]; jj++)
1767 for (
size_t k = 0; k < dim; k++) {
1768 P11colind[dim*jj+k] = dim*D0Pcolind[jj]+k;
1769 P11vals[dim*jj+k] = SC_ZERO;
1772 RCP<const Map> P_nodal_imported_colmap = P_nodal_imported->getColMap();
1773 RCP<const Map> D0_P_nodal_colmap = D0_P_nodal->getColMap();
1775 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
1779 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
1781 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
1782 for (
size_t i = 0; i < numLocalRows; i++) {
1783 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1784 LO l = D0colind[ll];
1786 if (Teuchos::ScalarTraits<Scalar>::magnitude(p) < tol)
1788 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1790 j = D0_P_nodal_colmap->getLocalElement(P_nodal_imported_colmap->getGlobalElement(j));
1792 for (
size_t k = 0; k < dim; k++) {
1794 SC n = nullspace[k][i];
1796 for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
1797 if (P11colind[m] == jNew)
1799 #ifdef HAVE_MUELU_DEBUG 1800 TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
1802 P11vals[m] += half * v * n;
1809 for (
size_t i = 0; i < numLocalRows; i++) {
1810 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1811 LO l = D0colind[ll];
1812 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1814 j = D0_P_nodal_colmap->getLocalElement(P_nodal_imported_colmap->getGlobalElement(j));
1816 for (
size_t k = 0; k < dim; k++) {
1818 SC n = nullspace[k][i];
1820 for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
1821 if (P11colind[m] == jNew)
1823 #ifdef HAVE_MUELU_DEBUG 1824 TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
1826 P11vals[m] += half * v * n;
1833 P11Crs->setAllValues(P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1834 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
1836 }
else if (algo ==
"gustavson") {
1837 ArrayRCP<const size_t> D0rowptr_RCP;
1838 ArrayRCP<const LO> D0colind_RCP;
1839 ArrayRCP<const SC> D0vals_RCP;
1840 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->getAllValues(D0rowptr_RCP, D0colind_RCP, D0vals_RCP);
1844 ArrayView<const size_t> D0rowptr;
1845 ArrayView<const LO> D0colind;
1846 ArrayView<const SC> D0vals;
1847 D0rowptr = D0rowptr_RCP(); D0colind = D0colind_RCP(); D0vals = D0vals_RCP();
1850 ArrayRCP<const size_t> Prowptr_RCP;
1851 ArrayRCP<const LO> Pcolind_RCP;
1852 ArrayRCP<const SC> Pvals_RCP;
1853 P_nodal_imported->getAllValues(Prowptr_RCP, Pcolind_RCP, Pvals_RCP);
1854 ArrayView<const size_t> Prowptr;
1855 ArrayView<const LO> Pcolind;
1856 ArrayView<const SC> Pvals;
1857 Prowptr = Prowptr_RCP(); Pcolind = Pcolind_RCP(); Pvals = Pvals_RCP();
1859 LO maxP11col = dim * P_nodal_imported->getColMap()->getMaxLocalIndex();
1860 const size_t ST_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
1861 Array<size_t> P11_status(dim*maxP11col, ST_INVALID);
1863 size_t nnz_alloc = dim*D0vals_RCP.size();
1866 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
1867 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
1868 P11_ = rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0));
1869 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
1870 P11Crs->allocateAllValues(nnz_alloc, P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1872 ArrayView<size_t> P11rowptr = P11rowptr_RCP();
1873 ArrayView<LO> P11colind = P11colind_RCP();
1874 ArrayView<SC> P11vals = P11vals_RCP();
1877 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
1881 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
1883 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
1886 for (
size_t i = 0; i < numLocalRows; i++) {
1888 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1889 LO l = D0colind[ll];
1891 if (Teuchos::ScalarTraits<Scalar>::magnitude(p) < tol)
1893 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1896 for (
size_t k = 0; k < dim; k++) {
1898 SC n = nullspace[k][i];
1900 if (P11_status[jNew] == ST_INVALID || P11_status[jNew] < nnz_old) {
1901 P11_status[jNew] = nnz;
1902 P11colind[nnz] = jNew;
1903 P11vals[nnz] = half * v * n;
1908 P11vals[P11_status[jNew]] += half * v * n;
1917 P11rowptr[numLocalRows] = nnz;
1921 for (
size_t i = 0; i < numLocalRows; i++) {
1923 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1924 LO l = D0colind[ll];
1925 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1928 for (
size_t k = 0; k < dim; k++) {
1930 SC n = nullspace[k][i];
1932 if (P11_status[jNew] == ST_INVALID || P11_status[jNew] < nnz_old) {
1933 P11_status[jNew] = nnz;
1934 P11colind[nnz] = jNew;
1935 P11vals[nnz] = half * v * n;
1940 P11vals[P11_status[jNew]] += half * v * n;
1949 P11rowptr[numLocalRows] = nnz;
1952 if (blockDomainMap->lib() == Xpetra::UseTpetra) {
1956 P11vals_RCP.resize(nnz);
1957 P11colind_RCP.resize(nnz);
1960 P11Crs->setAllValues(P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1961 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
1963 TEUCHOS_TEST_FOR_EXCEPTION(
false,std::invalid_argument,algo <<
" is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
1965 NullspaceH_ = MultiVectorFactory::Build(P11_->getDomainMap(), dim);
1967 ArrayRCP<const Scalar> ns_rcp = Nullspace_nodal->getData(0);
1968 ArrayView<const Scalar> ns_view = ns_rcp();
1969 for (
size_t i = 0; i < Nullspace_nodal->getLocalLength(); i++) {
1971 for (
size_t j = 0; j < dim; j++)
1972 NullspaceH_->replaceLocalValue(dim*i+j, j, val);
1977 ArrayRCP<const size_t> D0rowptr_RCP;
1978 ArrayRCP<const LO> D0colind_RCP;
1979 ArrayRCP<const SC> D0vals_RCP;
1980 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->getAllValues(D0rowptr_RCP, D0colind_RCP, D0vals_RCP);
1984 ArrayView<const size_t> D0rowptr;
1985 ArrayView<const LO> D0colind;
1986 ArrayView<const SC> D0vals;
1987 D0rowptr = D0rowptr_RCP(); D0colind = D0colind_RCP(); D0vals = D0vals_RCP();
1990 if (algo ==
"mat-mat") {
1993 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(D0_Matrix_->getColMap(), dim);
1994 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(D0_Matrix_->getDomainMap(), dim);
1995 P11_ = rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0));
1996 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
1997 size_t nnzEstimate = dim*D0rowptr[numLocalRows];
1998 P11Crs->allocateAllValues(nnzEstimate, P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
2000 ArrayView<size_t> P11rowptr = P11rowptr_RCP();
2001 ArrayView<LO> P11colind = P11colind_RCP();
2002 ArrayView<SC> P11vals = P11vals_RCP();
2005 for (
size_t i = 0; i < numLocalRows+1; i++) {
2006 P11rowptr[i] = dim*D0rowptr[i];
2010 for (
size_t jj = 0; jj < (size_t) D0rowptr[numLocalRows]; jj++)
2011 for (
size_t k = 0; k < dim; k++) {
2012 P11colind[dim*jj+k] = dim*D0colind[jj]+k;
2013 P11vals[dim*jj+k] = SC_ZERO;
2017 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
2021 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
2023 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
2024 for (
size_t i = 0; i < numLocalRows; i++) {
2025 for (
size_t jj = D0rowptr[i]; jj < D0rowptr[i+1]; jj++) {
2026 LO j = D0colind[jj];
2028 if (Teuchos::ScalarTraits<Scalar>::magnitude(p) < tol)
2030 for (
size_t k = 0; k < dim; k++) {
2032 SC n = nullspace[k][i];
2034 for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
2035 if (P11colind[m] == jNew)
2037 #ifdef HAVE_MUELU_DEBUG 2038 TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
2040 P11vals[m] += half * n;
2046 for (
size_t i = 0; i < numLocalRows; i++) {
2047 for (
size_t jj = D0rowptr[i]; jj < D0rowptr[i+1]; jj++) {
2048 LO j = D0colind[jj];
2050 for (
size_t k = 0; k < dim; k++) {
2052 SC n = nullspace[k][i];
2054 for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
2055 if (P11colind[m] == jNew)
2057 #ifdef HAVE_MUELU_DEBUG 2058 TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
2060 P11vals[m] += half * n;
2066 P11Crs->setAllValues(P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
2067 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
2070 TEUCHOS_TEST_FOR_EXCEPTION(
false,std::invalid_argument,algo <<
" is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
2075 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2077 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu RefMaxwell: Build coarse (1,1) matrix");
2079 const Scalar one = Teuchos::ScalarTraits<Scalar>::one();
2083 temp = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*SM_Matrix_,
false,*P11_,
false,temp,GetOStream(
Runtime0),
true,
true);
2084 if (ImporterH_.is_null())
2085 AH_ = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*P11_,
true,*temp,
false,AH_,GetOStream(
Runtime0),
true,
true);
2089 temp2 = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*P11_,
true,*temp,
false,temp2,GetOStream(
Runtime0),
true,
true);
2091 RCP<const Map> map = ImporterH_->getTargetMap()->removeEmptyProcesses();
2092 temp2->removeEmptyProcessesInPlace(map);
2093 if (!temp2.is_null() && temp2->getRowMap().is_null())
2094 temp2 = Teuchos::null;
2098 if (!disable_addon_) {
2102 if (!AH_.is_null() && Addon_Matrix_.is_null()) {
2104 RCP<Teuchos::TimeMonitor> tmAddon = getTimer(
"MueLu RefMaxwell: Build coarse addon matrix");
2106 TEUCHOS_TEST_FOR_EXCEPTION(M0inv_Matrix_==Teuchos::null,std::invalid_argument,
2107 "MueLu::RefMaxwell::formCoarseMatrix(): Inverse of " 2108 "lumped mass matrix required for add-on (i.e. M0inv_Matrix is null)");
2111 RCP<Matrix> Zaux, Z;
2114 Zaux = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M1_Matrix_,
false,*P11_,
false,Zaux,GetOStream(
Runtime0),
true,
true);
2116 Z = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
true,*Zaux,
false,Z,GetOStream(
Runtime0),
true,
true);
2119 if (M0inv_Matrix_->getGlobalMaxNumRowEntries()<=1) {
2122 RCP<Vector> diag = VectorFactory::Build(M0inv_Matrix_->getRowMap());
2123 M0inv_Matrix_->getLocalDiagCopy(*diag);
2125 ArrayRCP<Scalar> diagVals = diag->getDataNonConst(0);
2126 for (
size_t j=0; j < diag->getMap()->getNodeNumElements(); j++) {
2127 diagVals[j] = Teuchos::ScalarTraits<Scalar>::squareroot(diagVals[j]);
2130 if (Z->getRowMap()->isSameAs(*(diag->getMap())))
2131 Z->leftScale(*diag);
2133 RCP<Import> importer = ImportFactory::Build(diag->getMap(),Z->getRowMap());
2134 RCP<Vector> diag2 = VectorFactory::Build(Z->getRowMap());
2135 diag2->doImport(*diag,*importer,Xpetra::INSERT);
2136 Z->leftScale(*diag2);
2138 addon = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*Z,
true,*Z,
false,addon,GetOStream(
Runtime0),
true,
true);
2139 }
else if (parameterList_.get<
bool>(
"rap: triple product",
false) ==
false) {
2142 C2 = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M0inv_Matrix_,
false,*Z,
false,C2,GetOStream(
Runtime0),
true,
true);
2144 addon = Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*Z,
true,*C2,
false,addon,GetOStream(
Runtime0),
true,
true);
2146 addon = MatrixFactory::Build(Z->getDomainMap());
2148 Xpetra::TripleMatrixMultiply<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
2149 MultiplyRAP(*Z,
true, *M0inv_Matrix_,
false, *Z,
false, *addon,
true,
true);
2153 Addon_Matrix_ = addon;
2155 addon = Addon_Matrix_;
2157 if (!AH_.is_null()) {
2160 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*AH_,
false,one,*addon,
false,one,newAH,GetOStream(
Runtime0));
2161 newAH->fillComplete();
2166 if (!AH_.is_null() && !skipFirstLevel_) {
2167 ArrayRCP<bool> AHBCrows;
2168 AHBCrows.resize(AH_->getRowMap()->getNodeNumElements());
2169 size_t dim = Nullspace_->getNumVectors();
2170 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 2172 for (
size_t i = 0; i < BCdomainKokkos_.size(); i++)
2173 for (
size_t k = 0; k < dim; k++)
2174 AHBCrows[i*dim+k] = BCdomainKokkos_(i);
2177 for (
size_t i = 0; i < static_cast<size_t>(BCdomain_.size()); i++)
2178 for (
size_t k = 0; k < dim; k++)
2179 AHBCrows[i*dim+k] = BCdomain_[i];
2180 magnitudeType rowSumTol = parameterList_.get(
"refmaxwell: row sum drop tol (1,1)",-1.0);
2187 if (!AH_.is_null()) {
2193 bool fixZeroDiagonal = !applyBCsToAnodal_;
2194 if (precList11_.isParameter(
"rap: fix zero diagonals"))
2195 fixZeroDiagonal = precList11_.get<
bool>(
"rap: fix zero diagonals");
2197 if (fixZeroDiagonal) {
2199 Scalar replacement = 1.0;
2200 if (precList11_.isType<
magnitudeType>(
"rap: fix zero diagonals threshold"))
2201 threshold = precList11_.get<
magnitudeType>(
"rap: fix zero diagonals threshold");
2202 else if (precList11_.isType<
double>(
"rap: fix zero diagonals threshold"))
2203 threshold = Teuchos::as<magnitudeType>(precList11_.get<
double>(
"rap: fix zero diagonals threshold"));
2204 if (precList11_.isType<
double>(
"rap: fix zero diagonals replacement"))
2205 replacement = Teuchos::as<Scalar>(precList11_.get<
double>(
"rap: fix zero diagonals replacement"));
2206 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(AH_,
true, GetOStream(
Warnings1), threshold, replacement);
2210 size_t dim = Nullspace_->getNumVectors();
2211 AH_->SetFixedBlockSize(dim);
2216 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2218 bool reuse = !SM_Matrix_.is_null();
2219 SM_Matrix_ = SM_Matrix_new;
2220 dump(*SM_Matrix_,
"SM.m");
2226 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2229 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
2233 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu RefMaxwell: residual calculation");
2239 if (implicitTranspose_) {
2241 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu RefMaxwell: restriction coarse (1,1) (implicit)");
2242 P11_->apply(*residual_,*P11res_,Teuchos::TRANS);
2244 if (!allNodesBoundary_) {
2245 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: restriction (2,2) (implicit)");
2246 D0_Matrix_->apply(*residual_,*D0res_,Teuchos::TRANS);
2249 #ifdef MUELU_HAVE_TPETRA 2250 if (D0_T_R11_colMapsMatch_) {
2253 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: restrictions import");
2254 D0TR11Tmp_->doImport(*residual_, *rcp_dynamic_cast<CrsMatrixWrap>(R11_)->getCrsMatrix()->getCrsGraph()->getImporter(), Xpetra::INSERT);
2256 if (!allNodesBoundary_) {
2257 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: restriction (2,2) (explicit)");
2258 rcp_dynamic_cast<TpetraCrsMatrix>(rcp_dynamic_cast<CrsMatrixWrap>(D0_T_Matrix_)->getCrsMatrix())->getTpetra_CrsMatrix()->localApply(toTpetra(*D0TR11Tmp_),toTpetra(*D0res_),Teuchos::NO_TRANS);
2261 RCP<Teuchos::TimeMonitor> tmP11 = getTimer(
"MueLu RefMaxwell: restriction coarse (1,1) (explicit)");
2262 rcp_dynamic_cast<TpetraCrsMatrix>(rcp_dynamic_cast<CrsMatrixWrap>(R11_)->getCrsMatrix())->getTpetra_CrsMatrix()->localApply(toTpetra(*D0TR11Tmp_),toTpetra(*P11res_),Teuchos::NO_TRANS);
2268 RCP<Teuchos::TimeMonitor> tmP11 = getTimer(
"MueLu RefMaxwell: restriction coarse (1,1) (explicit)");
2269 R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
2271 if (!allNodesBoundary_) {
2272 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: restriction (2,2) (explicit)");
2273 D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
2280 RCP<Teuchos::TimeMonitor> tmSubSolves = getTimer(
"MueLu RefMaxwell: subsolves");
2284 if (!ImporterH_.is_null() && !implicitTranspose_) {
2285 RCP<Teuchos::TimeMonitor> tmH = getTimer(
"MueLu RefMaxwell: import coarse (1,1)");
2286 P11resTmp_->beginImport(*P11res_, *ImporterH_, Xpetra::INSERT);
2288 if (!allNodesBoundary_ && !Importer22_.is_null() && !implicitTranspose_) {
2289 RCP<Teuchos::TimeMonitor> tm22 = getTimer(
"MueLu RefMaxwell: import (2,2)");
2290 D0resTmp_->beginImport(*D0res_, *Importer22_, Xpetra::INSERT);
2294 if (!AH_.is_null()) {
2295 if (!ImporterH_.is_null() && !implicitTranspose_)
2296 P11resTmp_->endImport(*P11res_, *ImporterH_, Xpetra::INSERT);
2298 RCP<Teuchos::TimeMonitor> tmH = getTimer(
"MueLu RefMaxwell: solve coarse (1,1)", AH_->getRowMap()->getComm());
2300 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 2301 if (!thyraPrecH_.is_null()) {
2302 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
2303 RCP<Thyra::MultiVectorBase<Scalar> > thyraP11x = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(P11xSubComm_));
2304 RCP<const Thyra::MultiVectorBase<Scalar> > thyraP11res = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(P11resSubComm_);
2305 thyraPrecH_->getUnspecifiedPrecOp()->apply(Thyra::NOTRANS, *thyraP11res, thyraP11x.ptr(), one, zero);
2306 RCP<MultiVector> thyXpP11x = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraP11x, P11x_->getMap()->getComm());
2307 *P11xSubComm_ = *thyXpP11x;
2311 HierarchyH_->Iterate(*P11resSubComm_, *P11xSubComm_, numItersH_,
true);
2315 if (!A22_.is_null()) {
2316 if (!allNodesBoundary_ && !Importer22_.is_null() && !implicitTranspose_)
2317 D0resTmp_->endImport(*D0res_, *Importer22_, Xpetra::INSERT);
2319 RCP<Teuchos::TimeMonitor> tm22 = getTimer(
"MueLu RefMaxwell: solve (2,2)", A22_->getRowMap()->getComm());
2321 #if defined(HAVE_MUELU_STRATIMIKOS) && defined(HAVE_MUELU_THYRA) 2322 if (!thyraPrec22_.is_null()) {
2323 Scalar zero = Teuchos::ScalarTraits<Scalar>::zero();
2324 RCP<Thyra::MultiVectorBase<Scalar> > thyraD0x = Teuchos::rcp_const_cast<Thyra::MultiVectorBase<Scalar> >(Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(D0xSubComm_));
2325 RCP<const Thyra::MultiVectorBase<Scalar> > thyraD0res = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toThyraMultiVector(D0resSubComm_);
2326 thyraPrec22_->getUnspecifiedPrecOp()->apply(Thyra::NOTRANS, *thyraD0res, thyraD0x.ptr(), one, zero);
2327 RCP<MultiVector> thyXpD0x = Xpetra::ThyraUtils<Scalar,LocalOrdinal,GlobalOrdinal,Node>::toXpetra(thyraD0x, D0x_->getMap()->getComm());
2328 *D0xSubComm_ = *thyXpD0x;
2332 Hierarchy22_->Iterate(*D0resSubComm_, *D0xSubComm_, numIters22_,
true);
2335 if (AH_.is_null() && !ImporterH_.is_null() && !implicitTranspose_)
2336 P11resTmp_->endImport(*P11res_, *ImporterH_, Xpetra::INSERT);
2337 if (A22_.is_null() && !allNodesBoundary_ && !Importer22_.is_null() && !implicitTranspose_)
2338 D0resTmp_->endImport(*D0res_, *Importer22_, Xpetra::INSERT);
2341 if (fuseProlongationAndUpdate_) {
2343 RCP<Teuchos::TimeMonitor> tmP11 = getTimer(
"MueLu RefMaxwell: prolongation coarse (1,1) (fused)");
2344 P11_->apply(*P11x_,X,Teuchos::NO_TRANS,one,one);
2347 if (!allNodesBoundary_) {
2348 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: prolongation (2,2) (fused)");
2349 D0_Matrix_->apply(*D0x_,X,Teuchos::NO_TRANS,one,one);
2353 RCP<Teuchos::TimeMonitor> tmP11 = getTimer(
"MueLu RefMaxwell: prolongation coarse (1,1) (unfused)");
2354 P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
2357 if (!allNodesBoundary_) {
2358 RCP<Teuchos::TimeMonitor> tmD0 = getTimer(
"MueLu RefMaxwell: prolongation (2,2) (unfused)");
2359 D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS,one,one);
2363 RCP<Teuchos::TimeMonitor> tmUpdate = getTimer(
"MueLu RefMaxwell: update");
2364 X.update(one, *residual_, one);
2371 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2374 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
2377 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu RefMaxwell: residual calculation");
2379 if (implicitTranspose_)
2380 P11_->apply(*residual_,*P11res_,Teuchos::TRANS);
2382 R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
2386 if (!ImporterH_.is_null() && !implicitTranspose_) {
2387 RCP<Teuchos::TimeMonitor> tmH = getTimer(
"MueLu RefMaxwell: import coarse (1,1)");
2388 P11resTmp_->doImport(*P11res_, *ImporterH_, Xpetra::INSERT);
2390 if (!AH_.is_null()) {
2391 RCP<Teuchos::TimeMonitor> tmH = getTimer(
"MueLu RefMaxwell: solve coarse (1,1)", AH_->getRowMap()->getComm());
2392 HierarchyH_->Iterate(*P11resSubComm_, *P11xSubComm_, numItersH_,
true);
2397 RCP<Teuchos::TimeMonitor> tmUp = getTimer(
"MueLu RefMaxwell: update");
2398 P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
2399 X.update(one, *residual_, one);
2405 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2408 if (allNodesBoundary_)
2411 Scalar one = Teuchos::ScalarTraits<Scalar>::one();
2414 RCP<Teuchos::TimeMonitor> tmRes = getTimer(
"MueLu RefMaxwell: residual calculation");
2416 if (implicitTranspose_)
2417 D0_Matrix_->apply(*residual_,*D0res_,Teuchos::TRANS);
2419 D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
2423 if (!Importer22_.is_null() && !implicitTranspose_) {
2424 RCP<Teuchos::TimeMonitor> tm22 = getTimer(
"MueLu RefMaxwell: import (2,2)");
2425 D0resTmp_->doImport(*D0res_, *Importer22_, Xpetra::INSERT);
2427 if (!A22_.is_null()) {
2428 RCP<Teuchos::TimeMonitor> tm22 = getTimer(
"MueLu RefMaxwell: solve (2,2)", A22_->getRowMap()->getComm());
2429 Hierarchy22_->Iterate(*D0resSubComm_, *D0xSubComm_, numIters22_,
true);
2434 RCP<Teuchos::TimeMonitor> tmUp = getTimer(
"MueLu RefMaxwell: update");
2435 D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS);
2436 X.update(one, *residual_, one);
2442 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2448 RCP<Teuchos::TimeMonitor> tm = getTimer(
"MueLu RefMaxwell: solve");
2451 if (!allEdgesBoundary_ && X.getNumVectors() != P11res_->getNumVectors())
2452 allocateMemory(X.getNumVectors());
2456 RCP<Teuchos::TimeMonitor> tmSm = getTimer(
"MueLu RefMaxwell: smoothing");
2458 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 2459 if (useHiptmairSmoothing_) {
2462 hiptmairPreSmoother_->apply(tRHS, tX);
2466 PreSmoother_->Apply(X, RHS, use_as_preconditioner_);
2470 if(mode_==
"additive")
2471 applyInverseAdditive(RHS,X);
2472 else if(mode_==
"121") {
2476 }
else if(mode_==
"212") {
2480 }
else if(mode_==
"1")
2484 else if(mode_==
"7") {
2488 RCP<Teuchos::TimeMonitor> tmSm = getTimer(
"MueLu RefMaxwell: smoothing");
2490 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 2491 if (useHiptmairSmoothing_) {
2494 hiptmairPreSmoother_->apply(tRHS, tX);
2498 PreSmoother_->Apply(X, RHS,
false);
2503 RCP<Teuchos::TimeMonitor> tmSm = getTimer(
"MueLu RefMaxwell: smoothing");
2505 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 2506 if (useHiptmairSmoothing_)
2510 hiptmairPostSmoother_->apply(tRHS, tX);
2514 PostSmoother_->Apply(X, RHS,
false);
2517 }
else if(mode_==
"none") {
2521 applyInverseAdditive(RHS,X);
2525 RCP<Teuchos::TimeMonitor> tmSm = getTimer(
"MueLu RefMaxwell: smoothing");
2527 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 2528 if (useHiptmairSmoothing_)
2532 hiptmairPostSmoother_->apply(tRHS, tX);
2536 PostSmoother_->Apply(X, RHS,
false);
2542 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2548 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2551 const Teuchos::RCP<Matrix> & Ms_Matrix,
2552 const Teuchos::RCP<Matrix> & M0inv_Matrix,
2553 const Teuchos::RCP<Matrix> & M1_Matrix,
2554 const Teuchos::RCP<MultiVector> & Nullspace,
2555 const Teuchos::RCP<RealValuedMultiVector> & Coords,
2556 Teuchos::ParameterList& List)
2559 TEUCHOS_ASSERT(D0_Matrix!=Teuchos::null);
2560 TEUCHOS_ASSERT(Ms_Matrix!=Teuchos::null);
2561 TEUCHOS_ASSERT(M1_Matrix!=Teuchos::null);
2563 #ifdef HAVE_MUELU_DEBUG 2564 TEUCHOS_ASSERT(M1_Matrix->getDomainMap()->isSameAs(*M1_Matrix->getRangeMap()));
2565 TEUCHOS_ASSERT(M1_Matrix->getDomainMap()->isSameAs(*M1_Matrix->getRowMap()));
2566 TEUCHOS_ASSERT(M1_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getRangeMap()));
2568 TEUCHOS_ASSERT(Ms_Matrix->getDomainMap()->isSameAs(*Ms_Matrix->getRangeMap()));
2569 TEUCHOS_ASSERT(Ms_Matrix->getDomainMap()->isSameAs(*Ms_Matrix->getRowMap()));
2570 TEUCHOS_ASSERT(Ms_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getRangeMap()));
2572 TEUCHOS_ASSERT(D0_Matrix->getRangeMap()->isSameAs(*D0_Matrix->getRowMap()));
2574 if (!M0inv_Matrix) {
2575 TEUCHOS_ASSERT(M0inv_Matrix->getDomainMap()->isSameAs(*M0inv_Matrix->getRowMap()));
2576 TEUCHOS_ASSERT(M0inv_Matrix->getDomainMap()->isSameAs(*M0inv_Matrix->getRangeMap()));
2577 TEUCHOS_ASSERT(M0inv_Matrix->getDomainMap()->isSameAs(*D0_Matrix->getDomainMap()));
2581 HierarchyH_ = Teuchos::null;
2582 Hierarchy22_ = Teuchos::null;
2583 PreSmoother_ = Teuchos::null;
2584 PostSmoother_ = Teuchos::null;
2585 disable_addon_ =
false;
2589 setParameters(List);
2591 if (D0_Matrix->getRowMap()->lib() == Xpetra::UseTpetra) {
2596 RCP<Matrix> D0copy = MatrixFactory::Build(D0_Matrix->getRowMap(), D0_Matrix->getColMap(), 0);
2597 RCP<CrsMatrix> D0copyCrs = rcp_dynamic_cast<CrsMatrixWrap>(D0copy,
true)->getCrsMatrix();
2598 ArrayRCP<const size_t> D0rowptr_RCP;
2599 ArrayRCP<const LO> D0colind_RCP;
2600 ArrayRCP<const SC> D0vals_RCP;
2601 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getAllValues(D0rowptr_RCP,
2605 ArrayRCP<size_t> D0copyrowptr_RCP;
2606 ArrayRCP<LO> D0copycolind_RCP;
2607 ArrayRCP<SC> D0copyvals_RCP;
2608 D0copyCrs->allocateAllValues(D0vals_RCP.size(),D0copyrowptr_RCP,D0copycolind_RCP,D0copyvals_RCP);
2609 D0copyrowptr_RCP.deepCopy(D0rowptr_RCP());
2610 D0copycolind_RCP.deepCopy(D0colind_RCP());
2611 D0copyvals_RCP.deepCopy(D0vals_RCP());
2612 D0copyCrs->setAllValues(D0copyrowptr_RCP,
2615 D0copyCrs->expertStaticFillComplete(D0_Matrix->getDomainMap(), D0_Matrix->getRangeMap(),
2616 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getImporter(),
2617 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix,
true)->getCrsMatrix()->getCrsGraph()->getExporter());
2618 D0_Matrix_ = D0copy;
2620 D0_Matrix_ = MatrixFactory::BuildCopy(D0_Matrix);
2623 M0inv_Matrix_ = M0inv_Matrix;
2624 Ms_Matrix_ = Ms_Matrix;
2625 M1_Matrix_ = M1_Matrix;
2627 Nullspace_ = Nullspace;
2629 dump(*D0_Matrix_,
"D0_clean.m");
2630 dump(*Ms_Matrix_,
"Ms.m");
2631 dump(*M1_Matrix_,
"M1.m");
2632 if (!M0inv_Matrix_.is_null()) dump(*M0inv_Matrix_,
"M0inv.m");
2633 if (!Coords_.is_null()) dumpCoords(*Coords_,
"coords.m");
2638 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
2640 describe(Teuchos::FancyOStream& out,
const Teuchos::EVerbosityLevel )
const {
2642 std::ostringstream oss;
2644 RCP<const Teuchos::Comm<int> > comm = SM_Matrix_->getDomainMap()->getComm();
2649 root = comm->getRank();
2654 reduceAll(*comm, Teuchos::REDUCE_MAX, root, Teuchos::ptr(&actualRoot));
2659 oss <<
"\n--------------------------------------------------------------------------------\n" <<
2660 "--- RefMaxwell Summary ---\n" 2661 "--------------------------------------------------------------------------------" << std::endl;
2667 SM_Matrix_->getRowMap()->getComm()->barrier();
2669 numRows = SM_Matrix_->getGlobalNumRows();
2670 nnz = SM_Matrix_->getGlobalNumEntries();
2672 Xpetra::global_size_t tt = numRows;
2673 int rowspacer = 3;
while (tt != 0) { tt /= 10; rowspacer++; }
2675 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
2677 oss <<
"block " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " << std::setw(9) <<
" nnz/row" << std::endl;
2678 oss <<
"(1, 1)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
2680 if (!A22_.is_null()) {
2682 numRows = A22_->getGlobalNumRows();
2683 nnz = A22_->getGlobalNumEntries();
2685 oss <<
"(2, 2)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
2690 #if defined(MUELU_REFMAXWELL_CAN_USE_HIPTMAIR) 2691 if (useHiptmairSmoothing_) {
2692 if (hiptmairPreSmoother_ != null && hiptmairPreSmoother_ == hiptmairPostSmoother_)
2693 oss <<
"Smoother both : " << hiptmairPreSmoother_->description() << std::endl;
2695 oss <<
"Smoother pre : " 2696 << (hiptmairPreSmoother_ != null ? hiptmairPreSmoother_->description() :
"no smoother") << std::endl;
2697 oss <<
"Smoother post : " 2698 << (hiptmairPostSmoother_ != null ? hiptmairPostSmoother_->description() :
"no smoother") << std::endl;
2704 if (PreSmoother_ != null && PreSmoother_ == PostSmoother_)
2705 oss <<
"Smoother both : " << PreSmoother_->description() << std::endl;
2707 oss <<
"Smoother pre : " 2708 << (PreSmoother_ != null ? PreSmoother_->description() :
"no smoother") << std::endl;
2709 oss <<
"Smoother post : " 2710 << (PostSmoother_ != null ? PostSmoother_->description() :
"no smoother") << std::endl;
2715 std::string outstr = oss.str();
2718 RCP<const Teuchos::MpiComm<int> > mpiComm = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
2719 MPI_Comm rawComm = (*mpiComm->getRawMpiComm())();
2721 int strLength = outstr.size();
2722 MPI_Bcast(&strLength, 1, MPI_INT, root, rawComm);
2723 if (comm->getRank() != root)
2724 outstr.resize(strLength);
2725 MPI_Bcast(&outstr[0], strLength, MPI_CHAR, root, rawComm);
2730 if (!HierarchyH_.is_null())
2731 HierarchyH_->describe(out, GetVerbLevel());
2733 if (!Hierarchy22_.is_null())
2734 Hierarchy22_->describe(out, GetVerbLevel());
2738 std::ostringstream oss2;
2740 oss2 <<
"Sub-solver distribution over ranks" << std::endl;
2741 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;
2743 int numProcs = comm->getSize();
2745 RCP<const Teuchos::MpiComm<int> > tmpic = rcp_dynamic_cast<
const Teuchos::MpiComm<int> >(comm);
2746 TEUCHOS_TEST_FOR_EXCEPTION(tmpic == Teuchos::null,
Exceptions::RuntimeError,
"Cannot cast base Teuchos::Comm to Teuchos::MpiComm object.");
2747 RCP<const Teuchos::OpaqueWrapper<MPI_Comm> > rawMpiComm = tmpic->getRawMpiComm();
2753 if (!A22_.is_null())
2755 std::vector<char> states(numProcs, 0);
2757 MPI_Gather(&status, 1, MPI_CHAR, &states[0], 1, MPI_CHAR, 0, *rawMpiComm);
2759 states.push_back(status);
2762 int rowWidth = std::min(Teuchos::as<int>(ceil(sqrt(numProcs))), 100);
2763 for (
int proc = 0; proc < numProcs; proc += rowWidth) {
2764 for (
int j = 0; j < rowWidth; j++)
2765 if (proc + j < numProcs)
2766 if (states[proc+j] == 0)
2768 else if (states[proc+j] == 1)
2770 else if (states[proc+j] == 2)
2777 oss2 <<
" " << proc <<
":" << std::min(proc + rowWidth, numProcs) - 1 << std::endl;
2789 #define MUELU_REFMAXWELL_SHORT 2790 #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 HAVE_MUELU_KOKKOS_REFACTOR
#define MueLu_maxAll(rcpComm, in, out)
static void ApplyRowSumCriterion(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude rowSumTol, Teuchos::ArrayRCP< bool > &dirichletRows)
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.
static void SetDefaultVerbLevel(const VerbLevel defaultVerbLevel)
Set the default (global) verbosity level.
Class that encapsulates external library smoothers.
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 setFineLevelSmoother()
Set the fine level smoother.
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)
Utility functions for Maxwell.
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 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 void setMatvecParams(Matrix &A, RCP< ParameterList > matvecParams)
Sets matvec params on a matrix.
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
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.
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.
Factory for creating a graph based 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 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
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.
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
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.