46 #ifndef MUELU_PGPFACTORY_DEF_HPP 47 #define MUELU_PGPFACTORY_DEF_HPP 51 #include <Xpetra_Vector.hpp> 52 #include <Xpetra_VectorFactory.hpp> 53 #include <Xpetra_Import.hpp> 54 #include <Xpetra_ImportFactory.hpp> 55 #include <Xpetra_Export.hpp> 56 #include <Xpetra_ExportFactory.hpp> 57 #include <Xpetra_Matrix.hpp> 58 #include <Xpetra_MatrixMatrix.hpp> 62 #include "MueLu_PerfUtils.hpp" 65 #include "MueLu_SmootherFactory.hpp" 66 #include "MueLu_TentativePFactory.hpp" 67 #include "MueLu_Utilities.hpp" 71 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
75 validParamList->
set<
RCP<const FactoryBase> >(
"A", Teuchos::null,
"Generating factory of the matrix A used during the prolongator smoothing process");
78 validParamList->
set<
bool > (
"ReUseRowBasedOmegas",
false,
"Reuse omegas for prolongator for restrictor");
80 return validParamList;
83 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
88 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
94 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
96 Input(fineLevel,
"A");
102 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.
GetFactoryManager()->GetFactory(
"Ptent"); }
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);
131 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
137 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.
GetFactoryManager()->GetFactory(
"Ptent"); }
141 if(restrictionMode_) {
147 bool doFillComplete=
true;
148 bool optimizeStorage=
true;
149 RCP<Matrix> DinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *Ptent,
false, GetOStream(
Statistics2), doFillComplete, optimizeStorage);
152 optimizeStorage=
false;
161 bool bReUseRowBasedOmegas = pL.
get<
bool>(
"ReUseRowBasedOmegas");
162 if(restrictionMode_ ==
false || bReUseRowBasedOmegas ==
false) {
165 ComputeRowBasedOmega(fineLevel, coarseLevel, A, Ptent, DinvAP0, RowBasedOmega);
178 ExportFactory::Build(RowBasedOmega->getMap(), A->getRangeMap());
181 VectorFactory::Build(A->getRangeMap());
183 noRowBasedOmega->doExport(*RowBasedOmega, *exporter, Xpetra::INSERT);
189 ImportFactory::Build(A->getRangeMap(), A->getRowMap());
192 RowBasedOmega->doImport(*noRowBasedOmega, *importer, Xpetra::INSERT);
204 P_smoothed->fillComplete(Ptent->getDomainMap(), Ptent->getRangeMap());
209 params->
set(
"printLoadBalancingInfo",
true);
212 if (!restrictionMode_) {
214 Set(coarseLevel,
"P", P_smoothed);
220 if (Ptent->IsView(
"stridedMaps"))
221 P_smoothed->CreateView(
"stridedMaps", Ptent);
226 Set(coarseLevel,
"R", R);
232 if (Ptent->IsView(
"stridedMaps"))
233 R->CreateView(
"stridedMaps", Ptent,
true);
238 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
240 FactoryMonitor m(*
this,
"PgPFactory::ComputeRowBasedOmega", coarseLevel);
266 bool doFillComplete=
true;
267 bool optimizeStorage=
false;
268 RCP<Matrix> AP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *P0,
false, GetOStream(
Statistics2), doFillComplete, optimizeStorage);
271 RCP<Matrix> ADinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *DinvAP0,
false, GetOStream(
Statistics2), doFillComplete, optimizeStorage);
273 Numerator = VectorFactory::Build(ADinvAP0->getColMap(),
true);
274 Denominator = VectorFactory::Build(ADinvAP0->getColMap(),
true);
275 MultiplyAll(AP0, ADinvAP0, Numerator);
276 MultiplySelfAll(ADinvAP0, Denominator);
287 Numerator = VectorFactory::Build(DinvAP0->getColMap(),
true);
288 Denominator = VectorFactory::Build(DinvAP0->getColMap(),
true);
289 MultiplyAll(P0, DinvAP0, Numerator);
290 MultiplySelfAll(DinvAP0, Denominator);
304 bool doFillComplete=
true;
305 bool optimizeStorage=
false;
307 RCP<Matrix> DinvADinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *DinvAP0,
false, GetOStream(
Statistics2), doFillComplete, optimizeStorage);
310 Numerator = VectorFactory::Build(DinvADinvAP0->getColMap(),
true);
311 Denominator = VectorFactory::Build(DinvADinvAP0->getColMap(),
true);
312 MultiplyAll(DinvAP0, DinvADinvAP0, Numerator);
313 MultiplySelfAll(DinvADinvAP0, Denominator);
324 VectorFactory::Build(Numerator->getMap(),
true);
326 ColBasedOmega->putScalar(-666);
331 GlobalOrdinal zero_local = 0;
332 GlobalOrdinal nan_local = 0;
333 Magnitude min_local = 1000000.0;
334 Magnitude max_local = 0.0;
335 for(LocalOrdinal i = 0; i < Teuchos::as<LocalOrdinal>(Numerator->getLocalLength()); i++) {
338 ColBasedOmega_local[i] = 0.0;
343 ColBasedOmega_local[i] = Numerator_local[i] / Denominator_local[i];
356 ColBasedOmega_local[i] = 0.0;
366 GlobalOrdinal zero_all;
367 GlobalOrdinal nan_all;
370 MueLu_sumAll(A->getRowMap()->getComm(), zero_local, zero_all);
371 MueLu_sumAll(A->getRowMap()->getComm(), nan_local, nan_all);
372 MueLu_minAll(A->getRowMap()->getComm(), min_local, min_all);
373 MueLu_maxAll(A->getRowMap()->getComm(), max_local, max_all);
380 default: GetOStream(
Statistics1) <<
"unknown)" << std::endl;
break;
382 GetOStream(
Statistics1) <<
"Damping parameter: min = " << min_all <<
", max = " << max_all << std::endl;
383 GetOStream(
Statistics) <<
"# negative omegas: " << zero_all <<
" out of " << ColBasedOmega->getGlobalLength() <<
" column-based omegas" << std::endl;
384 GetOStream(
Statistics) <<
"# NaNs: " << nan_all <<
" out of " << ColBasedOmega->getGlobalLength() <<
" column-based omegas" << std::endl;
387 if(coarseLevel.
IsRequested(
"ColBasedOmega",
this)) {
388 coarseLevel.
Set(
"ColBasedOmega", ColBasedOmega,
this);
394 VectorFactory::Build(DinvAP0->getRowMap(),
true);
396 RowBasedOmega->putScalar(-666);
398 bool bAtLeastOneDefined =
false;
400 for(LocalOrdinal row = 0; row<Teuchos::as<LocalOrdinal>(A->getNodeNumRows()); row++) {
403 DinvAP0->getLocalRowView(row, lindices, lvals);
404 bAtLeastOneDefined =
false;
405 for(
size_t j=0; j<Teuchos::as<size_t>(lindices.
size()); j++) {
406 Scalar omega = ColBasedOmega_local[lindices[j]];
408 bAtLeastOneDefined =
true;
413 if(bAtLeastOneDefined ==
true) {
415 RowBasedOmega_local[row] = sZero;
419 if(coarseLevel.
IsRequested(
"RowBasedOmega",
this)) {
420 Set(coarseLevel,
"RowBasedOmega", RowBasedOmega);
424 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
435 for(
size_t n=0; n<Op->getNodeNumRows(); n++) {
436 Op->getLocalRowView(n, lindices, lvals);
437 for(
size_t i=0; i<Teuchos::as<size_t>(lindices.
size()); i++) {
438 InnerProd_local[lindices[i]] += lvals[i]*lvals[i];
444 ExportFactory::Build(Op->getColMap(), Op->getDomainMap());
447 VectorFactory::Build(Op->getDomainMap());
449 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
455 ImportFactory::Build(Op->getDomainMap(), Op->getColMap());
458 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
462 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
467 #if 1 // 1=new "fast code, 0=old "slow", but safe code 468 #if 0 // not necessary - remove me 469 if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
472 std::vector<LocalOrdinal> NewRightLocal(right->getColMap()->getNodeNumElements(), Teuchos::as<LocalOrdinal>(left->getColMap()->getNodeNumElements()+1));
475 for (
size_t j=0; j < right->getColMap()->getNodeNumElements(); j++) {
476 while ( (i < Teuchos::as<LocalOrdinal>(left->getColMap()->getNodeNumElements())) &&
477 (left->getColMap()->getGlobalElement(i) < right->getColMap()->getGlobalElement(j)) ) i++;
478 if (left->getColMap()->getGlobalElement(i) == right->getColMap()->getGlobalElement(j)) {
479 NewRightLocal[j] = i;
484 std::vector<Scalar> temp_array(left->getColMap()->getNodeNumElements()+1, 0.0);
486 for(
size_t n=0; n<right->getNodeNumRows(); n++) {
492 left->getLocalRowView (n, lindices_left, lvals_left);
493 right->getLocalRowView(n, lindices_right, lvals_right);
495 for(
size_t j=0; j<Teuchos::as<size_t>(lindices_right.
size()); j++) {
496 temp_array[NewRightLocal[lindices_right[j] ] ] = lvals_right[j];
498 for (
size_t j=0; j < Teuchos::as<size_t>(lindices_left.
size()); j++) {
499 InnerProd_local[lindices_left[j]] += temp_array[lindices_left[j] ]*lvals_left[j];
501 for (
size_t j=0; j < Teuchos::as<size_t>(lindices_right.
size()); j++) {
502 temp_array[NewRightLocal[lindices_right[j] ] ] = 0.0;
507 ExportFactory::Build(left->getColMap(), left->getDomainMap());
510 VectorFactory::Build(left->getDomainMap());
512 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
518 ImportFactory::Build(left->getDomainMap(), left->getColMap());
521 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
525 #endif // end remove me 526 if(InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
527 size_t szNewLeftLocal = TEUCHOS_MAX(left->getColMap()->getNodeNumElements(), right->getColMap()->getNodeNumElements());
531 for (
size_t i=0; i < left->getColMap()->getNodeNumElements(); i++) {
532 while ( (j < Teuchos::as<LocalOrdinal>(right->getColMap()->getNodeNumElements())) &&
533 (right->getColMap()->getGlobalElement(j) < left->getColMap()->getGlobalElement(i)) ) j++;
534 if (right->getColMap()->getGlobalElement(j) == left->getColMap()->getGlobalElement(i)) {
535 (*NewLeftLocal)[i] = j;
546 for(
size_t n=0; n<left->getNodeNumRows(); n++) {
552 left->getLocalRowView (n, lindices_left, lvals_left);
553 right->getLocalRowView(n, lindices_right, lvals_right);
555 for(
size_t i=0; i<Teuchos::as<size_t>(lindices_left.
size()); i++) {
556 (*temp_array)[(*NewLeftLocal)[lindices_left[i] ] ] = lvals_left[i];
558 for (
size_t i=0; i < Teuchos::as<size_t>(lindices_right.
size()); i++) {
559 InnerProd_local[lindices_right[i]] += (*temp_array)[lindices_right[i] ] * lvals_right[i];
561 for (
size_t i=0; i < Teuchos::as<size_t>(lindices_left.
size()); i++) {
562 (*temp_array)[(*NewLeftLocal)[lindices_left[i] ] ] = 0.0;
568 ExportFactory::Build(right->getColMap(), right->getDomainMap());
571 VectorFactory::Build(right->getDomainMap());
573 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
579 ImportFactory::Build(right->getDomainMap(), right->getColMap());
581 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
586 #else // old "safe" code 587 if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
597 for(
size_t n=0; n<left->getNodeNumRows(); n++)
600 left->getLocalRowView (n, lindices_left, lvals_left);
601 right->getLocalRowView(n, lindices_right, lvals_right);
603 for(
size_t i=0; i<Teuchos::as<size_t>(lindices_left.
size()); i++)
605 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
606 for(
size_t j=0; j<Teuchos::as<size_t>(lindices_right.
size()); j++)
608 GlobalOrdinal right_gid= right->getColMap()->getGlobalElement(lindices_right[j]);
609 if(left_gid == right_gid)
611 InnerProd_local[lindices_left[i]] += lvals_left[i]*lvals_right[j];
620 ExportFactory::Build(left->getColMap(), left->getDomainMap());
623 VectorFactory::Build(left->getDomainMap());
625 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
631 ImportFactory::Build(left->getDomainMap(), left->getColMap());
634 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
636 else if(InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
644 for(
size_t n=0; n<left->getNodeNumRows(); n++)
646 left->getLocalRowView(n, lindices_left, lvals_left);
647 right->getLocalRowView(n, lindices_right, lvals_right);
649 for(
size_t i=0; i<Teuchos::as<size_t>(lindices_left.
size()); i++)
651 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
652 for(
size_t j=0; j<Teuchos::as<size_t>(lindices_right.
size()); j++)
654 GlobalOrdinal right_gid= right->getColMap()->getGlobalElement(lindices_right[j]);
655 if(left_gid == right_gid)
657 InnerProd_local[lindices_right[j]] += lvals_left[i]*lvals_right[j];
666 ExportFactory::Build(right->getColMap(), right->getDomainMap());
669 VectorFactory::Build(right->getDomainMap());
671 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
677 ImportFactory::Build(right->getDomainMap(), right->getColMap());
680 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
688 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
690 std::cout <<
"TODO: remove me" << std::endl;
693 template <
class Scalar,
class LocalOrdinal,
class GlobalOrdinal,
class Node>
static Teuchos::ArrayRCP< Scalar > GetMatrixDiagonal(const Xpetra::Matrix< Scalar, LocalOrdinal, GlobalOrdinal, Node > &A)
void DeclareInput(Level &fineLevel, Level &coarseLevel) const
Input.
#define MueLu_sumAll(rcpComm, in, out)
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
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)
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())
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
Timer to be used in factories. Similar to Monitor but with additional timers.
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Namespace for MueLu classes and methods.
#define MueLu_minAll(rcpComm, in, out)
Print even more statistics.
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 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)
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 std::string PrintMatrixInfo(const Matrix &A, const std::string &msgTag, RCP< const Teuchos::ParameterList > params=Teuchos::null)
RCP< const ParameterList > GetValidParameterList() const
Return a const parameter list of valid parameters that setParameterList() will accept.
MinimizationNorm GetMinimizationMode()
return minimization mode
void Set(const std::string &ename, const T &entry, const FactoryBase *factory=NoFactory::get())
static magnitudeType magnitude(T a)
void ReUseDampingParameters(bool bReuse)
void BuildP(Level &fineLevel, Level &coarseLevel) const
Abstract Build method.
const RCP< const FactoryManagerBase > GetFactoryManager()
returns the current factory manager
Exception throws to report errors in the internal logical of the program.
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