46 #ifndef MUELU_REFMAXWELL_DEF_HPP 47 #define MUELU_REFMAXWELL_DEF_HPP 51 #include "Xpetra_Map.hpp" 52 #include "Xpetra_MatrixMatrix.hpp" 53 #include "Xpetra_TripleMatrixMultiply.hpp" 54 #include "Xpetra_CrsMatrixUtils.hpp" 55 #include "Xpetra_MatrixUtils.hpp" 59 #include "MueLu_AmalgamationFactory.hpp" 60 #include "MueLu_RAPFactory.hpp" 61 #include "MueLu_ThresholdAFilterFactory.hpp" 62 #include "MueLu_TransPFactory.hpp" 63 #include "MueLu_SmootherFactory.hpp" 65 #include "MueLu_CoalesceDropFactory.hpp" 66 #include "MueLu_CoarseMapFactory.hpp" 67 #include "MueLu_CoordinatesTransferFactory.hpp" 68 #include "MueLu_UncoupledAggregationFactory.hpp" 69 #include "MueLu_TentativePFactory.hpp" 70 #include "MueLu_Utilities.hpp" 72 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 73 #include "MueLu_CoalesceDropFactory_kokkos.hpp" 74 #include "MueLu_CoarseMapFactory_kokkos.hpp" 75 #include "MueLu_CoordinatesTransferFactory_kokkos.hpp" 76 #include "MueLu_UncoupledAggregationFactory_kokkos.hpp" 77 #include "MueLu_TentativePFactory_kokkos.hpp" 78 #include "MueLu_Utilities_kokkos.hpp" 79 #include <Kokkos_Core.hpp> 80 #include <KokkosSparse_CrsMatrix.hpp> 85 #include "MueLu_MLParameterListInterpreter.hpp" 86 #include "MueLu_ParameterListInterpreter.hpp" 91 #ifdef HAVE_MUELU_CUDA 92 #include "cuda_profiler_api.h" 98 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
100 return SM_Matrix_->getDomainMap();
104 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
106 return SM_Matrix_->getRangeMap();
110 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
113 parameterList_ = list;
114 disable_addon_ = list.get(
"refmaxwell: disable addon",
true);
115 mode_ = list.get(
"refmaxwell: mode",
"additive");
116 use_as_preconditioner_ = list.get<
bool>(
"refmaxwell: use as preconditioner");
117 dump_matrices_ = list.get(
"refmaxwell: dump matrices",
false);
119 if(list.isSublist(
"refmaxwell: 11list"))
120 precList11_ = list.sublist(
"refmaxwell: 11list");
122 if(list.isSublist(
"refmaxwell: 22list"))
123 precList22_ = list.sublist(
"refmaxwell: 22list");
125 if(list.isSublist(
"smoother: params")) {
126 smootherList_ = list.sublist(
"smoother: params");
129 #if !defined(HAVE_MUELU_KOKKOS_REFACTOR) 131 #elif defined(HAVE_MUELU_KOKKOS_REFACTOR_USE_BY_DEFAULT) 132 useKokkos_ = list.get(
"use kokkos refactor",
true);
134 useKokkos_ = list.get(
"use kokkos refactor",
false);
140 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
144 #ifdef HAVE_MUELU_CUDA 145 if (parameterList_.get<
bool>(
"refmaxwell: cuda profile setup",
false)) cudaProfilerStart();
148 Teuchos::TimeMonitor tmCompute(*Teuchos::TimeMonitor::getNewTimer(
"MueLu RefMaxwell: compute"));
150 bool defaultFilter =
false;
155 if (parameterList_.get<
bool>(
"refmaxwell: filter D0",
true) && D0_Matrix_->getNodeMaxNumRowEntries()>2) {
159 fineLevel.
Set(
"A",D0_Matrix_);
160 fineLevel.
setlib(D0_Matrix_->getDomainMap()->lib());
163 fineLevel.
Request(
"A",ThreshFact.get());
164 ThreshFact->Build(fineLevel);
165 D0_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
168 defaultFilter =
true;
171 if (parameterList_.get<
bool>(
"refmaxwell: filter SM", defaultFilter)) {
172 RCP<Vector> diag = VectorFactory::Build(SM_Matrix_->getRowMap());
174 SM_Matrix_->getLocalDiagCopy(*diag);
180 fineLevel.
Set(
"A",SM_Matrix_);
181 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
183 fineLevel.
Request(
"A",ThreshFact.get());
184 ThreshFact->Build(fineLevel);
185 SM_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
188 if (parameterList_.get<
bool>(
"refmaxwell: filter M1", defaultFilter)) {
189 RCP<Vector> diag = VectorFactory::Build(M1_Matrix_->getRowMap());
191 M1_Matrix_->getLocalDiagCopy(*diag);
197 fineLevel.
Set(
"A",M1_Matrix_);
198 fineLevel.
setlib(M1_Matrix_->getDomainMap()->lib());
200 fineLevel.
Request(
"A",ThreshFact.get());
201 ThreshFact->Build(fineLevel);
202 M1_Matrix_ = fineLevel.
Get< RCP<Matrix> >(
"A",ThreshFact.get());
209 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 223 if(Nullspace_ != Teuchos::null) {
226 else if(Nullspace_ == Teuchos::null && Coords_ != Teuchos::null) {
228 typedef typename RealValuedMultiVector::scalar_type realScalarType;
229 typedef typename Teuchos::ScalarTraits<realScalarType>::magnitudeType realMagnitudeType;
230 Teuchos::Array<realMagnitudeType> norms(Coords_->getNumVectors());
231 Coords_->norm2(norms);
232 for (
size_t i=0;i<Coords_->getNumVectors();i++)
233 norms[i] = ((realMagnitudeType)1.0)/norms[i];
234 Coords_->scale(norms());
235 Nullspace_ = MultiVectorFactory::Build(SM_Matrix_->getRowMap(),Coords_->getNumVectors());
238 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 239 RCP<MultiVector> CoordsSC;
241 CoordsSC = Utilities_kokkos::RealValuedToScalarMultiVector(Coords_);
247 D0_Matrix_->apply(*CoordsSC,*Nullspace_);
250 GetOStream(
Errors) <<
"MueLu::RefMaxwell::compute(): either the nullspace or the nodal coordinates must be provided." << std::endl;
254 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 264 Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string(
"D0_clean.mat"), *D0_Matrix_);
267 if(P11_==Teuchos::null) {
269 Level fineLevel, coarseLevel;
275 fineLevel.
Set(
"A",M1_Matrix_);
276 coarseLevel.
Set(
"P",D0_Matrix_);
277 coarseLevel.
setlib(M1_Matrix_->getDomainMap()->lib());
278 fineLevel.
setlib(M1_Matrix_->getDomainMap()->lib());
279 coarseLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
280 fineLevel.setObjectLabel(
"RefMaxwell (1,1) A_nodal");
282 RCP<RAPFactory> rapFact = rcp(
new RAPFactory());
283 Teuchos::ParameterList rapList = *(rapFact->GetValidParameterList());
284 rapList.set(
"transpose: use implicit", parameterList_.get<
bool>(
"transpose: use implicit",
false));
285 rapList.set(
"rap: fix zero diagonals", parameterList_.get<
bool>(
"rap: fix zero diagonals",
true));
286 rapList.set(
"rap: triple product", parameterList_.get<
bool>(
"rap: triple product",
false));
287 rapFact->SetParameterList(rapList);
289 RCP<TransPFactory> transPFactory;
290 if (!parameterList_.get<
bool>(
"transpose: use implicit",
false)) {
292 rapFact->SetFactory(
"R", transPFactory);
295 coarseLevel.
Request(
"A", rapFact.get());
297 A_nodal_Matrix_ = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
300 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): building special prolongator" << std::endl;
303 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 305 R11_ = Utilities_kokkos::Transpose(*P11_);
317 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 321 if (useKokkos_ && precList11_.isParameter(
"aggregation: drop tol") && precList11_.get<
double>(
"aggregation: drop tol") != 0.0) {
322 GetOStream(
Warnings0) <<
"RefMaxwell::compute(): Setting \"aggregation: drop tol\". to 0.0, since CoalesceDropFactory_kokkos does not " 323 <<
"support BlockSize > 1 and drop tol != 0.0" << std::endl;
324 precList11_.set(
"aggregation: drop tol", 0.0);
331 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): nuking BC edges of D0" << std::endl;
333 D0_Matrix_->resumeFill();
335 Scalar replaceWith = Teuchos::ScalarTraits<SC>::zero();
336 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 348 D0_Matrix_->fillComplete(D0_Matrix_->getDomainMap(),D0_Matrix_->getRangeMap());
352 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): building MG for (2,2)-block" << std::endl;
355 Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer(
"MueLu RefMaxwell: Build A22"));
357 Level fineLevel, coarseLevel;
363 fineLevel.
Set(
"A",SM_Matrix_);
364 coarseLevel.
Set(
"P",D0_Matrix_);
365 coarseLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
366 fineLevel.
setlib(SM_Matrix_->getDomainMap()->lib());
367 coarseLevel.setObjectLabel(
"RefMaxwell (2,2)");
368 fineLevel.setObjectLabel(
"RefMaxwell (2,2)");
370 RCP<RAPFactory> rapFact = rcp(
new RAPFactory());
371 Teuchos::ParameterList rapList = *(rapFact->GetValidParameterList());
372 rapList.set(
"transpose: use implicit",
false);
373 rapList.set(
"rap: fix zero diagonals", parameterList_.get<
bool>(
"rap: fix zero diagonals",
true));
374 rapList.set(
"rap: triple product", parameterList_.get<
bool>(
"rap: triple product",
false));
375 rapFact->SetParameterList(rapList);
377 RCP<TransPFactory> transPFactory;
379 rapFact->SetFactory(
"R", transPFactory);
381 coarseLevel.
Request(
"R", transPFactory.get());
382 coarseLevel.
Request(
"A", rapFact.get());
383 A22_ = coarseLevel.
Get< RCP<Matrix> >(
"A", rapFact.get());
384 A22_->setObjectLabel(
"RefMaxwell (2,2)");
385 D0_T_Matrix_ = coarseLevel.
Get< RCP<Matrix> >(
"R", transPFactory.get());
392 if (parameterList_.isType<std::string>(
"smoother: type") &&
393 parameterList_.get<std::string>(
"smoother: type") ==
"hiptmair" &&
394 SM_Matrix_->getDomainMap()->lib() == Xpetra::UseTpetra &&
395 A22_->getDomainMap()->lib() == Xpetra::UseTpetra &&
396 D0_Matrix_->getDomainMap()->lib() == Xpetra::UseTpetra) {
397 #if defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || (defined(HAVE_MUELU_INST_DOUBLE_INT_INT))) 398 Teuchos::ParameterList hiptmairPreList, hiptmairPostList, smootherPreList, smootherPostList;
400 if (smootherList_.isSublist(
"smoother: pre params"))
401 smootherPreList = smootherList_.sublist(
"smoother: pre params");
402 else if (smootherList_.isSublist(
"smoother: params"))
403 smootherPreList = smootherList_.sublist(
"smoother: params");
404 hiptmairPreList.set(
"hiptmair: smoother type 1",
405 smootherPreList.get<std::string>(
"hiptmair: smoother type 1",
"CHEBYSHEV"));
406 hiptmairPreList.set(
"hiptmair: smoother type 2",
407 smootherPreList.get<std::string>(
"hiptmair: smoother type 2",
"CHEBYSHEV"));
408 if(smootherPreList.isSublist(
"hiptmair: smoother list 1"))
409 hiptmairPreList.set(
"hiptmair: smoother list 1", smootherPreList.sublist(
"hiptmair: smoother list 1"));
410 if(smootherPreList.isSublist(
"hiptmair: smoother list 2"))
411 hiptmairPreList.set(
"hiptmair: smoother list 2", smootherPreList.sublist(
"hiptmair: smoother list 2"));
412 hiptmairPreList.set(
"hiptmair: pre or post",
413 smootherPreList.get<std::string>(
"hiptmair: pre or post",
"pre"));
414 hiptmairPreList.set(
"hiptmair: zero starting solution",
415 smootherPreList.get<
bool>(
"hiptmair: zero starting solution",
true));
417 if (smootherList_.isSublist(
"smoother: post params"))
418 smootherPostList = smootherList_.sublist(
"smoother: post params");
419 else if (smootherList_.isSublist(
"smoother: params"))
420 smootherPostList = smootherList_.sublist(
"smoother: params");
421 hiptmairPostList.set(
"hiptmair: smoother type 1",
422 smootherPostList.get<std::string>(
"hiptmair: smoother type 1",
"CHEBYSHEV"));
423 hiptmairPostList.set(
"hiptmair: smoother type 2",
424 smootherPostList.get<std::string>(
"hiptmair: smoother type 2",
"CHEBYSHEV"));
425 if(smootherPostList.isSublist(
"hiptmair: smoother list 1"))
426 hiptmairPostList.set(
"hiptmair: smoother list 1", smootherPostList.sublist(
"hiptmair: smoother list 1"));
427 if(smootherPostList.isSublist(
"hiptmair: smoother list 2"))
428 hiptmairPostList.set(
"hiptmair: smoother list 2", smootherPostList.sublist(
"hiptmair: smoother list 2"));
429 hiptmairPostList.set(
"hiptmair: pre or post",
430 smootherPostList.get<std::string>(
"hiptmair: pre or post",
"post"));
431 hiptmairPostList.set(
"hiptmair: zero starting solution",
432 smootherPostList.get<
bool>(
"hiptmair: zero starting solution",
false));
434 typedef Tpetra::RowMatrix<SC, LO, GO, NO> TROW;
440 hiptmairPreSmoother_ -> setParameters(hiptmairPreList);
441 hiptmairPreSmoother_ -> initialize();
442 hiptmairPreSmoother_ -> compute();
444 hiptmairPostSmoother_ -> setParameters(hiptmairPostList);
445 hiptmairPostSmoother_ -> initialize();
446 hiptmairPostSmoother_ -> compute();
447 useHiptmairSmoothing_ =
true;
449 throw(Xpetra::Exceptions::RuntimeError(
"MueLu must be compiled with Ifpack2 for Hiptmair smoothing."));
450 #endif // defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || defined(HAVE_MUELU_INST_DOUBLE_INT_INT)) 452 if (parameterList_.isType<std::string>(
"smoother: pre type") && parameterList_.isType<std::string>(
"smoother: post type")) {
453 std::string preSmootherType = parameterList_.get<std::string>(
"smoother: pre type");
454 std::string postSmootherType = parameterList_.get<std::string>(
"smoother: post type");
456 Teuchos::ParameterList preSmootherList, postSmootherList;
457 if (parameterList_.isSublist(
"smoother: pre params"))
458 preSmootherList = parameterList_.sublist(
"smoother: pre params");
459 if (parameterList_.isSublist(
"smoother: post params"))
460 postSmootherList = parameterList_.sublist(
"smoother: post params");
463 RCP<MueLu::FactoryManagerBase> factoryHandler = rcp(
new FactoryManager());
466 level.setObjectLabel(
"RefMaxwell (1,1)");
467 level.
Set(
"A",SM_Matrix_);
468 level.
setlib(SM_Matrix_->getDomainMap()->lib());
470 RCP<SmootherPrototype> preSmootherPrototype = rcp(
new TrilinosSmoother(preSmootherType, preSmootherList));
471 RCP<SmootherFactory> preSmootherFact = rcp(
new SmootherFactory(preSmootherPrototype));
473 RCP<SmootherPrototype> postSmootherPrototype = rcp(
new TrilinosSmoother(postSmootherType, postSmootherList));
474 RCP<SmootherFactory> postSmootherFact = rcp(
new SmootherFactory(postSmootherPrototype));
476 level.
Request(
"PreSmoother",preSmootherFact.get());
477 preSmootherFact->Build(level);
478 PreSmoother_ = level.
Get<RCP<SmootherBase> >(
"PreSmoother",preSmootherFact.get());
480 level.
Request(
"PostSmoother",postSmootherFact.get());
481 postSmootherFact->Build(level);
482 PostSmoother_ = level.
Get<RCP<SmootherBase> >(
"PostSmoother",postSmootherFact.get());
484 std::string smootherType = parameterList_.get<std::string>(
"smoother: type",
"CHEBYSHEV");
486 RCP<MueLu::FactoryManagerBase> factoryHandler = rcp(
new FactoryManager());
487 level.SetFactoryManager(factoryHandler);
489 level.setObjectLabel(
"RefMaxwell (1,1)");
490 level.Set(
"A",SM_Matrix_);
491 level.setlib(SM_Matrix_->getDomainMap()->lib());
492 RCP<SmootherPrototype> smootherPrototype = rcp(
new TrilinosSmoother(smootherType, smootherList_));
493 RCP<SmootherFactory> SmootherFact = rcp(
new SmootherFactory(smootherPrototype));
494 level.Request(
"PreSmoother",SmootherFact.get());
495 SmootherFact->Build(level);
496 PreSmoother_ = level.Get<RCP<SmootherBase> >(
"PreSmoother",SmootherFact.get());
497 PostSmoother_ = PreSmoother_;
499 useHiptmairSmoothing_ =
false;
504 P11res_ = MultiVectorFactory::Build(P11_->getDomainMap(),1);
505 P11x_ = MultiVectorFactory::Build(P11_->getDomainMap(),1);
506 D0res_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),1);
507 D0x_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),1);
508 residual_ = MultiVectorFactory::Build(SM_Matrix_->getDomainMap(),1);
510 #ifdef HAVE_MUELU_CUDA 511 if (parameterList_.get<
bool>(
"refmaxwell: cuda profile setup",
false)) cudaProfilerStop();
516 if (dump_matrices_) {
517 GetOStream(
Runtime0) <<
"RefMaxwell::compute(): dumping data" << std::endl;
518 Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string(
"SM.mat"), *SM_Matrix_);
519 Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string(
"M1.mat"), *M1_Matrix_);
520 Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string(
"M0inv.mat"), *M0inv_Matrix_);
521 #ifndef HAVE_MUELU_KOKKOS_REFACTOR 522 std::ofstream outBCrows(
"BCrows.mat");
523 std::copy(BCrows_.begin(), BCrows_.end(), std::ostream_iterator<LO>(outBCrows,
"\n"));
524 std::ofstream outBCcols(
"BCcols.mat");
525 std::copy(BCcols_.begin(), BCcols_.end(), std::ostream_iterator<LO>(outBCcols,
"\n"));
528 std::ofstream outBCrows(
"BCrows.mat");
529 auto BCrows = Kokkos::create_mirror_view (BCrowsKokkos_);
530 Kokkos::deep_copy(BCrows , BCrowsKokkos_);
531 for (
size_t i = 0; i < BCrows.size(); i++)
532 outBCrows << BCrows[i] <<
"\n";
534 std::ofstream outBCcols(
"BCcols.mat");
535 auto BCcols = Kokkos::create_mirror_view (BCcolsKokkos_);
536 Kokkos::deep_copy(BCcols , BCcolsKokkos_);
537 for (
size_t i = 0; i < BCcols.size(); i++)
538 outBCcols << BCcols[i] <<
"\n";
540 std::ofstream outBCrows(
"BCrows.mat");
541 std::copy(BCrows_.begin(), BCrows_.end(), std::ostream_iterator<LO>(outBCrows,
"\n"));
542 std::ofstream outBCcols(
"BCcols.mat");
543 std::copy(BCcols_.begin(), BCcols_.end(), std::ostream_iterator<LO>(outBCcols,
"\n"));
546 Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string(
"nullspace.mat"), *Nullspace_);
547 if (Coords_ != Teuchos::null)
548 Xpetra::IO<double, LO, GlobalOrdinal, Node>::Write(std::string(
"coords.mat"), *Coords_);
549 Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string(
"D0_nuked.mat"), *D0_Matrix_);
550 Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string(
"A_nodal.mat"), *A_nodal_Matrix_);
551 Xpetra::IO<SC, LO, GO, NO>::Write(std::string(
"P11.mat"), *P11_);
552 Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string(
"AH.mat"), *AH_);
553 Xpetra::IO<SC, LO, GlobalOrdinal, Node>::Write(std::string(
"A22.mat"), *A22_);
558 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
571 const SC SC_ZERO = Teuchos::ScalarTraits<SC>::zero();
572 const SC SC_ONE = Teuchos::ScalarTraits<SC>::one();
573 size_t dim = Nullspace_->getNumVectors();
574 size_t numLocalRows = SM_Matrix_->getNodeNumRows();
579 bool read_P_from_file = parameterList_.get(
"refmaxwell: read_P_from_file",
false);
580 if (read_P_from_file) {
583 std::string P_filename = parameterList_.get(
"refmaxwell: P_filename",std::string(
"P.mat"));
584 P_nodal = Xpetra::IO<SC, LO, GO, NO>::Read(P_filename, A_nodal_Matrix_->getDomainMap());
586 Level fineLevel, coarseLevel;
592 fineLevel.
Set(
"A",A_nodal_Matrix_);
593 fineLevel.
Set(
"Coordinates",Coords_);
594 fineLevel.
Set(
"DofsPerNode",1);
595 coarseLevel.
setlib(A_nodal_Matrix_->getDomainMap()->lib());
596 fineLevel.
setlib(A_nodal_Matrix_->getDomainMap()->lib());
597 coarseLevel.setObjectLabel(
"RefMaxwell (1,1)");
598 fineLevel.setObjectLabel(
"RefMaxwell (1,1)");
600 LocalOrdinal NSdim = 1;
601 RCP<MultiVector> nullSpace = MultiVectorFactory::Build(A_nodal_Matrix_->getRowMap(),NSdim);
602 nullSpace->putScalar(SC_ONE);
603 fineLevel.
Set(
"Nullspace",nullSpace);
606 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 607 RCP<Factory> dropFact, UncoupledAggFact, coarseMapFact, TentativePFact, Tfact;
609 dropFact = rcp(
new CoalesceDropFactory_kokkos());
610 UncoupledAggFact = rcp(
new UncoupledAggregationFactory_kokkos());
611 coarseMapFact = rcp(
new CoarseMapFactory_kokkos());
612 TentativePFact = rcp(
new TentativePFactory_kokkos());
613 Tfact = rcp(
new CoordinatesTransferFactory_kokkos());
628 dropFact->SetFactory(
"UnAmalgamationInfo", amalgFact);
629 double dropTol = parameterList_.get(
"aggregation: drop tol",0.0);
630 dropFact->SetParameter(
"aggregation: drop tol",Teuchos::ParameterEntry(dropTol));
632 UncoupledAggFact->SetFactory(
"Graph", dropFact);
634 coarseMapFact->SetFactory(
"Aggregates", UncoupledAggFact);
636 TentativePFact->SetFactory(
"Aggregates", UncoupledAggFact);
637 TentativePFact->SetFactory(
"UnAmalgamationInfo", amalgFact);
638 TentativePFact->SetFactory(
"CoarseMap", coarseMapFact);
640 Tfact->SetFactory(
"Aggregates", UncoupledAggFact);
641 Tfact->SetFactory(
"CoarseMap", coarseMapFact);
643 coarseLevel.
Request(
"P",TentativePFact.get());
644 coarseLevel.
Request(
"Coordinates",Tfact.get());
646 coarseLevel.
Get(
"P",P_nodal,TentativePFact.get());
647 coarseLevel.
Get(
"Coordinates",CoordsH_,Tfact.get());
650 Xpetra::IO<SC, LO, GO, NO>::Write(std::string(
"P_nodal.mat"), *P_nodal);
652 RCP<CrsMatrix> D0Crs = rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix();
655 RCP<CrsMatrix> P_nodal_imported;
656 int numProcs = P_nodal->getDomainMap()->getComm()->getSize();
658 RCP<CrsMatrixWrap> P_nodal_temp;
659 RCP<const Map> targetMap = D0Crs->getColMap();
660 P_nodal_temp = Teuchos::rcp(
new CrsMatrixWrap(targetMap, 0));
661 RCP<const Import> importer = D0Crs->getCrsGraph()->getImporter();
662 P_nodal_temp->doImport(*P_nodal, *importer, Xpetra::INSERT);
663 P_nodal_temp->fillComplete(rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix()->getDomainMap(),
664 rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix()->getRangeMap());
665 P_nodal_imported = P_nodal_temp->getCrsMatrix();
667 Xpetra::IO<SC, LO, GO, NO>::Write(std::string(
"P_nodal_imported.mat"), *P_nodal_temp);
669 P_nodal_imported = rcp_dynamic_cast<CrsMatrixWrap>(P_nodal)->getCrsMatrix();
671 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 674 using ATS = Kokkos::ArithTraits<SC>;
675 using range_type = Kokkos::RangePolicy<LO, typename NO::execution_space>;
677 typedef typename Matrix::local_matrix_type KCRS;
678 typedef typename KCRS::device_type device_t;
679 typedef typename KCRS::StaticCrsGraphType graph_t;
680 typedef typename graph_t::row_map_type::non_const_type lno_view_t;
681 typedef typename graph_t::entries_type::non_const_type lno_nnz_view_t;
682 typedef typename KCRS::values_type::non_const_type scalar_view_t;
685 auto localP = P_nodal_imported->getLocalMatrix();
686 auto localD0 = D0_Matrix_->getLocalMatrix();
691 std::string defaultAlgo =
"mat-mat";
692 std::string algo = parameterList_.get(
"refmaxwell: prolongator compute algorithm",defaultAlgo);
694 if (algo ==
"mat-mat") {
695 Teuchos::RCP<Matrix> D0_P_nodal = MatrixFactory::Build(SM_Matrix_->getRowMap(),0);
696 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
false,*P_nodal,
false,*D0_P_nodal,
true,
true);
699 auto localD0P = D0_P_nodal->getLocalMatrix();
702 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
703 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
705 lno_view_t P11rowptr(
"P11_rowptr", numLocalRows+1);
706 lno_nnz_view_t P11colind(
"P11_colind",dim*localD0P.graph.entries.size());
707 scalar_view_t P11vals(
"P11_vals",dim*localD0P.graph.entries.size());
710 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_adjustRowptr", range_type(0,numLocalRows+1),
711 KOKKOS_LAMBDA(
const size_t i) {
712 P11rowptr(i) = dim*localD0P.graph.row_map(i);
716 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_adjustColind", range_type(0,localD0P.graph.entries.size()),
717 KOKKOS_LAMBDA(
const size_t jj) {
718 for (
size_t k = 0; k < dim; k++) {
719 P11colind(dim*jj+k) = dim*localD0P.graph.entries(jj)+k;
720 P11vals(dim*jj+k) = SC_ZERO;
724 auto localNullspace = Nullspace_->template getLocalView<device_t>();
727 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
731 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
733 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
735 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_enterValues_D0wZeros", range_type(0,numLocalRows),
736 KOKKOS_LAMBDA(
const size_t i) {
737 for (
size_t ll = localD0.graph.row_map(i); ll < localD0.graph.row_map(i+1); ll++) {
738 LO l = localD0.graph.entries(ll);
739 SC p = localD0.values(ll);
740 if (ATS::magnitude(p) < tol)
742 for (
size_t jj = localP.graph.row_map(l); jj < localP.graph.row_map(l+1); jj++) {
743 LO j = localP.graph.entries(jj);
744 SC v = localP.values(jj);
745 for (
size_t k = 0; k < dim; k++) {
747 SC n = localNullspace(i,k);
749 for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
750 if (P11colind(m) == jNew)
752 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA) 753 TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
755 P11vals(m) += 0.5 * v * n;
762 Kokkos::parallel_for(
"MueLu:RefMaxwell::buildProlongator_enterValues", range_type(0,numLocalRows),
763 KOKKOS_LAMBDA(
const size_t i) {
764 for (
size_t ll = localD0.graph.row_map(i); ll < localD0.graph.row_map(i+1); ll++) {
765 LO l = localD0.graph.entries(ll);
766 for (
size_t jj = localP.graph.row_map(l); jj < localP.graph.row_map(l+1); jj++) {
767 LO j = localP.graph.entries(jj);
768 SC v = localP.values(jj);
769 for (
size_t k = 0; k < dim; k++) {
771 SC n = localNullspace(i,k);
773 for (m = P11rowptr(i); m < P11rowptr(i+1); m++)
774 if (P11colind(m) == jNew)
776 #if defined(HAVE_MUELU_DEBUG) && !defined(HAVE_MUELU_CUDA) 777 TEUCHOS_ASSERT_EQUALITY(P11colind(m),jNew);
779 P11vals(m) += 0.5 * v * n;
786 P11_ = Teuchos::rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0, Xpetra::StaticProfile));
787 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
788 P11Crs->setAllValues(P11rowptr, P11colind, P11vals);
789 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
792 TEUCHOS_TEST_FOR_EXCEPTION(
false,std::invalid_argument,algo <<
" is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
794 #endif // ifdef(HAVE_MUELU_KOKKOS_REFACTOR) 797 ArrayRCP<ArrayRCP<const SC> > nullspaceRCP(dim);
798 ArrayRCP<ArrayView<const SC> > nullspace(dim);
799 for(
size_t i=0; i<dim; i++) {
800 nullspaceRCP[i] = Nullspace_->getData(i);
801 nullspace[i] = nullspaceRCP[i]();
805 ArrayRCP<const size_t> Prowptr_RCP, D0rowptr_RCP;
806 ArrayRCP<const LO> Pcolind_RCP, D0colind_RCP;
807 ArrayRCP<const SC> Pvals_RCP, D0vals_RCP;
808 ArrayRCP<size_t> P11rowptr_RCP;
809 ArrayRCP<LO> P11colind_RCP;
810 ArrayRCP<SC> P11vals_RCP;
812 P_nodal_imported->getAllValues(Prowptr_RCP, Pcolind_RCP, Pvals_RCP);
813 rcp_dynamic_cast<CrsMatrixWrap>(D0_Matrix_)->getCrsMatrix()->getAllValues(D0rowptr_RCP, D0colind_RCP, D0vals_RCP);
818 ArrayView<const size_t> Prowptr, D0rowptr;
819 ArrayView<const LO> Pcolind, D0colind;
820 ArrayView<const SC> Pvals, D0vals;
821 Prowptr = Prowptr_RCP(); Pcolind = Pcolind_RCP(); Pvals = Pvals_RCP();
822 D0rowptr = D0rowptr_RCP(); D0colind = D0colind_RCP(); D0vals = D0vals_RCP();
830 std::string defaultAlgo =
"gustavson";
831 std::string algo = parameterList_.get(
"refmaxwell: prolongator compute algorithm",defaultAlgo);
833 if (algo ==
"mat-mat") {
834 Teuchos::RCP<Matrix> D0_P_nodal = MatrixFactory::Build(SM_Matrix_->getRowMap(),0);
835 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
false,*P_nodal,
false,*D0_P_nodal,
true,
true);
838 ArrayRCP<const size_t> D0Prowptr_RCP;
839 ArrayRCP<const LO> D0Pcolind_RCP;
840 ArrayRCP<const SC> D0Pvals_RCP;
841 rcp_dynamic_cast<CrsMatrixWrap>(D0_P_nodal)->getCrsMatrix()->getAllValues(D0Prowptr_RCP, D0Pcolind_RCP, D0Pvals_RCP);
846 ArrayView<const size_t> D0Prowptr;
847 ArrayView<const LO> D0Pcolind;
848 D0Prowptr = D0Prowptr_RCP(); D0Pcolind = D0Pcolind_RCP();
851 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
852 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
853 P11_ = Teuchos::rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0, Xpetra::StaticProfile));
854 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
855 P11Crs->allocateAllValues(dim*D0Pcolind.size(), P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
857 ArrayView<size_t> P11rowptr = P11rowptr_RCP();
858 ArrayView<LO> P11colind = P11colind_RCP();
859 ArrayView<SC> P11vals = P11vals_RCP();
862 for (
size_t i = 0; i < numLocalRows+1; i++) {
863 P11rowptr[i] = dim*D0Prowptr[i];
868 for (
size_t jj = 0; jj < (size_t) D0Pcolind.size(); jj++)
869 for (
size_t k = 0; k < dim; k++) {
870 P11colind[nnz] = dim*D0Pcolind[jj]+k;
871 P11vals[nnz] = SC_ZERO;
876 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
880 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
882 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
883 for (
size_t i = 0; i < numLocalRows; i++) {
884 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
887 if (Teuchos::ScalarTraits<Scalar>::magnitude(p) < tol)
889 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
892 for (
size_t k = 0; k < dim; k++) {
894 SC n = nullspace[k][i];
896 for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
897 if (P11colind[m] == jNew)
899 #ifdef HAVE_MUELU_DEBUG 900 TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
902 P11vals[m] += 0.5 * v * n;
909 for (
size_t i = 0; i < numLocalRows; i++) {
910 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
912 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
915 for (
size_t k = 0; k < dim; k++) {
917 SC n = nullspace[k][i];
919 for (m = P11rowptr[i]; m < P11rowptr[i+1]; m++)
920 if (P11colind[m] == jNew)
922 #ifdef HAVE_MUELU_DEBUG 923 TEUCHOS_ASSERT_EQUALITY(P11colind[m],jNew);
925 P11vals[m] += 0.5 * v * n;
932 P11Crs->setAllValues(P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
933 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
935 }
else if (algo ==
"gustavson") {
937 LO maxP11col = dim * P_nodal_imported->getColMap()->getMaxLocalIndex();
938 const size_t ST_INVALID = Teuchos::OrdinalTraits<LO>::invalid();
939 Array<size_t> P11_status(dim*maxP11col, ST_INVALID);
941 size_t nnz_alloc = dim*D0vals_RCP.size();
944 RCP<Map> blockColMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal_imported->getColMap(), dim);
945 RCP<Map> blockDomainMap = Xpetra::MapFactory<LO,GO,NO>::Build(P_nodal->getDomainMap(), dim);
946 P11_ = Teuchos::rcp(
new CrsMatrixWrap(SM_Matrix_->getRowMap(), blockColMap, 0, Xpetra::StaticProfile));
947 RCP<CrsMatrix> P11Crs = rcp_dynamic_cast<CrsMatrixWrap>(P11_)->getCrsMatrix();
948 P11Crs->allocateAllValues(nnz_alloc, P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
950 ArrayView<size_t> P11rowptr = P11rowptr_RCP();
951 ArrayView<LO> P11colind = P11colind_RCP();
952 ArrayView<SC> P11vals = P11vals_RCP();
955 if (D0_Matrix_->getNodeMaxNumRowEntries()>2) {
959 GetOStream(
Warnings0) <<
"RefMaxwell::buildProlongator(): D0 matrix has more than 2 entries per row. Taking inefficient code path." << std::endl;
961 magnitudeType tol = Teuchos::ScalarTraits<magnitudeType>::eps();
964 for (
size_t i = 0; i < numLocalRows; i++) {
966 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
969 if (Teuchos::ScalarTraits<Scalar>::magnitude(p) < tol)
971 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
974 for (
size_t k = 0; k < dim; k++) {
976 SC n = nullspace[k][i];
978 if (P11_status[jNew] == ST_INVALID || P11_status[jNew] < nnz_old) {
979 P11_status[jNew] = nnz;
980 P11colind[nnz] = jNew;
981 P11vals[nnz] = 0.5 * v * n;
986 P11vals[P11_status[jNew]] += 0.5 * v * n;
995 P11rowptr[numLocalRows] = nnz;
999 for (
size_t i = 0; i < numLocalRows; i++) {
1001 for (
size_t ll = D0rowptr[i]; ll < D0rowptr[i+1]; ll++) {
1002 LO l = D0colind[ll];
1003 for (
size_t jj = Prowptr[l]; jj < Prowptr[l+1]; jj++) {
1006 for (
size_t k = 0; k < dim; k++) {
1008 SC n = nullspace[k][i];
1010 if (P11_status[jNew] == ST_INVALID || P11_status[jNew] < nnz_old) {
1011 P11_status[jNew] = nnz;
1012 P11colind[nnz] = jNew;
1013 P11vals[nnz] = 0.5 * v * n;
1018 P11vals[P11_status[jNew]] += 0.5 * v * n;
1027 P11rowptr[numLocalRows] = nnz;
1030 if (blockDomainMap->lib() == Xpetra::UseTpetra) {
1034 P11vals_RCP.resize(nnz);
1035 P11colind_RCP.resize(nnz);
1038 P11Crs->setAllValues(P11rowptr_RCP, P11colind_RCP, P11vals_RCP);
1039 P11Crs->expertStaticFillComplete(blockDomainMap, SM_Matrix_->getRangeMap());
1041 TEUCHOS_TEST_FOR_EXCEPTION(
false,std::invalid_argument,algo <<
" is not a valid option for \"refmaxwell: prolongator compute algorithm\"");
1042 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1047 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1049 Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer(
"MueLu RefMaxwell: Build coarse (1,1) matrix"));
1052 Teuchos::RCP<Matrix> Matrix1 = MatrixFactory::Build(P11_->getDomainMap(),0);
1053 if (parameterList_.get<
bool>(
"rap: triple product",
false) ==
false) {
1054 Teuchos::RCP<Matrix> C = MatrixFactory::Build(SM_Matrix_->getRowMap(),0);
1056 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*SM_Matrix_,
false,*P11_,
false,*C,
true,
true);
1058 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*R11_,
false,*C,
false,*Matrix1,
true,
true);
1060 Xpetra::TripleMatrixMultiply<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
1061 MultiplyRAP(*P11_,
true, *SM_Matrix_,
false, *P11_,
false, *Matrix1,
true,
true);
1063 if (parameterList_.get<
bool>(
"rap: fix zero diagonals",
true))
1064 Xpetra::MatrixUtils<SC,LO,GO,NO>::CheckRepairMainDiagonal(Matrix1,
true, GetOStream(
Warnings1));
1066 if(disable_addon_==
true) {
1071 Teuchos::TimeMonitor tmAddon(*Teuchos::TimeMonitor::getNewTimer(
"MueLu RefMaxwell: Build coarse addon matrix"));
1073 TEUCHOS_TEST_FOR_EXCEPTION(M0inv_Matrix_==Teuchos::null,std::invalid_argument,
1074 "MueLu::RefMaxwell::formCoarseMatrix(): Inverse of " 1075 "lumped mass matrix required for add-on (i.e. M0inv_Matrix is null)");
1078 Teuchos::RCP<Matrix> Zaux = MatrixFactory::Build(M1_Matrix_->getRowMap(),0);
1079 Teuchos::RCP<Matrix> Z = MatrixFactory::Build(D0_Matrix_->getDomainMap(),0);
1082 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M1_Matrix_,
false,*P11_,
false,*Zaux,
true,
true);
1084 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*D0_Matrix_,
true,*Zaux,
false,*Z,
true,
true);
1087 Teuchos::RCP<Matrix> Matrix2 = MatrixFactory::Build(Z->getDomainMap(),0);
1088 if (M0inv_Matrix_->getGlobalMaxNumRowEntries()<=1) {
1091 #ifdef HAVE_MUELU_KOKKOS_REFACTOR 1094 ZT = Utilities_kokkos::Transpose(*Z);
1100 RCP<Vector> diag = VectorFactory::Build(M0inv_Matrix_->getRowMap());
1101 M0inv_Matrix_->getLocalDiagCopy(*diag);
1102 if (Z->getRowMap()->isSameAs(*(diag->getMap())))
1103 Z->leftScale(*diag);
1105 RCP<Import> importer = ImportFactory::Build(diag->getMap(),Z->getRowMap());
1106 RCP<Vector> diag2 = VectorFactory::Build(Z->getRowMap());
1107 diag2->doImport(*diag,*importer,Xpetra::INSERT);
1108 Z->leftScale(*diag2);
1110 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*ZT,
false,*Z,
false,*Matrix2,
true,
true);
1111 }
else if (parameterList_.get<
bool>(
"rap: triple product",
false) ==
false) {
1112 Teuchos::RCP<Matrix> C2 = MatrixFactory::Build(M0inv_Matrix_->getRowMap(),0);
1114 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*M0inv_Matrix_,
false,*Z,
false,*C2,
true,
true);
1116 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::Multiply(*Z,
true,*C2,
false,*Matrix2,
true,
true);
1119 Xpetra::TripleMatrixMultiply<Scalar,LocalOrdinal,GlobalOrdinal,Node>::
1120 MultiplyRAP(*Z,
true, *M0inv_Matrix_,
false, *Z,
false, *Matrix2,
true,
true);
1123 Xpetra::MatrixMatrix<Scalar,LocalOrdinal,GlobalOrdinal,Node>::TwoMatrixAdd(*Matrix1,
false,(Scalar)1.0,*Matrix2,
false,(Scalar)1.0,AH_,GetOStream(
Runtime0));
1124 AH_->fillComplete();
1128 size_t dim = Nullspace_->getNumVectors();
1129 AH_->SetFixedBlockSize(dim);
1130 AH_->setObjectLabel(
"RefMaxwell (1,1)");
1135 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1137 SM_Matrix_ = SM_Matrix_new;
1141 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1145 Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
1146 SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1147 residual_->update(one, RHS, negone);
1148 R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1149 D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
1152 HierarchyH_->Iterate(*P11res_, *P11x_, 1,
true);
1153 Hierarchy22_->Iterate(*D0res_, *D0x_, 1,
true);
1156 P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1157 D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS,one,one);
1158 X.update(one, *residual_, one);
1163 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1167 Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
1168 SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1169 residual_->update(one, RHS, negone);
1170 R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1171 HierarchyH_->Iterate(*P11res_, *P11x_, 1,
true);
1172 P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1173 X.update(one, *residual_, one);
1176 SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1177 residual_->update(one, RHS, negone);
1178 D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
1179 Hierarchy22_->Iterate(*D0res_, *D0x_, 1,
true);
1180 D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS);
1181 X.update(one, *residual_, one);
1184 SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1185 residual_->update(one, RHS, negone);
1186 R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1187 HierarchyH_->Iterate(*P11res_, *P11x_, 1,
true);
1188 P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1189 X.update(one, *residual_, one);
1194 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1198 Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
1199 SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1200 residual_->update(one, RHS, negone);
1201 D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
1202 Hierarchy22_->Iterate(*D0res_, *D0x_, 1,
true);
1203 D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS);
1204 X.update(one, *residual_, one);
1207 SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1208 residual_->update(one, RHS, negone);
1209 R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1210 HierarchyH_->Iterate(*P11res_, *P11x_, 1,
true);
1211 P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1212 X.update(one, *residual_, one);
1215 SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1216 residual_->update(one, RHS, negone);
1217 D0_T_Matrix_->apply(*residual_,*D0res_,Teuchos::NO_TRANS);
1218 Hierarchy22_->Iterate(*D0res_, *D0x_, 1,
true);
1219 D0_Matrix_->apply(*D0x_,*residual_,Teuchos::NO_TRANS);
1220 X.update(one, *residual_, one);
1224 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1228 Scalar one = Teuchos::ScalarTraits<Scalar>::one(), negone = -one, zero = Teuchos::ScalarTraits<Scalar>::zero();
1229 SM_Matrix_->apply(X, *residual_, Teuchos::NO_TRANS, one, zero);
1230 residual_->update(one, RHS, negone);
1231 R11_->apply(*residual_,*P11res_,Teuchos::NO_TRANS);
1234 HierarchyH_->Iterate(*P11res_, *P11x_, 1,
true);
1237 P11_->apply(*P11x_,*residual_,Teuchos::NO_TRANS);
1238 X.update(one, *residual_, one);
1243 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1245 Teuchos::ETransp mode,
1247 Scalar beta)
const {
1249 Teuchos::TimeMonitor tm(*Teuchos::TimeMonitor::getNewTimer(
"MueLuRefmaxwell: Solve"));
1252 if (X.getNumVectors() != P11res_->getNumVectors()) {
1253 P11res_ = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
1254 P11x_ = MultiVectorFactory::Build(P11_->getDomainMap(),X.getNumVectors());
1255 D0res_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
1256 D0x_ = MultiVectorFactory::Build(D0_Matrix_->getDomainMap(),X.getNumVectors());
1257 residual_ = MultiVectorFactory::Build(SM_Matrix_->getDomainMap(),1);
1261 #if defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || defined(HAVE_MUELU_INST_DOUBLE_INT_INT)) 1262 if (useHiptmairSmoothing_) {
1265 hiptmairPreSmoother_->apply(tRHS, tX);
1269 PreSmoother_->Apply(X, RHS, use_as_preconditioner_);
1272 if(mode_==
"additive")
1273 applyInverseAdditive(RHS,X);
1274 else if(mode_==
"121")
1275 applyInverse121(RHS,X);
1276 else if(mode_==
"212")
1277 applyInverse212(RHS,X);
1278 else if(mode_==
"11only")
1279 applyInverse11only(RHS,X);
1280 else if(mode_==
"none") {
1284 applyInverseAdditive(RHS,X);
1287 #if defined(HAVE_MUELU_IFPACK2) && (!defined(HAVE_MUELU_EPETRA) || defined(HAVE_MUELU_INST_DOUBLE_INT_INT)) 1288 if (useHiptmairSmoothing_)
1292 hiptmairPostSmoother_->apply(tRHS, tX);
1296 PostSmoother_->Apply(X, RHS,
false);
1301 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1307 template<
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1310 const Teuchos::RCP<Matrix> & M0inv_Matrix,
1311 const Teuchos::RCP<Matrix> & M1_Matrix,
1312 const Teuchos::RCP<MultiVector> & Nullspace,
1313 const Teuchos::RCP<RealValuedMultiVector> & Coords,
1314 Teuchos::ParameterList& List)
1317 TEUCHOS_ASSERT(D0_Matrix!=Teuchos::null);
1318 TEUCHOS_ASSERT(M1_Matrix!=Teuchos::null);
1320 HierarchyH_ = Teuchos::null;
1321 Hierarchy22_ = Teuchos::null;
1322 PreSmoother_ = Teuchos::null;
1323 PostSmoother_ = Teuchos::null;
1324 disable_addon_ =
false;
1328 setParameters(List);
1330 D0_Matrix_ = D0_Matrix;
1331 M0inv_Matrix_ = M0inv_Matrix;
1332 M1_Matrix_ = M1_Matrix;
1334 Nullspace_ = Nullspace;
1337 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
1339 describe(Teuchos::FancyOStream& out,
const Teuchos::EVerbosityLevel verbLevel)
const {
1340 out <<
"\n--------------------------------------------------------------------------------\n" <<
1341 "--- RefMaxwell Summary ---\n" 1342 "--------------------------------------------------------------------------------" << std::endl;
1345 GlobalOrdinal numRows;
1348 numRows = SM_Matrix_->getGlobalNumRows();
1349 nnz = SM_Matrix_->getGlobalNumEntries();
1351 Xpetra::global_size_t tt = numRows;
1352 int rowspacer = 3;
while (tt != 0) { tt /= 10; rowspacer++; }
1354 int nnzspacer = 2;
while (tt != 0) { tt /= 10; nnzspacer++; }
1356 out <<
"block " << std::setw(rowspacer) <<
" rows " << std::setw(nnzspacer) <<
" nnz " << std::setw(9) <<
" nnz/row" << std::endl;
1357 out <<
"(1, 1)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
1359 numRows = A22_->getGlobalNumRows();
1360 nnz = A22_->getGlobalNumEntries();
1362 out <<
"(2, 2)" << std::setw(rowspacer) << numRows << std::setw(nnzspacer) << nnz << std::setw(9) << as<double>(nnz) / numRows << std::endl;
1370 #define MUELU_REFMAXWELL_SHORT 1371 #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...
Various adapters that will create a MueLu preconditioner that is an Xpetra::Matrix.
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
static Teuchos::ArrayRCP< const bool > DetectDirichletCols(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Teuchos::ArrayRCP< const bool > &dirichletRows)
Detect Dirichlet columns based on Dirichlet rows.
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.
Class that encapsulates external library smoothers.
static Teuchos::ArrayRCP< const bool > DetectDirichletRows(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A, const Magnitude &tol=Teuchos::ScalarTraits< Scalar >::magnitude(0.), const bool count_twos_as_dirichlet=false)
void ZeroDirichletCols(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const Kokkos::View< const bool *, typename Node::device_type > &dirichletCols, Scalar replaceWith)
Teuchos::RCP< MueLu::Hierarchy< Scalar, LocalOrdinal, GlobalOrdinal, Node > > CreateXpetraPreconditioner(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > op, const Teuchos::ParameterList &inParamList, Teuchos::RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > coords=Teuchos::null, Teuchos::RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > nullspace=Teuchos::null)
Helper function to create a MueLu preconditioner that can be used by Xpetra.Given an Xpetra::Matrix...
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.
Namespace for MueLu classes and methods.
void SetPreviousLevel(const RCP< Level > &previousLevel)
void SetFactoryManager(const RCP< const FactoryManagerBase > &factoryManager)
Set default factories (used internally by Hierarchy::SetLevel()).
void setParameters(Teuchos::ParameterList &list)
Set parameters.
static RCP< Xpetra::MultiVector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > RealValuedToScalarMultiVector(RCP< Xpetra::MultiVector< double, LocalOrdinal, GlobalOrdinal, Node > > X)
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)
Teuchos::ScalarTraits< Scalar >::magnitudeType magnitudeType
Factory for building tentative prolongator.
void setlib(Xpetra::UnderlyingLib lib2)
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::VERB_HIGH) const
void applyInverse121(const MultiVector &RHS, MultiVector &X) const
apply 1-2-1 algorithm for 2x2 solve
Class that holds all level-specific information.
void compute()
Setup the preconditioner.
Teuchos::RCP< const Map > getRangeMap() const
Returns the Xpetra::Map object associated with the range of this operator.
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)
AmalgamationFactory for subblocks of strided map based amalgamation data.
void initialize(const Teuchos::RCP< Matrix > &D0_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)
static void ZeroDirichletRows(Teuchos::RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &A, const std::vector< LocalOrdinal > &dirichletRows, Scalar replaceWith=Teuchos::ScalarTraits< Scalar >::zero())
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
bool hasTransposeApply() const
Indicates whether this operator supports applying the adjoint operator.
void applyInverse212(const MultiVector &RHS, MultiVector &X) const
apply 2-1-2 algorithm for 2x2 solve
Factory for creating a graph base on a given matrix.
void SetLevelID(int levelID)
Set level number.
Class for transferring coordinates from a finer level to a coarser one.
Factory for building restriction operators.
Kokkos::View< const bool *, typename NO::device_type > DetectDirichletCols(const Xpetra::Matrix< SC, LO, GO, NO > &A, const Kokkos::View< const bool *, typename NO::device_type > &dirichletRows)
void formCoarseMatrix()
Compute P11^{T}*A*P11 efficiently.
void applyInverse11only(const MultiVector &RHS, MultiVector &X) const
apply solve to 1-1 block only
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
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.
static RCP< Tpetra::RowMatrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op2NonConstTpetraRow(RCP< Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > > Op)
Factory for building uncoupled aggregates.
Factory for building a thresholded operator.
Kokkos::View< const bool *, typename NO::device_type > DetectDirichletRows(const Xpetra::Matrix< SC, LO, GO, NO > &A, const typename Teuchos::ScalarTraits< SC >::magnitudeType &tol, const bool count_twos_as_dirichlet)
void resetMatrix(Teuchos::RCP< Matrix > SM_Matrix_new)
Reset system matrix.
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.