10 #ifndef MUELU_PGPFACTORY_DEF_HPP
11 #define MUELU_PGPFACTORY_DEF_HPP
16 #include <Xpetra_MultiVectorFactory.hpp>
27 #include "MueLu_PerfUtils.hpp"
29 #include "MueLu_Utilities.hpp"
33 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
37 validParamList->
set<
RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
40 validParamList->
set<
bool>(
"ReUseRowBasedOmegas",
false,
"Reuse omegas for prolongator for restrictor");
42 return validParamList;
45 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
50 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
56 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
58 Input(fineLevel,
"A");
64 if (initialPFact == Teuchos::null) {
84 bool bReUseRowBasedOmegas = pL.
get<
bool>(
"ReUseRowBasedOmegas");
85 if (bReUseRowBasedOmegas ==
true && restrictionMode_ ==
true) {
90 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 FactoryMonitor m(*
this,
"Prolongator smoothing (PG-AMG)", coarseLevel);
101 if (initialPFact == Teuchos::null) {
107 if (restrictionMode_) {
113 bool doFillComplete =
true;
114 bool optimizeStorage =
true;
117 const auto rowMap = A->getRowMap();
119 A->getLocalDiagCopy(*diag);
120 diag->reciprocal(*diag);
121 DinvAP0->leftScale(*diag);
128 bool bReUseRowBasedOmegas = pL.
get<
bool>(
"ReUseRowBasedOmegas");
129 if (restrictionMode_ ==
false || bReUseRowBasedOmegas ==
false) {
132 ComputeRowBasedOmega(fineLevel, coarseLevel, A, Ptent, DinvAP0, RowBasedOmega);
145 ExportFactory::Build(RowBasedOmega->getMap(), A->getRangeMap());
148 VectorFactory::Build(A->getRangeMap());
150 noRowBasedOmega->doExport(*RowBasedOmega, *exporter,
Xpetra::INSERT);
156 ImportFactory::Build(A->getRangeMap(), A->getRowMap());
159 RowBasedOmega->doImport(*noRowBasedOmega, *importer,
Xpetra::INSERT);
162 DinvAP0->leftScale(*RowBasedOmega);
169 P_smoothed->fillComplete(Ptent->getDomainMap(), Ptent->getRangeMap());
174 params->
set(
"printLoadBalancingInfo",
true);
177 if (!restrictionMode_) {
179 Set(coarseLevel,
"P", P_smoothed);
187 Set(coarseLevel,
"RfromPfactory", dummy);
193 if (Ptent->IsView(
"stridedMaps"))
194 P_smoothed->CreateView(
"stridedMaps", Ptent);
199 Set(coarseLevel,
"R", R);
205 if (Ptent->IsView(
"stridedMaps"))
206 R->CreateView(
"stridedMaps", Ptent,
true);
210 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
212 FactoryMonitor m(*
this,
"PgPFactory::ComputeRowBasedOmega", coarseLevel);
237 bool doFillComplete =
true;
238 bool optimizeStorage =
false;
244 Numerator = VectorFactory::Build(ADinvAP0->getColMap(),
true);
245 Denominator = VectorFactory::Build(ADinvAP0->getColMap(),
true);
246 MultiplyAll(AP0, ADinvAP0, Numerator);
247 MultiplySelfAll(ADinvAP0, Denominator);
256 Numerator = VectorFactory::Build(DinvAP0->getColMap(),
true);
257 Denominator = VectorFactory::Build(DinvAP0->getColMap(),
true);
258 MultiplyAll(P0, DinvAP0, Numerator);
259 MultiplySelfAll(DinvAP0, Denominator);
272 bool doFillComplete =
true;
273 bool optimizeStorage =
true;
277 const auto rowMap = A->getRowMap();
279 A->getLocalDiagCopy(*diag);
280 diag->reciprocal(*diag);
281 DinvADinvAP0->leftScale(*diag);
283 Numerator = VectorFactory::Build(DinvADinvAP0->getColMap(),
true);
284 Denominator = VectorFactory::Build(DinvADinvAP0->getColMap(),
true);
285 MultiplyAll(DinvAP0, DinvADinvAP0, Numerator);
286 MultiplySelfAll(DinvADinvAP0, Denominator);
294 VectorFactory::Build(Numerator->getMap() ,
true);
296 ColBasedOmega->putScalar(-666 );
303 Magnitude min_local = 1000000.0;
304 Magnitude max_local = 0.0;
305 for (
LocalOrdinal i = 0; i < Teuchos::as<LocalOrdinal>(Numerator->getLocalLength()); i++) {
307 ColBasedOmega_local[i] = 0.0;
310 ColBasedOmega_local[i] = Numerator_local[i] / Denominator_local[i];
323 ColBasedOmega_local[i] = 0.0;
339 MueLu_sumAll(A->getRowMap()->getComm(), zero_local, zero_all);
340 MueLu_sumAll(A->getRowMap()->getComm(), nan_local, nan_all);
341 MueLu_minAll(A->getRowMap()->getComm(), min_local, min_all);
342 MueLu_maxAll(A->getRowMap()->getComm(), max_local, max_all);
349 default: GetOStream(
Statistics1) <<
"unknown)" << std::endl;
break;
351 GetOStream(
Statistics1) <<
"Damping parameter: min = " << min_all <<
", max = " << max_all << std::endl;
352 GetOStream(
Statistics) <<
"# negative omegas: " << zero_all <<
" out of " << ColBasedOmega->getGlobalLength() <<
" column-based omegas" << std::endl;
353 GetOStream(
Statistics) <<
"# NaNs: " << nan_all <<
" out of " << ColBasedOmega->getGlobalLength() <<
" column-based omegas" << std::endl;
356 if (coarseLevel.
IsRequested(
"ColBasedOmega",
this)) {
357 coarseLevel.
Set(
"ColBasedOmega", ColBasedOmega,
this);
363 VectorFactory::Build(DinvAP0->getRowMap(),
true);
365 RowBasedOmega->putScalar(-666);
367 bool bAtLeastOneDefined =
false;
369 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(A->getLocalNumRows()); row++) {
372 DinvAP0->getLocalRowView(row, lindices, lvals);
373 bAtLeastOneDefined =
false;
374 for (
size_t j = 0; j < Teuchos::as<size_t>(lindices.
size()); j++) {
375 Scalar omega = ColBasedOmega_local[lindices[j]];
377 bAtLeastOneDefined =
true;
379 RowBasedOmega_local[row] = omega;
381 RowBasedOmega_local[row] = omega;
384 if (bAtLeastOneDefined ==
true) {
386 RowBasedOmega_local[row] = sZero;
390 if (coarseLevel.
IsRequested(
"RowBasedOmega",
this)) {
391 Set(coarseLevel,
"RowBasedOmega", RowBasedOmega);
395 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
405 for (
size_t n = 0; n < Op->getLocalNumRows(); n++) {
406 Op->getLocalRowView(n, lindices, lvals);
407 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices.
size()); i++) {
408 InnerProd_local[lindices[i]] += lvals[i] * lvals[i];
415 ExportFactory::Build(Op->getColMap(), Op->getDomainMap());
418 VectorFactory::Build(Op->getDomainMap());
420 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
426 ImportFactory::Build(Op->getDomainMap(), Op->getColMap());
432 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
436 #if 1 // 1=new "fast code, 0=old "slow", but safe code
437 #if 0 // not necessary - remove me
438 if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
441 std::vector<LocalOrdinal> NewRightLocal(right->getColMap()->getLocalNumElements(), Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements()+1));
444 for (
size_t j=0; j < right->getColMap()->getLocalNumElements(); j++) {
445 while ( (i < Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements())) &&
446 (left->getColMap()->getGlobalElement(i) < right->getColMap()->getGlobalElement(j)) ) i++;
447 if (left->getColMap()->getGlobalElement(i) == right->getColMap()->getGlobalElement(j)) {
448 NewRightLocal[j] = i;
453 std::vector<Scalar> temp_array(left->getColMap()->getLocalNumElements()+1, 0.0);
455 for(
size_t n=0; n<right->getLocalNumRows(); n++) {
461 left->getLocalRowView (n, lindices_left, lvals_left);
462 right->getLocalRowView(n, lindices_right, lvals_right);
464 for(
size_t j=0; j<Teuchos::as<size_t>(lindices_right.
size()); j++) {
465 temp_array[NewRightLocal[lindices_right[j] ] ] = lvals_right[j];
467 for (
size_t j=0; j < Teuchos::as<size_t>(lindices_left.
size()); j++) {
468 InnerProd_local[lindices_left[j]] += temp_array[lindices_left[j] ]*lvals_left[j];
470 for (
size_t j=0; j < Teuchos::as<size_t>(lindices_right.
size()); j++) {
471 temp_array[NewRightLocal[lindices_right[j] ] ] = 0.0;
476 ExportFactory::Build(left->getColMap(), left->getDomainMap());
479 VectorFactory::Build(left->getDomainMap());
481 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
487 ImportFactory::Build(left->getDomainMap(), left->getColMap());
494 #endif // end remove me
495 if (InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
496 size_t szNewLeftLocal = TEUCHOS_MAX(left->getColMap()->getLocalNumElements(), right->getColMap()->getLocalNumElements());
500 for (
size_t i = 0; i < left->getColMap()->getLocalNumElements(); i++) {
501 while ((j < Teuchos::as<LocalOrdinal>(right->getColMap()->getLocalNumElements())) &&
502 (right->getColMap()->getGlobalElement(j) < left->getColMap()->getGlobalElement(i))) j++;
503 if (right->getColMap()->getGlobalElement(j) == left->getColMap()->getGlobalElement(i)) {
504 (*NewLeftLocal)[i] = j;
515 for (
size_t n = 0; n < left->getLocalNumRows(); n++) {
521 left->getLocalRowView(n, lindices_left, lvals_left);
522 right->getLocalRowView(n, lindices_right, lvals_right);
524 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_left.
size()); i++) {
525 (*temp_array)[(*NewLeftLocal)[lindices_left[i]]] = lvals_left[i];
527 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_right.
size()); i++) {
528 InnerProd_local[lindices_right[i]] += (*temp_array)[lindices_right[i]] * lvals_right[i];
530 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_left.
size()); i++) {
531 (*temp_array)[(*NewLeftLocal)[lindices_left[i]]] = 0.0;
537 ExportFactory::Build(right->getColMap(), right->getDomainMap());
540 VectorFactory::Build(right->getDomainMap());
542 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
548 ImportFactory::Build(right->getDomainMap(), right->getColMap());
555 #else // old "safe" code
556 if (InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
565 for (
size_t n = 0; n < left->getLocalNumRows(); n++) {
566 left->getLocalRowView(n, lindices_left, lvals_left);
567 right->getLocalRowView(n, lindices_right, lvals_right);
569 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_left.
size()); i++) {
570 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
571 for (
size_t j = 0; j < Teuchos::as<size_t>(lindices_right.
size()); j++) {
572 GlobalOrdinal right_gid = right->getColMap()->getGlobalElement(lindices_right[j]);
573 if (left_gid == right_gid) {
574 InnerProd_local[lindices_left[i]] += lvals_left[i] * lvals_right[j];
583 ExportFactory::Build(left->getColMap(), left->getDomainMap());
586 VectorFactory::Build(left->getDomainMap());
588 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
594 ImportFactory::Build(left->getDomainMap(), left->getColMap());
598 }
else if (InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
606 for (
size_t n = 0; n < left->getLocalNumRows(); n++) {
607 left->getLocalRowView(n, lindices_left, lvals_left);
608 right->getLocalRowView(n, lindices_right, lvals_right);
610 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_left.
size()); i++) {
611 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
612 for (
size_t j = 0; j < Teuchos::as<size_t>(lindices_right.
size()); j++) {
613 GlobalOrdinal right_gid = right->getColMap()->getGlobalElement(lindices_right[j]);
614 if (left_gid == right_gid) {
615 InnerProd_local[lindices_right[j]] += lvals_left[i] * lvals_right[j];
624 ExportFactory::Build(right->getColMap(), right->getDomainMap());
627 VectorFactory::Build(right->getDomainMap());
629 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
635 ImportFactory::Build(right->getDomainMap(), right->getColMap());
645 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
647 std::cout <<
"TODO: remove me" << std::endl;
650 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
#define MueLu_sumAll(rcpComm, in, out)
MueLu::DefaultLocalOrdinal LocalOrdinal
T & Get(const std::string &ename, const FactoryBase *factory=NoFactory::get())
Get data without decrementing associated storage counter (i.e., read-only access). Usage: Level->Get< RCP<Matrix> >("A", factory) if factory == NULL => use default factory.
#define MueLu_maxAll(rcpComm, in, out)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
T & get(const std::string &name, T def_value)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
void ComputeRowBasedOmega(Level &fineLevel, Level &coarseLevel, const RCP< Matrix > &A, const RCP< Matrix > &P0, const RCP< Matrix > &DinvAP0, RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &RowBasedOmega) const
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define MueLu_minAll(rcpComm, in, out)
Print even more statistics.
void SetMinimizationMode(MinimizationNorm minnorm)
Set minimization mode (L2NORM for cheapest, ANORM more expensive, DINVANORM = default) ...
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
MueLu::DefaultScalar Scalar
MueLu::DefaultGlobalOrdinal GlobalOrdinal
void Build(Level &fineLevel, Level &coarseLevel) const
Build method.
Class that holds all level-specific information.
Timer to be used in factories. Similar to SubMonitor but adds a timer level by level.
void MultiplySelfAll(const RCP< Matrix > &Op, Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &InnerProdVec) const
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal_arcp(const Matrix &A)
Extract Matrix Diagonal.
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)
static std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
MinimizationNorm GetMinimizationMode()
return minimization mode
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static RCP< Vector > Build(const Teuchos::RCP< const Map > &map, bool zeroOut=true)
static void TwoMatrixAdd(const Matrix &A, bool transposeA, SC alpha, Matrix &B, SC beta)
static magnitudeType magnitude(T a)
void ReUseDampingParameters(bool bReuse)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
static void Multiply(const Matrix &A, bool transposeA, const Matrix &B, bool transposeB, Matrix &C, bool call_FillComplete_on_result=true, bool doOptimizeStorage=true, const std::string &label=std::string(), const RCP< ParameterList > ¶ms=null)
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Exception throws to report errors in the internal logical of the program.
bool IsRequested(const std::string &ename, const FactoryBase *factory=NoFactory::get()) const
Test whether a need has been requested. Note: this tells nothing about whether the need's value exist...
void DeclareInput(const std::string &ename, const FactoryBase *factory, const FactoryBase *requestedBy=NoFactory::get())
Callback from FactoryBase::CallDeclareInput() and FactoryBase::DeclareInput()
void MultiplyAll(const RCP< Matrix > &left, const RCP< Matrix > &right, Teuchos::RCP< Xpetra::Vector< Scalar, LocalOrdinal, GlobalOrdinal, Node > > &InnerProdVec) const