350template<
class ScalarType,
class MV>
357 using Teuchos::rcp_const_cast;
369 const ST defaultLambdaMax = STS::nan ();
370 const ST defaultLambdaMin = STS::nan ();
377 const ST defaultEigRatio = Teuchos::as<ST> (30);
378 const MT defaultBoostFactor =
static_cast<MT> (1.1);
379 const ST defaultMinDiagVal = STS::eps ();
380 const int defaultNumIters = 1;
381 const int defaultEigMaxIters = 10;
382 const MT defaultEigRelTolerance = Teuchos::ScalarTraits<MT>::zero ();
383 const bool defaultEigKeepVectors =
false;
384 const int defaultEigNormalizationFreq = 1;
385 const bool defaultZeroStartingSolution =
true;
386 const bool defaultAssumeMatrixUnchanged =
false;
387 const bool defaultTextbookAlgorithm =
false;
388 const bool defaultComputeMaxResNorm =
false;
389 const bool defaultComputeSpectralRadius =
true;
390 const bool defaultCkUseNativeSpMV =
false;
391 const bool defaultDebug =
false;
397 RCP<const V> userInvDiagCopy;
398 ST lambdaMax = defaultLambdaMax;
399 ST lambdaMin = defaultLambdaMin;
400 ST eigRatio = defaultEigRatio;
401 MT boostFactor = defaultBoostFactor;
402 ST minDiagVal = defaultMinDiagVal;
403 int numIters = defaultNumIters;
404 int eigMaxIters = defaultEigMaxIters;
405 MT eigRelTolerance = defaultEigRelTolerance;
406 bool eigKeepVectors = defaultEigKeepVectors;
407 int eigNormalizationFreq = defaultEigNormalizationFreq;
408 bool zeroStartingSolution = defaultZeroStartingSolution;
409 bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
410 bool textbookAlgorithm = defaultTextbookAlgorithm;
411 bool computeMaxResNorm = defaultComputeMaxResNorm;
412 bool computeSpectralRadius = defaultComputeSpectralRadius;
413 bool ckUseNativeSpMV = defaultCkUseNativeSpMV;
414 bool debug = defaultDebug;
421 if (plist.isType<
bool> (
"debug")) {
422 debug = plist.get<
bool> (
"debug");
424 else if (plist.isType<
int> (
"debug")) {
425 const int debugInt = plist.get<
bool> (
"debug");
426 debug = debugInt != 0;
439 const char opInvDiagLabel[] =
"chebyshev: operator inv diagonal";
440 if (plist.isParameter (opInvDiagLabel)) {
442 RCP<const V> userInvDiag;
444 if (plist.isType<
const V*> (opInvDiagLabel)) {
445 const V* rawUserInvDiag =
446 plist.get<
const V*> (opInvDiagLabel);
448 userInvDiag = rcp (rawUserInvDiag,
false);
450 else if (plist.isType<
const V*> (opInvDiagLabel)) {
451 V* rawUserInvDiag = plist.get<
V*> (opInvDiagLabel);
453 userInvDiag = rcp (
const_cast<const V*
> (rawUserInvDiag),
false);
455 else if (plist.isType<RCP<const V>> (opInvDiagLabel)) {
456 userInvDiag = plist.get<RCP<const V> > (opInvDiagLabel);
458 else if (plist.isType<RCP<V>> (opInvDiagLabel)) {
459 RCP<V> userInvDiagNonConst =
460 plist.get<RCP<V> > (opInvDiagLabel);
461 userInvDiag = rcp_const_cast<const V> (userInvDiagNonConst);
463 else if (plist.isType<
const V> (opInvDiagLabel)) {
464 const V& userInvDiagRef = plist.get<
const V> (opInvDiagLabel);
465 userInvDiagCopy = rcp (
new V (userInvDiagRef, Teuchos::Copy));
466 userInvDiag = userInvDiagCopy;
468 else if (plist.isType<
V> (opInvDiagLabel)) {
469 V& userInvDiagNonConstRef = plist.get<
V> (opInvDiagLabel);
470 const V& userInvDiagRef =
const_cast<const V&
> (userInvDiagNonConstRef);
471 userInvDiagCopy = rcp (
new V (userInvDiagRef, Teuchos::Copy));
472 userInvDiag = userInvDiagCopy;
482 if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
483 userInvDiagCopy = rcp (
new V (*userInvDiag, Teuchos::Copy));
493 if (plist.isParameter (
"chebyshev: use native spmv"))
494 ckUseNativeSpMV = plist.get(
"chebyshev: use native spmv", ckUseNativeSpMV);
499 if (plist.isParameter (
"chebyshev: max eigenvalue")) {
500 if (plist.isType<
double>(
"chebyshev: max eigenvalue"))
501 lambdaMax = plist.get<
double> (
"chebyshev: max eigenvalue");
503 lambdaMax = plist.get<
ST> (
"chebyshev: max eigenvalue");
504 TEUCHOS_TEST_FOR_EXCEPTION(
505 STS::isnaninf (lambdaMax), std::invalid_argument,
506 "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" "
507 "parameter is NaN or Inf. This parameter is optional, but if you "
508 "choose to supply it, it must have a finite value.");
510 if (plist.isParameter (
"chebyshev: min eigenvalue")) {
511 if (plist.isType<
double>(
"chebyshev: min eigenvalue"))
512 lambdaMin = plist.get<
double> (
"chebyshev: min eigenvalue");
514 lambdaMin = plist.get<
ST> (
"chebyshev: min eigenvalue");
515 TEUCHOS_TEST_FOR_EXCEPTION(
516 STS::isnaninf (lambdaMin), std::invalid_argument,
517 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" "
518 "parameter is NaN or Inf. This parameter is optional, but if you "
519 "choose to supply it, it must have a finite value.");
523 if (plist.isParameter (
"smoother: Chebyshev alpha")) {
524 if (plist.isType<
double>(
"smoother: Chebyshev alpha"))
525 eigRatio = plist.get<
double> (
"smoother: Chebyshev alpha");
527 eigRatio = plist.get<
ST> (
"smoother: Chebyshev alpha");
530 eigRatio = plist.get (
"chebyshev: ratio eigenvalue", eigRatio);
531 TEUCHOS_TEST_FOR_EXCEPTION(
532 STS::isnaninf (eigRatio), std::invalid_argument,
533 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" "
534 "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. "
535 "This parameter is optional, but if you choose to supply it, it must have "
542 TEUCHOS_TEST_FOR_EXCEPTION(
543 STS::real (eigRatio) < STS::real (STS::one ()),
544 std::invalid_argument,
545 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\""
546 "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, "
547 "but you supplied the value " << eigRatio <<
".");
552 const char paramName[] =
"chebyshev: boost factor";
554 if (plist.isParameter (paramName)) {
555 if (plist.isType<
MT> (paramName)) {
556 boostFactor = plist.get<
MT> (paramName);
558 else if (! std::is_same<double, MT>::value &&
559 plist.isType<
double> (paramName)) {
560 const double dblBF = plist.get<
double> (paramName);
561 boostFactor =
static_cast<MT> (dblBF);
564 TEUCHOS_TEST_FOR_EXCEPTION
565 (
true, std::invalid_argument,
566 "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\""
567 "parameter must have type magnitude_type (MT) or double.");
577 plist.set (paramName, defaultBoostFactor);
579 TEUCHOS_TEST_FOR_EXCEPTION
580 (boostFactor < Teuchos::ScalarTraits<MT>::one (), std::invalid_argument,
581 "Ifpack2::Chebyshev::setParameters: \"" << paramName <<
"\" parameter "
582 "must be >= 1, but you supplied the value " << boostFactor <<
".");
586 minDiagVal = plist.get (
"chebyshev: min diagonal value", minDiagVal);
587 TEUCHOS_TEST_FOR_EXCEPTION(
588 STS::isnaninf (minDiagVal), std::invalid_argument,
589 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" "
590 "parameter is NaN or Inf. This parameter is optional, but if you choose "
591 "to supply it, it must have a finite value.");
594 if (plist.isParameter (
"smoother: sweeps")) {
595 numIters = plist.get<
int> (
"smoother: sweeps");
597 if (plist.isParameter (
"relaxation: sweeps")) {
598 numIters = plist.get<
int> (
"relaxation: sweeps");
600 numIters = plist.get (
"chebyshev: degree", numIters);
601 TEUCHOS_TEST_FOR_EXCEPTION(
602 numIters < 0, std::invalid_argument,
603 "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also "
604 "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a "
605 "nonnegative integer. You gave a value of " << numIters <<
".");
608 if (plist.isParameter (
"eigen-analysis: iterations")) {
609 eigMaxIters = plist.get<
int> (
"eigen-analysis: iterations");
611 eigMaxIters = plist.get (
"chebyshev: eigenvalue max iterations", eigMaxIters);
612 TEUCHOS_TEST_FOR_EXCEPTION(
613 eigMaxIters < 0, std::invalid_argument,
614 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations"
615 "\" parameter (also called \"eigen-analysis: iterations\") must be a "
616 "nonnegative integer. You gave a value of " << eigMaxIters <<
".");
618 if (plist.isType<
double>(
"chebyshev: eigenvalue relative tolerance"))
619 eigRelTolerance = Teuchos::as<MT>(plist.get<
double> (
"chebyshev: eigenvalue relative tolerance"));
620 else if (plist.isType<
MT>(
"chebyshev: eigenvalue relative tolerance"))
621 eigRelTolerance = plist.get<
MT> (
"chebyshev: eigenvalue relative tolerance");
622 else if (plist.isType<
ST>(
"chebyshev: eigenvalue relative tolerance"))
623 eigRelTolerance = Teuchos::ScalarTraits<ST>::magnitude(plist.get<
ST> (
"chebyshev: eigenvalue relative tolerance"));
625 eigKeepVectors = plist.get (
"chebyshev: eigenvalue keep vectors", eigKeepVectors);
627 eigNormalizationFreq = plist.get (
"chebyshev: eigenvalue normalization frequency", eigNormalizationFreq);
628 TEUCHOS_TEST_FOR_EXCEPTION(
629 eigNormalizationFreq < 0, std::invalid_argument,
630 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue normalization frequency"
631 "\" parameter must be a "
632 "nonnegative integer. You gave a value of " << eigNormalizationFreq <<
".")
634 zeroStartingSolution = plist.get (
"chebyshev: zero starting solution",
635 zeroStartingSolution);
636 assumeMatrixUnchanged = plist.get (
"chebyshev: assume matrix does not change",
637 assumeMatrixUnchanged);
641 if (plist.isParameter (
"chebyshev: textbook algorithm")) {
642 textbookAlgorithm = plist.get<
bool> (
"chebyshev: textbook algorithm");
644 if (plist.isParameter (
"chebyshev: compute max residual norm")) {
645 computeMaxResNorm = plist.get<
bool> (
"chebyshev: compute max residual norm");
647 if (plist.isParameter (
"chebyshev: compute spectral radius")) {
648 computeSpectralRadius = plist.get<
bool> (
"chebyshev: compute spectral radius");
654 TEUCHOS_TEST_FOR_EXCEPTION
655 (plist.isType<
bool> (
"chebyshev: use block mode") &&
656 ! plist.get<
bool> (
"chebyshev: use block mode"),
657 std::invalid_argument,
658 "Ifpack2::Chebyshev requires that if you set \"chebyshev: use "
659 "block mode\" at all, you must set it to false. "
660 "Ifpack2::Chebyshev does not implement Ifpack's block mode.");
661 TEUCHOS_TEST_FOR_EXCEPTION
662 (plist.isType<
bool> (
"chebyshev: solve normal equations") &&
663 ! plist.get<
bool> (
"chebyshev: solve normal equations"),
664 std::invalid_argument,
665 "Ifpack2::Chebyshev does not and will never implement the Ifpack "
666 "parameter \"chebyshev: solve normal equations\". If you want to "
667 "solve the normal equations, construct a Tpetra::Operator that "
668 "implements A^* A, and use Chebyshev to solve A^* A x = A^* b.");
676 std::string eigenAnalysisType (
"power-method");
677 if (plist.isParameter (
"eigen-analysis: type")) {
678 eigenAnalysisType = plist.get<std::string> (
"eigen-analysis: type");
679 TEUCHOS_TEST_FOR_EXCEPTION(
680 eigenAnalysisType !=
"power-method" &&
681 eigenAnalysisType !=
"power method" &&
682 eigenAnalysisType !=
"cg",
683 std::invalid_argument,
684 "Ifpack2::Chebyshev: Ifpack2 only supports \"power method\" and \"cg\" for \"eigen-analysis: type\".");
688 userInvDiag_ = userInvDiagCopy;
689 userLambdaMax_ = lambdaMax;
690 userLambdaMin_ = lambdaMin;
691 userEigRatio_ = eigRatio;
692 boostFactor_ =
static_cast<MT> (boostFactor);
693 minDiagVal_ = minDiagVal;
694 numIters_ = numIters;
695 eigMaxIters_ = eigMaxIters;
696 eigRelTolerance_ = eigRelTolerance;
697 eigKeepVectors_ = eigKeepVectors;
698 eigNormalizationFreq_ = eigNormalizationFreq;
699 eigenAnalysisType_ = eigenAnalysisType;
700 zeroStartingSolution_ = zeroStartingSolution;
701 assumeMatrixUnchanged_ = assumeMatrixUnchanged;
702 textbookAlgorithm_ = textbookAlgorithm;
703 computeMaxResNorm_ = computeMaxResNorm;
704 computeSpectralRadius_ = computeSpectralRadius;
705 ckUseNativeSpMV_ = ckUseNativeSpMV;
711 if (A_.is_null () || A_->getComm ().is_null ()) {
717 myRank = A_->getComm ()->getRank ();
721 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
724 using Teuchos::oblackholestream;
725 RCP<oblackholestream> blackHole (
new oblackholestream ());
726 out_ = Teuchos::getFancyOStream (blackHole);
731 out_ = Teuchos::null;
736template<
class ScalarType,
class MV>
742 diagOffsets_ = offsets_type ();
743 savedDiagOffsets_ =
false;
745 computedLambdaMax_ = STS::nan ();
746 computedLambdaMin_ = STS::nan ();
747 eigVector_ = Teuchos::null;
748 eigVector2_ = Teuchos::null;
752template<
class ScalarType,
class MV>
755setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
757 if (A.getRawPtr () != A_.getRawPtr ()) {
758 if (! assumeMatrixUnchanged_) {
770 if (A.is_null () || A->getComm ().is_null ()) {
776 myRank = A->getComm ()->getRank ();
780 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
783 Teuchos::RCP<Teuchos::oblackholestream> blackHole (
new Teuchos::oblackholestream ());
784 out_ = Teuchos::getFancyOStream (blackHole);
789 out_ = Teuchos::null;
794template<
class ScalarType,
class MV>
803 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
804 typename MV::local_ordinal_type,
805 typename MV::global_ordinal_type,
806 typename MV::node_type> crs_matrix_type;
808 TEUCHOS_TEST_FOR_EXCEPTION(
809 A_.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::compute: The input "
810 "matrix A is null. Please call setMatrix() with a nonnull input matrix "
811 "before calling this method.");
826 if (userInvDiag_.is_null ()) {
827 Teuchos::RCP<const crs_matrix_type> A_crsMat =
828 Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
831 if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
833 const size_t lclNumRows = A_crsMat->getLocalNumRows ();
834 if (diagOffsets_.extent (0) < lclNumRows) {
835 diagOffsets_ = offsets_type ();
836 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
838 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
839 savedDiagOffsets_ =
true;
840 D_ = makeInverseDiagonal (*A_,
true);
843 D_ = makeInverseDiagonal (*A_);
846 else if (! assumeMatrixUnchanged_) {
847 if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
850 if (! savedDiagOffsets_) {
851 const size_t lclNumRows = A_crsMat->getLocalNumRows ();
852 if (diagOffsets_.extent (0) < lclNumRows) {
853 diagOffsets_ = offsets_type ();
854 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
856 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
857 savedDiagOffsets_ =
true;
860 D_ = makeInverseDiagonal (*A_,
true);
863 D_ = makeInverseDiagonal (*A_);
868 D_ = makeRangeMapVectorConst (userInvDiag_);
872 const bool computedEigenvalueEstimates =
873 STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
883 if (! assumeMatrixUnchanged_ ||
884 (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
885 ST computedLambdaMax;
886 if ((eigenAnalysisType_ ==
"power method") || (eigenAnalysisType_ ==
"power-method")) {
888 if (eigVector_.is_null()) {
889 x = Teuchos::rcp(
new V(A_->getDomainMap ()));
892 PowerMethod::computeInitialGuessForPowerMethod (*x,
false);
897 if (eigVector2_.is_null()) {
898 y = rcp(
new V(A_->getRangeMap ()));
904 Teuchos::RCP<Teuchos::FancyOStream> stream = (debug_ ? out_ : Teuchos::null);
905 computedLambdaMax = PowerMethod::powerMethodWithInitGuess (*A_, *D_, eigMaxIters_, x, y,
906 eigRelTolerance_, eigNormalizationFreq_, stream,
907 computeSpectralRadius_);
910 computedLambdaMax = cgMethod (*A_, *D_, eigMaxIters_);
911 TEUCHOS_TEST_FOR_EXCEPTION(
912 STS::isnaninf (computedLambdaMax),
914 "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue "
915 "of D^{-1} A failed, by producing Inf or NaN. This probably means that "
916 "the matrix contains Inf or NaN values, or that it is badly scaled.");
917 TEUCHOS_TEST_FOR_EXCEPTION(
918 STS::isnaninf (userEigRatio_),
920 "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
921 << endl <<
"This should be impossible." << endl <<
922 "Please report this bug to the Ifpack2 developers.");
928 const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
931 computedLambdaMax_ = computedLambdaMax;
932 computedLambdaMin_ = computedLambdaMin;
934 TEUCHOS_TEST_FOR_EXCEPTION(
935 STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
937 "Ifpack2::Chebyshev::compute: " << endl <<
938 "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
939 << endl <<
"This should be impossible." << endl <<
940 "Please report this bug to the Ifpack2 developers.");
948 lambdaMaxForApply_ = STS::isnaninf (userLambdaMax_) ? computedLambdaMax_ : userLambdaMax_;
961 lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
962 eigRatioForApply_ = userEigRatio_;
964 if (! textbookAlgorithm_) {
967 const ST one = Teuchos::as<ST> (1);
970 if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
971 lambdaMinForApply_ = one;
972 lambdaMaxForApply_ = lambdaMinForApply_;
973 eigRatioForApply_ = one;
979template<
class ScalarType,
class MV>
983 return lambdaMaxForApply_;
987template<
class ScalarType,
class MV>
991 const char prefix[] =
"Ifpack2::Chebyshev::apply: ";
994 *out_ <<
"apply: " << std::endl;
996 TEUCHOS_TEST_FOR_EXCEPTION
997 (A_.is_null (), std::runtime_error, prefix <<
"The input matrix A is null. "
998 " Please call setMatrix() with a nonnull input matrix, and then call "
999 "compute(), before calling this method.");
1000 TEUCHOS_TEST_FOR_EXCEPTION
1001 (STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
1002 prefix <<
"There is no estimate for the max eigenvalue."
1003 << std::endl << computeBeforeApplyReminder);
1004 TEUCHOS_TEST_FOR_EXCEPTION
1005 (STS::isnaninf (lambdaMinForApply_), std::runtime_error,
1006 prefix <<
"There is no estimate for the min eigenvalue."
1007 << std::endl << computeBeforeApplyReminder);
1008 TEUCHOS_TEST_FOR_EXCEPTION
1009 (STS::isnaninf (eigRatioForApply_), std::runtime_error,
1010 prefix <<
"There is no estimate for the ratio of the max "
1011 "eigenvalue to the min eigenvalue."
1012 << std::endl << computeBeforeApplyReminder);
1013 TEUCHOS_TEST_FOR_EXCEPTION
1014 (D_.is_null (), std::runtime_error, prefix <<
"The vector of inverse "
1015 "diagonal entries of the matrix has not yet been computed."
1016 << std::endl << computeBeforeApplyReminder);
1018 if (textbookAlgorithm_) {
1019 textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1020 lambdaMinForApply_, eigRatioForApply_, *D_);
1023 ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1024 lambdaMinForApply_, eigRatioForApply_, *D_);
1027 if (computeMaxResNorm_ && B.getNumVectors () > 0) {
1028 MV R (B.getMap (), B.getNumVectors ());
1029 computeResidual (R, B, *A_, X);
1030 Teuchos::Array<MT> norms (B.getNumVectors ());
1032 return *std::max_element (norms.begin (), norms.end ());
1035 return Teuchos::ScalarTraits<MT>::zero ();
1039template<
class ScalarType,
class MV>
1042print (std::ostream& out)
1044 using Teuchos::rcpFromRef;
1045 this->describe (* (Teuchos::getFancyOStream (rcpFromRef (out))),
1046 Teuchos::VERB_MEDIUM);
1049template<
class ScalarType,
class MV>
1053 const ScalarType& alpha,
1058 solve (W, alpha, D_inv, B);
1059 Tpetra::deep_copy (X, W);
1062template<
class ScalarType,
class MV>
1066 const Teuchos::ETransp mode)
1068 Tpetra::Details::residual(A,X,B,R);
1071template <
class ScalarType,
class MV>
1073Chebyshev<ScalarType, MV>::
1074solve (MV& Z,
const V& D_inv,
const MV& R) {
1075 Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
1078template<
class ScalarType,
class MV>
1080Chebyshev<ScalarType, MV>::
1081solve (MV& Z,
const ST alpha,
const V& D_inv,
const MV& R) {
1082 Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
1085template<
class ScalarType,
class MV>
1086Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
1087Chebyshev<ScalarType, MV>::
1088makeInverseDiagonal (
const row_matrix_type& A,
const bool useDiagOffsets)
const
1091 using Teuchos::rcpFromRef;
1092 using Teuchos::rcp_dynamic_cast;
1095 if (!D_.is_null() &&
1096 D_->getMap()->isSameAs(*(A.getRowMap ()))) {
1098 *out_ <<
"Reusing pre-existing vector for diagonal extraction" << std::endl;
1099 D_rowMap = Teuchos::rcp_const_cast<V>(D_);
1101 D_rowMap = Teuchos::rcp(
new V (A.getRowMap (),
false));
1103 *out_ <<
"Allocated new vector for diagonal extraction" << std::endl;
1105 if (useDiagOffsets) {
1109 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
1110 typename MV::local_ordinal_type,
1111 typename MV::global_ordinal_type,
1112 typename MV::node_type> crs_matrix_type;
1113 RCP<const crs_matrix_type> A_crsMat =
1114 rcp_dynamic_cast<const crs_matrix_type> (rcpFromRef (A));
1115 if (! A_crsMat.is_null ()) {
1116 TEUCHOS_TEST_FOR_EXCEPTION(
1117 ! savedDiagOffsets_, std::logic_error,
1118 "Ifpack2::Details::Chebyshev::makeInverseDiagonal: "
1119 "It is not allowed to call this method with useDiagOffsets=true, "
1120 "if you have not previously saved offsets of diagonal entries. "
1121 "This situation should never arise if this class is used properly. "
1122 "Please report this bug to the Ifpack2 developers.");
1123 A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_);
1129 A.getLocalDiagCopy (*D_rowMap);
1131 RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
1137 bool foundNonpositiveValue =
false;
1139 auto D_lcl = D_rangeMap->getLocalViewHost (Tpetra::Access::ReadOnly);
1140 auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1142 typedef typename MV::impl_scalar_type IST;
1143 typedef typename MV::local_ordinal_type LO;
1144 typedef Kokkos::Details::ArithTraits<IST> ATS;
1145 typedef Kokkos::Details::ArithTraits<typename ATS::mag_type> STM;
1147 const LO lclNumRows =
static_cast<LO
> (D_rangeMap->getLocalLength ());
1148 for (LO i = 0; i < lclNumRows; ++i) {
1149 if (STS::real (D_lcl_1d(i)) <= STM::zero ()) {
1150 foundNonpositiveValue =
true;
1156 using Teuchos::outArg;
1157 using Teuchos::REDUCE_MIN;
1158 using Teuchos::reduceAll;
1160 const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1161 int gblSuccess = lclSuccess;
1162 if (! D_rangeMap->getMap ().is_null () && D_rangeMap->getMap ()->getComm ().is_null ()) {
1163 const Teuchos::Comm<int>& comm = * (D_rangeMap->getMap ()->getComm ());
1164 reduceAll<int, int> (comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
1166 if (gblSuccess == 1) {
1167 *out_ <<
"makeInverseDiagonal: The matrix's diagonal entries all have "
1168 "positive real part (this is good for Chebyshev)." << std::endl;
1171 *out_ <<
"makeInverseDiagonal: The matrix's diagonal has at least one "
1172 "entry with nonpositive real part, on at least one process in the "
1173 "matrix's communicator. This is bad for Chebyshev." << std::endl;
1179 reciprocal_threshold (*D_rangeMap, minDiagVal_);
1180 return Teuchos::rcp_const_cast<const V> (D_rangeMap);
1184template<
class ScalarType,
class MV>
1185Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
1186Chebyshev<ScalarType, MV>::
1187makeRangeMapVectorConst (
const Teuchos::RCP<const V>& D)
const
1191 typedef Tpetra::Export<
typename MV::local_ordinal_type,
1192 typename MV::global_ordinal_type,
1193 typename MV::node_type> export_type;
1197 TEUCHOS_TEST_FOR_EXCEPTION(
1198 A_.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1199 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() "
1200 "with a nonnull input matrix before calling this method. This is probably "
1201 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1202 TEUCHOS_TEST_FOR_EXCEPTION(
1203 D.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1204 "makeRangeMapVector: The input Vector D is null. This is probably "
1205 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1207 RCP<const map_type> sourceMap = D->getMap ();
1208 RCP<const map_type> rangeMap = A_->getRangeMap ();
1209 RCP<const map_type> rowMap = A_->getRowMap ();
1211 if (rangeMap->isSameAs (*sourceMap)) {
1217 RCP<const export_type> exporter;
1221 if (sourceMap->isSameAs (*rowMap)) {
1223 exporter = A_->getGraph ()->getExporter ();
1226 exporter = rcp (
new export_type (sourceMap, rangeMap));
1229 if (exporter.is_null ()) {
1233 RCP<V> D_out = rcp (
new V (*D, Teuchos::Copy));
1234 D_out->doExport (*D, *exporter, Tpetra::ADD);
1235 return Teuchos::rcp_const_cast<const V> (D_out);
1241template<
class ScalarType,
class MV>
1242Teuchos::RCP<typename Chebyshev<ScalarType, MV>::V>
1243Chebyshev<ScalarType, MV>::
1244makeRangeMapVector (
const Teuchos::RCP<V>& D)
const
1246 using Teuchos::rcp_const_cast;
1247 return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1251template<
class ScalarType,
class MV>
1253Chebyshev<ScalarType, MV>::
1254textbookApplyImpl (
const op_type& A,
1261 const V& D_inv)
const
1264 const ST myLambdaMin = lambdaMax / eigRatio;
1266 const ST zero = Teuchos::as<ST> (0);
1267 const ST one = Teuchos::as<ST> (1);
1268 const ST two = Teuchos::as<ST> (2);
1269 const ST d = (lambdaMax + myLambdaMin) / two;
1270 const ST c = (lambdaMax - myLambdaMin) / two;
1272 if (zeroStartingSolution_ && numIters > 0) {
1276 MV R (B.getMap (), B.getNumVectors (),
false);
1277 MV P (B.getMap (), B.getNumVectors (),
false);
1278 MV Z (B.getMap (), B.getNumVectors (),
false);
1280 for (
int i = 0; i < numIters; ++i) {
1281 computeResidual (R, B, A, X);
1282 solve (Z, D_inv, R);
1290 beta = alpha * (c/two) * (c/two);
1291 alpha = one / (d - beta);
1292 P.update (one, Z, beta);
1294 X.update (alpha, P, one);
1301template<
class ScalarType,
class MV>
1303Chebyshev<ScalarType, MV>::maxNormInf (
const MV& X) {
1304 Teuchos::Array<MT> norms (X.getNumVectors ());
1305 X.normInf (norms());
1306 return *std::max_element (norms.begin (), norms.end ());
1309template<
class ScalarType,
class MV>
1311Chebyshev<ScalarType, MV>::
1312ifpackApplyImpl (
const op_type& A,
1322#ifdef HAVE_IFPACK2_DEBUG
1323 const bool debug = debug_;
1325 const bool debug =
false;
1329 *out_ <<
" \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1330 *out_ <<
" \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1333 if (numIters <= 0) {
1336 const ST zero =
static_cast<ST
> (0.0);
1337 const ST one =
static_cast<ST
> (1.0);
1338 const ST two =
static_cast<ST
> (2.0);
1341 if (lambdaMin == one && lambdaMax == lambdaMin) {
1342 solve (X, D_inv, B);
1347 const ST alpha = lambdaMax / eigRatio;
1348 const ST beta = boostFactor_ * lambdaMax;
1349 const ST delta = two / (beta - alpha);
1350 const ST theta = (beta + alpha) / two;
1351 const ST s1 = theta * delta;
1354 *out_ <<
" alpha = " << alpha << endl
1355 <<
" beta = " << beta << endl
1356 <<
" delta = " << delta << endl
1357 <<
" theta = " << theta << endl
1358 <<
" s1 = " << s1 << endl;
1362 Teuchos::RCP<MV> W_ptr = makeTempMultiVector (B);
1366 *out_ <<
" Iteration " << 1 <<
":" << endl
1367 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1371 if (! zeroStartingSolution_) {
1374 if (ck_.is_null ()) {
1375 Teuchos::RCP<const op_type> A_op = A_;
1376 ck_ = Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1380 ck_->compute (W, one/theta,
const_cast<V&
> (D_inv),
1381 const_cast<MV&
> (B), X, zero);
1385 firstIterationWithZeroStartingSolution (W, one/theta, D_inv, B, X);
1389 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1390 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1393 if (numIters > 1 && ck_.is_null ()) {
1394 Teuchos::RCP<const op_type> A_op = A_;
1395 ck_ = Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1400 ST rhokp1, dtemp1, dtemp2;
1401 for (
int deg = 1; deg < numIters; ++deg) {
1403 *out_ <<
" Iteration " << deg+1 <<
":" << endl
1404 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl
1405 <<
" - \\|B\\|_{\\infty} = " << maxNormInf (B) << endl
1406 <<
" - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm ()
1407 << endl <<
" - rhok = " << rhok << endl;
1410 rhokp1 = one / (two * s1 - rhok);
1411 dtemp1 = rhokp1 * rhok;
1412 dtemp2 = two * rhokp1 * delta;
1416 *out_ <<
" - dtemp1 = " << dtemp1 << endl
1417 <<
" - dtemp2 = " << dtemp2 << endl;
1422 ck_->compute (W, dtemp2,
const_cast<V&
> (D_inv),
1423 const_cast<MV&
> (B), (X), dtemp1);
1426 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1427 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1433template<
class ScalarType,
class MV>
1435Chebyshev<ScalarType, MV>::
1436cgMethodWithInitGuess (
const op_type& A,
1442 using MagnitudeType =
typename STS::magnitudeType;
1444 *out_ <<
" cgMethodWithInitGuess:" << endl;
1447 const ST one = STS::one();
1448 ST beta, betaOld = one, pAp, pApOld = one, alpha, rHz, rHzOld, rHzOld2 = one, lambdaMax;
1450 Teuchos::ArrayRCP<MagnitudeType> diag, offdiag;
1451 Teuchos::RCP<V> p, z, Ap;
1452 diag.resize(numIters);
1453 offdiag.resize(numIters-1);
1455 p = rcp(
new V(A.getRangeMap ()));
1456 z = rcp(
new V(A.getRangeMap ()));
1457 Ap = rcp(
new V(A.getRangeMap ()));
1460 solve (*p, D_inv, r);
1463 for (
int iter = 0; iter < numIters; ++iter) {
1465 *out_ <<
" Iteration " << iter << endl;
1470 r.update (-alpha, *Ap, one);
1472 solve (*z, D_inv, r);
1474 beta = rHz / rHzOld;
1475 p->update (one, *z, beta);
1477 diag[iter] = STS::real((betaOld*betaOld * pApOld + pAp) / rHzOld);
1478 offdiag[iter-1] = -STS::real((betaOld * pApOld / (sqrt(rHzOld * rHzOld2))));
1480 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1481 *out_ <<
" offdiag["<< iter-1 <<
"] = " << offdiag[iter-1] << endl;
1484 diag[iter] = STS::real(pAp/rHzOld);
1486 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1494 lambdaMax = LapackHelper<ST>::tri_diag_spectral_radius(diag, offdiag);
1500template<
class ScalarType,
class MV>
1502Chebyshev<ScalarType, MV>::
1503cgMethod (
const op_type& A,
const V& D_inv,
const int numIters)
1507 *out_ <<
"cgMethod:" << endl;
1511 if (eigVector_.is_null()) {
1512 r = rcp(
new V(A.getDomainMap ()));
1513 if (eigKeepVectors_)
1518 PowerMethod::computeInitialGuessForPowerMethod (*r,
false);
1522 ST lambdaMax = cgMethodWithInitGuess (A, D_inv, numIters, *r);
1527template<
class ScalarType,
class MV>
1528Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::row_matrix_type>
1533template<
class ScalarType,
class MV>
1541template<
class ScalarType,
class MV>
1550 const size_t B_numVecs = B.getNumVectors ();
1551 if (W_.is_null () || W_->getNumVectors () != B_numVecs) {
1552 W_ = Teuchos::rcp (
new MV (B.getMap (), B_numVecs,
false));
1557template<
class ScalarType,
class MV>
1561 std::ostringstream oss;
1564 oss <<
"\"Ifpack2::Details::Chebyshev\":"
1566 <<
"degree: " << numIters_
1567 <<
", lambdaMax: " << lambdaMaxForApply_
1568 <<
", alpha: " << eigRatioForApply_
1569 <<
", lambdaMin: " << lambdaMinForApply_
1570 <<
", boost factor: " << boostFactor_;
1571 if (!userInvDiag_.is_null())
1572 oss <<
", diagonal: user-supplied";
1577template<
class ScalarType,
class MV>
1580describe (Teuchos::FancyOStream& out,
1581 const Teuchos::EVerbosityLevel verbLevel)
const
1583 using Teuchos::TypeNameTraits;
1586 const Teuchos::EVerbosityLevel vl =
1587 (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1588 if (vl == Teuchos::VERB_NONE) {
1598 Teuchos::OSTab tab0 (out);
1604 if (A_.is_null () || A_->getComm ().is_null ()) {
1608 myRank = A_->getComm ()->getRank ();
1613 out <<
"\"Ifpack2::Details::Chebyshev\":" << endl;
1615 Teuchos::OSTab tab1 (out);
1617 if (vl == Teuchos::VERB_LOW) {
1619 out <<
"degree: " << numIters_ << endl
1620 <<
"lambdaMax: " << lambdaMaxForApply_ << endl
1621 <<
"alpha: " << eigRatioForApply_ << endl
1622 <<
"lambdaMin: " << lambdaMinForApply_ << endl
1623 <<
"boost factor: " << boostFactor_ << endl;
1631 out <<
"Template parameters:" << endl;
1633 Teuchos::OSTab tab2 (out);
1634 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1635 <<
"MV: " << TypeNameTraits<MV>::name () << endl;
1641 out <<
"Computed parameters:" << endl;
1647 Teuchos::OSTab tab2 (out);
1653 if (D_.is_null ()) {
1655 out <<
"unset" << endl;
1658 else if (vl <= Teuchos::VERB_HIGH) {
1660 out <<
"set" << endl;
1669 D_->describe (out, vl);
1674 out <<
"W_: " << (W_.is_null () ?
"unset" :
"set") << endl
1675 <<
"computedLambdaMax_: " << computedLambdaMax_ << endl
1676 <<
"computedLambdaMin_: " << computedLambdaMin_ << endl
1677 <<
"lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1678 <<
"lambdaMinForApply_: " << lambdaMinForApply_ << endl
1679 <<
"eigRatioForApply_: " << eigRatioForApply_ << endl;
1684 out <<
"User parameters:" << endl;
1689 Teuchos::OSTab tab2 (out);
1690 out <<
"userInvDiag_: ";
1691 if (userInvDiag_.is_null ()) {
1692 out <<
"unset" << endl;
1694 else if (vl <= Teuchos::VERB_HIGH) {
1695 out <<
"set" << endl;
1701 userInvDiag_->describe (out, vl);
1704 out <<
"userLambdaMax_: " << userLambdaMax_ << endl
1705 <<
"userLambdaMin_: " << userLambdaMin_ << endl
1706 <<
"userEigRatio_: " << userEigRatio_ << endl
1707 <<
"numIters_: " << numIters_ << endl
1708 <<
"eigMaxIters_: " << eigMaxIters_ << endl
1709 <<
"eigRelTolerance_: " << eigRelTolerance_ << endl
1710 <<
"eigNormalizationFreq_: " << eigNormalizationFreq_ << endl
1711 <<
"zeroStartingSolution_: " << zeroStartingSolution_ << endl
1712 <<
"assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1713 <<
"textbookAlgorithm_: " << textbookAlgorithm_ << endl
1714 <<
"computeMaxResNorm_: " << computeMaxResNorm_ << endl;