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.extent (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 ()),
286 userEigRatio_ (Teuchos::
as<
ST> (30)),
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 ()),
312 userEigRatio_ (Teuchos::
as<
ST> (30)),
313 minDiagVal_ (STS::eps ()),
316 zeroStartingSolution_ (true),
317 assumeMatrixUnchanged_ (false),
318 textbookAlgorithm_ (false),
319 computeMaxResNorm_ (false),
322 checkConstructorInput ();
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");
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);
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 ()) {
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_.extent (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_.extent (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);
971 TEUCHOS_TEST_FOR_EXCEPTION
972 (STS::isnaninf (lambdaMinForApply_), std::runtime_error,
973 prefix <<
"There is no estimate for the min eigenvalue."
974 << std::endl << computeBeforeApplyReminder);
975 TEUCHOS_TEST_FOR_EXCEPTION
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);
980 TEUCHOS_TEST_FOR_EXCEPTION
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 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
1081 D_rangeMap->template sync<Kokkos::HostSpace> ();
1082 auto D_lcl = D_rangeMap->template getLocalView<Kokkos::HostSpace> ();
1084 D_rangeMap->sync_host ();
1085 auto D_lcl = D_rangeMap->getLocalViewHost ();
1087 auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1089 typedef typename MV::impl_scalar_type IST;
1090 typedef typename MV::local_ordinal_type LO;
1091 typedef Kokkos::Details::ArithTraits<IST> STS;
1092 typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
1094 const LO lclNumRows =
static_cast<LO
> (D_rangeMap->getLocalLength ());
1095 bool foundNonpositiveValue =
false;
1096 for (LO i = 0; i < lclNumRows; ++i) {
1097 if (STS::real (D_lcl_1d(i)) <= STM::zero ()) {
1098 foundNonpositiveValue =
true;
1103 using Teuchos::outArg;
1107 const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1108 int gblSuccess = lclSuccess;
1109 if (! D_rangeMap->getMap ().is_null () && D_rangeMap->getMap ()->getComm ().is_null ()) {
1111 reduceAll<int, int> (comm,
REDUCE_MIN, lclSuccess, outArg (gblSuccess));
1113 if (gblSuccess == 1) {
1114 *out_ <<
"makeInverseDiagonal: The matrix's diagonal entries all have "
1115 "positive real part (this is good for Chebyshev)." << std::endl;
1118 *out_ <<
"makeInverseDiagonal: The matrix's diagonal has at least one "
1119 "entry with nonpositive real part, on at least one process in the "
1120 "matrix's communicator. This is bad for Chebyshev." << std::endl;
1126 reciprocal_threshold (*D_rangeMap, minDiagVal_);
1127 return Teuchos::rcp_const_cast<
const V> (D_rangeMap);
1131 template<
class ScalarType,
class MV>
1133 Chebyshev<ScalarType, MV>::
1138 typedef Tpetra::Export<
typename MV::local_ordinal_type,
1139 typename MV::global_ordinal_type,
1140 typename MV::node_type> export_type;
1145 A_.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1146 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() "
1147 "with a nonnull input matrix before calling this method. This is probably "
1148 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1150 D.
is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1151 "makeRangeMapVector: The input Vector D is null. This is probably "
1152 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1154 RCP<const map_type> sourceMap = D->getMap ();
1155 RCP<const map_type> rangeMap = A_->getRangeMap ();
1156 RCP<const map_type> rowMap = A_->getRowMap ();
1158 if (rangeMap->isSameAs (*sourceMap)) {
1164 RCP<const export_type> exporter;
1168 if (sourceMap->isSameAs (*rowMap)) {
1170 exporter = A_->getGraph ()->getExporter ();
1173 exporter =
rcp (
new export_type (sourceMap, rangeMap));
1176 if (exporter.is_null ()) {
1181 D_out->doExport (*D, *exporter, Tpetra::ADD);
1182 return Teuchos::rcp_const_cast<
const V> (D_out);
1188 template<
class ScalarType,
class MV>
1190 Chebyshev<ScalarType, MV>::
1193 using Teuchos::rcp_const_cast;
1194 return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1198 template<
class ScalarType,
class MV>
1200 Chebyshev<ScalarType, MV>::
1201 textbookApplyImpl (
const op_type& A,
1208 const V& D_inv)
const
1211 const ST myLambdaMin = lambdaMax / eigRatio;
1213 const ST zero = Teuchos::as<ST> (0);
1214 const ST one = Teuchos::as<ST> (1);
1215 const ST two = Teuchos::as<ST> (2);
1216 const ST d = (lambdaMax + myLambdaMin) / two;
1217 const ST c = (lambdaMax - myLambdaMin) / two;
1219 if (zeroStartingSolution_ && numIters > 0) {
1223 MV R (B.getMap (), B.getNumVectors (),
false);
1224 MV P (B.getMap (), B.getNumVectors (),
false);
1225 MV Z (B.getMap (), B.getNumVectors (),
false);
1227 for (
int i = 0; i < numIters; ++i) {
1228 computeResidual (R, B, A, X);
1229 solve (Z, D_inv, R);
1237 beta = alpha * (c/two) * (c/two);
1238 alpha = one / (d - beta);
1239 P.update (one, Z, beta);
1241 X.update (alpha, P, one);
1248 template<
class ScalarType,
class MV>
1249 typename Chebyshev<ScalarType, MV>::MT
1250 Chebyshev<ScalarType, MV>::maxNormInf (
const MV& X) {
1252 X.normInf (norms());
1253 return *std::max_element (norms.begin (), norms.end ());
1256 template<
class ScalarType,
class MV>
1258 Chebyshev<ScalarType, MV>::
1259 ifpackApplyImpl (
const op_type& A,
1269 #ifdef HAVE_IFPACK2_DEBUG
1270 const bool debug = debug_;
1272 const bool debug =
false;
1276 *out_ <<
" \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1277 *out_ <<
" \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1280 if (numIters <= 0) {
1283 const ST one =
static_cast<ST
> (1.0);
1284 const ST two =
static_cast<ST
> (2.0);
1287 if (lambdaMin == one && lambdaMax == lambdaMin) {
1288 solve (X, D_inv, B);
1293 const ST alpha = lambdaMax / eigRatio;
1294 const ST beta = boostFactor_ * lambdaMax;
1295 const ST delta = two / (beta - alpha);
1296 const ST theta = (beta + alpha) / two;
1297 const ST s1 = theta * delta;
1300 *out_ <<
" alpha = " << alpha << endl
1301 <<
" beta = " << beta << endl
1302 <<
" delta = " << delta << endl
1303 <<
" theta = " << theta << endl
1304 <<
" s1 = " << s1 << endl;
1309 makeTempMultiVectors (V_ptr, W_ptr, B);
1319 *out_ <<
" Iteration " << 1 <<
":" << endl
1320 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1324 if (! zeroStartingSolution_) {
1325 computeResidual (V1, B, A, X);
1327 *out_ <<
" - \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1330 solve (W, one/theta, D_inv, V1);
1332 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1334 X.update (one, W, one);
1337 solve (W, one/theta, D_inv, B);
1339 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1341 Tpetra::deep_copy (X, W);
1344 *out_ <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1349 ST rhokp1, dtemp1, dtemp2;
1350 for (
int deg = 1; deg < numIters; ++deg) {
1352 *out_ <<
" Iteration " << deg+1 <<
":" << endl
1353 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl
1354 <<
" - \\|B\\|_{\\infty} = " << maxNormInf (B) << endl
1355 <<
" - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm () << endl
1356 <<
" - rhok = " << rhok << endl;
1357 V1.putScalar (STS::zero ());
1360 computeResidual (V1, B, A, X);
1362 *out_ <<
" - \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1365 rhokp1 = one / (two * s1 - rhok);
1366 dtemp1 = rhokp1 * rhok;
1367 dtemp2 = two * rhokp1 * delta;
1371 *out_ <<
" - dtemp1 = " << dtemp1 << endl
1372 <<
" - dtemp2 = " << dtemp2 << endl;
1376 W.elementWiseMultiply (dtemp2, D_inv, V1, dtemp1);
1377 X.update (one, W, one);
1380 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1381 *out_ <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1386 template<
class ScalarType,
class MV>
1387 typename Chebyshev<ScalarType, MV>::ST
1388 Chebyshev<ScalarType, MV>::
1389 powerMethodWithInitGuess (
const op_type& A,
1396 *out_ <<
" powerMethodWithInitGuess:" << endl;
1399 const ST zero =
static_cast<ST
> (0.0);
1400 const ST one =
static_cast<ST
> (1.0);
1401 ST lambdaMax = zero;
1402 ST RQ_top, RQ_bottom, norm;
1404 V y (A.getRangeMap ());
1407 (norm == zero, std::runtime_error,
1408 "Ifpack2::Chebyshev::powerMethodWithInitGuess: "
1409 "The initial guess has zero norm. This could be either because Tpetra::"
1410 "Vector::randomize() filled the vector with zeros (if that was used to "
1411 "compute the initial guess), or because the norm2 method has a bug. The "
1412 "first is not impossible, but unlikely.");
1415 *out_ <<
" Original norm1(x): " << x.norm1 () <<
", norm2(x): " << norm << endl;
1418 x.scale (one / norm);
1421 *out_ <<
" norm1(x.scale(one/norm)): " << x.norm1 () << endl;
1424 for (
int iter = 0; iter < numIters; ++iter) {
1426 *out_ <<
" Iteration " << iter << endl;
1429 solve (y, D_inv, y);
1431 RQ_bottom = x.dot (x);
1433 *out_ <<
" RQ_top: " << RQ_top <<
", RQ_bottom: " << RQ_bottom << endl;
1435 lambdaMax = RQ_top / RQ_bottom;
1439 *out_ <<
" norm is zero; returning zero" << endl;
1443 x.update (one / norm, y, zero);
1446 *out_ <<
" lambdaMax: " << lambdaMax << endl;
1451 template<
class ScalarType,
class MV>
1453 Chebyshev<ScalarType, MV>::
1454 computeInitialGuessForPowerMethod (V& x,
const bool nonnegativeRealParts)
const
1456 typedef typename MV::device_type::execution_space dev_execution_space;
1457 typedef typename MV::device_type::memory_space dev_memory_space;
1458 typedef typename MV::local_ordinal_type LO;
1462 if (nonnegativeRealParts) {
1463 x.template modify<dev_memory_space> ();
1464 auto x_lcl = x.template getLocalView<dev_memory_space> ();
1465 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
1467 const LO lclNumRows =
static_cast<LO
> (x.getLocalLength ());
1468 Kokkos::RangePolicy<dev_execution_space, LO> range (0, lclNumRows);
1469 PositivizeVector<decltype (x_lcl_1d), LO> functor (x_lcl_1d);
1470 Kokkos::parallel_for (range, functor);
1474 template<
class ScalarType,
class MV>
1475 typename Chebyshev<ScalarType, MV>::ST
1476 Chebyshev<ScalarType, MV>::
1477 powerMethod (
const op_type& A,
const V& D_inv,
const int numIters)
1481 *out_ <<
"powerMethod:" << endl;
1484 const ST zero =
static_cast<ST
> (0.0);
1485 V x (A.getDomainMap ());
1489 computeInitialGuessForPowerMethod (x,
false);
1491 ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1500 if (STS::real (lambdaMax) < STS::real (zero)) {
1502 *out_ <<
"real(lambdaMax) = " << STS::real (lambdaMax) <<
" < 0; try "
1503 "again with a different random initial guess" << endl;
1512 computeInitialGuessForPowerMethod (x,
true);
1513 lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1518 template<
class ScalarType,
class MV>
1524 template<
class ScalarType,
class MV>
1532 template<
class ScalarType,
class MV>
1542 if (V_.is_null () || V_->getNumVectors () != X.getNumVectors ()) {
1543 V_ =
Teuchos::rcp (
new MV (X.getMap (), X.getNumVectors (),
false));
1546 if (W_.is_null () || W_->getNumVectors () != X.getNumVectors ()) {
1547 W_ =
Teuchos::rcp (
new MV (X.getMap (), X.getNumVectors (),
true));
1553 template<
class ScalarType,
class MV>
1557 std::ostringstream oss;
1560 oss <<
"\"Ifpack2::Details::Chebyshev\":"
1562 <<
"degree: " << numIters_
1563 <<
", lambdaMax: " << lambdaMaxForApply_
1564 <<
", alpha: " << eigRatioForApply_
1565 <<
", lambdaMin: " << lambdaMinForApply_
1566 <<
", boost factor: " << boostFactor_
1571 template<
class ScalarType,
class MV>
1598 if (A_.is_null () || A_->getComm ().is_null ()) {
1602 myRank = A_->getComm ()->getRank ();
1607 out <<
"\"Ifpack2::Details::Chebyshev\":" << endl;
1613 out <<
"degree: " << numIters_ << endl
1614 <<
"lambdaMax: " << lambdaMaxForApply_ << endl
1615 <<
"alpha: " << eigRatioForApply_ << endl
1616 <<
"lambdaMin: " << lambdaMinForApply_ << endl
1617 <<
"boost factor: " << boostFactor_ << endl;
1625 out <<
"Template parameters:" << endl;
1628 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1629 <<
"MV: " << TypeNameTraits<MV>::name () << endl;
1635 out <<
"Computed parameters:" << endl;
1647 if (D_.is_null ()) {
1649 out <<
"unset" << endl;
1654 out <<
"set" << endl;
1663 D_->describe (out, vl);
1668 out <<
"V_: " << (V_.is_null () ?
"unset" :
"set") << endl
1669 <<
"W_: " << (W_.is_null () ?
"unset" :
"set") << endl
1670 <<
"computedLambdaMax_: " << computedLambdaMax_ << endl
1671 <<
"computedLambdaMin_: " << computedLambdaMin_ << endl
1672 <<
"lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1673 <<
"lambdaMinForApply_: " << lambdaMinForApply_ << endl
1674 <<
"eigRatioForApply_: " << eigRatioForApply_ << endl;
1679 out <<
"User parameters:" << endl;
1685 out <<
"userInvDiag_: ";
1686 if (userInvDiag_.is_null ()) {
1687 out <<
"unset" << endl;
1690 out <<
"set" << endl;
1696 userInvDiag_->describe (out, vl);
1699 out <<
"userLambdaMax_: " << userLambdaMax_ << endl
1700 <<
"userLambdaMin_: " << userLambdaMin_ << endl
1701 <<
"userEigRatio_: " << userEigRatio_ << endl
1702 <<
"numIters_: " << numIters_ << endl
1703 <<
"eigMaxIters_: " << eigMaxIters_ << endl
1704 <<
"zeroStartingSolution_: " << zeroStartingSolution_ << endl
1705 <<
"assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1706 <<
"textbookAlgorithm_: " << textbookAlgorithm_ << endl
1707 <<
"computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1715 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \
1716 template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >;
1718 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:329
MT apply(const MV &B, MV &X)
Solve Ax=b for x with Chebyshev iteration with left diagonal scaling.
Definition: Ifpack2_Details_Chebyshev_def.hpp:956
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print a description of the Chebyshev solver to out.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1574
T & get(const std::string &name, T def_value)
void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set the matrix.
Definition: Ifpack2_Details_Chebyshev_def.hpp:745
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)
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the matrix given to the constructor.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1520
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:276
bool isParameter(const std::string &name) const
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
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
void print(std::ostream &out)
Print instance data to the given output stream.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1009
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:111
ScalarType 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.
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1556
TypeTo as(const TypeFrom &t)
bool isType(const std::string &name) const
bool hasTransposeApply() const
Whether it's possible to apply the transpose of this operator.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1527
void compute()
(Re)compute the left scaling D_inv, and estimate min and max eigenvalues of D_inv * A...
Definition: Ifpack2_Details_Chebyshev_def.hpp:786