44 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
45 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
58 #include "Ifpack2_Details_ChebyshevKernel.hpp"
59 #include "Kokkos_ArithTraits.hpp"
60 #include "Teuchos_FancyOStream.hpp"
61 #include "Teuchos_oblackholestream.hpp"
62 #include "Tpetra_Details_residual.hpp"
64 #include "Ifpack2_Details_LapackSupportsScalar.hpp"
74 const char computeBeforeApplyReminder[] =
75 "This means one of the following:\n"
76 " - you have not yet called compute() on this instance, or \n"
77 " - you didn't call compute() after calling setParameters().\n\n"
78 "After creating an Ifpack2::Chebyshev instance,\n"
79 "you must _always_ call compute() at least once before calling apply().\n"
80 "After calling compute() once, you do not need to call it again,\n"
81 "unless the matrix has changed or you have changed parameters\n"
82 "(by calling setParameters()).";
88 template<
class XV,
class SizeType =
typename XV::
size_type>
89 struct V_ReciprocalThresholdSelfFunctor
91 typedef typename XV::execution_space execution_space;
92 typedef typename XV::non_const_value_type value_type;
93 typedef SizeType size_type;
94 typedef Kokkos::ArithTraits<value_type> KAT;
95 typedef typename KAT::mag_type mag_type;
98 const value_type minVal_;
99 const mag_type minValMag_;
101 V_ReciprocalThresholdSelfFunctor (
const XV& X,
102 const value_type& min_val) :
105 minValMag_ (KAT::abs (min_val))
108 KOKKOS_INLINE_FUNCTION
109 void operator() (
const size_type& i)
const
111 const mag_type X_i_abs = KAT::abs (X_[i]);
113 if (X_i_abs < minValMag_) {
117 X_[i] = KAT::one () / X_[i];
122 template<
class XV,
class SizeType =
typename XV::
size_type>
123 struct LocalReciprocalThreshold {
125 compute (
const XV& X,
126 const typename XV::non_const_value_type& minVal)
128 typedef typename XV::execution_space execution_space;
129 Kokkos::RangePolicy<execution_space, SizeType> policy (0, X.extent (0));
130 V_ReciprocalThresholdSelfFunctor<XV, SizeType> op (X, minVal);
131 Kokkos::parallel_for (policy, op);
135 template <
class TpetraVectorType,
136 const bool classic = TpetraVectorType::node_type::classic>
137 struct GlobalReciprocalThreshold {};
139 template <
class TpetraVectorType>
140 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
142 compute (TpetraVectorType& V,
143 const typename TpetraVectorType::scalar_type& min_val)
145 typedef typename TpetraVectorType::scalar_type scalar_type;
146 typedef typename TpetraVectorType::mag_type mag_type;
147 typedef Kokkos::ArithTraits<scalar_type> STS;
149 const scalar_type ONE = STS::one ();
150 const mag_type min_val_abs = STS::abs (min_val);
153 const size_t lclNumRows = V.getLocalLength ();
155 for (
size_t i = 0; i < lclNumRows; ++i) {
156 const scalar_type V_0i = V_0[i];
157 if (STS::abs (V_0i) < min_val_abs) {
166 template <
class TpetraVectorType>
167 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
169 compute (TpetraVectorType& X,
170 const typename TpetraVectorType::scalar_type& minVal)
172 typedef typename TpetraVectorType::impl_scalar_type value_type;
174 const value_type minValS =
static_cast<value_type
> (minVal);
175 auto X_0 = Kokkos::subview (X.getLocalViewDevice (Tpetra::Access::ReadWrite),
177 LocalReciprocalThreshold<decltype (X_0) >::compute (X_0, minValS);
182 template <
typename S,
typename L,
typename G,
typename N>
184 reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V,
const S& minVal)
186 GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
190 template<class ScalarType, const bool lapackSupportsScalarType = LapackSupportsScalar<ScalarType>::value>
196 throw std::runtime_error(
"LAPACK does not support the scalar type.");
203 computeInitialGuessForCG (
const V& diagonal, V& x) {
204 using device_type =
typename V::node_type::device_type;
205 using range_policy = Kokkos::RangePolicy<typename device_type::execution_space>;
213 size_t N = x.getLocalLength();
214 auto d_view = diagonal.template getLocalView<device_type>(Tpetra::Access::ReadOnly);
215 auto x_view = x.template getLocalView<device_type>(Tpetra::Access::ReadWrite);
220 Kokkos::parallel_for(
"computeInitialGuessforCG::zero_bcs", range_policy(0,N), KOKKOS_LAMBDA(
const size_t & i) {
221 if(d_view(i,0) == ONE)
230 template<
class ScalarType>
231 struct LapackHelper<ScalarType,true> {
237 using MagnitudeType =
typename STS::magnitudeType;
239 const int N = diag.size ();
240 ScalarType scalar_dummy;
241 std::vector<MagnitudeType> mag_dummy(4*N);
245 ScalarType lambdaMax = STS::one();
248 lapack.
PTEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
249 &scalar_dummy, 1, &mag_dummy[0], &info);
251 (info < 0, std::logic_error,
"Ifpack2::Details::LapackHelper::tri_diag_spectral_radius:"
252 "LAPACK's _PTEQR failed with info = "
253 << info <<
" < 0. This suggests there might be a bug in the way Ifpack2 "
254 "is calling LAPACK. Please report this to the Ifpack2 developers.");
256 lambdaMax = Teuchos::as<ScalarType> (diag[0]);
263 template<
class ScalarType,
class MV>
264 void Chebyshev<ScalarType, MV>::checkInputMatrix ()
const
267 ! A_.
is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
268 std::invalid_argument,
269 "Ifpack2::Chebyshev: The input matrix A must be square. "
270 "A has " << A_->getGlobalNumRows () <<
" rows and "
271 << A_->getGlobalNumCols () <<
" columns.");
275 if (debug_ && ! A_.
is_null ()) {
283 ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
284 "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be "
285 "the same (in the sense of isSameAs())." << std::endl <<
"We only check "
286 "for this in debug mode.");
291 template<
class ScalarType,
class MV>
293 Chebyshev<ScalarType, MV>::
294 checkConstructorInput ()
const
299 if (STS::isComplex) {
301 (
true, std::logic_error,
"Ifpack2::Chebyshev: This class' implementation "
302 "of Chebyshev iteration only works for real-valued, symmetric positive "
303 "definite matrices. However, you instantiated this class for ScalarType"
305 "complex-valued type. While this may be algorithmically correct if all "
306 "of the complex numbers in the matrix have zero imaginary part, we "
307 "forbid using complex ScalarType altogether in order to remind you of "
308 "the limitations of our implementation (and of the algorithm itself).");
315 template<
class ScalarType,
class MV>
319 savedDiagOffsets_ (false),
320 computedLambdaMax_ (STS::nan ()),
321 computedLambdaMin_ (STS::nan ()),
322 lambdaMaxForApply_ (STS::nan ()),
323 lambdaMinForApply_ (STS::nan ()),
324 eigRatioForApply_ (STS::nan ()),
325 userLambdaMax_ (STS::nan ()),
326 userLambdaMin_ (STS::nan ()),
327 userEigRatio_ (Teuchos::
as<
ST> (30)),
328 minDiagVal_ (STS::eps ()),
331 eigRelTolerance_(Teuchos::ScalarTraits<
MT>::zero ()),
332 eigKeepVectors_(false),
333 eigenAnalysisType_(
"power method"),
334 eigNormalizationFreq_(1),
335 zeroStartingSolution_ (true),
336 assumeMatrixUnchanged_ (false),
337 chebyshevAlgorithm_(
"first"),
338 computeMaxResNorm_ (false),
339 computeSpectralRadius_(true),
340 ckUseNativeSpMV_(false),
343 checkConstructorInput ();
346 template<
class ScalarType,
class MV>
351 savedDiagOffsets_ (false),
352 computedLambdaMax_ (STS::nan ()),
353 computedLambdaMin_ (STS::nan ()),
354 lambdaMaxForApply_ (STS::nan ()),
355 boostFactor_ (static_cast<
MT> (1.1)),
356 lambdaMinForApply_ (STS::nan ()),
357 eigRatioForApply_ (STS::nan ()),
358 userLambdaMax_ (STS::nan ()),
359 userLambdaMin_ (STS::nan ()),
360 userEigRatio_ (Teuchos::
as<
ST> (30)),
361 minDiagVal_ (STS::eps ()),
364 eigRelTolerance_(Teuchos::ScalarTraits<
MT>::zero ()),
365 eigKeepVectors_(false),
366 eigenAnalysisType_(
"power method"),
367 eigNormalizationFreq_(1),
368 zeroStartingSolution_ (true),
369 assumeMatrixUnchanged_ (false),
370 chebyshevAlgorithm_(
"first"),
371 computeMaxResNorm_ (false),
372 computeSpectralRadius_(true),
373 ckUseNativeSpMV_(false),
376 checkConstructorInput ();
380 template<
class ScalarType,
class MV>
387 using Teuchos::rcp_const_cast;
399 const ST defaultLambdaMax = STS::nan ();
400 const ST defaultLambdaMin = STS::nan ();
407 const ST defaultEigRatio = Teuchos::as<ST> (30);
408 const MT defaultBoostFactor =
static_cast<MT> (1.1);
409 const ST defaultMinDiagVal = STS::eps ();
410 const int defaultNumIters = 1;
411 const int defaultEigMaxIters = 10;
413 const bool defaultEigKeepVectors =
false;
414 const int defaultEigNormalizationFreq = 1;
415 const bool defaultZeroStartingSolution =
true;
416 const bool defaultAssumeMatrixUnchanged =
false;
417 const std::string defaultChebyshevAlgorithm =
"first";
418 const bool defaultComputeMaxResNorm =
false;
419 const bool defaultComputeSpectralRadius =
true;
420 const bool defaultCkUseNativeSpMV =
false;
421 const bool defaultDebug =
false;
427 RCP<const V> userInvDiagCopy;
428 ST lambdaMax = defaultLambdaMax;
429 ST lambdaMin = defaultLambdaMin;
430 ST eigRatio = defaultEigRatio;
431 MT boostFactor = defaultBoostFactor;
432 ST minDiagVal = defaultMinDiagVal;
433 int numIters = defaultNumIters;
434 int eigMaxIters = defaultEigMaxIters;
435 MT eigRelTolerance = defaultEigRelTolerance;
436 bool eigKeepVectors = defaultEigKeepVectors;
437 int eigNormalizationFreq = defaultEigNormalizationFreq;
438 bool zeroStartingSolution = defaultZeroStartingSolution;
439 bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
440 std::string chebyshevAlgorithm = defaultChebyshevAlgorithm;
441 bool computeMaxResNorm = defaultComputeMaxResNorm;
442 bool computeSpectralRadius = defaultComputeSpectralRadius;
443 bool ckUseNativeSpMV = defaultCkUseNativeSpMV;
444 bool debug = defaultDebug;
451 if (plist.
isType<
bool> (
"debug")) {
452 debug = plist.
get<
bool> (
"debug");
454 else if (plist.
isType<
int> (
"debug")) {
455 const int debugInt = plist.
get<
bool> (
"debug");
456 debug = debugInt != 0;
469 const char opInvDiagLabel[] =
"chebyshev: operator inv diagonal";
472 RCP<const V> userInvDiag;
474 if (plist.
isType<
const V*> (opInvDiagLabel)) {
475 const V* rawUserInvDiag =
476 plist.
get<
const V*> (opInvDiagLabel);
478 userInvDiag =
rcp (rawUserInvDiag,
false);
480 else if (plist.
isType<
const V*> (opInvDiagLabel)) {
481 V* rawUserInvDiag = plist.
get<V*> (opInvDiagLabel);
483 userInvDiag =
rcp (const_cast<const V*> (rawUserInvDiag),
false);
485 else if (plist.
isType<RCP<const V>> (opInvDiagLabel)) {
486 userInvDiag = plist.
get<RCP<const V> > (opInvDiagLabel);
488 else if (plist.
isType<RCP<V>> (opInvDiagLabel)) {
489 RCP<V> userInvDiagNonConst =
490 plist.
get<RCP<V> > (opInvDiagLabel);
491 userInvDiag = rcp_const_cast<
const V> (userInvDiagNonConst);
493 else if (plist.
isType<
const V> (opInvDiagLabel)) {
494 const V& userInvDiagRef = plist.
get<
const V> (opInvDiagLabel);
496 userInvDiag = userInvDiagCopy;
498 else if (plist.
isType<V> (opInvDiagLabel)) {
499 V& userInvDiagNonConstRef = plist.
get<V> (opInvDiagLabel);
500 const V& userInvDiagRef =
const_cast<const V&
> (userInvDiagNonConstRef);
502 userInvDiag = userInvDiagCopy;
512 if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
523 if (plist.
isParameter (
"chebyshev: use native spmv"))
524 ckUseNativeSpMV = plist.
get(
"chebyshev: use native spmv", ckUseNativeSpMV);
529 if (plist.
isParameter (
"chebyshev: max eigenvalue")) {
530 if (plist.
isType<
double>(
"chebyshev: max eigenvalue"))
531 lambdaMax = plist.
get<
double> (
"chebyshev: max eigenvalue");
533 lambdaMax = plist.
get<
ST> (
"chebyshev: max eigenvalue");
535 STS::isnaninf (lambdaMax), std::invalid_argument,
536 "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" "
537 "parameter is NaN or Inf. This parameter is optional, but if you "
538 "choose to supply it, it must have a finite value.");
540 if (plist.
isParameter (
"chebyshev: min eigenvalue")) {
541 if (plist.
isType<
double>(
"chebyshev: min eigenvalue"))
542 lambdaMin = plist.
get<
double> (
"chebyshev: min eigenvalue");
544 lambdaMin = plist.
get<
ST> (
"chebyshev: min eigenvalue");
546 STS::isnaninf (lambdaMin), std::invalid_argument,
547 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" "
548 "parameter is NaN or Inf. This parameter is optional, but if you "
549 "choose to supply it, it must have a finite value.");
553 if (plist.
isParameter (
"smoother: Chebyshev alpha")) {
554 if (plist.
isType<
double>(
"smoother: Chebyshev alpha"))
555 eigRatio = plist.
get<
double> (
"smoother: Chebyshev alpha");
557 eigRatio = plist.
get<
ST> (
"smoother: Chebyshev alpha");
560 eigRatio = plist.
get (
"chebyshev: ratio eigenvalue", eigRatio);
562 STS::isnaninf (eigRatio), std::invalid_argument,
563 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" "
564 "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. "
565 "This parameter is optional, but if you choose to supply it, it must have "
573 STS::real (eigRatio) < STS::real (STS::one ()),
574 std::invalid_argument,
575 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\""
576 "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, "
577 "but you supplied the value " << eigRatio <<
".");
582 const char paramName[] =
"chebyshev: boost factor";
586 boostFactor = plist.
get<
MT> (paramName);
588 else if (! std::is_same<double, MT>::value &&
589 plist.
isType<
double> (paramName)) {
590 const double dblBF = plist.
get<
double> (paramName);
591 boostFactor =
static_cast<MT> (dblBF);
595 (
true, std::invalid_argument,
596 "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\""
597 "parameter must have type magnitude_type (MT) or double.");
607 plist.
set (paramName, defaultBoostFactor);
611 "Ifpack2::Chebyshev::setParameters: \"" << paramName <<
"\" parameter "
612 "must be >= 1, but you supplied the value " << boostFactor <<
".");
616 minDiagVal = plist.
get (
"chebyshev: min diagonal value", minDiagVal);
618 STS::isnaninf (minDiagVal), std::invalid_argument,
619 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" "
620 "parameter is NaN or Inf. This parameter is optional, but if you choose "
621 "to supply it, it must have a finite value.");
625 numIters = plist.
get<
int> (
"smoother: sweeps");
628 numIters = plist.
get<
int> (
"relaxation: sweeps");
630 numIters = plist.
get (
"chebyshev: degree", numIters);
632 numIters < 0, std::invalid_argument,
633 "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also "
634 "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a "
635 "nonnegative integer. You gave a value of " << numIters <<
".");
638 if (plist.
isParameter (
"eigen-analysis: iterations")) {
639 eigMaxIters = plist.
get<
int> (
"eigen-analysis: iterations");
641 eigMaxIters = plist.
get (
"chebyshev: eigenvalue max iterations", eigMaxIters);
643 eigMaxIters < 0, std::invalid_argument,
644 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations"
645 "\" parameter (also called \"eigen-analysis: iterations\") must be a "
646 "nonnegative integer. You gave a value of " << eigMaxIters <<
".");
648 if (plist.
isType<
double>(
"chebyshev: eigenvalue relative tolerance"))
649 eigRelTolerance = Teuchos::as<MT>(plist.
get<
double> (
"chebyshev: eigenvalue relative tolerance"));
650 else if (plist.
isType<
MT>(
"chebyshev: eigenvalue relative tolerance"))
651 eigRelTolerance = plist.
get<
MT> (
"chebyshev: eigenvalue relative tolerance");
652 else if (plist.
isType<
ST>(
"chebyshev: eigenvalue relative tolerance"))
655 eigKeepVectors = plist.
get (
"chebyshev: eigenvalue keep vectors", eigKeepVectors);
657 eigNormalizationFreq = plist.
get (
"chebyshev: eigenvalue normalization frequency", eigNormalizationFreq);
659 eigNormalizationFreq < 0, std::invalid_argument,
660 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue normalization frequency"
661 "\" parameter must be a "
662 "nonnegative integer. You gave a value of " << eigNormalizationFreq <<
".")
664 zeroStartingSolution = plist.
get (
"chebyshev: zero starting solution",
665 zeroStartingSolution);
666 assumeMatrixUnchanged = plist.
get (
"chebyshev: assume matrix does not change",
667 assumeMatrixUnchanged);
672 chebyshevAlgorithm = plist.
get<std::string> (
"chebyshev: algorithm");
674 chebyshevAlgorithm !=
"first" &&
675 chebyshevAlgorithm !=
"textbook" &&
676 chebyshevAlgorithm !=
"fourth" &&
677 chebyshevAlgorithm !=
"opt_fourth",
678 std::invalid_argument,
679 "Ifpack2::Chebyshev: Ifpack2 only supports \"first\", \"textbook\", \"fourth\", and \"opt_fourth\", for \"chebyshev: algorithm\".");
682 #ifdef IFPACK2_ENABLE_DEPRECATED_CODE
686 if (plist.
isParameter (
"chebyshev: textbook algorithm")) {
687 const bool textbookAlgorithm = plist.
get<
bool> (
"chebyshev: textbook algorithm");
688 if(textbookAlgorithm){
689 chebyshevAlgorithm =
"textbook";
691 chebyshevAlgorithm =
"first";
697 if (plist.
isParameter (
"chebyshev: compute max residual norm")) {
698 computeMaxResNorm = plist.
get<
bool> (
"chebyshev: compute max residual norm");
700 if (plist.
isParameter (
"chebyshev: compute spectral radius")) {
701 computeSpectralRadius = plist.
get<
bool> (
"chebyshev: compute spectral radius");
710 (plist.
isType<
bool> (
"chebyshev: use block mode") &&
711 ! plist.
get<
bool> (
"chebyshev: use block mode"),
712 std::invalid_argument,
713 "Ifpack2::Chebyshev requires that if you set \"chebyshev: use "
714 "block mode\" at all, you must set it to false. "
715 "Ifpack2::Chebyshev does not implement Ifpack's block mode.");
717 (plist.
isType<
bool> (
"chebyshev: solve normal equations") &&
718 ! plist.
get<
bool> (
"chebyshev: solve normal equations"),
719 std::invalid_argument,
720 "Ifpack2::Chebyshev does not and will never implement the Ifpack "
721 "parameter \"chebyshev: solve normal equations\". If you want to "
722 "solve the normal equations, construct a Tpetra::Operator that "
723 "implements A^* A, and use Chebyshev to solve A^* A x = A^* b.");
731 std::string eigenAnalysisType (
"power-method");
733 eigenAnalysisType = plist.
get<std::string> (
"eigen-analysis: type");
735 eigenAnalysisType !=
"power-method" &&
736 eigenAnalysisType !=
"power method" &&
737 eigenAnalysisType !=
"cg",
738 std::invalid_argument,
739 "Ifpack2::Chebyshev: Ifpack2 only supports \"power method\" and \"cg\" for \"eigen-analysis: type\".");
743 userInvDiag_ = userInvDiagCopy;
744 userLambdaMax_ = lambdaMax;
745 userLambdaMin_ = lambdaMin;
746 userEigRatio_ = eigRatio;
747 boostFactor_ =
static_cast<MT> (boostFactor);
748 minDiagVal_ = minDiagVal;
749 numIters_ = numIters;
750 eigMaxIters_ = eigMaxIters;
751 eigRelTolerance_ = eigRelTolerance;
752 eigKeepVectors_ = eigKeepVectors;
753 eigNormalizationFreq_ = eigNormalizationFreq;
754 eigenAnalysisType_ = eigenAnalysisType;
755 zeroStartingSolution_ = zeroStartingSolution;
756 assumeMatrixUnchanged_ = assumeMatrixUnchanged;
757 chebyshevAlgorithm_ = chebyshevAlgorithm;
758 computeMaxResNorm_ = computeMaxResNorm;
759 computeSpectralRadius_ = computeSpectralRadius;
760 ckUseNativeSpMV_ = ckUseNativeSpMV;
772 myRank = A_->getComm ()->getRank ();
776 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
781 out_ = Teuchos::getFancyOStream (blackHole);
786 out_ = Teuchos::null;
791 template<
class ScalarType,
class MV>
797 diagOffsets_ = offsets_type ();
798 savedDiagOffsets_ =
false;
800 computedLambdaMax_ = STS::nan ();
801 computedLambdaMin_ = STS::nan ();
802 eigVector_ = Teuchos::null;
803 eigVector2_ = Teuchos::null;
807 template<
class ScalarType,
class MV>
813 if (! assumeMatrixUnchanged_) {
831 myRank = A->getComm ()->getRank ();
835 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
839 out_ = Teuchos::getFancyOStream (blackHole);
844 out_ = Teuchos::null;
849 template<
class ScalarType,
class MV>
858 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
859 typename MV::local_ordinal_type,
860 typename MV::global_ordinal_type,
861 typename MV::node_type> crs_matrix_type;
864 A_.
is_null (), std::runtime_error,
"Ifpack2::Chebyshev::compute: The input "
865 "matrix A is null. Please call setMatrix() with a nonnull input matrix "
866 "before calling this method.");
881 if (userInvDiag_.is_null ()) {
883 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A_);
885 if (! A_crsMat.
is_null () && A_crsMat->isFillComplete ()) {
887 const size_t lclNumRows = A_crsMat->getLocalNumRows ();
888 if (diagOffsets_.extent (0) < lclNumRows) {
889 diagOffsets_ = offsets_type ();
890 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
892 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
893 savedDiagOffsets_ =
true;
894 D_ = makeInverseDiagonal (*A_,
true);
897 D_ = makeInverseDiagonal (*A_);
900 else if (! assumeMatrixUnchanged_) {
901 if (! A_crsMat.
is_null () && A_crsMat->isFillComplete ()) {
904 if (! savedDiagOffsets_) {
905 const size_t lclNumRows = A_crsMat->getLocalNumRows ();
906 if (diagOffsets_.extent (0) < lclNumRows) {
907 diagOffsets_ = offsets_type ();
908 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
910 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
911 savedDiagOffsets_ =
true;
914 D_ = makeInverseDiagonal (*A_,
true);
917 D_ = makeInverseDiagonal (*A_);
922 D_ = makeRangeMapVectorConst (userInvDiag_);
926 const bool computedEigenvalueEstimates =
927 STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
937 if (! assumeMatrixUnchanged_ ||
938 (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
939 ST computedLambdaMax;
940 if ((eigenAnalysisType_ ==
"power method") || (eigenAnalysisType_ ==
"power-method")) {
942 if (eigVector_.is_null()) {
946 PowerMethod::computeInitialGuessForPowerMethod (*x,
false);
951 if (eigVector2_.is_null()) {
952 y =
rcp(
new V(A_->getRangeMap ()));
959 computedLambdaMax = PowerMethod::powerMethodWithInitGuess (*A_, *D_, eigMaxIters_, x, y,
960 eigRelTolerance_, eigNormalizationFreq_, stream,
961 computeSpectralRadius_);
964 computedLambdaMax = cgMethod (*A_, *D_, eigMaxIters_);
967 STS::isnaninf (computedLambdaMax),
969 "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue "
970 "of D^{-1} A failed, by producing Inf or NaN. This probably means that "
971 "the matrix contains Inf or NaN values, or that it is badly scaled.");
973 STS::isnaninf (userEigRatio_),
975 "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
976 << endl <<
"This should be impossible." << endl <<
977 "Please report this bug to the Ifpack2 developers.");
983 const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
986 computedLambdaMax_ = computedLambdaMax;
987 computedLambdaMin_ = computedLambdaMin;
990 STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
992 "Ifpack2::Chebyshev::compute: " << endl <<
993 "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
994 << endl <<
"This should be impossible." << endl <<
995 "Please report this bug to the Ifpack2 developers.");
1003 lambdaMaxForApply_ = STS::isnaninf (userLambdaMax_) ? computedLambdaMax_ : userLambdaMax_;
1016 lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
1017 eigRatioForApply_ = userEigRatio_;
1019 if (chebyshevAlgorithm_ ==
"first") {
1022 const ST one = Teuchos::as<ST> (1);
1025 if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
1026 lambdaMinForApply_ = one;
1027 lambdaMaxForApply_ = lambdaMinForApply_;
1028 eigRatioForApply_ = one;
1034 template<
class ScalarType,
class MV>
1038 return lambdaMaxForApply_;
1042 template<
class ScalarType,
class MV>
1046 const char prefix[] =
"Ifpack2::Chebyshev::apply: ";
1049 *out_ <<
"apply: " << std::endl;
1052 (A_.
is_null (), std::runtime_error, prefix <<
"The input matrix A is null. "
1053 " Please call setMatrix() with a nonnull input matrix, and then call "
1054 "compute(), before calling this method.");
1056 (STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
1057 prefix <<
"There is no estimate for the max eigenvalue."
1058 << std::endl << computeBeforeApplyReminder);
1059 TEUCHOS_TEST_FOR_EXCEPTION
1060 (STS::isnaninf (lambdaMinForApply_), std::runtime_error,
1061 prefix <<
"There is no estimate for the min eigenvalue."
1062 << std::endl << computeBeforeApplyReminder);
1063 TEUCHOS_TEST_FOR_EXCEPTION
1064 (STS::isnaninf (eigRatioForApply_), std::runtime_error,
1065 prefix <<
"There is no estimate for the ratio of the max "
1066 "eigenvalue to the min eigenvalue."
1067 << std::endl << computeBeforeApplyReminder);
1068 TEUCHOS_TEST_FOR_EXCEPTION
1069 (D_.is_null (), std::runtime_error, prefix <<
"The vector of inverse "
1070 "diagonal entries of the matrix has not yet been computed."
1071 << std::endl << computeBeforeApplyReminder);
1073 if (chebyshevAlgorithm_ ==
"fourth" || chebyshevAlgorithm_ ==
"opt_fourth") {
1074 fourthKindApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_, *D_);
1076 else if (chebyshevAlgorithm_ ==
"textbook") {
1077 textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1078 lambdaMinForApply_, eigRatioForApply_, *D_);
1081 ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1082 lambdaMinForApply_, eigRatioForApply_, *D_);
1085 if (computeMaxResNorm_ && B.getNumVectors () > 0) {
1086 MV R (B.getMap (), B.getNumVectors ());
1087 computeResidual (R, B, *A_, X);
1090 return *std::max_element (norms.begin (), norms.end ());
1097 template<
class ScalarType,
class MV>
1102 using Teuchos::rcpFromRef;
1103 this->describe (* (Teuchos::getFancyOStream (rcpFromRef (out))),
1107 template<
class ScalarType,
class MV>
1111 const ScalarType& alpha,
1116 solve (W, alpha, D_inv, B);
1117 Tpetra::deep_copy (X, W);
1120 template<
class ScalarType,
class MV>
1126 Tpetra::Details::residual(A,X,B,R);
1129 template <
class ScalarType,
class MV>
1131 Chebyshev<ScalarType, MV>::
1132 solve (MV& Z,
const V& D_inv,
const MV& R) {
1133 Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
1136 template<
class ScalarType,
class MV>
1138 Chebyshev<ScalarType, MV>::
1139 solve (MV& Z,
const ST alpha,
const V& D_inv,
const MV& R) {
1140 Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
1143 template<
class ScalarType,
class MV>
1145 Chebyshev<ScalarType, MV>::
1146 makeInverseDiagonal (
const row_matrix_type& A,
const bool useDiagOffsets)
const
1149 using Teuchos::rcpFromRef;
1150 using Teuchos::rcp_dynamic_cast;
1153 if (!D_.is_null() &&
1154 D_->getMap()->isSameAs(*(A.getRowMap ()))) {
1156 *out_ <<
"Reusing pre-existing vector for diagonal extraction" << std::endl;
1157 D_rowMap = Teuchos::rcp_const_cast<V>(D_);
1159 D_rowMap =
Teuchos::rcp(
new V (A.getRowMap (),
false));
1161 *out_ <<
"Allocated new vector for diagonal extraction" << std::endl;
1163 if (useDiagOffsets) {
1167 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
1168 typename MV::local_ordinal_type,
1169 typename MV::global_ordinal_type,
1170 typename MV::node_type> crs_matrix_type;
1171 RCP<const crs_matrix_type> A_crsMat =
1172 rcp_dynamic_cast<
const crs_matrix_type> (rcpFromRef (A));
1173 if (! A_crsMat.is_null ()) {
1175 ! savedDiagOffsets_, std::logic_error,
1176 "Ifpack2::Details::Chebyshev::makeInverseDiagonal: "
1177 "It is not allowed to call this method with useDiagOffsets=true, "
1178 "if you have not previously saved offsets of diagonal entries. "
1179 "This situation should never arise if this class is used properly. "
1180 "Please report this bug to the Ifpack2 developers.");
1181 A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_);
1187 A.getLocalDiagCopy (*D_rowMap);
1189 RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
1195 bool foundNonpositiveValue =
false;
1197 auto D_lcl = D_rangeMap->getLocalViewHost (Tpetra::Access::ReadOnly);
1198 auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1200 typedef typename MV::impl_scalar_type IST;
1201 typedef typename MV::local_ordinal_type LO;
1202 typedef Kokkos::ArithTraits<IST> ATS;
1203 typedef Kokkos::ArithTraits<typename ATS::mag_type> STM;
1205 const LO lclNumRows =
static_cast<LO
> (D_rangeMap->getLocalLength ());
1206 for (LO i = 0; i < lclNumRows; ++i) {
1207 if (STS::real (D_lcl_1d(i)) <= STM::zero ()) {
1208 foundNonpositiveValue =
true;
1214 using Teuchos::outArg;
1216 using Teuchos::reduceAll;
1218 const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1219 int gblSuccess = lclSuccess;
1220 if (! D_rangeMap->getMap ().is_null () && D_rangeMap->getMap ()->getComm ().is_null ()) {
1222 reduceAll<int, int> (comm,
REDUCE_MIN, lclSuccess, outArg (gblSuccess));
1224 if (gblSuccess == 1) {
1225 *out_ <<
"makeInverseDiagonal: The matrix's diagonal entries all have "
1226 "positive real part (this is good for Chebyshev)." << std::endl;
1229 *out_ <<
"makeInverseDiagonal: The matrix's diagonal has at least one "
1230 "entry with nonpositive real part, on at least one process in the "
1231 "matrix's communicator. This is bad for Chebyshev." << std::endl;
1237 reciprocal_threshold (*D_rangeMap, minDiagVal_);
1238 return Teuchos::rcp_const_cast<
const V> (D_rangeMap);
1242 template<
class ScalarType,
class MV>
1244 Chebyshev<ScalarType, MV>::
1249 typedef Tpetra::Export<
typename MV::local_ordinal_type,
1250 typename MV::global_ordinal_type,
1251 typename MV::node_type> export_type;
1256 A_.
is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1257 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() "
1258 "with a nonnull input matrix before calling this method. This is probably "
1259 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1261 D.
is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1262 "makeRangeMapVector: The input Vector D is null. This is probably "
1263 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1265 RCP<const map_type> sourceMap = D->getMap ();
1266 RCP<const map_type> rangeMap = A_->getRangeMap ();
1267 RCP<const map_type> rowMap = A_->getRowMap ();
1269 if (rangeMap->isSameAs (*sourceMap)) {
1275 RCP<const export_type> exporter;
1279 if (sourceMap->isSameAs (*rowMap)) {
1281 exporter = A_->getGraph ()->getExporter ();
1284 exporter =
rcp (
new export_type (sourceMap, rangeMap));
1287 if (exporter.is_null ()) {
1292 D_out->doExport (*D, *exporter, Tpetra::ADD);
1293 return Teuchos::rcp_const_cast<
const V> (D_out);
1299 template<
class ScalarType,
class MV>
1301 Chebyshev<ScalarType, MV>::
1304 using Teuchos::rcp_const_cast;
1305 return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1309 template<
class ScalarType,
class MV>
1311 Chebyshev<ScalarType, MV>::
1312 textbookApplyImpl (
const op_type& A,
1319 const V& D_inv)
const
1322 const ST myLambdaMin = lambdaMax / eigRatio;
1324 const ST zero = Teuchos::as<ST> (0);
1325 const ST one = Teuchos::as<ST> (1);
1326 const ST two = Teuchos::as<ST> (2);
1327 const ST d = (lambdaMax + myLambdaMin) / two;
1328 const ST c = (lambdaMax - myLambdaMin) / two;
1330 if (zeroStartingSolution_ && numIters > 0) {
1334 MV R (B.getMap (), B.getNumVectors (),
false);
1335 MV P (B.getMap (), B.getNumVectors (),
false);
1336 MV Z (B.getMap (), B.getNumVectors (),
false);
1338 for (
int i = 0; i < numIters; ++i) {
1339 computeResidual (R, B, A, X);
1340 solve (Z, D_inv, R);
1348 beta = alpha * (c/two) * (c/two);
1349 alpha = one / (d - beta);
1350 P.update (one, Z, beta);
1352 X.update (alpha, P, one);
1359 template<
class ScalarType,
class MV>
1361 Chebyshev<ScalarType, MV>::
1362 fourthKindApplyImpl (
const op_type& A,
1370 std::vector<ScalarType> betas(numIters, 1.0);
1371 if(chebyshevAlgorithm_ ==
"opt_fourth"){
1372 betas = optimalWeightsImpl<ScalarType>(numIters);
1375 const ST invEig = MT(1) / (lambdaMax * boostFactor_);
1387 if (! zeroStartingSolution_) {
1390 Tpetra::deep_copy (X4, X);
1392 if (ck_.is_null ()) {
1394 ck_ =
Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1398 ck_->compute (Z, MT(4.0/3.0) * invEig, const_cast<V&> (D_inv),
1399 const_cast<MV&> (B), X4, STS::zero());
1402 X.update (betas[0], Z, STS::one());
1406 firstIterationWithZeroStartingSolution (Z, MT(4.0/3.0) * invEig, D_inv, B, X4);
1409 X.update (betas[0], Z, STS::zero());
1412 if (numIters > 1 && ck_.is_null ()) {
1414 ck_ =
Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1417 for (
int i = 1; i < numIters; ++i) {
1418 const ST zScale = (2.0 * i - 1.0) / (2.0 * i + 3.0);
1419 const ST rScale = MT((8.0 * i + 4.0) / (2.0 * i + 3.0)) * invEig;
1423 ck_->compute (Z, rScale, const_cast<V&> (D_inv),
1424 const_cast<MV&> (B), (X4), zScale);
1427 X.update (betas[i], Z, STS::one());
1431 template<
class ScalarType,
class MV>
1432 typename Chebyshev<ScalarType, MV>::MT
1433 Chebyshev<ScalarType, MV>::maxNormInf (
const MV& X) {
1435 X.normInf (norms());
1436 return *std::max_element (norms.begin (), norms.end ());
1439 template<
class ScalarType,
class MV>
1441 Chebyshev<ScalarType, MV>::
1442 ifpackApplyImpl (
const op_type& A,
1452 #ifdef HAVE_IFPACK2_DEBUG
1453 const bool debug = debug_;
1455 const bool debug =
false;
1459 *out_ <<
" \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1460 *out_ <<
" \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1463 if (numIters <= 0) {
1466 const ST zero =
static_cast<ST
> (0.0);
1467 const ST one =
static_cast<ST
> (1.0);
1468 const ST two =
static_cast<ST
> (2.0);
1471 if (lambdaMin == one && lambdaMax == lambdaMin) {
1472 solve (X, D_inv, B);
1477 const ST alpha = lambdaMax / eigRatio;
1478 const ST beta = boostFactor_ * lambdaMax;
1479 const ST delta = two / (beta - alpha);
1480 const ST theta = (beta + alpha) / two;
1481 const ST s1 = theta * delta;
1484 *out_ <<
" alpha = " << alpha << endl
1485 <<
" beta = " << beta << endl
1486 <<
" delta = " << delta << endl
1487 <<
" theta = " << theta << endl
1488 <<
" s1 = " << s1 << endl;
1496 *out_ <<
" Iteration " << 1 <<
":" << endl
1497 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1501 if (! zeroStartingSolution_) {
1504 if (ck_.is_null ()) {
1506 ck_ =
Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1510 ck_->compute (W, one/theta, const_cast<V&> (D_inv),
1511 const_cast<MV&> (B), X, zero);
1515 firstIterationWithZeroStartingSolution (W, one/theta, D_inv, B, X);
1519 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1520 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1523 if (numIters > 1 && ck_.is_null ()) {
1525 ck_ =
Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1530 ST rhokp1, dtemp1, dtemp2;
1531 for (
int deg = 1; deg < numIters; ++deg) {
1533 *out_ <<
" Iteration " << deg+1 <<
":" << endl
1534 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl
1535 <<
" - \\|B\\|_{\\infty} = " << maxNormInf (B) << endl
1536 <<
" - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm ()
1537 << endl <<
" - rhok = " << rhok << endl;
1540 rhokp1 = one / (two * s1 - rhok);
1541 dtemp1 = rhokp1 * rhok;
1542 dtemp2 = two * rhokp1 * delta;
1546 *out_ <<
" - dtemp1 = " << dtemp1 << endl
1547 <<
" - dtemp2 = " << dtemp2 << endl;
1552 ck_->compute (W, dtemp2, const_cast<V&> (D_inv),
1553 const_cast<MV&> (B), (X), dtemp1);
1556 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1557 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1565 template<
class ScalarType,
class MV>
1566 typename Chebyshev<ScalarType, MV>::ST
1567 Chebyshev<ScalarType, MV>::
1568 cgMethodWithInitGuess (
const op_type& A,
1574 using MagnitudeType =
typename STS::magnitudeType;
1576 *out_ <<
" cgMethodWithInitGuess:" << endl;
1581 const ST one = STS::one();
1582 ST beta, betaOld = one, pAp, pApOld = one, alpha, rHz, rHzOld, rHzOld2 = one, lambdaMax;
1587 offdiag.
resize(numIters-1);
1589 p =
rcp(
new V(A.getRangeMap ()));
1590 z =
rcp(
new V(A.getRangeMap ()));
1591 Ap =
rcp(
new V(A.getRangeMap ()));
1594 solve (*p, D_inv, r);
1597 for (
int iter = 0; iter < numIters; ++iter) {
1599 *out_ <<
" Iteration " << iter << endl;
1604 r.update (-alpha, *Ap, one);
1606 solve (*z, D_inv, r);
1608 beta = rHz / rHzOld;
1609 p->update (one, *z, beta);
1611 diag[iter] = STS::real((betaOld*betaOld * pApOld + pAp) / rHzOld);
1612 offdiag[iter-1] = -STS::real((betaOld * pApOld / (sqrt(rHzOld * rHzOld2))));
1614 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1615 *out_ <<
" offdiag["<< iter-1 <<
"] = " << offdiag[iter-1] << endl;
1616 *out_ <<
" rHz = "<<rHz <<endl;
1617 *out_ <<
" alpha = "<<alpha<<endl;
1618 *out_ <<
" beta = "<<beta<<endl;
1621 diag[iter] = STS::real(pAp/rHzOld);
1623 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1624 *out_ <<
" rHz = "<<rHz <<endl;
1625 *out_ <<
" alpha = "<<alpha<<endl;
1626 *out_ <<
" beta = "<<beta<<endl;
1634 lambdaMax = LapackHelper<ST>::tri_diag_spectral_radius(diag, offdiag);
1640 template<
class ScalarType,
class MV>
1641 typename Chebyshev<ScalarType, MV>::ST
1642 Chebyshev<ScalarType, MV>::
1643 cgMethod (
const op_type& A,
const V& D_inv,
const int numIters)
1648 *out_ <<
"cgMethod:" << endl;
1652 if (eigVector_.is_null()) {
1653 r =
rcp(
new V(A.getDomainMap ()));
1654 if (eigKeepVectors_)
1657 Details::computeInitialGuessForCG (D_inv,*r);
1661 ST lambdaMax = cgMethodWithInitGuess (A, D_inv, numIters, *r);
1666 template<
class ScalarType,
class MV>
1672 template<
class ScalarType,
class MV>
1680 template<
class ScalarType,
class MV>
1689 const size_t B_numVecs = B.getNumVectors ();
1690 if (W_.is_null () || W_->getNumVectors () != B_numVecs) {
1691 W_ =
Teuchos::rcp (
new MV (B.getMap (), B_numVecs,
false));
1696 template<
class ScalarType,
class MV>
1698 Chebyshev<ScalarType, MV>::
1699 makeSecondTempMultiVector (
const MV& B)
1705 const size_t B_numVecs = B.getNumVectors ();
1706 if (W2_.is_null () || W2_->getNumVectors () != B_numVecs) {
1707 W2_ =
Teuchos::rcp (
new MV (B.getMap (), B_numVecs,
false));
1713 template<
class ScalarType,
class MV>
1717 std::ostringstream oss;
1720 oss <<
"\"Ifpack2::Details::Chebyshev\":"
1722 <<
"degree: " << numIters_
1723 <<
", lambdaMax: " << lambdaMaxForApply_
1724 <<
", alpha: " << eigRatioForApply_
1725 <<
", lambdaMin: " << lambdaMinForApply_
1726 <<
", boost factor: " << boostFactor_
1727 <<
", algorithm: " << chebyshevAlgorithm_;
1728 if (!userInvDiag_.is_null())
1729 oss <<
", diagonal: user-supplied";
1734 template<
class ScalarType,
class MV>
1765 myRank = A_->getComm ()->getRank ();
1770 out <<
"\"Ifpack2::Details::Chebyshev\":" << endl;
1776 out <<
"degree: " << numIters_ << endl
1777 <<
"lambdaMax: " << lambdaMaxForApply_ << endl
1778 <<
"alpha: " << eigRatioForApply_ << endl
1779 <<
"lambdaMin: " << lambdaMinForApply_ << endl
1780 <<
"boost factor: " << boostFactor_ << endl;
1788 out <<
"Template parameters:" << endl;
1791 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1792 <<
"MV: " << TypeNameTraits<MV>::name () << endl;
1798 out <<
"Computed parameters:" << endl;
1810 if (D_.is_null ()) {
1812 out <<
"unset" << endl;
1817 out <<
"set" << endl;
1826 D_->describe (out, vl);
1831 out <<
"W_: " << (W_.is_null () ?
"unset" :
"set") << endl
1832 <<
"computedLambdaMax_: " << computedLambdaMax_ << endl
1833 <<
"computedLambdaMin_: " << computedLambdaMin_ << endl
1834 <<
"lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1835 <<
"lambdaMinForApply_: " << lambdaMinForApply_ << endl
1836 <<
"eigRatioForApply_: " << eigRatioForApply_ << endl;
1841 out <<
"User parameters:" << endl;
1847 out <<
"userInvDiag_: ";
1848 if (userInvDiag_.is_null ()) {
1849 out <<
"unset" << endl;
1852 out <<
"set" << endl;
1858 userInvDiag_->describe (out, vl);
1861 out <<
"userLambdaMax_: " << userLambdaMax_ << endl
1862 <<
"userLambdaMin_: " << userLambdaMin_ << endl
1863 <<
"userEigRatio_: " << userEigRatio_ << endl
1864 <<
"numIters_: " << numIters_ << endl
1865 <<
"eigMaxIters_: " << eigMaxIters_ << endl
1866 <<
"eigRelTolerance_: " << eigRelTolerance_ << endl
1867 <<
"eigNormalizationFreq_: " << eigNormalizationFreq_ << endl
1868 <<
"zeroStartingSolution_: " << zeroStartingSolution_ << endl
1869 <<
"assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1870 <<
"chebyshevAlgorithm_: " << chebyshevAlgorithm_ << endl
1871 <<
"computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1879 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \
1880 template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >;
1882 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:383
Definition of power methods.
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:1044
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:1737
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:810
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:1668
void PTEQR(const char &COMPZ, const OrdinalType &n, MagnitudeType *D, MagnitudeType *E, ScalarType *Z, const OrdinalType &ldz, MagnitudeType *WORK, OrdinalType *info) const
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:317
void resize(const size_type n, const T &val=T())
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
void print(std::ostream &out)
Print instance data to the given output stream.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1100
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:119
ScalarType ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:115
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:109
Declaration of Chebyshev implementation.
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1716
TypeTo as(const TypeFrom &t)
bool isType(const std::string &name) const
Definition of Chebyshev implementation.
bool hasTransposeApply() const
Whether it's possible to apply the transpose of this operator.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1675
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:851
basic_oblackholestream< char, std::char_traits< char > > oblackholestream