44 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
45 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
55 #include "Ifpack2_Details_ScaledDampedResidual.hpp"
56 #include "Kokkos_ArithTraits.hpp"
57 #include "Teuchos_FancyOStream.hpp"
58 #include "Teuchos_oblackholestream.hpp"
68 const char computeBeforeApplyReminder[] =
69 "This means one of the following:\n"
70 " - you have not yet called compute() on this instance, or \n"
71 " - you didn't call compute() after calling setParameters().\n\n"
72 "After creating an Ifpack2::Chebyshev instance,\n"
73 "you must _always_ call compute() at least once before calling apply().\n"
74 "After calling compute() once, you do not need to call it again,\n"
75 "unless the matrix has changed or you have changed parameters\n"
76 "(by calling setParameters()).";
82 template<
class XV,
class SizeType =
typename XV::
size_type>
83 struct V_ReciprocalThresholdSelfFunctor
85 typedef typename XV::execution_space execution_space;
86 typedef typename XV::non_const_value_type value_type;
87 typedef SizeType size_type;
88 typedef Kokkos::Details::ArithTraits<value_type> KAT;
89 typedef typename KAT::mag_type mag_type;
92 const value_type minVal_;
93 const mag_type minValMag_;
95 V_ReciprocalThresholdSelfFunctor (
const XV& X,
96 const value_type& min_val) :
99 minValMag_ (KAT::abs (min_val))
102 KOKKOS_INLINE_FUNCTION
103 void operator() (
const size_type& i)
const
105 const mag_type X_i_abs = KAT::abs (X_[i]);
107 if (X_i_abs < minValMag_) {
111 X_[i] = KAT::one () / X_[i];
116 template<
class XV,
class SizeType =
typename XV::
size_type>
117 struct LocalReciprocalThreshold {
119 compute (
const XV& X,
120 const typename XV::non_const_value_type& minVal)
122 typedef typename XV::execution_space execution_space;
123 Kokkos::RangePolicy<execution_space, SizeType> policy (0, X.extent (0));
124 V_ReciprocalThresholdSelfFunctor<XV, SizeType> op (X, minVal);
125 Kokkos::parallel_for (policy, op);
129 template <
class TpetraVectorType,
130 const bool classic = TpetraVectorType::node_type::classic>
131 struct GlobalReciprocalThreshold {};
133 template <
class TpetraVectorType>
134 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
136 compute (TpetraVectorType& V,
137 const typename TpetraVectorType::scalar_type& min_val)
139 typedef typename TpetraVectorType::scalar_type scalar_type;
140 typedef typename TpetraVectorType::mag_type mag_type;
141 typedef Kokkos::Details::ArithTraits<scalar_type> STS;
143 const scalar_type ONE = STS::one ();
144 const mag_type min_val_abs = STS::abs (min_val);
147 const size_t lclNumRows = V.getLocalLength ();
149 for (
size_t i = 0; i < lclNumRows; ++i) {
150 const scalar_type V_0i = V_0[i];
151 if (STS::abs (V_0i) < min_val_abs) {
160 template <
class TpetraVectorType>
161 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
163 compute (TpetraVectorType& X,
164 const typename TpetraVectorType::scalar_type& minVal)
166 typedef typename TpetraVectorType::impl_scalar_type value_type;
167 typedef typename TpetraVectorType::device_type::memory_space memory_space;
169 X.template sync<memory_space> ();
170 X.template modify<memory_space> ();
172 const value_type minValS =
static_cast<value_type
> (minVal);
173 auto X_0 = Kokkos::subview (X.template getLocalView<memory_space> (),
175 LocalReciprocalThreshold<decltype (X_0) >::compute (X_0, minValS);
180 template <
typename S,
typename L,
typename G,
typename N>
182 reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V,
const S& minVal)
184 GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
192 template<
class OneDViewType,
193 class LocalOrdinal =
typename OneDViewType::size_type>
194 class PositivizeVector {
195 static_assert (Kokkos::Impl::is_view<OneDViewType>::value,
196 "OneDViewType must be a 1-D Kokkos::View.");
197 static_assert (static_cast<int> (OneDViewType::rank) == 1,
198 "This functor only works with a 1-D View.");
199 static_assert (std::is_integral<LocalOrdinal>::value,
200 "The second template parameter, LocalOrdinal, "
201 "must be an integer type.");
203 PositivizeVector (
const OneDViewType& x) : x_ (x) {}
205 KOKKOS_INLINE_FUNCTION
void
206 operator () (
const LocalOrdinal& i)
const
208 typedef typename OneDViewType::non_const_value_type IST;
209 typedef Kokkos::Details::ArithTraits<IST> STS;
210 typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
212 if (STS::real (x_(i)) < STM::zero ()) {
223 template<
class ScalarType,
class MV>
224 void Chebyshev<ScalarType, MV>::checkInputMatrix ()
const
227 ! A_.is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
228 std::invalid_argument,
229 "Ifpack2::Chebyshev: The input matrix A must be square. "
230 "A has " << A_->getGlobalNumRows () <<
" rows and "
231 << A_->getGlobalNumCols () <<
" columns.");
235 if (debug_ && ! A_.is_null ()) {
243 ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
244 "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be "
245 "the same (in the sense of isSameAs())." << std::endl <<
"We only check "
246 "for this in debug mode.");
251 template<
class ScalarType,
class MV>
253 Chebyshev<ScalarType, MV>::
254 checkConstructorInput ()
const
259 if (STS::isComplex) {
261 (
true, std::logic_error,
"Ifpack2::Chebyshev: This class' implementation "
262 "of Chebyshev iteration only works for real-valued, symmetric positive "
263 "definite matrices. However, you instantiated this class for ScalarType"
265 "complex-valued type. While this may be algorithmically correct if all "
266 "of the complex numbers in the matrix have zero imaginary part, we "
267 "forbid using complex ScalarType altogether in order to remind you of "
268 "the limitations of our implementation (and of the algorithm itself).");
275 template<
class ScalarType,
class MV>
279 savedDiagOffsets_ (false),
280 computedLambdaMax_ (STS::nan ()),
281 computedLambdaMin_ (STS::nan ()),
282 lambdaMaxForApply_ (STS::nan ()),
283 lambdaMinForApply_ (STS::nan ()),
284 eigRatioForApply_ (STS::nan ()),
285 userLambdaMax_ (STS::nan ()),
286 userLambdaMin_ (STS::nan ()),
287 userEigRatio_ (Teuchos::
as<
ST> (30)),
288 minDiagVal_ (STS::eps ()),
291 zeroStartingSolution_ (true),
292 assumeMatrixUnchanged_ (false),
293 textbookAlgorithm_ (false),
294 computeMaxResNorm_ (false),
297 checkConstructorInput ();
300 template<
class ScalarType,
class MV>
305 savedDiagOffsets_ (false),
306 computedLambdaMax_ (STS::nan ()),
307 computedLambdaMin_ (STS::nan ()),
308 lambdaMaxForApply_ (STS::nan ()),
309 boostFactor_ (static_cast<
MT> (1.1)),
310 lambdaMinForApply_ (STS::nan ()),
311 eigRatioForApply_ (STS::nan ()),
312 userLambdaMax_ (STS::nan ()),
313 userLambdaMin_ (STS::nan ()),
314 userEigRatio_ (Teuchos::
as<
ST> (30)),
315 minDiagVal_ (STS::eps ()),
318 zeroStartingSolution_ (true),
319 assumeMatrixUnchanged_ (false),
320 textbookAlgorithm_ (false),
321 computeMaxResNorm_ (false),
324 checkConstructorInput ();
328 template<
class ScalarType,
class MV>
335 using Teuchos::rcp_const_cast;
347 const ST defaultLambdaMax = STS::nan ();
348 const ST defaultLambdaMin = STS::nan ();
355 const ST defaultEigRatio = Teuchos::as<ST> (30);
356 const MT defaultBoostFactor =
static_cast<MT> (1.1);
357 const ST defaultMinDiagVal = STS::eps ();
358 const int defaultNumIters = 1;
359 const int defaultEigMaxIters = 10;
360 const bool defaultZeroStartingSolution =
true;
361 const bool defaultAssumeMatrixUnchanged =
false;
362 const bool defaultTextbookAlgorithm =
false;
363 const bool defaultComputeMaxResNorm =
false;
364 const bool defaultDebug =
false;
370 RCP<const V> userInvDiagCopy;
371 ST lambdaMax = defaultLambdaMax;
372 ST lambdaMin = defaultLambdaMin;
373 ST eigRatio = defaultEigRatio;
374 MT boostFactor = defaultBoostFactor;
375 ST minDiagVal = defaultMinDiagVal;
376 int numIters = defaultNumIters;
377 int eigMaxIters = defaultEigMaxIters;
378 bool zeroStartingSolution = defaultZeroStartingSolution;
379 bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
380 bool textbookAlgorithm = defaultTextbookAlgorithm;
381 bool computeMaxResNorm = defaultComputeMaxResNorm;
382 bool debug = defaultDebug;
389 if (plist.
isType<
bool> (
"debug")) {
390 debug = plist.
get<
bool> (
"debug");
392 else if (plist.
isType<
int> (
"debug")) {
393 const int debugInt = plist.
get<
bool> (
"debug");
394 debug = debugInt != 0;
407 const char opInvDiagLabel[] =
"chebyshev: operator inv diagonal";
410 RCP<const V> userInvDiag;
412 if (plist.
isType<
const V*> (opInvDiagLabel)) {
413 const V* rawUserInvDiag =
414 plist.
get<
const V*> (opInvDiagLabel);
416 userInvDiag =
rcp (rawUserInvDiag,
false);
418 else if (plist.
isType<
const V*> (opInvDiagLabel)) {
419 V* rawUserInvDiag = plist.
get<V*> (opInvDiagLabel);
421 userInvDiag =
rcp (const_cast<const V*> (rawUserInvDiag),
false);
423 else if (plist.
isType<RCP<const V>> (opInvDiagLabel)) {
424 userInvDiag = plist.
get<RCP<const V> > (opInvDiagLabel);
426 else if (plist.
isType<RCP<V>> (opInvDiagLabel)) {
427 RCP<V> userInvDiagNonConst =
428 plist.
get<RCP<V> > (opInvDiagLabel);
429 userInvDiag = rcp_const_cast<
const V> (userInvDiagNonConst);
431 else if (plist.
isType<
const V> (opInvDiagLabel)) {
432 const V& userInvDiagRef = plist.
get<
const V> (opInvDiagLabel);
434 userInvDiag = userInvDiagCopy;
436 else if (plist.
isType<V> (opInvDiagLabel)) {
437 V& userInvDiagNonConstRef = plist.
get<V> (opInvDiagLabel);
438 const V& userInvDiagRef =
const_cast<const V&
> (userInvDiagNonConstRef);
440 userInvDiag = userInvDiagCopy;
450 if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
463 if (plist.
isParameter (
"chebyshev: max eigenvalue")) {
464 if (plist.
isType<
double>(
"chebyshev: max eigenvalue"))
465 lambdaMax = plist.
get<
double> (
"chebyshev: max eigenvalue");
467 lambdaMax = plist.
get<
ST> (
"chebyshev: max eigenvalue");
469 STS::isnaninf (lambdaMax), std::invalid_argument,
470 "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" "
471 "parameter is NaN or Inf. This parameter is optional, but if you "
472 "choose to supply it, it must have a finite value.");
474 if (plist.
isParameter (
"chebyshev: min eigenvalue")) {
475 if (plist.
isType<
double>(
"chebyshev: min eigenvalue"))
476 lambdaMin = plist.
get<
double> (
"chebyshev: min eigenvalue");
478 lambdaMin = plist.
get<
ST> (
"chebyshev: min eigenvalue");
480 STS::isnaninf (lambdaMin), std::invalid_argument,
481 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" "
482 "parameter is NaN or Inf. This parameter is optional, but if you "
483 "choose to supply it, it must have a finite value.");
487 if (plist.
isParameter (
"smoother: Chebyshev alpha")) {
488 if (plist.
isType<
double>(
"smoother: Chebyshev alpha"))
489 eigRatio = plist.
get<
double> (
"smoother: Chebyshev alpha");
491 eigRatio = plist.
get<
ST> (
"smoother: Chebyshev alpha");
494 eigRatio = plist.
get (
"chebyshev: ratio eigenvalue", eigRatio);
496 STS::isnaninf (eigRatio), std::invalid_argument,
497 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" "
498 "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. "
499 "This parameter is optional, but if you choose to supply it, it must have "
507 STS::real (eigRatio) < STS::real (STS::one ()),
508 std::invalid_argument,
509 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\""
510 "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, "
511 "but you supplied the value " << eigRatio <<
".");
516 const char paramName[] =
"chebyshev: boost factor";
520 boostFactor = plist.
get<
MT> (paramName);
522 else if (! std::is_same<double, MT>::value &&
523 plist.
isType<
double> (paramName)) {
524 const double dblBF = plist.
get<
double> (paramName);
525 boostFactor =
static_cast<MT> (dblBF);
529 (
true, std::invalid_argument,
530 "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\""
531 "parameter must have type magnitude_type (MT) or double.");
541 plist.
set (paramName, defaultBoostFactor);
545 "Ifpack2::Chebyshev::setParameters: \"" << paramName <<
"\" parameter "
546 "must be >= 1, but you supplied the value " << boostFactor <<
".");
550 minDiagVal = plist.
get (
"chebyshev: min diagonal value", minDiagVal);
552 STS::isnaninf (minDiagVal), std::invalid_argument,
553 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" "
554 "parameter is NaN or Inf. This parameter is optional, but if you choose "
555 "to supply it, it must have a finite value.");
559 numIters = plist.
get<
int> (
"smoother: sweeps");
562 numIters = plist.
get<
int> (
"relaxation: sweeps");
564 numIters = plist.
get (
"chebyshev: degree", numIters);
566 numIters < 0, std::invalid_argument,
567 "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also "
568 "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a "
569 "nonnegative integer. You gave a value of " << numIters <<
".");
572 if (plist.
isParameter (
"eigen-analysis: iterations")) {
573 eigMaxIters = plist.
get<
int> (
"eigen-analysis: iterations");
575 eigMaxIters = plist.
get (
"chebyshev: eigenvalue max iterations", eigMaxIters);
577 eigMaxIters < 0, std::invalid_argument,
578 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations"
579 "\" parameter (also called \"eigen-analysis: iterations\") must be a "
580 "nonnegative integer. You gave a value of " << eigMaxIters <<
".");
582 zeroStartingSolution = plist.
get (
"chebyshev: zero starting solution",
583 zeroStartingSolution);
584 assumeMatrixUnchanged = plist.
get (
"chebyshev: assume matrix does not change",
585 assumeMatrixUnchanged);
589 if (plist.
isParameter (
"chebyshev: textbook algorithm")) {
590 textbookAlgorithm = plist.
get<
bool> (
"chebyshev: textbook algorithm");
592 if (plist.
isParameter (
"chebyshev: compute max residual norm")) {
593 computeMaxResNorm = plist.
get<
bool> (
"chebyshev: compute max residual norm");
600 (plist.
isType<
bool> (
"chebyshev: use block mode") &&
601 ! plist.
get<
bool> (
"chebyshev: use block mode"),
602 std::invalid_argument,
603 "Ifpack2::Chebyshev requires that if you set \"chebyshev: use "
604 "block mode\" at all, you must set it to false. "
605 "Ifpack2::Chebyshev does not implement Ifpack's block mode.");
607 (plist.
isType<
bool> (
"chebyshev: solve normal equations") &&
608 ! plist.
get<
bool> (
"chebyshev: solve normal equations"),
609 std::invalid_argument,
610 "Ifpack2::Chebyshev does not and will never implement the Ifpack "
611 "parameter \"chebyshev: solve normal equations\". If you want to "
612 "solve the normal equations, construct a Tpetra::Operator that "
613 "implements A^* A, and use Chebyshev to solve A^* A x = A^* b.");
621 std::string eigenAnalysisType (
"power-method");
623 eigenAnalysisType = plist.
get<std::string> (
"eigen-analysis: type");
625 eigenAnalysisType !=
"power-method" &&
626 eigenAnalysisType !=
"power method",
627 std::invalid_argument,
628 "Ifpack2::Chebyshev: This class supports the ML parameter \"eigen-"
629 "analysis: type\" for backwards compatibility. However, Ifpack2 "
630 "currently only supports the \"power-method\" option for this "
631 "parameter. This imitates Ifpack, which only implements the power "
632 "method for eigenanalysis.");
636 userInvDiag_ = userInvDiagCopy;
637 userLambdaMax_ = lambdaMax;
638 userLambdaMin_ = lambdaMin;
639 userEigRatio_ = eigRatio;
640 boostFactor_ =
static_cast<MT> (boostFactor);
641 minDiagVal_ = minDiagVal;
642 numIters_ = numIters;
643 eigMaxIters_ = eigMaxIters;
644 zeroStartingSolution_ = zeroStartingSolution;
645 assumeMatrixUnchanged_ = assumeMatrixUnchanged;
646 textbookAlgorithm_ = textbookAlgorithm;
647 computeMaxResNorm_ = computeMaxResNorm;
653 if (A_.is_null () || A_->getComm ().is_null ()) {
659 myRank = A_->getComm ()->getRank ();
663 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
668 out_ = Teuchos::getFancyOStream (blackHole);
673 out_ = Teuchos::null;
678 template<
class ScalarType,
class MV>
682 sdr_ = Teuchos::null;
684 diagOffsets_ = offsets_type ();
685 savedDiagOffsets_ =
false;
687 computedLambdaMax_ = STS::nan ();
688 computedLambdaMin_ = STS::nan ();
692 template<
class ScalarType,
class MV>
698 if (! assumeMatrixUnchanged_) {
702 sdr_ = Teuchos::null;
716 myRank = A->getComm ()->getRank ();
720 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
724 out_ = Teuchos::getFancyOStream (blackHole);
729 out_ = Teuchos::null;
735 template<
class ScalarType,
class MV>
744 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
745 typename MV::local_ordinal_type,
746 typename MV::global_ordinal_type,
747 typename MV::node_type> crs_matrix_type;
750 A_.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::compute: The input "
751 "matrix A is null. Please call setMatrix() with a nonnull input matrix "
752 "before calling this method.");
767 if (userInvDiag_.is_null ()) {
769 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A_);
772 if (! A_crsMat.
is_null () && A_crsMat->isStaticGraph ()) {
774 const size_t lclNumRows = A_crsMat->getNodeNumRows ();
775 if (diagOffsets_.extent (0) < lclNumRows) {
776 diagOffsets_ = offsets_type ();
777 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
779 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
780 savedDiagOffsets_ =
true;
781 D_ = makeInverseDiagonal (*A_,
true);
784 D_ = makeInverseDiagonal (*A_);
787 else if (! assumeMatrixUnchanged_) {
788 if (! A_crsMat.
is_null () && A_crsMat->isStaticGraph ()) {
791 if (! savedDiagOffsets_) {
792 const size_t lclNumRows = A_crsMat->getNodeNumRows ();
793 if (diagOffsets_.extent (0) < lclNumRows) {
794 diagOffsets_ = offsets_type ();
795 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
797 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
798 savedDiagOffsets_ =
true;
801 D_ = makeInverseDiagonal (*A_,
true);
804 D_ = makeInverseDiagonal (*A_);
809 D_ = makeRangeMapVectorConst (userInvDiag_);
813 const bool computedEigenvalueEstimates =
814 STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
823 if (! assumeMatrixUnchanged_ ||
824 (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
825 const ST computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
827 STS::isnaninf (computedLambdaMax),
829 "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue "
830 "of D^{-1} A failed, by producing Inf or NaN. This probably means that "
831 "the matrix contains Inf or NaN values, or that it is badly scaled.");
833 STS::isnaninf (userEigRatio_),
835 "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
836 << endl <<
"This should be impossible." << endl <<
837 "Please report this bug to the Ifpack2 developers.");
843 const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
846 computedLambdaMax_ = computedLambdaMax;
847 computedLambdaMin_ = computedLambdaMin;
850 STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
852 "Ifpack2::Chebyshev::compute: " << endl <<
853 "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
854 << endl <<
"This should be impossible." << endl <<
855 "Please report this bug to the Ifpack2 developers.");
863 if (STS::isnaninf (userLambdaMax_)) {
864 lambdaMaxForApply_ = computedLambdaMax_;
866 lambdaMaxForApply_ = userLambdaMax_;
879 lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
880 eigRatioForApply_ = userEigRatio_;
882 if (! textbookAlgorithm_) {
885 const ST one = Teuchos::as<ST> (1);
888 if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
889 lambdaMinForApply_ = one;
890 lambdaMaxForApply_ = lambdaMinForApply_;
891 eigRatioForApply_ = one;
897 template<
class ScalarType,
class MV>
901 return lambdaMaxForApply_;
905 template<
class ScalarType,
class MV>
909 const char prefix[] =
"Ifpack2::Chebyshev::apply: ";
912 *out_ <<
"apply: " << std::endl;
915 (A_.is_null (), std::runtime_error, prefix <<
"The input matrix A is null. "
916 " Please call setMatrix() with a nonnull input matrix, and then call "
917 "compute(), before calling this method.");
919 (STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
920 prefix <<
"There is no estimate for the max eigenvalue."
921 << std::endl << computeBeforeApplyReminder);
922 TEUCHOS_TEST_FOR_EXCEPTION
923 (STS::isnaninf (lambdaMinForApply_), std::runtime_error,
924 prefix <<
"There is no estimate for the min eigenvalue."
925 << std::endl << computeBeforeApplyReminder);
926 TEUCHOS_TEST_FOR_EXCEPTION
927 (STS::isnaninf (eigRatioForApply_), std::runtime_error,
928 prefix <<
"There is no estimate for the ratio of the max "
929 "eigenvalue to the min eigenvalue."
930 << std::endl << computeBeforeApplyReminder);
931 TEUCHOS_TEST_FOR_EXCEPTION
932 (D_.is_null (), std::runtime_error, prefix <<
"The vector of inverse "
933 "diagonal entries of the matrix has not yet been computed."
934 << std::endl << computeBeforeApplyReminder);
936 if (textbookAlgorithm_) {
937 textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
938 lambdaMinForApply_, eigRatioForApply_, *D_);
941 ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
942 lambdaMinForApply_, eigRatioForApply_, *D_);
945 if (computeMaxResNorm_ && B.getNumVectors () > 0) {
946 MV R (B.getMap (), B.getNumVectors ());
947 computeResidual (R, B, *A_, X);
950 return *std::max_element (norms.begin (), norms.end ());
957 template<
class ScalarType,
class MV>
962 using Teuchos::rcpFromRef;
963 this->describe (* (Teuchos::getFancyOStream (rcpFromRef (out))),
967 template<
class ScalarType,
class MV>
971 const ScalarType& alpha,
976 solve (W, alpha, D_inv, B);
977 Tpetra::deep_copy (X, W);
980 template<
class ScalarType,
class MV>
986 Tpetra::deep_copy(R, B);
987 A.apply (X, R, mode, -STS::one(), STS::one());
990 template<
class ScalarType,
class MV>
992 Chebyshev<ScalarType, MV>::
993 solve (MV& Z,
const V& D_inv,
const MV& R) {
994 Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
997 template<
class ScalarType,
class MV>
999 Chebyshev<ScalarType, MV>::
1000 solve (MV& Z,
const ST alpha,
const V& D_inv,
const MV& R) {
1001 Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
1004 template<
class ScalarType,
class MV>
1006 Chebyshev<ScalarType, MV>::
1007 makeInverseDiagonal (
const row_matrix_type& A,
const bool useDiagOffsets)
const
1010 using Teuchos::rcpFromRef;
1011 using Teuchos::rcp_dynamic_cast;
1013 RCP<V> D_rowMap (
new V (A.getGraph ()->getRowMap ()));
1014 if (useDiagOffsets) {
1018 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
1019 typename MV::local_ordinal_type,
1020 typename MV::global_ordinal_type,
1021 typename MV::node_type> crs_matrix_type;
1022 RCP<const crs_matrix_type> A_crsMat =
1023 rcp_dynamic_cast<
const crs_matrix_type> (rcpFromRef (A));
1024 if (! A_crsMat.is_null ()) {
1026 ! savedDiagOffsets_, std::logic_error,
1027 "Ifpack2::Details::Chebyshev::makeInverseDiagonal: "
1028 "It is not allowed to call this method with useDiagOffsets=true, "
1029 "if you have not previously saved offsets of diagonal entries. "
1030 "This situation should never arise if this class is used properly. "
1031 "Please report this bug to the Ifpack2 developers.");
1032 A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_);
1038 A.getLocalDiagCopy (*D_rowMap);
1040 RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
1046 D_rangeMap->sync_host ();
1047 auto D_lcl = D_rangeMap->getLocalViewHost ();
1048 auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1050 typedef typename MV::impl_scalar_type IST;
1051 typedef typename MV::local_ordinal_type LO;
1052 typedef Kokkos::Details::ArithTraits<IST> STS;
1053 typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
1055 const LO lclNumRows =
static_cast<LO
> (D_rangeMap->getLocalLength ());
1056 bool foundNonpositiveValue =
false;
1057 for (LO i = 0; i < lclNumRows; ++i) {
1058 if (STS::real (D_lcl_1d(i)) <= STM::zero ()) {
1059 foundNonpositiveValue =
true;
1064 using Teuchos::outArg;
1068 const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1069 int gblSuccess = lclSuccess;
1070 if (! D_rangeMap->getMap ().is_null () && D_rangeMap->getMap ()->getComm ().is_null ()) {
1072 reduceAll<int, int> (comm,
REDUCE_MIN, lclSuccess, outArg (gblSuccess));
1074 if (gblSuccess == 1) {
1075 *out_ <<
"makeInverseDiagonal: The matrix's diagonal entries all have "
1076 "positive real part (this is good for Chebyshev)." << std::endl;
1079 *out_ <<
"makeInverseDiagonal: The matrix's diagonal has at least one "
1080 "entry with nonpositive real part, on at least one process in the "
1081 "matrix's communicator. This is bad for Chebyshev." << std::endl;
1087 reciprocal_threshold (*D_rangeMap, minDiagVal_);
1088 return Teuchos::rcp_const_cast<
const V> (D_rangeMap);
1092 template<
class ScalarType,
class MV>
1094 Chebyshev<ScalarType, MV>::
1099 typedef Tpetra::Export<
typename MV::local_ordinal_type,
1100 typename MV::global_ordinal_type,
1101 typename MV::node_type> export_type;
1106 A_.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1107 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() "
1108 "with a nonnull input matrix before calling this method. This is probably "
1109 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1111 D.
is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1112 "makeRangeMapVector: The input Vector D is null. This is probably "
1113 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1115 RCP<const map_type> sourceMap = D->getMap ();
1116 RCP<const map_type> rangeMap = A_->getRangeMap ();
1117 RCP<const map_type> rowMap = A_->getRowMap ();
1119 if (rangeMap->isSameAs (*sourceMap)) {
1125 RCP<const export_type> exporter;
1129 if (sourceMap->isSameAs (*rowMap)) {
1131 exporter = A_->getGraph ()->getExporter ();
1134 exporter =
rcp (
new export_type (sourceMap, rangeMap));
1137 if (exporter.is_null ()) {
1142 D_out->doExport (*D, *exporter, Tpetra::ADD);
1143 return Teuchos::rcp_const_cast<
const V> (D_out);
1149 template<
class ScalarType,
class MV>
1151 Chebyshev<ScalarType, MV>::
1154 using Teuchos::rcp_const_cast;
1155 return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1159 template<
class ScalarType,
class MV>
1161 Chebyshev<ScalarType, MV>::
1162 textbookApplyImpl (
const op_type& A,
1169 const V& D_inv)
const
1172 const ST myLambdaMin = lambdaMax / eigRatio;
1174 const ST zero = Teuchos::as<ST> (0);
1175 const ST one = Teuchos::as<ST> (1);
1176 const ST two = Teuchos::as<ST> (2);
1177 const ST d = (lambdaMax + myLambdaMin) / two;
1178 const ST c = (lambdaMax - myLambdaMin) / two;
1180 if (zeroStartingSolution_ && numIters > 0) {
1184 MV R (B.getMap (), B.getNumVectors (),
false);
1185 MV P (B.getMap (), B.getNumVectors (),
false);
1186 MV Z (B.getMap (), B.getNumVectors (),
false);
1188 for (
int i = 0; i < numIters; ++i) {
1189 computeResidual (R, B, A, X);
1190 solve (Z, D_inv, R);
1198 beta = alpha * (c/two) * (c/two);
1199 alpha = one / (d - beta);
1200 P.update (one, Z, beta);
1202 X.update (alpha, P, one);
1209 template<
class ScalarType,
class MV>
1210 typename Chebyshev<ScalarType, MV>::MT
1211 Chebyshev<ScalarType, MV>::maxNormInf (
const MV& X) {
1213 X.normInf (norms());
1214 return *std::max_element (norms.begin (), norms.end ());
1217 template<
class ScalarType,
class MV>
1219 Chebyshev<ScalarType, MV>::
1220 ifpackApplyImpl (
const op_type& A,
1230 #ifdef HAVE_IFPACK2_DEBUG
1231 const bool debug = debug_;
1233 const bool debug =
false;
1237 *out_ <<
" \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1238 *out_ <<
" \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1241 if (numIters <= 0) {
1244 const ST zero =
static_cast<ST
> (0.0);
1245 const ST one =
static_cast<ST
> (1.0);
1246 const ST two =
static_cast<ST
> (2.0);
1249 if (lambdaMin == one && lambdaMax == lambdaMin) {
1250 solve (X, D_inv, B);
1255 const ST alpha = lambdaMax / eigRatio;
1256 const ST beta = boostFactor_ * lambdaMax;
1257 const ST delta = two / (beta - alpha);
1258 const ST theta = (beta + alpha) / two;
1259 const ST s1 = theta * delta;
1262 *out_ <<
" alpha = " << alpha << endl
1263 <<
" beta = " << beta << endl
1264 <<
" delta = " << delta << endl
1265 <<
" theta = " << theta << endl
1266 <<
" s1 = " << s1 << endl;
1274 *out_ <<
" Iteration " << 1 <<
":" << endl
1275 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1279 if (! zeroStartingSolution_) {
1282 if (sdr_.is_null ()) {
1284 sdr_ =
Teuchos::rcp (
new ScaledDampedResidual<op_type> (A_op));
1287 sdr_->compute (W, one/theta, const_cast<V&> (D_inv),
1288 const_cast<MV&> (B), X, zero);
1289 X.update (one, W, one);
1293 firstIterationWithZeroStartingSolution (W, one/theta, D_inv, B, X);
1297 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1298 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1301 if (numIters > 1 && sdr_.is_null ()) {
1303 sdr_ =
Teuchos::rcp (
new ScaledDampedResidual<op_type> (A_op));
1308 ST rhokp1, dtemp1, dtemp2;
1309 for (
int deg = 1; deg < numIters; ++deg) {
1311 *out_ <<
" Iteration " << deg+1 <<
":" << endl
1312 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl
1313 <<
" - \\|B\\|_{\\infty} = " << maxNormInf (B) << endl
1314 <<
" - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm ()
1315 << endl <<
" - rhok = " << rhok << endl;
1318 rhokp1 = one / (two * s1 - rhok);
1319 dtemp1 = rhokp1 * rhok;
1320 dtemp2 = two * rhokp1 * delta;
1324 *out_ <<
" - dtemp1 = " << dtemp1 << endl
1325 <<
" - dtemp2 = " << dtemp2 << endl;
1329 sdr_->compute (W, dtemp2, const_cast<V&> (D_inv),
1330 const_cast<MV&> (B), (X), dtemp1);
1331 X.update (one, W, one);
1334 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1335 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1340 template<
class ScalarType,
class MV>
1341 typename Chebyshev<ScalarType, MV>::ST
1342 Chebyshev<ScalarType, MV>::
1343 powerMethodWithInitGuess (
const op_type& A,
1350 *out_ <<
" powerMethodWithInitGuess:" << endl;
1353 const ST zero =
static_cast<ST
> (0.0);
1354 const ST one =
static_cast<ST
> (1.0);
1355 ST lambdaMax = zero;
1356 ST RQ_top, RQ_bottom, norm;
1358 V y (A.getRangeMap ());
1361 (norm == zero, std::runtime_error,
1362 "Ifpack2::Chebyshev::powerMethodWithInitGuess: The initial guess "
1363 "has zero norm. This could be either because Tpetra::Vector::"
1364 "randomize filled the vector with zeros (if that was used to "
1365 "compute the initial guess), or because the norm2 method has a "
1366 "bug. The first is not impossible, but unlikely.");
1369 *out_ <<
" Original norm1(x): " << x.norm1 ()
1370 <<
", norm2(x): " << norm << endl;
1373 x.scale (one / norm);
1376 *out_ <<
" norm1(x.scale(one/norm)): " << x.norm1 () << endl;
1379 for (
int iter = 0; iter < numIters; ++iter) {
1381 *out_ <<
" Iteration " << iter << endl;
1384 solve (y, D_inv, y);
1386 RQ_bottom = x.dot (x);
1388 *out_ <<
" RQ_top: " << RQ_top
1389 <<
", RQ_bottom: " << RQ_bottom << endl;
1391 lambdaMax = RQ_top / RQ_bottom;
1395 *out_ <<
" norm is zero; returning zero" << endl;
1399 x.update (one / norm, y, zero);
1402 *out_ <<
" lambdaMax: " << lambdaMax << endl;
1407 template<
class ScalarType,
class MV>
1409 Chebyshev<ScalarType, MV>::
1410 computeInitialGuessForPowerMethod (V& x,
const bool nonnegativeRealParts)
const
1412 typedef typename MV::device_type::execution_space dev_execution_space;
1413 typedef typename MV::device_type::memory_space dev_memory_space;
1414 typedef typename MV::local_ordinal_type LO;
1418 if (nonnegativeRealParts) {
1419 x.template modify<dev_memory_space> ();
1420 auto x_lcl = x.template getLocalView<dev_memory_space> ();
1421 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
1423 const LO lclNumRows =
static_cast<LO
> (x.getLocalLength ());
1424 Kokkos::RangePolicy<dev_execution_space, LO> range (0, lclNumRows);
1425 PositivizeVector<decltype (x_lcl_1d), LO> functor (x_lcl_1d);
1426 Kokkos::parallel_for (range, functor);
1430 template<
class ScalarType,
class MV>
1431 typename Chebyshev<ScalarType, MV>::ST
1432 Chebyshev<ScalarType, MV>::
1433 powerMethod (
const op_type& A,
const V& D_inv,
const int numIters)
1437 *out_ <<
"powerMethod:" << endl;
1440 const ST zero =
static_cast<ST
> (0.0);
1441 V x (A.getDomainMap ());
1445 computeInitialGuessForPowerMethod (x,
false);
1447 ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1456 if (STS::real (lambdaMax) < STS::real (zero)) {
1458 *out_ <<
"real(lambdaMax) = " << STS::real (lambdaMax) <<
" < 0; "
1459 "try again with a different random initial guess" << endl;
1468 computeInitialGuessForPowerMethod (x,
true);
1469 lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1474 template<
class ScalarType,
class MV>
1480 template<
class ScalarType,
class MV>
1488 template<
class ScalarType,
class MV>
1498 const size_t B_numVecs = B.getNumVectors ();
1499 if (W_.is_null () || W_->getNumVectors () != B_numVecs) {
1500 W_ =
Teuchos::rcp (
new MV (B.getMap (), B_numVecs,
true));
1505 template<
class ScalarType,
class MV>
1509 std::ostringstream oss;
1512 oss <<
"\"Ifpack2::Details::Chebyshev\":"
1514 <<
"degree: " << numIters_
1515 <<
", lambdaMax: " << lambdaMaxForApply_
1516 <<
", alpha: " << eigRatioForApply_
1517 <<
", lambdaMin: " << lambdaMinForApply_
1518 <<
", boost factor: " << boostFactor_
1523 template<
class ScalarType,
class MV>
1550 if (A_.is_null () || A_->getComm ().is_null ()) {
1554 myRank = A_->getComm ()->getRank ();
1559 out <<
"\"Ifpack2::Details::Chebyshev\":" << endl;
1565 out <<
"degree: " << numIters_ << endl
1566 <<
"lambdaMax: " << lambdaMaxForApply_ << endl
1567 <<
"alpha: " << eigRatioForApply_ << endl
1568 <<
"lambdaMin: " << lambdaMinForApply_ << endl
1569 <<
"boost factor: " << boostFactor_ << endl;
1577 out <<
"Template parameters:" << endl;
1580 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1581 <<
"MV: " << TypeNameTraits<MV>::name () << endl;
1587 out <<
"Computed parameters:" << endl;
1599 if (D_.is_null ()) {
1601 out <<
"unset" << endl;
1606 out <<
"set" << endl;
1615 D_->describe (out, vl);
1620 out <<
"W_: " << (W_.is_null () ?
"unset" :
"set") << endl
1621 <<
"computedLambdaMax_: " << computedLambdaMax_ << endl
1622 <<
"computedLambdaMin_: " << computedLambdaMin_ << endl
1623 <<
"lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1624 <<
"lambdaMinForApply_: " << lambdaMinForApply_ << endl
1625 <<
"eigRatioForApply_: " << eigRatioForApply_ << endl;
1630 out <<
"User parameters:" << endl;
1636 out <<
"userInvDiag_: ";
1637 if (userInvDiag_.is_null ()) {
1638 out <<
"unset" << endl;
1641 out <<
"set" << endl;
1647 userInvDiag_->describe (out, vl);
1650 out <<
"userLambdaMax_: " << userLambdaMax_ << endl
1651 <<
"userLambdaMin_: " << userLambdaMin_ << endl
1652 <<
"userEigRatio_: " << userEigRatio_ << endl
1653 <<
"numIters_: " << numIters_ << endl
1654 <<
"eigMaxIters_: " << eigMaxIters_ << endl
1655 <<
"zeroStartingSolution_: " << zeroStartingSolution_ << endl
1656 <<
"assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1657 <<
"textbookAlgorithm_: " << textbookAlgorithm_ << endl
1658 <<
"computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1666 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \
1667 template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >;
1669 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:331
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:907
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:1526
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:695
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:1476
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:277
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:960
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:116
ScalarType ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:112
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:106
Declaration of Chebyshev implementation.
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1508
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:1483
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:737
basic_oblackholestream< char, std::char_traits< char > > oblackholestream