44 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP 45 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP 55 #include "Kokkos_ArithTraits.hpp" 56 #include "Teuchos_FancyOStream.hpp" 57 #include "Teuchos_oblackholestream.hpp" 67 const char computeBeforeApplyReminder[] =
68 "This means one of the following:\n" 69 " - you have not yet called compute() on this instance, or \n" 70 " - you didn't call compute() after calling setParameters().\n\n" 71 "After creating an Ifpack2::Chebyshev instance,\n" 72 "you must _always_ call compute() at least once before calling apply().\n" 73 "After calling compute() once, you do not need to call it again,\n" 74 "unless the matrix has changed or you have changed parameters\n" 75 "(by calling setParameters()).";
81 template<
class XV,
class SizeType =
typename XV::
size_type>
82 struct V_ReciprocalThresholdSelfFunctor
84 typedef typename XV::execution_space execution_space;
85 typedef typename XV::non_const_value_type value_type;
86 typedef SizeType size_type;
87 typedef Kokkos::Details::ArithTraits<value_type> KAT;
88 typedef typename KAT::mag_type mag_type;
91 const value_type minVal_;
92 const mag_type minValMag_;
94 V_ReciprocalThresholdSelfFunctor (
const XV& X,
95 const value_type& min_val) :
98 minValMag_ (KAT::abs (min_val))
101 KOKKOS_INLINE_FUNCTION
102 void operator() (
const size_type& i)
const 104 const mag_type X_i_abs = KAT::abs (X_[i]);
106 if (X_i_abs < minValMag_) {
110 X_[i] = KAT::one () / X_[i];
115 template<
class XV,
class SizeType =
typename XV::
size_type>
116 struct LocalReciprocalThreshold {
118 compute (
const XV& X,
119 const typename XV::non_const_value_type& minVal)
121 typedef typename XV::execution_space execution_space;
122 Kokkos::RangePolicy<execution_space, SizeType> policy (0, X.dimension_0 ());
123 V_ReciprocalThresholdSelfFunctor<XV, SizeType> op (X, minVal);
124 Kokkos::parallel_for (policy, op);
128 template <
class TpetraVectorType,
129 const bool classic = TpetraVectorType::node_type::classic>
130 struct GlobalReciprocalThreshold {};
132 template <
class TpetraVectorType>
133 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
135 compute (TpetraVectorType& V,
136 const typename TpetraVectorType::scalar_type& min_val)
138 typedef typename TpetraVectorType::scalar_type scalar_type;
139 typedef typename TpetraVectorType::mag_type mag_type;
140 typedef Kokkos::Details::ArithTraits<scalar_type> STS;
142 const scalar_type ONE = STS::one ();
143 const mag_type min_val_abs = STS::abs (min_val);
146 const size_t lclNumRows = V.getLocalLength ();
148 for (
size_t i = 0; i < lclNumRows; ++i) {
149 const scalar_type V_0i = V_0[i];
150 if (STS::abs (V_0i) < min_val_abs) {
159 template <
class TpetraVectorType>
160 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
162 compute (TpetraVectorType& X,
163 const typename TpetraVectorType::scalar_type& minVal)
165 typedef typename TpetraVectorType::impl_scalar_type value_type;
166 typedef typename TpetraVectorType::device_type::memory_space memory_space;
168 X.template sync<memory_space> ();
169 X.template modify<memory_space> ();
171 const value_type minValS =
static_cast<value_type
> (minVal);
172 auto X_0 = Kokkos::subview (X.template getLocalView<memory_space> (),
174 LocalReciprocalThreshold<decltype (X_0) >::compute (X_0, minValS);
179 template <
typename S,
typename L,
typename G,
typename N>
181 reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V,
const S& minVal)
183 GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
191 template<
class OneDViewType,
192 class LocalOrdinal =
typename OneDViewType::size_type>
193 class PositivizeVector {
194 static_assert (Kokkos::Impl::is_view<OneDViewType>::value,
195 "OneDViewType must be a 1-D Kokkos::View.");
196 static_assert (static_cast<int> (OneDViewType::rank) == 1,
197 "This functor only works with a 1-D View.");
198 static_assert (std::is_integral<LocalOrdinal>::value,
199 "The second template parameter, LocalOrdinal, " 200 "must be an integer type.");
202 PositivizeVector (
const OneDViewType& x) : x_ (x) {}
204 KOKKOS_INLINE_FUNCTION
void 205 operator () (
const LocalOrdinal& i)
const 207 typedef typename OneDViewType::non_const_value_type IST;
208 typedef Kokkos::Details::ArithTraits<IST> STS;
209 typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
211 if (STS::real (x_(i)) < STM::zero ()) {
222 template<
class ScalarType,
class MV>
223 void Chebyshev<ScalarType, MV>::checkInputMatrix ()
const 226 ! A_.is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
227 std::invalid_argument,
228 "Ifpack2::Chebyshev: The input matrix A must be square. " 229 "A has " << A_->getGlobalNumRows () <<
" rows and " 230 << A_->getGlobalNumCols () <<
" columns.");
234 if (debug_ && ! A_.is_null ()) {
242 ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
243 "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be " 244 "the same (in the sense of isSameAs())." << std::endl <<
"We only check " 245 "for this in debug mode.");
250 template<
class ScalarType,
class MV>
252 Chebyshev<ScalarType, MV>::
253 checkConstructorInput ()
const 258 if (STS::isComplex) {
260 (
true, std::logic_error,
"Ifpack2::Chebyshev: This class' implementation " 261 "of Chebyshev iteration only works for real-valued, symmetric positive " 262 "definite matrices. However, you instantiated this class for ScalarType" 264 "complex-valued type. While this may be algorithmically correct if all " 265 "of the complex numbers in the matrix have zero imaginary part, we " 266 "forbid using complex ScalarType altogether in order to remind you of " 267 "the limitations of our implementation (and of the algorithm itself).");
274 template<
class ScalarType,
class MV>
278 savedDiagOffsets_ (false),
279 computedLambdaMax_ (
STS::nan ()),
280 computedLambdaMin_ (
STS::nan ()),
281 lambdaMaxForApply_ (
STS::nan ()),
282 lambdaMinForApply_ (
STS::nan ()),
283 eigRatioForApply_ (
STS::nan ()),
284 userLambdaMax_ (
STS::nan ()),
285 userLambdaMin_ (
STS::nan ()),
287 minDiagVal_ (
STS::eps ()),
290 zeroStartingSolution_ (true),
291 assumeMatrixUnchanged_ (false),
292 textbookAlgorithm_ (false),
293 computeMaxResNorm_ (false),
296 checkConstructorInput ();
299 template<
class ScalarType,
class MV>
303 savedDiagOffsets_ (false),
304 computedLambdaMax_ (
STS::nan ()),
305 computedLambdaMin_ (
STS::nan ()),
306 lambdaMaxForApply_ (
STS::nan ()),
307 boostFactor_ (static_cast<
MT> (1.1)),
308 lambdaMinForApply_ (
STS::nan ()),
309 eigRatioForApply_ (
STS::nan ()),
310 userLambdaMax_ (
STS::nan ()),
311 userLambdaMin_ (
STS::nan ()),
313 minDiagVal_ (
STS::eps ()),
316 zeroStartingSolution_ (true),
317 assumeMatrixUnchanged_ (false),
318 textbookAlgorithm_ (false),
319 computeMaxResNorm_ (false),
322 checkConstructorInput ();
323 setParameters (params);
326 template<
class ScalarType,
class MV>
333 using Teuchos::rcp_const_cast;
345 const ST defaultLambdaMax = STS::nan ();
346 const ST defaultLambdaMin = STS::nan ();
353 const ST defaultEigRatio = Teuchos::as<ST> (30);
354 const MT defaultBoostFactor =
static_cast<MT> (1.1);
355 const ST defaultMinDiagVal = STS::eps ();
356 const int defaultNumIters = 1;
357 const int defaultEigMaxIters = 10;
358 const bool defaultZeroStartingSolution =
true;
359 const bool defaultAssumeMatrixUnchanged =
false;
360 const bool defaultTextbookAlgorithm =
false;
361 const bool defaultComputeMaxResNorm =
false;
362 const bool defaultDebug =
false;
368 RCP<const V> userInvDiagCopy;
369 ST lambdaMax = defaultLambdaMax;
370 ST lambdaMin = defaultLambdaMin;
371 ST eigRatio = defaultEigRatio;
372 MT boostFactor = defaultBoostFactor;
373 ST minDiagVal = defaultMinDiagVal;
374 int numIters = defaultNumIters;
375 int eigMaxIters = defaultEigMaxIters;
376 bool zeroStartingSolution = defaultZeroStartingSolution;
377 bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
378 bool textbookAlgorithm = defaultTextbookAlgorithm;
379 bool computeMaxResNorm = defaultComputeMaxResNorm;
380 bool debug = defaultDebug;
389 debug = plist.
get<
bool> (
"debug");
394 int debugInt = plist.
get<
bool> (
"debug");
395 debug = debugInt != 0;
409 if (plist.
isParameter (
"chebyshev: operator inv diagonal")) {
411 RCP<const V> userInvDiag;
414 const V* rawUserInvDiag =
415 plist.
get<
const V*> (
"chebyshev: operator inv diagonal");
417 userInvDiag =
rcp (rawUserInvDiag,
false);
420 if (userInvDiag.is_null ()) {
422 V* rawUserInvDiag = plist.
get<
V*> (
"chebyshev: operator inv diagonal");
424 userInvDiag =
rcp (const_cast<const V*> (rawUserInvDiag),
false);
428 if (userInvDiag.is_null ()) {
431 plist.
get<RCP<const V> > (
"chebyshev: operator inv diagonal");
435 if (userInvDiag.is_null ()) {
437 RCP<V> userInvDiagNonConst =
438 plist.
get<RCP<V> > (
"chebyshev: operator inv diagonal");
439 userInvDiag = rcp_const_cast<
const V> (userInvDiagNonConst);
443 if (userInvDiag.is_null ()) {
452 const V& userInvDiagRef =
453 plist.
get<
const V> (
"chebyshev: operator inv diagonal");
454 userInvDiagCopy =
rcp (
new V (userInvDiagRef, Teuchos::Copy));
456 userInvDiag = userInvDiagCopy;
461 true, std::runtime_error,
462 "Ifpack2::Chebyshev::setParameters: \"chebyshev: operator inv diagonal\" " 463 "plist.get<const V> does not compile around return held == other_held " 464 "in Teuchos::any in Visual Studio. Can't fix it now, so throwing " 465 "in case someone builds there.");
468 if (userInvDiag.is_null ()) {
477 V& userInvDiagNonConstRef =
478 plist.
get<
V> (
"chebyshev: operator inv diagonal");
479 const V& userInvDiagRef =
const_cast<const V&
> (userInvDiagNonConstRef);
480 userInvDiagCopy =
rcp (
new V (userInvDiagRef, Teuchos::Copy));
482 userInvDiag = userInvDiagCopy;
487 true, std::runtime_error,
488 "Ifpack2::Chebyshev::setParameters: \"chebyshev: operator inv diagonal\" " 489 "plist.get<V> does not compile around return held == other_held " 490 "in Teuchos::any in Visual Studio. Can't fix it now, so throwing " 491 "in case someone builds there.");
502 if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
503 userInvDiagCopy =
rcp (
new V (*userInvDiag, Teuchos::Copy));
515 if (plist.
isParameter (
"chebyshev: max eigenvalue")) {
516 if (plist.
isType<
double>(
"chebyshev: max eigenvalue"))
517 lambdaMax = plist.
get<
double> (
"chebyshev: max eigenvalue");
519 lambdaMax = plist.
get<
ST> (
"chebyshev: max eigenvalue");
521 STS::isnaninf (lambdaMax), std::invalid_argument,
522 "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" " 523 "parameter is NaN or Inf. This parameter is optional, but if you " 524 "choose to supply it, it must have a finite value.");
526 if (plist.
isParameter (
"chebyshev: min eigenvalue")) {
527 if (plist.
isType<
double>(
"chebyshev: min eigenvalue"))
528 lambdaMin = plist.
get<
double> (
"chebyshev: min eigenvalue");
530 lambdaMin = plist.
get<
ST> (
"chebyshev: min eigenvalue");
532 STS::isnaninf (lambdaMin), std::invalid_argument,
533 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" " 534 "parameter is NaN or Inf. This parameter is optional, but if you " 535 "choose to supply it, it must have a finite value.");
539 if (plist.
isParameter (
"smoother: Chebyshev alpha")) {
540 if (plist.
isType<
double>(
"smoother: Chebyshev alpha"))
541 eigRatio = plist.
get<
double> (
"smoother: Chebyshev alpha");
543 eigRatio = plist.
get<
ST> (
"smoother: Chebyshev alpha");
546 eigRatio = plist.
get (
"chebyshev: ratio eigenvalue", eigRatio);
548 STS::isnaninf (eigRatio), std::invalid_argument,
549 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" " 550 "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. " 551 "This parameter is optional, but if you choose to supply it, it must have " 559 STS::real (eigRatio) < STS::real (STS::one ()),
560 std::invalid_argument,
561 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\"" 562 "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, " 563 "but you supplied the value " << eigRatio <<
".");
568 const char paramName[] =
"chebyshev: boost factor";
572 boostFactor = plist.
get<
MT> (paramName);
574 else if (! std::is_same<double, MT>::value &&
575 plist.
isType<
double> (paramName)) {
576 const double dblBF = plist.
get<
double> (paramName);
577 boostFactor =
static_cast<MT> (dblBF);
581 (
true, std::invalid_argument,
582 "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\"" 583 "parameter must have type magnitude_type (MT) or double.");
593 plist.
set (paramName, defaultBoostFactor);
597 "Ifpack2::Chebyshev::setParameters: \"" << paramName <<
"\" parameter " 598 "must be >= 1, but you supplied the value " << boostFactor <<
".");
602 minDiagVal = plist.
get (
"chebyshev: min diagonal value", minDiagVal);
604 STS::isnaninf (minDiagVal), std::invalid_argument,
605 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" " 606 "parameter is NaN or Inf. This parameter is optional, but if you choose " 607 "to supply it, it must have a finite value.");
611 numIters = plist.
get<
int> (
"smoother: sweeps");
614 numIters = plist.
get<
int> (
"relaxation: sweeps");
616 numIters = plist.
get (
"chebyshev: degree", numIters);
618 numIters < 0, std::invalid_argument,
619 "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also " 620 "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a " 621 "nonnegative integer. You gave a value of " << numIters <<
".");
624 if (plist.
isParameter (
"eigen-analysis: iterations")) {
625 eigMaxIters = plist.
get<
int> (
"eigen-analysis: iterations");
627 eigMaxIters = plist.
get (
"chebyshev: eigenvalue max iterations", eigMaxIters);
629 eigMaxIters < 0, std::invalid_argument,
630 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations" 631 "\" parameter (also called \"eigen-analysis: iterations\") must be a " 632 "nonnegative integer. You gave a value of " << eigMaxIters <<
".");
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");
653 ! plist.
get<
bool> (
"chebyshev: use block mode"),
654 std::invalid_argument,
655 "Ifpack2::Chebyshev requires that if you supply the Ifpack parameter " 656 "\"chebyshev: use block mode\", it must be set to false. Ifpack2's " 657 "Chebyshev does not implement Ifpack's block mode.");
659 plist.
isParameter (
"chebyshev: solve normal equations") &&
660 ! plist.
get<
bool> (
"chebyshev: solve normal equations"),
661 std::invalid_argument,
662 "Ifpack2::Chebyshev does not and will never implement the Ifpack " 663 "parameter \"chebyshev: solve normal equations\". If you want to solve " 664 "the normal equations, construct a Tpetra::Operator that implements " 665 "A^* A, and use Chebyshev to solve A^* A x = A^* b.");
673 std::string eigenAnalysisType (
"power-method");
675 eigenAnalysisType = plist.
get<std::string> (
"eigen-analysis: type");
677 eigenAnalysisType !=
"power-method" &&
678 eigenAnalysisType !=
"power method",
679 std::invalid_argument,
680 "Ifpack2::Chebyshev: This class supports the ML parameter \"eigen-" 681 "analysis: type\" for backwards compatibility. However, Ifpack2 " 682 "currently only supports the \"power-method\" option for this " 683 "parameter. This imitates Ifpack, which only implements the power " 684 "method for eigenanalysis.");
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 zeroStartingSolution_ = zeroStartingSolution;
697 assumeMatrixUnchanged_ = assumeMatrixUnchanged;
698 textbookAlgorithm_ = textbookAlgorithm;
699 computeMaxResNorm_ = computeMaxResNorm;
705 if (A_.is_null () || A_->getComm ().is_null ()) {
711 myRank = A_->getComm ()->getRank ();
715 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
719 out_ = Teuchos::getFancyOStream (blackHole);
724 out_ = Teuchos::null;
729 template<
class ScalarType,
class MV>
734 diagOffsets_ = offsets_type ();
735 savedDiagOffsets_ =
false;
738 computedLambdaMax_ = STS::nan ();
739 computedLambdaMin_ = STS::nan ();
743 template<
class ScalarType,
class MV>
748 if (! assumeMatrixUnchanged_) {
765 myRank = A->getComm ()->getRank ();
769 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
773 out_ = Teuchos::getFancyOStream (blackHole);
778 out_ = Teuchos::null;
784 template<
class ScalarType,
class MV>
793 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
794 typename MV::local_ordinal_type,
795 typename MV::global_ordinal_type,
796 typename MV::node_type> crs_matrix_type;
799 A_.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::compute: The input " 800 "matrix A is null. Please call setMatrix() with a nonnull input matrix " 801 "before calling this method.");
816 if (userInvDiag_.is_null ()) {
818 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A_);
821 if (! A_crsMat.
is_null () && A_crsMat->isStaticGraph ()) {
823 const size_t lclNumRows = A_crsMat->getNodeNumRows ();
824 if (diagOffsets_.dimension_0 () < lclNumRows) {
825 diagOffsets_ = offsets_type ();
826 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
828 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
829 savedDiagOffsets_ =
true;
830 D_ = makeInverseDiagonal (*A_,
true);
833 D_ = makeInverseDiagonal (*A_);
836 else if (! assumeMatrixUnchanged_) {
837 if (! A_crsMat.
is_null () && A_crsMat->isStaticGraph ()) {
840 if (! savedDiagOffsets_) {
841 const size_t lclNumRows = A_crsMat->getNodeNumRows ();
842 if (diagOffsets_.dimension_0 () < lclNumRows) {
843 diagOffsets_ = offsets_type ();
844 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
846 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
847 savedDiagOffsets_ =
true;
850 D_ = makeInverseDiagonal (*A_,
true);
853 D_ = makeInverseDiagonal (*A_);
858 D_ = makeRangeMapVectorConst (userInvDiag_);
862 const bool computedEigenvalueEstimates =
863 STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
872 if (! assumeMatrixUnchanged_ ||
873 (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
874 const ST computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
876 STS::isnaninf (computedLambdaMax),
878 "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue " 879 "of D^{-1} A failed, by producing Inf or NaN. This probably means that " 880 "the matrix contains Inf or NaN values, or that it is badly scaled.");
882 STS::isnaninf (userEigRatio_),
884 "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN." 885 << endl <<
"This should be impossible." << endl <<
886 "Please report this bug to the Ifpack2 developers.");
892 const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
895 computedLambdaMax_ = computedLambdaMax;
896 computedLambdaMin_ = computedLambdaMin;
899 STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
901 "Ifpack2::Chebyshev::compute: " << endl <<
902 "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN." 903 << endl <<
"This should be impossible." << endl <<
904 "Please report this bug to the Ifpack2 developers.");
912 if (STS::isnaninf (userLambdaMax_)) {
913 lambdaMaxForApply_ = computedLambdaMax_;
915 lambdaMaxForApply_ = userLambdaMax_;
928 lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
929 eigRatioForApply_ = userEigRatio_;
931 if (! textbookAlgorithm_) {
934 const ST one = Teuchos::as<ST> (1);
937 if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
938 lambdaMinForApply_ = one;
939 lambdaMaxForApply_ = lambdaMinForApply_;
940 eigRatioForApply_ = one;
946 template<
class ScalarType,
class MV>
950 return lambdaMaxForApply_;
954 template<
class ScalarType,
class MV>
958 const char prefix[] =
"Ifpack2::Chebyshev::apply: ";
961 *out_ <<
"apply: " << std::endl;
964 (A_.is_null (), std::runtime_error, prefix <<
"The input matrix A is null. " 965 " Please call setMatrix() with a nonnull input matrix, and then call " 966 "compute(), before calling this method.");
968 (STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
969 prefix <<
"There is no estimate for the max eigenvalue." 970 << std::endl << computeBeforeApplyReminder);
972 (STS::isnaninf (lambdaMinForApply_), std::runtime_error,
973 prefix <<
"There is no estimate for the min eigenvalue." 974 << std::endl << computeBeforeApplyReminder);
976 (STS::isnaninf (eigRatioForApply_), std::runtime_error,
977 prefix <<
"There is no estimate for the ratio of the max " 978 "eigenvalue to the min eigenvalue." 979 << std::endl << computeBeforeApplyReminder);
981 (D_.is_null (), std::runtime_error, prefix <<
"The vector of inverse " 982 "diagonal entries of the matrix has not yet been computed." 983 << std::endl << computeBeforeApplyReminder);
985 if (textbookAlgorithm_) {
986 textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
987 lambdaMinForApply_, eigRatioForApply_, *D_);
990 ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
991 lambdaMinForApply_, eigRatioForApply_, *D_);
994 if (computeMaxResNorm_ && B.getNumVectors () > 0) {
995 MV R (B.getMap (), B.getNumVectors ());
996 computeResidual (R, B, *A_, X);
999 return *std::max_element (norms.begin (), norms.end ());
1006 template<
class ScalarType,
class MV>
1010 this->describe (* (Teuchos::getFancyOStream (Teuchos::rcpFromRef (out))),
1014 template<
class ScalarType,
class MV>
1020 Tpetra::deep_copy(R, B);
1021 A.apply (X, R, mode, -STS::one(), STS::one());
1024 template<
class ScalarType,
class MV>
1027 solve (MV& Z,
const V& D_inv,
const MV& R) {
1028 Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
1031 template<
class ScalarType,
class MV>
1033 Chebyshev<ScalarType, MV>::
1034 solve (MV& Z,
const ST alpha,
const V& D_inv,
const MV& R) {
1035 Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
1038 template<
class ScalarType,
class MV>
1040 Chebyshev<ScalarType, MV>::
1041 makeInverseDiagonal (
const row_matrix_type& A,
const bool useDiagOffsets)
const 1044 using Teuchos::rcpFromRef;
1045 using Teuchos::rcp_dynamic_cast;
1047 RCP<V> D_rowMap (
new V (A.getGraph ()->getRowMap ()));
1048 if (useDiagOffsets) {
1052 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
1053 typename MV::local_ordinal_type,
1054 typename MV::global_ordinal_type,
1055 typename MV::node_type> crs_matrix_type;
1056 RCP<const crs_matrix_type> A_crsMat =
1057 rcp_dynamic_cast<
const crs_matrix_type> (rcpFromRef (A));
1058 if (! A_crsMat.is_null ()) {
1060 ! savedDiagOffsets_, std::logic_error,
1061 "Ifpack2::Details::Chebyshev::makeInverseDiagonal: " 1062 "It is not allowed to call this method with useDiagOffsets=true, " 1063 "if you have not previously saved offsets of diagonal entries. " 1064 "This situation should never arise if this class is used properly. " 1065 "Please report this bug to the Ifpack2 developers.");
1066 A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_);
1072 A.getLocalDiagCopy (*D_rowMap);
1074 RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
1080 D_rangeMap->template sync<Kokkos::HostSpace> ();
1081 auto D_lcl = D_rangeMap->template getLocalView<Kokkos::HostSpace> ();
1082 auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1084 typedef typename MV::impl_scalar_type IST;
1085 typedef typename MV::local_ordinal_type LO;
1086 typedef Kokkos::Details::ArithTraits<IST> STS;
1087 typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
1089 const LO lclNumRows =
static_cast<LO
> (D_rangeMap->getLocalLength ());
1090 bool foundNonpositiveValue =
false;
1091 for (LO i = 0; i < lclNumRows; ++i) {
1092 if (STS::real (D_lcl_1d(i)) <= STM::zero ()) {
1093 foundNonpositiveValue =
true;
1098 using Teuchos::outArg;
1099 using Teuchos::REDUCE_MIN;
1102 const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1103 int gblSuccess = lclSuccess;
1104 if (! D_rangeMap->getMap ().is_null () && D_rangeMap->getMap ()->getComm ().is_null ()) {
1106 reduceAll<int, int> (comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
1108 if (gblSuccess == 1) {
1109 *out_ <<
"makeInverseDiagonal: The matrix's diagonal entries all have " 1110 "positive real part (this is good for Chebyshev)." << std::endl;
1113 *out_ <<
"makeInverseDiagonal: The matrix's diagonal has at least one " 1114 "entry with nonpositive real part, on at least one process in the " 1115 "matrix's communicator. This is bad for Chebyshev." << std::endl;
1121 reciprocal_threshold (*D_rangeMap, minDiagVal_);
1122 return Teuchos::rcp_const_cast<
const V> (D_rangeMap);
1126 template<
class ScalarType,
class MV>
1128 Chebyshev<ScalarType, MV>::
1133 typedef Tpetra::Export<
typename MV::local_ordinal_type,
1134 typename MV::global_ordinal_type,
1135 typename MV::node_type> export_type;
1140 A_.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::" 1141 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() " 1142 "with a nonnull input matrix before calling this method. This is probably " 1143 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1145 D.
is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::" 1146 "makeRangeMapVector: The input Vector D is null. This is probably " 1147 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1149 RCP<const map_type> sourceMap = D->getMap ();
1150 RCP<const map_type> rangeMap = A_->getRangeMap ();
1151 RCP<const map_type> rowMap = A_->getRowMap ();
1153 if (rangeMap->isSameAs (*sourceMap)) {
1159 RCP<const export_type> exporter;
1163 if (sourceMap->isSameAs (*rowMap)) {
1165 exporter = A_->getGraph ()->getExporter ();
1168 exporter =
rcp (
new export_type (sourceMap, rangeMap));
1171 if (exporter.is_null ()) {
1175 RCP<V> D_out =
rcp (
new V (*D, Teuchos::Copy));
1176 D_out->doExport (*D, *exporter, Tpetra::ADD);
1177 return Teuchos::rcp_const_cast<
const V> (D_out);
1183 template<
class ScalarType,
class MV>
1185 Chebyshev<ScalarType, MV>::
1188 using Teuchos::rcp_const_cast;
1189 return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1193 template<
class ScalarType,
class MV>
1195 Chebyshev<ScalarType, MV>::
1196 textbookApplyImpl (
const op_type& A,
1203 const V& D_inv)
const 1206 const ST myLambdaMin = lambdaMax / eigRatio;
1208 const ST
zero = Teuchos::as<ST> (0);
1209 const ST one = Teuchos::as<ST> (1);
1210 const ST two = Teuchos::as<ST> (2);
1211 const ST d = (lambdaMax + myLambdaMin) / two;
1212 const ST c = (lambdaMax - myLambdaMin) / two;
1214 if (zeroStartingSolution_ && numIters > 0) {
1218 MV R (B.getMap (), B.getNumVectors (),
false);
1219 MV P (B.getMap (), B.getNumVectors (),
false);
1220 MV Z (B.getMap (), B.getNumVectors (),
false);
1222 for (
int i = 0; i < numIters; ++i) {
1223 computeResidual (R, B, A, X);
1224 solve (Z, D_inv, R);
1232 beta = alpha * (c/two) * (c/two);
1233 alpha = one / (d - beta);
1234 P.update (one, Z, beta);
1236 X.update (alpha, P, one);
1243 template<
class ScalarType,
class MV>
1244 typename Chebyshev<ScalarType, MV>::MT
1245 Chebyshev<ScalarType, MV>::maxNormInf (
const MV& X) {
1247 X.normInf (norms());
1248 return *std::max_element (norms.begin (), norms.end ());
1251 template<
class ScalarType,
class MV>
1253 Chebyshev<ScalarType, MV>::
1254 ifpackApplyImpl (
const op_type& A,
1264 const bool debug = debug_;
1267 *out_ <<
" \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1268 *out_ <<
" \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1271 if (numIters <= 0) {
1274 const ST one =
static_cast<ST
> (1.0);
1275 const ST two =
static_cast<ST
> (2.0);
1278 if (lambdaMin == one && lambdaMax == lambdaMin) {
1279 solve (X, D_inv, B);
1284 const ST alpha = lambdaMax / eigRatio;
1285 const ST beta = boostFactor_ * lambdaMax;
1286 const ST delta = two / (beta - alpha);
1287 const ST theta = (beta + alpha) / two;
1288 const ST s1 = theta * delta;
1291 *out_ <<
" alpha = " << alpha << endl
1292 <<
" beta = " << beta << endl
1293 <<
" delta = " << delta << endl
1294 <<
" theta = " << theta << endl
1295 <<
" s1 = " << s1 << endl;
1300 makeTempMultiVectors (V_ptr, W_ptr, B);
1310 *out_ <<
" Iteration " << 1 <<
":" << endl
1311 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1315 if (! zeroStartingSolution_) {
1316 computeResidual (V1, B, A, X);
1319 *out_ <<
" - \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1322 solve (W, one/theta, D_inv, V1);
1325 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1328 X.update (one, W, one);
1331 solve (W, one/theta, D_inv, B);
1333 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1335 Tpetra::deep_copy (X, W);
1338 *out_ <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1343 ST rhokp1, dtemp1, dtemp2;
1344 for (
int deg = 1; deg < numIters; ++deg) {
1347 *out_ <<
" Iteration " << deg+1 <<
":" << endl
1348 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl
1349 <<
" - \\|B\\|_{\\infty} = " << maxNormInf (B) << endl
1350 <<
" - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm () << endl
1351 <<
" - rhok = " << rhok << endl;
1352 V1.putScalar (STS::zero ());
1355 computeResidual (V1, B, A, X);
1358 *out_ <<
" - \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1361 rhokp1 = one / (two * s1 - rhok);
1362 dtemp1 = rhokp1 * rhok;
1363 dtemp2 = two * rhokp1 * delta;
1367 *out_ <<
" - dtemp1 = " << dtemp1 << endl
1368 <<
" - dtemp2 = " << dtemp2 << endl;
1372 W.elementWiseMultiply (dtemp2, D_inv, V1, one);
1373 X.update (one, W, one);
1376 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1377 *out_ <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1382 template<
class ScalarType,
class MV>
1383 typename Chebyshev<ScalarType, MV>::ST
1384 Chebyshev<ScalarType, MV>::
1385 powerMethodWithInitGuess (
const op_type& A,
1392 *out_ <<
" powerMethodWithInitGuess:" << endl;
1395 const ST
zero =
static_cast<ST
> (0.0);
1396 const ST one =
static_cast<ST
> (1.0);
1397 ST lambdaMax =
zero;
1398 ST RQ_top, RQ_bottom, norm;
1400 V y (A.getRangeMap ());
1403 (norm == zero, std::runtime_error,
1404 "Ifpack2::Chebyshev::powerMethodWithInitGuess: " 1405 "The initial guess has zero norm. This could be either because Tpetra::" 1406 "Vector::randomize() filled the vector with zeros (if that was used to " 1407 "compute the initial guess), or because the norm2 method has a bug. The " 1408 "first is not impossible, but unlikely.");
1411 *out_ <<
" Original norm1(x): " << x.norm1 () <<
", norm2(x): " << norm << endl;
1414 x.scale (one / norm);
1417 *out_ <<
" norm1(x.scale(one/norm)): " << x.norm1 () << endl;
1420 for (
int iter = 0; iter < numIters; ++iter) {
1422 *out_ <<
" Iteration " << iter << endl;
1425 solve (y, D_inv, y);
1427 RQ_bottom = x.dot (x);
1429 *out_ <<
" RQ_top: " << RQ_top <<
", RQ_bottom: " << RQ_bottom << endl;
1431 lambdaMax = RQ_top / RQ_bottom;
1435 *out_ <<
" norm is zero; returning zero" << endl;
1439 x.update (one / norm, y, zero);
1442 *out_ <<
" lambdaMax: " << lambdaMax << endl;
1447 template<
class ScalarType,
class MV>
1449 Chebyshev<ScalarType, MV>::
1450 computeInitialGuessForPowerMethod (V& x,
const bool nonnegativeRealParts)
const 1452 typedef typename MV::device_type::execution_space dev_execution_space;
1453 typedef typename MV::device_type::memory_space dev_memory_space;
1454 typedef typename MV::local_ordinal_type LO;
1458 if (nonnegativeRealParts) {
1459 x.template modify<dev_memory_space> ();
1460 auto x_lcl = x.template getLocalView<dev_memory_space> ();
1461 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
1463 const LO lclNumRows =
static_cast<LO
> (x.getLocalLength ());
1464 Kokkos::RangePolicy<dev_execution_space, LO> range (0, lclNumRows);
1465 PositivizeVector<decltype (x_lcl_1d), LO> functor (x_lcl_1d);
1466 Kokkos::parallel_for (range, functor);
1470 template<
class ScalarType,
class MV>
1471 typename Chebyshev<ScalarType, MV>::ST
1472 Chebyshev<ScalarType, MV>::
1473 powerMethod (
const op_type& A,
const V& D_inv,
const int numIters)
1477 *out_ <<
"powerMethod:" << endl;
1480 const ST
zero =
static_cast<ST
> (0.0);
1481 V x (A.getDomainMap ());
1485 computeInitialGuessForPowerMethod (x,
false);
1487 ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1496 if (STS::real (lambdaMax) < STS::real (zero)) {
1498 *out_ <<
"real(lambdaMax) = " << STS::real (lambdaMax) <<
" < 0; try " 1499 "again with a different random initial guess" << endl;
1508 computeInitialGuessForPowerMethod (x,
true);
1509 lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1514 template<
class ScalarType,
class MV>
1520 template<
class ScalarType,
class MV>
1528 template<
class ScalarType,
class MV>
1535 if (V_.is_null ()) {
1536 V_ = Teuchos::rcp (
new MV (X.getMap (), X.getNumVectors (),
false));
1539 if (W_.is_null ()) {
1540 W_ = Teuchos::rcp (
new MV (X.getMap (), X.getNumVectors (),
true));
1546 template<
class ScalarType,
class MV>
1548 Chebyshev<ScalarType, MV>::
1549 description ()
const {
1550 std::ostringstream oss;
1553 oss <<
"\"Ifpack2::Details::Chebyshev\":" 1555 <<
"degree: " << numIters_
1556 <<
", lambdaMax: " << lambdaMaxForApply_
1557 <<
", alpha: " << eigRatioForApply_
1558 <<
", lambdaMin: " << lambdaMinForApply_
1559 <<
", boost factor: " << boostFactor_
1564 template<
class ScalarType,
class MV>
1591 if (A_.is_null () || A_->getComm ().is_null ()) {
1595 myRank = A_->getComm ()->getRank ();
1600 out <<
"\"Ifpack2::Details::Chebyshev\":" << endl;
1606 out <<
"degree: " << numIters_ << endl
1607 <<
"lambdaMax: " << lambdaMaxForApply_ << endl
1608 <<
"alpha: " << eigRatioForApply_ << endl
1609 <<
"lambdaMin: " << lambdaMinForApply_ << endl
1610 <<
"boost factor: " << boostFactor_ << endl;
1618 out <<
"Template parameters:" << endl;
1621 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1622 <<
"MV: " << TypeNameTraits<MV>::name () << endl;
1628 out <<
"Computed parameters:" << endl;
1640 if (D_.is_null ()) {
1642 out <<
"unset" << endl;
1647 out <<
"set" << endl;
1656 D_->describe (out, vl);
1661 out <<
"V_: " << (V_.is_null () ?
"unset" :
"set") << endl
1662 <<
"W_: " << (W_.is_null () ?
"unset" :
"set") << endl
1663 <<
"computedLambdaMax_: " << computedLambdaMax_ << endl
1664 <<
"computedLambdaMin_: " << computedLambdaMin_ << endl
1665 <<
"lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1666 <<
"lambdaMinForApply_: " << lambdaMinForApply_ << endl
1667 <<
"eigRatioForApply_: " << eigRatioForApply_ << endl;
1672 out <<
"User parameters:" << endl;
1678 out <<
"userInvDiag_: ";
1679 if (userInvDiag_.is_null ()) {
1680 out <<
"unset" << endl;
1683 out <<
"set" << endl;
1689 userInvDiag_->describe (out, vl);
1692 out <<
"userLambdaMax_: " << userLambdaMax_ << endl
1693 <<
"userLambdaMin_: " << userLambdaMin_ << endl
1694 <<
"userEigRatio_: " << userEigRatio_ << endl
1695 <<
"numIters_: " << numIters_ << endl
1696 <<
"eigMaxIters_: " << eigMaxIters_ << endl
1697 <<
"zeroStartingSolution_: " << zeroStartingSolution_ << endl
1698 <<
"assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1699 <<
"textbookAlgorithm_: " << textbookAlgorithm_ << endl
1700 <<
"computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1708 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \ 1709 template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >; 1711 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
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)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:276
Ifpack2 implementation details.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:199
LinearOp zero(const VectorSpace &vs)
bool isType(const std::string &name) const
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:111
scalar_type ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:107
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:101
Declaration of Chebyshev implementation.
TypeTo as(const TypeFrom &t)
bool isParameter(const std::string &name) const
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
Tpetra::Vector< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > V
Specialization of Tpetra::Vector.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:126