129 FactoryMonitor m(*
this,
"Prolongator smoothing (PG-AMG)", coarseLevel);
132 RCP<Matrix> A = Get< RCP<Matrix> >(fineLevel,
"A");
137 RCP<const FactoryBase> initialPFact = GetFactory(
"P");
138 if (initialPFact == Teuchos::null) { initialPFact = coarseLevel.
GetFactoryManager()->GetFactory(
"Ptent"); }
139 RCP<Matrix> Ptent = coarseLevel.
Get< RCP<Matrix> >(
"P", initialPFact.get());
142 if(restrictionMode_) {
148 bool doFillComplete=
true;
149 bool optimizeStorage=
true;
150 RCP<Matrix> DinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *Ptent,
false, GetOStream(
Statistics2), doFillComplete, optimizeStorage);
153 optimizeStorage=
false;
159 Teuchos::RCP<Vector > RowBasedOmega = Teuchos::null;
161 const ParameterList & pL = GetParameterList();
162 bool bReUseRowBasedOmegas = pL.get<
bool>(
"ReUseRowBasedOmegas");
163 if(restrictionMode_ ==
false || bReUseRowBasedOmegas ==
false) {
166 ComputeRowBasedOmega(fineLevel, coarseLevel, A, Ptent, DinvAP0, RowBasedOmega);
170 RowBasedOmega = coarseLevel.
Get<Teuchos::RCP<Vector > >(
"RowBasedOmega",
this);
178 Teuchos::RCP<const Export> exporter =
179 ExportFactory::Build(RowBasedOmega->getMap(), A->getRangeMap());
181 Teuchos::RCP<Vector > noRowBasedOmega =
182 VectorFactory::Build(A->getRangeMap());
184 noRowBasedOmega->doExport(*RowBasedOmega, *exporter, Xpetra::INSERT);
189 Teuchos::RCP<const Import > importer =
190 ImportFactory::Build(A->getRangeMap(), A->getRowMap());
193 RowBasedOmega->doImport(*noRowBasedOmega, *importer, Xpetra::INSERT);
196 Teuchos::ArrayRCP< Scalar > RowBasedOmega_local = RowBasedOmega->getDataNonConst(0);
199 RCP<Matrix> P_smoothed = Teuchos::null;
202 Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::TwoMatrixAdd(*Ptent,
false, Teuchos::ScalarTraits<Scalar>::one(),
203 *DinvAP0,
false, -Teuchos::ScalarTraits<Scalar>::one(),
205 P_smoothed->fillComplete(Ptent->getDomainMap(), Ptent->getRangeMap());
209 RCP<ParameterList> params = rcp(
new ParameterList());
210 params->set(
"printLoadBalancingInfo",
true);
213 if (!restrictionMode_) {
215 Set(coarseLevel,
"P", P_smoothed);
223 Set(coarseLevel,
"RfromPfactory", dummy);
229 if (Ptent->IsView(
"stridedMaps"))
230 P_smoothed->CreateView(
"stridedMaps", Ptent);
235 Set(coarseLevel,
"R", R);
241 if (Ptent->IsView(
"stridedMaps"))
242 R->CreateView(
"stridedMaps", Ptent,
true);
249 FactoryMonitor m(*
this,
"PgPFactory::ComputeRowBasedOmega", coarseLevel);
251 typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType Magnitude;
252 Scalar sZero = Teuchos::ScalarTraits<Scalar>::zero();
253 Magnitude mZero = Teuchos::ScalarTraits<Scalar>::magnitude(sZero);
255 Teuchos::RCP<Vector > Numerator = Teuchos::null;
256 Teuchos::RCP<Vector > Denominator = Teuchos::null;
258 const ParameterList & pL = GetParameterList();
275 bool doFillComplete=
true;
276 bool optimizeStorage=
false;
277 RCP<Matrix> AP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *P0,
false, GetOStream(
Statistics2), doFillComplete, optimizeStorage);
280 RCP<Matrix> ADinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *DinvAP0,
false, GetOStream(
Statistics2), doFillComplete, optimizeStorage);
282 Numerator = VectorFactory::Build(ADinvAP0->getColMap(),
true);
283 Denominator = VectorFactory::Build(ADinvAP0->getColMap(),
true);
284 MultiplyAll(AP0, ADinvAP0, Numerator);
285 MultiplySelfAll(ADinvAP0, Denominator);
296 Numerator = VectorFactory::Build(DinvAP0->getColMap(),
true);
297 Denominator = VectorFactory::Build(DinvAP0->getColMap(),
true);
298 MultiplyAll(P0, DinvAP0, Numerator);
299 MultiplySelfAll(DinvAP0, Denominator);
313 bool doFillComplete=
true;
314 bool optimizeStorage=
true;
316 RCP<Matrix> DinvADinvAP0 = Xpetra::MatrixMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>::Multiply(*A,
false, *DinvAP0,
false, GetOStream(
Statistics2), doFillComplete, optimizeStorage);
318 diagA = Teuchos::ArrayRCP<Scalar>();
320 Numerator = VectorFactory::Build(DinvADinvAP0->getColMap(),
true);
321 Denominator = VectorFactory::Build(DinvADinvAP0->getColMap(),
true);
322 MultiplyAll(DinvAP0, DinvADinvAP0, Numerator);
323 MultiplySelfAll(DinvADinvAP0, Denominator);
327 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"MueLu::PgPFactory::Build: minimization mode not supported. error");
332 Teuchos::RCP<Vector > ColBasedOmega =
333 VectorFactory::Build(Numerator->getMap(),
true);
335 ColBasedOmega->putScalar(-666);
337 Teuchos::ArrayRCP< const Scalar > Numerator_local = Numerator->getData(0);
338 Teuchos::ArrayRCP< const Scalar > Denominator_local = Denominator->getData(0);
339 Teuchos::ArrayRCP< Scalar > ColBasedOmega_local = ColBasedOmega->getDataNonConst(0);
342 Magnitude min_local = 1000000.0;
343 Magnitude max_local = 0.0;
344 for(
LocalOrdinal i = 0; i < Teuchos::as<LocalOrdinal>(Numerator->getLocalLength()); i++) {
345 if(Teuchos::ScalarTraits<Scalar>::magnitude(Denominator_local[i]) == mZero)
347 ColBasedOmega_local[i] = 0.0;
352 ColBasedOmega_local[i] = Numerator_local[i] / Denominator_local[i];
355 if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) < mZero) {
356 ColBasedOmega_local[i] = Teuchos::ScalarTraits<Scalar>::zero();
364 if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) >= 0.8) {
365 ColBasedOmega_local[i] = 0.0;
368 if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) < min_local)
369 { min_local = Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]); }
370 if(Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]) > max_local)
371 { max_local = Teuchos::ScalarTraits<Scalar>::magnitude(ColBasedOmega_local[i]); }
379 MueLu_sumAll(A->getRowMap()->getComm(), zero_local, zero_all);
380 MueLu_sumAll(A->getRowMap()->getComm(), nan_local, nan_all);
381 MueLu_minAll(A->getRowMap()->getComm(), min_local, min_all);
382 MueLu_maxAll(A->getRowMap()->getComm(), max_local, max_all);
389 default: GetOStream(
Statistics1) <<
"unknown)" << std::endl;
break;
391 GetOStream(
Statistics1) <<
"Damping parameter: min = " << min_all <<
", max = " << max_all << std::endl;
392 GetOStream(
Statistics) <<
"# negative omegas: " << zero_all <<
" out of " << ColBasedOmega->getGlobalLength() <<
" column-based omegas" << std::endl;
393 GetOStream(
Statistics) <<
"# NaNs: " << nan_all <<
" out of " << ColBasedOmega->getGlobalLength() <<
" column-based omegas" << std::endl;
396 if(coarseLevel.
IsRequested(
"ColBasedOmega",
this)) {
397 coarseLevel.
Set(
"ColBasedOmega", ColBasedOmega,
this);
403 VectorFactory::Build(DinvAP0->getRowMap(),
true);
405 RowBasedOmega->putScalar(-666);
407 bool bAtLeastOneDefined =
false;
408 Teuchos::ArrayRCP< Scalar > RowBasedOmega_local = RowBasedOmega->getDataNonConst(0);
409 for(
LocalOrdinal row = 0; row<Teuchos::as<LocalOrdinal>(A->getLocalNumRows()); row++) {
410 Teuchos::ArrayView<const LocalOrdinal> lindices;
411 Teuchos::ArrayView<const Scalar> lvals;
412 DinvAP0->getLocalRowView(row, lindices, lvals);
413 bAtLeastOneDefined =
false;
414 for(
size_t j=0; j<Teuchos::as<size_t>(lindices.size()); j++) {
415 Scalar omega = ColBasedOmega_local[lindices[j]];
416 if (Teuchos::ScalarTraits<Scalar>::magnitude(omega) != -666) {
417 bAtLeastOneDefined =
true;
418 if(Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row]) == -666) RowBasedOmega_local[row] = omega;
419 else if(Teuchos::ScalarTraits<Scalar>::magnitude(omega) < Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row])) RowBasedOmega_local[row] = omega;
422 if(bAtLeastOneDefined ==
true) {
423 if(Teuchos::ScalarTraits<Scalar>::magnitude(RowBasedOmega_local[row]) < mZero)
424 RowBasedOmega_local[row] = sZero;
428 if(coarseLevel.
IsRequested(
"RowBasedOmega",
this)) {
429 Set(coarseLevel,
"RowBasedOmega", RowBasedOmega);
475 TEUCHOS_TEST_FOR_EXCEPTION(!left->getDomainMap()->isSameAs(*right->getDomainMap()),
Exceptions::RuntimeError,
"MueLu::PgPFactory::MultiplyAll: domain maps of left and right do not match. Error.");
476 TEUCHOS_TEST_FOR_EXCEPTION(!left->getRowMap()->isSameAs(*right->getRowMap()),
Exceptions::RuntimeError,
"MueLu::PgPFactory::MultiplyAll: row maps of left and right do not match. Error.");
479 if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
482 std::vector<LocalOrdinal> NewRightLocal(right->getColMap()->getLocalNumElements(), Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements()+1));
485 for (
size_t j=0; j < right->getColMap()->getLocalNumElements(); j++) {
486 while ( (i < Teuchos::as<LocalOrdinal>(left->getColMap()->getLocalNumElements())) &&
487 (left->getColMap()->getGlobalElement(i) < right->getColMap()->getGlobalElement(j)) ) i++;
488 if (left->getColMap()->getGlobalElement(i) == right->getColMap()->getGlobalElement(j)) {
489 NewRightLocal[j] = i;
493 Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
494 std::vector<Scalar> temp_array(left->getColMap()->getLocalNumElements()+1, 0.0);
496 for(
size_t n=0; n<right->getLocalNumRows(); n++) {
497 Teuchos::ArrayView<const LocalOrdinal> lindices_left;
498 Teuchos::ArrayView<const Scalar> lvals_left;
499 Teuchos::ArrayView<const LocalOrdinal> lindices_right;
500 Teuchos::ArrayView<const Scalar> lvals_right;
502 left->getLocalRowView (n, lindices_left, lvals_left);
503 right->getLocalRowView(n, lindices_right, lvals_right);
505 for(
size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++) {
506 temp_array[NewRightLocal[lindices_right[j] ] ] = lvals_right[j];
508 for (
size_t j=0; j < Teuchos::as<size_t>(lindices_left.size()); j++) {
509 InnerProd_local[lindices_left[j]] += temp_array[lindices_left[j] ]*lvals_left[j];
511 for (
size_t j=0; j < Teuchos::as<size_t>(lindices_right.size()); j++) {
512 temp_array[NewRightLocal[lindices_right[j] ] ] = 0.0;
516 Teuchos::RCP<const Export> exporter =
517 ExportFactory::Build(left->getColMap(), left->getDomainMap());
519 Teuchos::RCP<Vector > nonoverlap =
520 VectorFactory::Build(left->getDomainMap());
522 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
527 Teuchos::RCP<const Import > importer =
528 ImportFactory::Build(left->getDomainMap(), left->getColMap());
531 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
536 if(InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
537 size_t szNewLeftLocal = TEUCHOS_MAX(left->getColMap()->getLocalNumElements(), right->getColMap()->getLocalNumElements());
538 Teuchos::RCP<std::vector<LocalOrdinal> > NewLeftLocal = Teuchos::rcp(
new std::vector<LocalOrdinal>(szNewLeftLocal, Teuchos::as<LocalOrdinal>(right->getColMap()->getMaxLocalIndex()+1)));
541 for (
size_t i=0; i < left->getColMap()->getLocalNumElements(); i++) {
542 while ( (j < Teuchos::as<LocalOrdinal>(right->getColMap()->getLocalNumElements())) &&
543 (right->getColMap()->getGlobalElement(j) < left->getColMap()->getGlobalElement(i)) ) j++;
544 if (right->getColMap()->getGlobalElement(j) == left->getColMap()->getGlobalElement(i)) {
545 (*NewLeftLocal)[i] = j;
553 Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
554 Teuchos::RCP<std::vector<Scalar> > temp_array = Teuchos::rcp(
new std::vector<Scalar>(right->getColMap()->getMaxLocalIndex()+2, 0.0));
556 for(
size_t n=0; n<left->getLocalNumRows(); n++) {
557 Teuchos::ArrayView<const LocalOrdinal> lindices_left;
558 Teuchos::ArrayView<const Scalar> lvals_left;
559 Teuchos::ArrayView<const LocalOrdinal> lindices_right;
560 Teuchos::ArrayView<const Scalar> lvals_right;
562 left->getLocalRowView (n, lindices_left, lvals_left);
563 right->getLocalRowView(n, lindices_right, lvals_right);
565 for(
size_t i=0; i<Teuchos::as<size_t>(lindices_left.size()); i++) {
566 (*temp_array)[(*NewLeftLocal)[lindices_left[i] ] ] = lvals_left[i];
568 for (
size_t i=0; i < Teuchos::as<size_t>(lindices_right.size()); i++) {
569 InnerProd_local[lindices_right[i]] += (*temp_array)[lindices_right[i] ] * lvals_right[i];
571 for (
size_t i=0; i < Teuchos::as<size_t>(lindices_left.size()); i++) {
572 (*temp_array)[(*NewLeftLocal)[lindices_left[i] ] ] = 0.0;
575 InnerProd_local = Teuchos::ArrayRCP< Scalar >();
577 Teuchos::RCP<const Export> exporter =
578 ExportFactory::Build(right->getColMap(), right->getDomainMap());
580 Teuchos::RCP<Vector> nonoverlap =
581 VectorFactory::Build(right->getDomainMap());
583 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
588 Teuchos::RCP<const Import > importer =
589 ImportFactory::Build(right->getDomainMap(), right->getColMap());
591 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
593 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"MueLu::PgPFactory::MultiplyAll: map of InnerProdVec must be same as column map of left operator? Error.");
597 if(InnerProdVec->getMap()->isSameAs(*left->getColMap())) {
599 Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
602 Teuchos::ArrayView<const LocalOrdinal> lindices_left;
603 Teuchos::ArrayView<const Scalar> lvals_left;
604 Teuchos::ArrayView<const LocalOrdinal> lindices_right;
605 Teuchos::ArrayView<const Scalar> lvals_right;
607 for(
size_t n=0; n<left->getLocalNumRows(); n++)
610 left->getLocalRowView (n, lindices_left, lvals_left);
611 right->getLocalRowView(n, lindices_right, lvals_right);
613 for(
size_t i=0; i<Teuchos::as<size_t>(lindices_left.size()); i++)
615 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
616 for(
size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++)
618 GlobalOrdinal right_gid= right->getColMap()->getGlobalElement(lindices_right[j]);
619 if(left_gid == right_gid)
621 InnerProd_local[lindices_left[i]] += lvals_left[i]*lvals_right[j];
629 Teuchos::RCP<const Export> exporter =
630 ExportFactory::Build(left->getColMap(), left->getDomainMap());
632 Teuchos::RCP<Vector > nonoverlap =
633 VectorFactory::Build(left->getDomainMap());
635 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
640 Teuchos::RCP<const Import > importer =
641 ImportFactory::Build(left->getDomainMap(), left->getColMap());
644 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
646 else if(InnerProdVec->getMap()->isSameAs(*right->getColMap())) {
647 Teuchos::ArrayRCP< Scalar > InnerProd_local = InnerProdVec->getDataNonConst(0);
649 Teuchos::ArrayView<const LocalOrdinal> lindices_left;
650 Teuchos::ArrayView<const Scalar> lvals_left;
651 Teuchos::ArrayView<const LocalOrdinal> lindices_right;
652 Teuchos::ArrayView<const Scalar> lvals_right;
654 for(
size_t n=0; n<left->getLocalNumRows(); n++)
656 left->getLocalRowView(n, lindices_left, lvals_left);
657 right->getLocalRowView(n, lindices_right, lvals_right);
659 for(
size_t i=0; i<Teuchos::as<size_t>(lindices_left.size()); i++)
661 GlobalOrdinal left_gid = left->getColMap()->getGlobalElement(lindices_left[i]);
662 for(
size_t j=0; j<Teuchos::as<size_t>(lindices_right.size()); j++)
664 GlobalOrdinal right_gid= right->getColMap()->getGlobalElement(lindices_right[j]);
665 if(left_gid == right_gid)
667 InnerProd_local[lindices_right[j]] += lvals_left[i]*lvals_right[j];
675 Teuchos::RCP<const Export> exporter =
676 ExportFactory::Build(right->getColMap(), right->getDomainMap());
678 Teuchos::RCP<Vector> nonoverlap =
679 VectorFactory::Build(right->getDomainMap());
681 nonoverlap->doExport(*InnerProdVec, *exporter, Xpetra::ADD);
686 Teuchos::RCP<const Import > importer =
687 ImportFactory::Build(right->getDomainMap(), right->getColMap());
690 InnerProdVec->doImport(*nonoverlap, *importer, Xpetra::INSERT);
693 TEUCHOS_TEST_FOR_EXCEPTION(
true,
Exceptions::RuntimeError,
"MueLu::PgPFactory::MultiplyAll: map of InnerProdVec must be same as column map of left or right operator? Error.");