70 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP 71 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP 81 #include <Kokkos_ArithTraits.hpp> 97 const char computeBeforeApplyReminder[] =
98 "This means one of the following:\n" 99 " - you have not yet called compute() on this instance, or \n" 100 " - you didn't call compute() after calling setParameters().\n\n" 101 "After creating an Ifpack2::Chebyshev instance,\n" 102 "you must _always_ call compute() at least once before calling apply().\n" 103 "After calling compute() once, you do not need to call it again,\n" 104 "unless the matrix has changed or you have changed parameters\n" 105 "(by calling setParameters()).";
111 template<
class XV,
class SizeType =
typename XV::
size_type>
112 struct V_ReciprocalThresholdSelfFunctor
114 typedef typename XV::execution_space execution_space;
115 typedef typename XV::non_const_value_type value_type;
116 typedef SizeType size_type;
117 typedef Kokkos::Details::ArithTraits<value_type> KAT;
118 typedef typename KAT::mag_type mag_type;
121 const value_type m_min_val;
122 const mag_type m_min_val_mag;
124 V_ReciprocalThresholdSelfFunctor (
const XV& X,
125 const value_type& min_val) :
128 m_min_val_mag (KAT::abs (min_val))
131 KOKKOS_INLINE_FUNCTION
132 void operator() (
const size_type& i)
const 134 if (KAT::abs (X_(i)) < m_min_val_mag) {
138 X_[i] = KAT::one () / X_[i];
143 template<
class XV,
class SizeType =
typename XV::
size_type>
144 struct LocalReciprocalThreshold {
146 compute (
const XV& X,
147 const typename XV::non_const_value_type& minVal)
149 typedef typename XV::execution_space execution_space;
150 Kokkos::RangePolicy<execution_space, SizeType> policy (0, X.dimension_0 ());
151 V_ReciprocalThresholdSelfFunctor<XV, SizeType> op (X, minVal);
152 Kokkos::parallel_for (policy, op);
156 template <
class TpetraVectorType,
157 const bool classic = TpetraVectorType::node_type::classic>
158 struct GlobalReciprocalThreshold {};
160 template <
class TpetraVectorType>
161 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
163 compute (TpetraVectorType& V,
164 const typename TpetraVectorType::scalar_type& min_val)
166 typedef typename TpetraVectorType::scalar_type scalar_type;
167 typedef typename TpetraVectorType::mag_type mag_type;
168 typedef Kokkos::Details::ArithTraits<scalar_type> STS;
170 const scalar_type ONE = STS::one ();
171 const mag_type min_val_abs = STS::abs (min_val);
173 Teuchos::ArrayRCP<scalar_type> V_0 = V.getDataNonConst (0);
174 const size_t lclNumRows = V.getLocalLength ();
176 for (
size_t i = 0; i < lclNumRows; ++i) {
177 const scalar_type V_0i = V_0[i];
178 if (STS::abs (V_0i) < min_val_abs) {
187 template <
class TpetraVectorType>
188 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
190 compute (TpetraVectorType& X,
191 const typename TpetraVectorType::scalar_type& minVal)
193 typedef typename TpetraVectorType::dual_view_type::non_const_value_type value_type;
194 typedef typename TpetraVectorType::execution_space::memory_space memory_space;
196 auto X_lcl = X.getDualView ();
197 X_lcl.template sync<memory_space> ();
198 X_lcl.template modify<memory_space> ();
200 const value_type minValS =
static_cast<value_type
> (minVal);
202 auto X_0 = Kokkos::subview (X_lcl.d_view, Kokkos::ALL (), 0);
203 LocalReciprocalThreshold<decltype (X_0) >::compute (X_0, minValS);
208 template <
typename S,
typename L,
typename G,
typename N>
210 reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V,
const S& minVal)
212 GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
215 template<
class ScalarType,
class MV>
216 void Chebyshev<ScalarType, MV>::checkInputMatrix ()
const 218 TEUCHOS_TEST_FOR_EXCEPTION(
219 ! A_.is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
220 std::invalid_argument,
221 "Ifpack2::Chebyshev: The input matrix A must be square. " 222 "A has " << A_->getGlobalNumRows () <<
" rows and " 223 << A_->getGlobalNumCols () <<
" columns.");
227 #ifdef HAVE_TEUCHOS_DEBUG 228 if (! A_.is_null ()) {
229 Teuchos::RCP<const map_type> domainMap = A_->getDomainMap ();
230 Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap ();
235 TEUCHOS_TEST_FOR_EXCEPTION(
236 ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
237 "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be " 238 "the same (in the sense of isSameAs())." << std::endl <<
"We only check " 239 "for this if Trilinos was built with the CMake configuration option " 240 "Teuchos_ENABLE_DEBUG set to ON.");
242 #endif // HAVE_TEUCHOS_DEBUG 246 template<
class ScalarType,
class MV>
248 Chebyshev<ScalarType, MV>::
249 checkConstructorInput ()
const 251 TEUCHOS_TEST_FOR_EXCEPTION(STS::isComplex, std::logic_error,
252 "Ifpack2::Chebyshev: This class' implementation of Chebyshev iteration " 253 "only works for real-valued, symmetric positive definite matrices. " 254 "However, you instantiated this class for ScalarType = " 255 << Teuchos::TypeNameTraits<ScalarType>::name () <<
", which is a complex-" 256 "valued type. While this may be algorithmically correct if all of the " 257 "complex numbers in the matrix have zero imaginary part, we forbid using " 258 "complex ScalarType altogether in order to remind you of the limitations " 259 "of our implementation (and of the algorithm itself).");
264 template<
class ScalarType,
class MV>
268 savedDiagOffsets_ (false),
269 computedLambdaMax_ (STS::nan ()),
270 computedLambdaMin_ (STS::nan ()),
271 lambdaMaxForApply_ (STS::nan ()),
272 lambdaMinForApply_ (STS::nan ()),
273 eigRatioForApply_ (STS::nan ()),
274 userLambdaMax_ (STS::nan ()),
275 userLambdaMin_ (STS::nan ()),
277 minDiagVal_ (STS::eps ()),
280 zeroStartingSolution_ (true),
281 assumeMatrixUnchanged_ (false),
282 textbookAlgorithm_ (false),
283 computeMaxResNorm_ (false)
285 checkConstructorInput ();
288 template<
class ScalarType,
class MV>
290 Chebyshev (Teuchos::RCP<const row_matrix_type> A, Teuchos::ParameterList& params) :
292 savedDiagOffsets_ (false),
293 computedLambdaMax_ (STS::nan ()),
294 computedLambdaMin_ (STS::nan ()),
295 lambdaMaxForApply_ (STS::nan ()),
296 lambdaMinForApply_ (STS::nan ()),
297 eigRatioForApply_ (STS::nan ()),
298 userLambdaMax_ (STS::nan ()),
299 userLambdaMin_ (STS::nan ()),
301 minDiagVal_ (STS::eps ()),
304 zeroStartingSolution_ (true),
305 assumeMatrixUnchanged_ (false),
306 textbookAlgorithm_ (false),
307 computeMaxResNorm_ (false)
309 checkConstructorInput ();
313 template<
class ScalarType,
class MV>
320 using Teuchos::rcp_const_cast;
332 const ST defaultLambdaMax = STS::nan ();
333 const ST defaultLambdaMin = STS::nan ();
340 const ST defaultEigRatio = Teuchos::as<ST> (30);
341 const ST defaultMinDiagVal = STS::eps ();
342 const int defaultNumIters = 1;
343 const int defaultEigMaxIters = 10;
344 const bool defaultZeroStartingSolution =
true;
345 const bool defaultAssumeMatrixUnchanged =
false;
346 const bool defaultTextbookAlgorithm =
false;
347 const bool defaultComputeMaxResNorm =
false;
353 RCP<const V> userInvDiagCopy;
354 ST lambdaMax = defaultLambdaMax;
355 ST lambdaMin = defaultLambdaMin;
356 ST eigRatio = defaultEigRatio;
357 ST minDiagVal = defaultMinDiagVal;
358 int numIters = defaultNumIters;
359 int eigMaxIters = defaultEigMaxIters;
360 bool zeroStartingSolution = defaultZeroStartingSolution;
361 bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
362 bool textbookAlgorithm = defaultTextbookAlgorithm;
363 bool computeMaxResNorm = defaultComputeMaxResNorm;
379 if (plist.isParameter (
"chebyshev: operator inv diagonal")) {
381 RCP<const V> userInvDiag;
384 const V* rawUserInvDiag =
385 plist.get<
const V*> (
"chebyshev: operator inv diagonal");
387 userInvDiag = rcp (rawUserInvDiag,
false);
388 }
catch (Teuchos::Exceptions::InvalidParameterType&) {
390 if (userInvDiag.is_null ()) {
392 V* rawUserInvDiag = plist.get<V*> (
"chebyshev: operator inv diagonal");
394 userInvDiag = rcp (const_cast<const V*> (rawUserInvDiag),
false);
395 }
catch (Teuchos::Exceptions::InvalidParameterType&) {
398 if (userInvDiag.is_null ()) {
401 plist.get<RCP<const V> > (
"chebyshev: operator inv diagonal");
402 }
catch (Teuchos::Exceptions::InvalidParameterType&) {
405 if (userInvDiag.is_null ()) {
407 RCP<V> userInvDiagNonConst =
408 plist.get<RCP<V> > (
"chebyshev: operator inv diagonal");
409 userInvDiag = rcp_const_cast<
const V> (userInvDiagNonConst);
410 }
catch (Teuchos::Exceptions::InvalidParameterType&) {
413 if (userInvDiag.is_null ()) {
422 const V& userInvDiagRef =
423 plist.get<
const V> (
"chebyshev: operator inv diagonal");
424 userInvDiagCopy = rcp (
new V (userInvDiagRef, Teuchos::Copy));
426 userInvDiag = userInvDiagCopy;
427 }
catch (Teuchos::Exceptions::InvalidParameterType&) {
430 TEUCHOS_TEST_FOR_EXCEPTION(
431 true, std::runtime_error,
432 "Ifpack2::Chebyshev::setParameters: \"chebyshev: operator inv diagonal\" " 433 "plist.get<const V> does not compile around return held == other_held " 434 "in Teuchos::any in Visual Studio. Can't fix it now, so throwing " 435 "in case someone builds there.");
438 if (userInvDiag.is_null ()) {
447 V& userInvDiagNonConstRef =
448 plist.get<V> (
"chebyshev: operator inv diagonal");
449 const V& userInvDiagRef =
const_cast<const V&
> (userInvDiagNonConstRef);
450 userInvDiagCopy = rcp (
new V (userInvDiagRef, Teuchos::Copy));
452 userInvDiag = userInvDiagCopy;
453 }
catch (Teuchos::Exceptions::InvalidParameterType&) {
456 TEUCHOS_TEST_FOR_EXCEPTION(
457 true, std::runtime_error,
458 "Ifpack2::Chebyshev::setParameters: \"chebyshev: operator inv diagonal\" " 459 "plist.get<V> does not compile around return held == other_held " 460 "in Teuchos::any in Visual Studio. Can't fix it now, so throwing " 461 "in case someone builds there.");
472 if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
473 userInvDiagCopy = rcp (
new V (*userInvDiag, Teuchos::Copy));
485 if (plist.isParameter (
"chebyshev: max eigenvalue")) {
486 if (plist.isType<
double>(
"chebyshev: max eigenvalue"))
487 lambdaMax = plist.get<
double> (
"chebyshev: max eigenvalue");
489 lambdaMax = plist.get<
ST> (
"chebyshev: max eigenvalue");
490 TEUCHOS_TEST_FOR_EXCEPTION(
491 STS::isnaninf (lambdaMax), std::invalid_argument,
492 "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" " 493 "parameter is NaN or Inf. This parameter is optional, but if you " 494 "choose to supply it, it must have a finite value.");
496 if (plist.isParameter (
"chebyshev: min eigenvalue")) {
497 if (plist.isType<
double>(
"chebyshev: min eigenvalue"))
498 lambdaMin = plist.get<
double> (
"chebyshev: min eigenvalue");
500 lambdaMin = plist.get<
ST> (
"chebyshev: min eigenvalue");
501 TEUCHOS_TEST_FOR_EXCEPTION(
502 STS::isnaninf (lambdaMin), std::invalid_argument,
503 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" " 504 "parameter is NaN or Inf. This parameter is optional, but if you " 505 "choose to supply it, it must have a finite value.");
509 if (plist.isParameter (
"smoother: Chebyshev alpha")) {
510 if (plist.isType<
double>(
"smoother: Chebyshev alpha"))
511 eigRatio = plist.get<
double> (
"smoother: Chebyshev alpha");
513 eigRatio = plist.get<
ST> (
"smoother: Chebyshev alpha");
516 eigRatio = plist.get (
"chebyshev: ratio eigenvalue", eigRatio);
517 TEUCHOS_TEST_FOR_EXCEPTION(
518 STS::isnaninf (eigRatio), std::invalid_argument,
519 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" " 520 "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. " 521 "This parameter is optional, but if you choose to supply it, it must have " 528 TEUCHOS_TEST_FOR_EXCEPTION(
529 STS::real (eigRatio) < STS::real (STS::one ()),
530 std::invalid_argument,
531 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\"" 532 "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, " 533 "but you supplied the value " << eigRatio <<
".");
536 minDiagVal = plist.get (
"chebyshev: min diagonal value", minDiagVal);
537 TEUCHOS_TEST_FOR_EXCEPTION(
538 STS::isnaninf (minDiagVal), std::invalid_argument,
539 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" " 540 "parameter is NaN or Inf. This parameter is optional, but if you choose " 541 "to supply it, it must have a finite value.");
544 if (plist.isParameter (
"smoother: sweeps")) {
545 numIters = plist.get<
int> (
"smoother: sweeps");
547 if (plist.isParameter (
"relaxation: sweeps")) {
548 numIters = plist.get<
int> (
"relaxation: sweeps");
550 numIters = plist.get (
"chebyshev: degree", numIters);
551 TEUCHOS_TEST_FOR_EXCEPTION(
552 numIters < 0, std::invalid_argument,
553 "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also " 554 "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a " 555 "nonnegative integer. You gave a value of " << numIters <<
".");
558 if (plist.isParameter (
"eigen-analysis: iterations")) {
559 eigMaxIters = plist.get<
int> (
"eigen-analysis: iterations");
561 eigMaxIters = plist.get (
"chebyshev: eigenvalue max iterations", eigMaxIters);
562 TEUCHOS_TEST_FOR_EXCEPTION(
563 eigMaxIters < 0, std::invalid_argument,
564 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations" 565 "\" parameter (also called \"eigen-analysis: iterations\") must be a " 566 "nonnegative integer. You gave a value of " << eigMaxIters <<
".");
568 zeroStartingSolution = plist.get (
"chebyshev: zero starting solution",
569 zeroStartingSolution);
570 assumeMatrixUnchanged = plist.get (
"chebyshev: assume matrix does not change",
571 assumeMatrixUnchanged);
575 if (plist.isParameter (
"chebyshev: textbook algorithm")) {
576 textbookAlgorithm = plist.get<
bool> (
"chebyshev: textbook algorithm");
578 if (plist.isParameter (
"chebyshev: compute max residual norm")) {
579 computeMaxResNorm = plist.get<
bool> (
"chebyshev: compute max residual norm");
585 TEUCHOS_TEST_FOR_EXCEPTION(
586 plist.isParameter (
"chebyshev: use block mode") &&
587 ! plist.get<
bool> (
"chebyshev: use block mode"),
588 std::invalid_argument,
589 "Ifpack2::Chebyshev requires that if you supply the Ifpack parameter " 590 "\"chebyshev: use block mode\", it must be set to false. Ifpack2's " 591 "Chebyshev does not implement Ifpack's block mode.");
592 TEUCHOS_TEST_FOR_EXCEPTION(
593 plist.isParameter (
"chebyshev: solve normal equations") &&
594 ! plist.get<
bool> (
"chebyshev: solve normal equations"),
595 std::invalid_argument,
596 "Ifpack2::Chebyshev does not and will never implement the Ifpack " 597 "parameter \"chebyshev: solve normal equations\". If you want to solve " 598 "the normal equations, construct a Tpetra::Operator that implements " 599 "A^* A, and use Chebyshev to solve A^* A x = A^* b.");
607 std::string eigenAnalysisType (
"power-method");
608 if (plist.isParameter (
"eigen-analysis: type")) {
609 eigenAnalysisType = plist.get<std::string> (
"eigen-analysis: type");
610 TEUCHOS_TEST_FOR_EXCEPTION(
611 eigenAnalysisType !=
"power-method" &&
612 eigenAnalysisType !=
"power method",
613 std::invalid_argument,
614 "Ifpack2::Chebyshev: This class supports the ML parameter \"eigen-" 615 "analysis: type\" for backwards compatibility. However, Ifpack2 " 616 "currently only supports the \"power-method\" option for this " 617 "parameter. This imitates Ifpack, which only implements the power " 618 "method for eigenanalysis.");
622 userInvDiag_ = userInvDiagCopy;
623 userLambdaMax_ = lambdaMax;
624 userLambdaMin_ = lambdaMin;
625 userEigRatio_ = eigRatio;
626 minDiagVal_ = minDiagVal;
627 numIters_ = numIters;
628 eigMaxIters_ = eigMaxIters;
629 zeroStartingSolution_ = zeroStartingSolution;
630 assumeMatrixUnchanged_ = assumeMatrixUnchanged;
631 textbookAlgorithm_ = textbookAlgorithm;
632 computeMaxResNorm_ = computeMaxResNorm;
636 template<
class ScalarType,
class MV>
641 diagOffsets_ = Teuchos::null;
642 savedDiagOffsets_ =
false;
645 computedLambdaMax_ = STS::nan ();
646 computedLambdaMin_ = STS::nan ();
650 template<
class ScalarType,
class MV>
654 if (A.getRawPtr () != A_.getRawPtr ()) {
655 if (! assumeMatrixUnchanged_) {
663 template<
class ScalarType,
class MV>
672 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
673 typename MV::local_ordinal_type,
674 typename MV::global_ordinal_type,
675 typename MV::node_type> crs_matrix_type;
677 TEUCHOS_TEST_FOR_EXCEPTION(
678 A_.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::compute: The input " 679 "matrix A is null. Please call setMatrix() with a nonnull input matrix " 680 "before calling this method.");
695 if (userInvDiag_.is_null ()) {
696 Teuchos::RCP<const crs_matrix_type> A_crsMat =
697 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A_);
700 if (! A_crsMat.is_null () && A_crsMat->isStaticGraph ()) {
702 A_crsMat->getLocalDiagOffsets (diagOffsets_);
703 savedDiagOffsets_ =
true;
704 D_ = makeInverseDiagonal (*A_,
true);
707 D_ = makeInverseDiagonal (*A_);
710 else if (! assumeMatrixUnchanged_) {
711 if (! A_crsMat.is_null () && A_crsMat->isStaticGraph ()) {
714 if (! savedDiagOffsets_) {
715 A_crsMat->getLocalDiagOffsets (diagOffsets_);
716 savedDiagOffsets_ =
true;
719 D_ = makeInverseDiagonal (*A_,
true);
722 D_ = makeInverseDiagonal (*A_);
727 D_ = makeRangeMapVectorConst (userInvDiag_);
731 const bool computedEigenvalueEstimates =
732 STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
741 if (! assumeMatrixUnchanged_ ||
742 (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
743 const ST computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
744 TEUCHOS_TEST_FOR_EXCEPTION(
745 STS::isnaninf (computedLambdaMax),
747 "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue " 748 "of D^{-1} A failed, by producing Inf or NaN. This probably means that " 749 "the matrix contains Inf or NaN values, or that it is badly scaled.");
750 TEUCHOS_TEST_FOR_EXCEPTION(
751 STS::isnaninf (userEigRatio_),
753 "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN." 754 << endl <<
"This should be impossible." << endl <<
755 "Please report this bug to the Ifpack2 developers.");
761 const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
764 computedLambdaMax_ = computedLambdaMax;
765 computedLambdaMin_ = computedLambdaMin;
767 TEUCHOS_TEST_FOR_EXCEPTION(
768 STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
770 "Ifpack2::Chebyshev::compute: " << endl <<
771 "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN." 772 << endl <<
"This should be impossible." << endl <<
773 "Please report this bug to the Ifpack2 developers.");
781 if (STS::isnaninf (userLambdaMax_)) {
782 lambdaMaxForApply_ = computedLambdaMax_;
784 lambdaMaxForApply_ = userLambdaMax_;
797 lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
798 eigRatioForApply_ = userEigRatio_;
800 if (! textbookAlgorithm_) {
803 const ST one = Teuchos::as<ST> (1);
806 if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
807 lambdaMinForApply_ = one;
808 lambdaMaxForApply_ = lambdaMinForApply_;
809 eigRatioForApply_ = one;
815 template<
class ScalarType,
class MV>
819 return lambdaMaxForApply_;
823 template<
class ScalarType,
class MV>
827 TEUCHOS_TEST_FOR_EXCEPTION(
828 A_.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::apply: The input " 829 "matrix A is null. Please call setMatrix() with a nonnull input matrix, " 830 "and then call compute(), before calling this method.");
831 TEUCHOS_TEST_FOR_EXCEPTION(
832 STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
833 "Ifpack2::Chebyshev::apply: There is no estimate for the max eigenvalue." 834 << std::endl << computeBeforeApplyReminder);
835 TEUCHOS_TEST_FOR_EXCEPTION(
836 STS::isnaninf (lambdaMinForApply_), std::runtime_error,
837 "Ifpack2::Chebyshev::apply: There is no estimate for the min eigenvalue." 838 << std::endl << computeBeforeApplyReminder);
839 TEUCHOS_TEST_FOR_EXCEPTION(
840 STS::isnaninf (eigRatioForApply_), std::runtime_error,
841 "Ifpack2::Chebyshev::apply: There is no estimate for the ratio of the max " 842 "eigenvalue to the min eigenvalue." 843 << std::endl << computeBeforeApplyReminder);
844 TEUCHOS_TEST_FOR_EXCEPTION(
845 D_.is_null (), std::runtime_error,
846 "Ifpack2::Chebyshev::apply: " 847 "The vector of inverse diagonal entries of the matrix has not yet been " 848 "computed." << std::endl << computeBeforeApplyReminder);
850 if (textbookAlgorithm_) {
851 textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
852 lambdaMinForApply_, eigRatioForApply_, *D_);
854 ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
855 lambdaMinForApply_, eigRatioForApply_, *D_);
858 if (computeMaxResNorm_ && B.getNumVectors () > 0) {
859 MV R (B.getMap (), B.getNumVectors ());
860 computeResidual (R, B, *A_, X);
861 Teuchos::Array<MT> norms (B.getNumVectors ());
863 return *std::max_element (norms.begin (), norms.end ());
865 return Teuchos::ScalarTraits<MT>::zero ();
869 template<
class ScalarType,
class MV>
873 this->
describe (* (Teuchos::getFancyOStream (Teuchos::rcpFromRef (out))),
874 Teuchos::VERB_MEDIUM);
877 template<
class ScalarType,
class MV>
881 const Teuchos::ETransp mode)
883 Tpetra::deep_copy(R, B);
884 A.apply (X, R, mode, -STS::one(), STS::one());
887 template<
class ScalarType,
class MV>
890 solve (MV& Z,
const V& D_inv,
const MV& R) {
891 Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
894 template<
class ScalarType,
class MV>
897 solve (MV& Z,
const ST alpha,
const V& D_inv,
const MV& R) {
898 Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
901 template<
class ScalarType,
class MV>
902 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
907 using Teuchos::rcpFromRef;
908 using Teuchos::rcp_dynamic_cast;
910 RCP<V> D_rowMap (
new V (A.getGraph ()->getRowMap ()));
911 if (useDiagOffsets) {
915 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
916 typename MV::local_ordinal_type,
917 typename MV::global_ordinal_type,
918 typename MV::node_type> crs_matrix_type;
919 RCP<const crs_matrix_type> A_crsMat =
920 rcp_dynamic_cast<
const crs_matrix_type> (rcpFromRef (A));
921 if (! A_crsMat.is_null ()) {
922 TEUCHOS_TEST_FOR_EXCEPTION(
923 ! savedDiagOffsets_, std::logic_error,
924 "Ifpack2::Details::Chebyshev::makeInverseDiagonal: " 925 "It is not allowed to call this method with useDiagOffsets=true, " 926 "if you have not previously saved offsets of diagonal entries. " 927 "This situation should never arise if this class is used properly. " 928 "Please report this bug to the Ifpack2 developers.");
929 A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_ ());
935 A.getLocalDiagCopy (*D_rowMap);
937 RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
941 reciprocal_threshold (*D_rangeMap, minDiagVal_);
942 return Teuchos::rcp_const_cast<
const V> (D_rangeMap);
946 template<
class ScalarType,
class MV>
947 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
953 typedef Tpetra::Export<
typename MV::local_ordinal_type,
954 typename MV::global_ordinal_type,
955 typename MV::node_type> export_type;
959 TEUCHOS_TEST_FOR_EXCEPTION(
960 A_.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::" 961 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() " 962 "with a nonnull input matrix before calling this method. This is probably " 963 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
964 TEUCHOS_TEST_FOR_EXCEPTION(
965 D.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::" 966 "makeRangeMapVector: The input Vector D is null. This is probably " 967 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
969 RCP<const map_type> sourceMap = D->getMap ();
970 RCP<const map_type> rangeMap = A_->getRangeMap ();
971 RCP<const map_type> rowMap = A_->getRowMap ();
973 if (rangeMap->isSameAs (*sourceMap)) {
979 RCP<const export_type> exporter;
983 if (sourceMap->isSameAs (*rowMap)) {
985 exporter = A_->getGraph ()->getExporter ();
988 exporter = rcp (
new export_type (sourceMap, rangeMap));
991 if (exporter.is_null ()) {
995 RCP<V> D_out = rcp (
new V (*D, Teuchos::Copy));
996 D_out->doExport (*D, *exporter, Tpetra::ADD);
997 return Teuchos::rcp_const_cast<
const V> (D_out);
1003 template<
class ScalarType,
class MV>
1004 Teuchos::RCP<typename Chebyshev<ScalarType, MV>::V>
1008 using Teuchos::rcp_const_cast;
1009 return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1013 template<
class ScalarType,
class MV>
1023 const V& D_inv)
const 1026 const ST myLambdaMin = lambdaMax / eigRatio;
1028 const ST zero = Teuchos::as<ST> (0);
1029 const ST one = Teuchos::as<ST> (1);
1030 const ST two = Teuchos::as<ST> (2);
1031 const ST d = (lambdaMax + myLambdaMin) / two;
1032 const ST c = (lambdaMax - myLambdaMin) / two;
1034 if (zeroStartingSolution_ && numIters > 0) {
1038 MV R (B.getMap (), B.getNumVectors (),
false);
1039 MV P (B.getMap (), B.getNumVectors (),
false);
1040 MV Z (B.getMap (), B.getNumVectors (),
false);
1042 for (
int i = 0; i < numIters; ++i) {
1043 computeResidual (R, B, A, X);
1044 solve (Z, D_inv, R);
1052 beta = alpha * (c/two) * (c/two);
1053 alpha = one / (d - beta);
1054 P.update (one, Z, beta);
1056 X.update (alpha, P, one);
1063 template<
class ScalarType,
class MV>
1066 Teuchos::Array<MT> norms (X.getNumVectors ());
1067 X.normInf (norms());
1068 return *std::max_element (norms.begin (), norms.end ());
1071 template<
class ScalarType,
class MV>
1083 #ifdef HAVE_TEUCHOS_DEBUG 1084 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1087 cerr <<
"\\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1088 cerr <<
"\\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1089 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1090 #endif // HAVE_TEUCHOS_DEBUG 1092 if (numIters <= 0) {
1095 const ST one = Teuchos::as<ST> (1);
1096 const ST two = Teuchos::as<ST> (2);
1099 if (lambdaMin == one && lambdaMax == lambdaMin) {
1100 solve (X, D_inv, B);
1105 const ST alpha = lambdaMax / eigRatio;
1106 const ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
1107 const ST delta = two / (beta - alpha);
1108 const ST theta = (beta + alpha) / two;
1109 const ST s1 = theta * delta;
1111 #ifdef HAVE_TEUCHOS_DEBUG 1112 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1113 cerr <<
"alpha = " << alpha << endl
1114 <<
"beta = " << beta << endl
1115 <<
"delta = " << delta << endl
1116 <<
"theta = " << theta << endl
1117 <<
"s1 = " << s1 << endl;
1118 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1119 #endif // HAVE_TEUCHOS_DEBUG 1122 Teuchos::RCP<MV> V_ptr, W_ptr;
1123 makeTempMultiVectors (V_ptr, W_ptr, B);
1132 #ifdef HAVE_TEUCHOS_DEBUG 1133 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1134 cerr <<
"Iteration " << 1 <<
":" << endl
1135 <<
"- \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1136 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1137 #endif // HAVE_TEUCHOS_DEBUG 1140 if (! zeroStartingSolution_) {
1141 computeResidual (V1, B, A, X);
1143 #ifdef HAVE_TEUCHOS_DEBUG 1144 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1145 cerr <<
"- \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1146 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1147 #endif // HAVE_TEUCHOS_DEBUG 1149 solve (W, one/theta, D_inv, V1);
1151 #ifdef HAVE_TEUCHOS_DEBUG 1152 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1153 cerr <<
"- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1154 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1155 #endif // HAVE_TEUCHOS_DEBUG 1157 X.update (one, W, one);
1159 solve (W, one/theta, D_inv, B);
1161 #ifdef HAVE_TEUCHOS_DEBUG 1162 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1163 cerr <<
"- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1164 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1165 #endif // HAVE_TEUCHOS_DEBUG 1167 Tpetra::deep_copy(X, W);
1169 #ifdef HAVE_TEUCHOS_DEBUG 1170 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1171 cerr <<
"- \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1172 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1173 #endif // HAVE_TEUCHOS_DEBUG 1177 ST rhokp1, dtemp1, dtemp2;
1178 for (
int deg = 1; deg < numIters; ++deg) {
1180 #ifdef HAVE_TEUCHOS_DEBUG 1181 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1182 cerr <<
"Iteration " << deg+1 <<
":" << endl;
1183 cerr <<
"- \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1184 cerr <<
"- \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1185 cerr <<
"- \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm () << endl;
1186 cerr <<
"- rhok = " << rhok << endl;
1187 V1.putScalar (STS::zero ());
1188 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1189 #endif // HAVE_TEUCHOS_DEBUG 1191 computeResidual (V1, B, A, X);
1193 #ifdef HAVE_TEUCHOS_DEBUG 1194 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1195 cerr <<
"- \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1196 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1197 #endif // HAVE_TEUCHOS_DEBUG 1199 rhokp1 = one / (two * s1 - rhok);
1200 dtemp1 = rhokp1 * rhok;
1201 dtemp2 = two * rhokp1 * delta;
1204 #ifdef HAVE_TEUCHOS_DEBUG 1205 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1206 cerr <<
"- dtemp1 = " << dtemp1 << endl
1207 <<
"- dtemp2 = " << dtemp2 << endl;
1208 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1209 #endif // HAVE_TEUCHOS_DEBUG 1212 W.elementWiseMultiply (dtemp2, D_inv, V1, one);
1213 X.update (one, W, one);
1215 #ifdef HAVE_TEUCHOS_DEBUG 1216 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG 1217 cerr <<
"- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1218 cerr <<
"- \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1219 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG 1220 #endif // HAVE_TEUCHOS_DEBUG 1224 template<
class ScalarType,
class MV>
1232 const ST zero = Teuchos::as<ST> (0);
1233 const ST one = Teuchos::as<ST> (1);
1234 ST lambdaMax = zero;
1235 ST RQ_top, RQ_bottom, norm;
1237 V y (A.getRangeMap ());
1239 TEUCHOS_TEST_FOR_EXCEPTION(
1240 norm == zero, std::runtime_error,
1241 "Ifpack2::Chebyshev::powerMethodWithInitGuess: " 1242 "The initial guess has zero norm. This could be either because Tpetra::" 1243 "Vector::randomize() filled the vector with zeros (if that was used to " 1244 "compute the initial guess), or because the norm2 method has a bug. The " 1245 "first is not impossible, but unlikely.");
1248 x.scale (one / norm);
1250 for (
int iter = 0; iter < numIters; ++iter) {
1253 solve (y, D_inv, y);
1256 RQ_bottom = x.dot (x);
1258 lambdaMax = RQ_top / RQ_bottom;
1263 x.update (one / norm, y, zero);
1268 template<
class ScalarType,
class MV>
1273 const ST zero = Teuchos::as<ST> (0);
1274 V x (A.getDomainMap ());
1277 ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1285 if (STS::real (lambdaMax) < STS::real (zero)) {
1291 lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1296 template<
class ScalarType,
class MV>
1297 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::row_matrix_type>
1302 template<
class ScalarType,
class MV>
1310 template<
class ScalarType,
class MV>
1314 Teuchos::RCP<MV>& W,
1317 if (V_.is_null ()) {
1318 V_ = Teuchos::rcp (
new MV (X.getMap (), X.getNumVectors (),
false));
1321 if (W_.is_null ()) {
1322 W_ = Teuchos::rcp (
new MV (X.getMap (), X.getNumVectors (),
true));
1328 template<
class ScalarType,
class MV>
1332 std::ostringstream oss;
1335 oss <<
"\"Ifpack2::Details::Chebyshev\":" 1337 <<
"degree: " << numIters_
1338 <<
", lambdaMax: " << lambdaMaxForApply_
1339 <<
", alpha: " << eigRatioForApply_
1340 <<
", lambdaMin: " << lambdaMinForApply_
1345 template<
class ScalarType,
class MV>
1349 const Teuchos::EVerbosityLevel verbLevel)
const 1351 using Teuchos::TypeNameTraits;
1354 const Teuchos::EVerbosityLevel vl =
1355 (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1356 if (vl > Teuchos::VERB_NONE) {
1357 if (vl == Teuchos::VERB_LOW) {
1362 out <<
"\"Ifpack2::Details::Chebyshev\":" << endl;
1363 Teuchos::OSTab tab1 (Teuchos::rcpFromRef (out));
1364 out <<
"Template parameters:" << endl;
1366 Teuchos::OSTab tab2 (Teuchos::rcpFromRef (out));
1367 out <<
"ScalarType: \"" << TypeNameTraits<ScalarType>::name () <<
"\"" << endl
1368 <<
"MV: \"" << TypeNameTraits<MV>::name () <<
"\"" << endl;
1372 out << endl <<
"Computed parameters:" << endl;
1374 Teuchos::OSTab tab2 (Teuchos::rcpFromRef (out));
1378 if (D_.is_null ()) {
1379 out <<
"unset" << endl;
1380 }
else if (vl <= Teuchos::VERB_HIGH) {
1381 out <<
"set" << endl;
1385 Teuchos::OSTab tab3 (Teuchos::rcpFromRef (out));
1386 D_->describe (out, vl);
1392 out <<
"V_: " << (V_.is_null () ?
"unset" :
"set") << endl
1393 <<
"W_: " << (W_.is_null () ?
"unset" :
"set") << endl
1394 <<
"computedLambdaMax_: " << computedLambdaMax_ << endl
1395 <<
"computedLambdaMin_: " << computedLambdaMin_ << endl
1396 <<
"lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1397 <<
"lambdaMinForApply_: " << lambdaMinForApply_ << endl
1398 <<
"eigRatioForApply_: " << eigRatioForApply_ << endl;
1400 out <<
"User parameters:" << endl;
1402 Teuchos::OSTab tab2 (Teuchos::rcpFromRef (out));
1403 out <<
"userInvDiag_: ";
1404 if (userInvDiag_.is_null ()) {
1405 out <<
"unset" << endl;
1406 }
else if (vl <= Teuchos::VERB_HIGH) {
1407 out <<
"set" << endl;
1411 Teuchos::OSTab tab3 (Teuchos::rcpFromRef (out));
1412 userInvDiag_->describe (out, vl);
1416 out <<
"userLambdaMax_: " << userLambdaMax_ << endl
1417 <<
"userLambdaMin_: " << userLambdaMin_ << endl
1418 <<
"userEigRatio_: " << userEigRatio_ << endl
1419 <<
"numIters_: " << numIters_ << endl
1420 <<
"eigMaxIters_: " << eigMaxIters_ << endl
1421 <<
"zeroStartingSolution_: " << zeroStartingSolution_ << endl
1422 <<
"assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1423 <<
"textbookAlgorithm_: " << textbookAlgorithm_ << endl
1424 <<
"computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1434 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \ 1435 template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >; 1437 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:316
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:825
Tpetra::Operator< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > op_type
Specialization of Tpetra::Operator.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:142
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:1348
void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set the matrix.
Definition: Ifpack2_Details_Chebyshev_def.hpp:652
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the matrix given to the constructor.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1298
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:266
Ifpack2 implementation details.
void print(std::ostream &out)
Print instance data to the given output stream.
Definition: Ifpack2_Details_Chebyshev_def.hpp:872
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:137
ScalarType ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:133
Tpetra::RowMatrix< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > row_matrix_type
Specialization of Tpetra::RowMatrix.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:147
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:127
Declaration of Chebyshev implementation.
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:63
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1331
bool hasTransposeApply() const
Whether it's possible to apply the transpose of this operator.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1305
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
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:665