46 #ifndef MUELU_PGPFACTORY_DEF_HPP
47 #define MUELU_PGPFACTORY_DEF_HPP
52 #include <Xpetra_MultiVectorFactory.hpp>
63 #include "MueLu_PerfUtils.hpp"
65 #include "MueLu_Utilities.hpp"
69 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
73 validParamList->
set<
RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
76 validParamList->
set<
bool>(
"ReUseRowBasedOmegas",
false,
"Reuse omegas for prolongator for restrictor");
78 return validParamList;
81 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
86 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
92 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 Input(fineLevel,
"A");
100 if (initialPFact == Teuchos::null) {
120 bool bReUseRowBasedOmegas = pL.
get<
bool>(
"ReUseRowBasedOmegas");
121 if (bReUseRowBasedOmegas ==
true && restrictionMode_ ==
true) {
126 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
128 FactoryMonitor m(*
this,
"Prolongator smoothing (PG-AMG)", coarseLevel);
137 if (initialPFact == Teuchos::null) {
143 if (restrictionMode_) {
149 bool doFillComplete =
true;
150 bool optimizeStorage =
true;
153 doFillComplete =
true;
154 optimizeStorage =
false;
163 bool bReUseRowBasedOmegas = pL.
get<
bool>(
"ReUseRowBasedOmegas");
164 if (restrictionMode_ ==
false || bReUseRowBasedOmegas ==
false) {
167 ComputeRowBasedOmega(fineLevel, coarseLevel, A, Ptent, DinvAP0, RowBasedOmega);
180 ExportFactory::Build(RowBasedOmega->getMap(), A->getRangeMap());
183 VectorFactory::Build(A->getRangeMap());
185 noRowBasedOmega->doExport(*RowBasedOmega, *exporter,
Xpetra::INSERT);
191 ImportFactory::Build(A->getRangeMap(), A->getRowMap());
194 RowBasedOmega->doImport(*noRowBasedOmega, *importer,
Xpetra::INSERT);
206 P_smoothed->fillComplete(Ptent->getDomainMap(), Ptent->getRangeMap());
211 params->
set(
"printLoadBalancingInfo",
true);
214 if (!restrictionMode_) {
216 Set(coarseLevel,
"P", P_smoothed);
224 Set(coarseLevel,
"RfromPfactory", dummy);
230 if (Ptent->IsView(
"stridedMaps"))
231 P_smoothed->CreateView(
"stridedMaps", Ptent);
236 Set(coarseLevel,
"R", R);
242 if (Ptent->IsView(
"stridedMaps"))
243 R->CreateView(
"stridedMaps", Ptent,
true);
247 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
249 FactoryMonitor m(*
this,
"PgPFactory::ComputeRowBasedOmega", coarseLevel);
274 bool doFillComplete =
true;
275 bool optimizeStorage =
false;
281 Numerator = VectorFactory::Build(ADinvAP0->getColMap(),
true);
282 Denominator = VectorFactory::Build(ADinvAP0->getColMap(),
true);
283 MultiplyAll(AP0, ADinvAP0, Numerator);
284 MultiplySelfAll(ADinvAP0, Denominator);
293 Numerator = VectorFactory::Build(DinvAP0->getColMap(),
true);
294 Denominator = VectorFactory::Build(DinvAP0->getColMap(),
true);
295 MultiplyAll(P0, DinvAP0, Numerator);
296 MultiplySelfAll(DinvAP0, Denominator);
309 bool doFillComplete =
true;
310 bool optimizeStorage =
true;
316 Numerator = VectorFactory::Build(DinvADinvAP0->getColMap(),
true);
317 Denominator = VectorFactory::Build(DinvADinvAP0->getColMap(),
true);
318 MultiplyAll(DinvAP0, DinvADinvAP0, Numerator);
319 MultiplySelfAll(DinvADinvAP0, Denominator);
327 VectorFactory::Build(Numerator->getMap() ,
true);
329 ColBasedOmega->putScalar(-666 );
336 Magnitude min_local = 1000000.0;
337 Magnitude max_local = 0.0;
338 for (
LocalOrdinal i = 0; i < Teuchos::as<LocalOrdinal>(Numerator->getLocalLength()); i++) {
340 ColBasedOmega_local[i] = 0.0;
343 ColBasedOmega_local[i] = Numerator_local[i] / Denominator_local[i];
356 ColBasedOmega_local[i] = 0.0;
372 MueLu_sumAll(A->getRowMap()->getComm(), zero_local, zero_all);
373 MueLu_sumAll(A->getRowMap()->getComm(), nan_local, nan_all);
374 MueLu_minAll(A->getRowMap()->getComm(), min_local, min_all);
375 MueLu_maxAll(A->getRowMap()->getComm(), max_local, max_all);
382 default: GetOStream(
Statistics1) <<
"unknown)" << std::endl;
break;
384 GetOStream(
Statistics1) <<
"Damping parameter: min = " << min_all <<
", max = " << max_all << std::endl;
385 GetOStream(
Statistics) <<
"# negative omegas: " << zero_all <<
" out of " << ColBasedOmega->getGlobalLength() <<
" column-based omegas" << std::endl;
386 GetOStream(
Statistics) <<
"# NaNs: " << nan_all <<
" out of " << ColBasedOmega->getGlobalLength() <<
" column-based omegas" << std::endl;
389 if (coarseLevel.
IsRequested(
"ColBasedOmega",
this)) {
390 coarseLevel.
Set(
"ColBasedOmega", ColBasedOmega,
this);
396 VectorFactory::Build(DinvAP0->getRowMap(),
true);
398 RowBasedOmega->putScalar(-666);
400 bool bAtLeastOneDefined =
false;
402 for (
LocalOrdinal row = 0; row < Teuchos::as<LocalOrdinal>(A->getLocalNumRows()); row++) {
405 DinvAP0->getLocalRowView(row, lindices, lvals);
406 bAtLeastOneDefined =
false;
407 for (
size_t j = 0; j < Teuchos::as<size_t>(lindices.
size()); j++) {
408 Scalar omega = ColBasedOmega_local[lindices[j]];
410 bAtLeastOneDefined =
true;
412 RowBasedOmega_local[row] = omega;
414 RowBasedOmega_local[row] = omega;
417 if (bAtLeastOneDefined ==
true) {
419 RowBasedOmega_local[row] = sZero;
423 if (coarseLevel.
IsRequested(
"RowBasedOmega",
this)) {
424 Set(coarseLevel,
"RowBasedOmega", RowBasedOmega);
428 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
438 for (
size_t n = 0; n < Op->getLocalNumRows(); n++) {
439 Op->getLocalRowView(n, lindices, lvals);
440 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices.
size()); i++) {
441 InnerProd_local[lindices[i]] += lvals[i] * lvals[i];
448 ExportFactory::Build(Op->getColMap(), Op->getDomainMap());
451 VectorFactory::Build(Op->getDomainMap());
453 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
459 ImportFactory::Build(Op->getDomainMap(), Op->getColMap());
465 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
469 #if 1 // 1=new "fast code, 0=old "slow", but safe code
470 #if 0 // not necessary - remove me
471 if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
474 std::vector<LocalOrdinal> NewRightLocal(right->getColMap()->getLocalNumElements(), Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements()+1));
477 for (
size_t j=0; j < right->getColMap()->getLocalNumElements(); j++) {
478 while ( (i < Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements())) &&
479 (left->getColMap()->getGlobalElement(i) < right->getColMap()->getGlobalElement(j)) ) i++;
480 if (left->getColMap()->getGlobalElement(i) == right->getColMap()->getGlobalElement(j)) {
481 NewRightLocal[j] = i;
486 std::vector<Scalar> temp_array(left->getColMap()->getLocalNumElements()+1, 0.0);
488 for(
size_t n=0; n<right->getLocalNumRows(); n++) {
494 left->getLocalRowView (n, lindices_left, lvals_left);
495 right->getLocalRowView(n, lindices_right, lvals_right);
497 for(
size_t j=0; j<Teuchos::as<size_t>(lindices_right.
size()); j++) {
498 temp_array[NewRightLocal[lindices_right[j] ] ] = lvals_right[j];
500 for (
size_t j=0; j < Teuchos::as<size_t>(lindices_left.
size()); j++) {
501 InnerProd_local[lindices_left[j]] += temp_array[lindices_left[j] ]*lvals_left[j];
503 for (
size_t j=0; j < Teuchos::as<size_t>(lindices_right.
size()); j++) {
504 temp_array[NewRightLocal[lindices_right[j] ] ] = 0.0;
509 ExportFactory::Build(left->getColMap(), left->getDomainMap());
512 VectorFactory::Build(left->getDomainMap());
514 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
520 ImportFactory::Build(left->getDomainMap(), left->getColMap());
527 #endif // end remove me
528 if (InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
529 size_t szNewLeftLocal = TEUCHOS_MAX(left->getColMap()->getLocalNumElements(), right->getColMap()->getLocalNumElements());
533 for (
size_t i = 0; i < left->getColMap()->getLocalNumElements(); i++) {
534 while ((j < Teuchos::as<LocalOrdinal>(right->getColMap()->getLocalNumElements())) &&
535 (right->getColMap()->getGlobalElement(j) < left->getColMap()->getGlobalElement(i))) j++;
536 if (right->getColMap()->getGlobalElement(j) == left->getColMap()->getGlobalElement(i)) {
537 (*NewLeftLocal)[i] = j;
548 for (
size_t n = 0; n < left->getLocalNumRows(); n++) {
554 left->getLocalRowView(n, lindices_left, lvals_left);
555 right->getLocalRowView(n, lindices_right, lvals_right);
557 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_left.
size()); i++) {
558 (*temp_array)[(*NewLeftLocal)[lindices_left[i]]] = lvals_left[i];
560 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_right.
size()); i++) {
561 InnerProd_local[lindices_right[i]] += (*temp_array)[lindices_right[i]] * lvals_right[i];
563 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_left.
size()); i++) {
564 (*temp_array)[(*NewLeftLocal)[lindices_left[i]]] = 0.0;
570 ExportFactory::Build(right->getColMap(), right->getDomainMap());
573 VectorFactory::Build(right->getDomainMap());
575 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
581 ImportFactory::Build(right->getDomainMap(), right->getColMap());
588 #else // old "safe" code
589 if (InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
598 for (
size_t n = 0; n < left->getLocalNumRows(); n++) {
599 left->getLocalRowView(n, lindices_left, lvals_left);
600 right->getLocalRowView(n, lindices_right, lvals_right);
602 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_left.
size()); i++) {
603 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
604 for (
size_t j = 0; j < Teuchos::as<size_t>(lindices_right.
size()); j++) {
605 GlobalOrdinal right_gid = right->getColMap()->getGlobalElement(lindices_right[j]);
606 if (left_gid == right_gid) {
607 InnerProd_local[lindices_left[i]] += lvals_left[i] * lvals_right[j];
616 ExportFactory::Build(left->getColMap(), left->getDomainMap());
619 VectorFactory::Build(left->getDomainMap());
621 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
627 ImportFactory::Build(left->getDomainMap(), left->getColMap());
631 }
else if (InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
639 for (
size_t n = 0; n < left->getLocalNumRows(); n++) {
640 left->getLocalRowView(n, lindices_left, lvals_left);
641 right->getLocalRowView(n, lindices_right, lvals_right);
643 for (
size_t i = 0; i < Teuchos::as<size_t>(lindices_left.
size()); i++) {
644 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
645 for (
size_t j = 0; j < Teuchos::as<size_t>(lindices_right.
size()); j++) {
646 GlobalOrdinal right_gid = right->getColMap()->getGlobalElement(lindices_right[j]);
647 if (left_gid == right_gid) {
648 InnerProd_local[lindices_right[j]] += lvals_left[i] * lvals_right[j];
657 ExportFactory::Build(right->getColMap(), right->getDomainMap());
660 VectorFactory::Build(right->getDomainMap());
662 nonoverlap->doExport(*InnerProdVec, *exporter,
Xpetra::ADD);
668 ImportFactory::Build(right->getDomainMap(), right->getColMap());
678 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
680 std::cout <<
"TODO: remove me" << std::endl;
683 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.
static void MyOldScaleMatrix(Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &Op, const Teuchos::ArrayRCP< const Scalar > &scalingVector, bool doInverse=true, bool doFillComplete=true, bool doOptimizeStorage=true)
#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.
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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)
#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 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