10 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
11 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
24 #include "Ifpack2_Details_ChebyshevKernel.hpp"
25 #include "Kokkos_ArithTraits.hpp"
26 #include "Teuchos_FancyOStream.hpp"
27 #include "Teuchos_oblackholestream.hpp"
28 #include "Tpetra_Details_residual.hpp"
30 #include "Ifpack2_Details_LapackSupportsScalar.hpp"
40 const char computeBeforeApplyReminder[] =
41 "This means one of the following:\n"
42 " - you have not yet called compute() on this instance, or \n"
43 " - you didn't call compute() after calling setParameters().\n\n"
44 "After creating an Ifpack2::Chebyshev instance,\n"
45 "you must _always_ call compute() at least once before calling apply().\n"
46 "After calling compute() once, you do not need to call it again,\n"
47 "unless the matrix has changed or you have changed parameters\n"
48 "(by calling setParameters()).";
54 template <
class XV,
class SizeType =
typename XV::
size_type>
55 struct V_ReciprocalThresholdSelfFunctor {
56 typedef typename XV::execution_space execution_space;
57 typedef typename XV::non_const_value_type value_type;
58 typedef SizeType size_type;
59 typedef Kokkos::ArithTraits<value_type> KAT;
60 typedef typename KAT::mag_type mag_type;
63 const value_type minVal_;
64 const mag_type minValMag_;
66 V_ReciprocalThresholdSelfFunctor(
const XV& X,
67 const value_type& min_val)
70 , minValMag_(KAT::abs(min_val)) {}
72 KOKKOS_INLINE_FUNCTION
73 void operator()(
const size_type& i)
const {
74 const mag_type X_i_abs = KAT::abs(X_[i]);
76 if (X_i_abs < minValMag_) {
79 X_[i] = KAT::one() / X_[i];
84 template <
class XV,
class SizeType =
typename XV::
size_type>
85 struct LocalReciprocalThreshold {
88 const typename XV::non_const_value_type& minVal) {
89 typedef typename XV::execution_space execution_space;
90 Kokkos::RangePolicy<execution_space, SizeType> policy(0, X.extent(0));
91 V_ReciprocalThresholdSelfFunctor<XV, SizeType> op(X, minVal);
92 Kokkos::parallel_for(policy, op);
96 template <
class TpetraVectorType,
97 const bool classic = TpetraVectorType::node_type::classic>
98 struct GlobalReciprocalThreshold {};
100 template <
class TpetraVectorType>
101 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
103 compute(TpetraVectorType& V,
104 const typename TpetraVectorType::scalar_type& min_val) {
105 typedef typename TpetraVectorType::scalar_type scalar_type;
106 typedef typename TpetraVectorType::mag_type mag_type;
107 typedef Kokkos::ArithTraits<scalar_type> STS;
109 const scalar_type ONE = STS::one();
110 const mag_type min_val_abs = STS::abs(min_val);
113 const size_t lclNumRows = V.getLocalLength();
115 for (
size_t i = 0; i < lclNumRows; ++i) {
116 const scalar_type V_0i = V_0[i];
117 if (STS::abs(V_0i) < min_val_abs) {
126 template <
class TpetraVectorType>
127 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
129 compute(TpetraVectorType& X,
130 const typename TpetraVectorType::scalar_type& minVal) {
131 typedef typename TpetraVectorType::impl_scalar_type value_type;
133 const value_type minValS =
static_cast<value_type
>(minVal);
134 auto X_0 = Kokkos::subview(X.getLocalViewDevice(Tpetra::Access::ReadWrite),
136 LocalReciprocalThreshold<decltype(X_0)>::compute(X_0, minValS);
141 template <
typename S,
typename L,
typename G,
typename N>
142 void reciprocal_threshold(Tpetra::Vector<S, L, G, N>& V,
const S& minVal) {
143 GlobalReciprocalThreshold<Tpetra::Vector<S, L, G, N>>::compute(V, minVal);
146 template <class ScalarType, const bool lapackSupportsScalarType = LapackSupportsScalar<ScalarType>::value>
147 struct LapackHelper {
151 throw std::runtime_error(
"LAPACK does not support the scalar type.");
156 void computeInitialGuessForCG(
const V& diagonal, V& x) {
157 using device_type =
typename V::node_type::device_type;
158 using range_policy = Kokkos::RangePolicy<typename device_type::execution_space>;
165 size_t N = x.getLocalLength();
166 auto d_view = diagonal.template getLocalView<device_type>(Tpetra::Access::ReadOnly);
167 auto x_view = x.template getLocalView<device_type>(Tpetra::Access::ReadWrite);
172 Kokkos::parallel_for(
173 "computeInitialGuessforCG::zero_bcs", range_policy(0, N), KOKKOS_LAMBDA(
const size_t& i) {
174 if (d_view(i, 0) == ONE)
179 template <
class ScalarType>
180 struct LapackHelper<ScalarType, true> {
185 using MagnitudeType =
typename STS::magnitudeType;
187 const int N = diag.size();
188 ScalarType scalar_dummy;
189 std::vector<MagnitudeType> mag_dummy(4 * N);
193 ScalarType lambdaMax = STS::one();
196 lapack.
PTEQR(char_N, N, diag.getRawPtr(), offdiag.getRawPtr(),
197 &scalar_dummy, 1, &mag_dummy[0], &info);
199 "Ifpack2::Details::LapackHelper::tri_diag_spectral_radius:"
200 "LAPACK's _PTEQR failed with info = "
201 << info <<
" < 0. This suggests there might be a bug in the way Ifpack2 "
202 "is calling LAPACK. Please report this to the Ifpack2 developers.");
204 lambdaMax = Teuchos::as<ScalarType>(diag[0]);
210 template <
class ScalarType,
class MV>
211 void Chebyshev<ScalarType, MV>::checkInputMatrix()
const {
213 !A_.
is_null() && A_->getGlobalNumRows() != A_->getGlobalNumCols(),
214 std::invalid_argument,
215 "Ifpack2::Chebyshev: The input matrix A must be square. "
217 << A_->getGlobalNumRows() <<
" rows and "
218 << A_->getGlobalNumCols() <<
" columns.");
230 !domainMap->isSameAs(*rangeMap), std::invalid_argument,
231 "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be "
232 "the same (in the sense of isSameAs())."
235 "for this in debug mode.");
239 template <
class ScalarType,
class MV>
240 void Chebyshev<ScalarType, MV>::
241 checkConstructorInput()
const {
245 if (STS::isComplex) {
247 "Ifpack2::Chebyshev: This class' implementation "
248 "of Chebyshev iteration only works for real-valued, symmetric positive "
249 "definite matrices. However, you instantiated this class for ScalarType"
252 "complex-valued type. While this may be algorithmically correct if all "
253 "of the complex numbers in the matrix have zero imaginary part, we "
254 "forbid using complex ScalarType altogether in order to remind you of "
255 "the limitations of our implementation (and of the algorithm itself).");
261 template <
class ScalarType,
class MV>
265 , savedDiagOffsets_(false)
266 , computedLambdaMax_(STS::nan())
267 , computedLambdaMin_(STS::nan())
268 , lambdaMaxForApply_(STS::nan())
269 , lambdaMinForApply_(STS::nan())
270 , eigRatioForApply_(STS::nan())
271 , userLambdaMax_(STS::nan())
272 , userLambdaMin_(STS::nan())
273 , userEigRatio_(Teuchos::
as<
ST>(30))
274 , minDiagVal_(STS::eps())
277 , eigRelTolerance_(Teuchos::ScalarTraits<
MT>::zero())
278 , eigKeepVectors_(false)
279 , eigenAnalysisType_(
"power method")
280 , eigNormalizationFreq_(1)
281 , zeroStartingSolution_(true)
282 , assumeMatrixUnchanged_(false)
283 , chebyshevAlgorithm_(
"first")
284 , computeMaxResNorm_(false)
285 , computeSpectralRadius_(true)
286 , ckUseNativeSpMV_(MV::node_type::is_gpu)
287 , preAllocateTempVector_(true)
289 checkConstructorInput();
292 template <
class ScalarType,
class MV>
297 , savedDiagOffsets_(false)
298 , computedLambdaMax_(STS::nan())
299 , computedLambdaMin_(STS::nan())
300 , lambdaMaxForApply_(STS::nan())
301 , boostFactor_(static_cast<
MT>(1.1))
302 , lambdaMinForApply_(STS::nan())
303 , eigRatioForApply_(STS::nan())
304 , userLambdaMax_(STS::nan())
305 , userLambdaMin_(STS::nan())
306 , userEigRatio_(Teuchos::
as<
ST>(30))
307 , minDiagVal_(STS::eps())
310 , eigRelTolerance_(Teuchos::ScalarTraits<
MT>::zero())
311 , eigKeepVectors_(false)
312 , eigenAnalysisType_(
"power method")
313 , eigNormalizationFreq_(1)
314 , zeroStartingSolution_(true)
315 , assumeMatrixUnchanged_(false)
316 , chebyshevAlgorithm_(
"first")
317 , computeMaxResNorm_(false)
318 , computeSpectralRadius_(true)
319 , ckUseNativeSpMV_(MV::node_type::is_gpu)
320 , preAllocateTempVector_(true)
322 checkConstructorInput();
326 template <
class ScalarType,
class MV>
331 using Teuchos::rcp_const_cast;
343 const ST defaultLambdaMax = STS::nan();
344 const ST defaultLambdaMin = STS::nan();
351 const ST defaultEigRatio = Teuchos::as<ST>(30);
352 const MT defaultBoostFactor =
static_cast<MT>(1.1);
353 const ST defaultMinDiagVal = STS::eps();
354 const int defaultNumIters = 1;
355 const int defaultEigMaxIters = 10;
357 const bool defaultEigKeepVectors =
false;
358 const int defaultEigNormalizationFreq = 1;
359 const bool defaultZeroStartingSolution =
true;
360 const bool defaultAssumeMatrixUnchanged =
false;
361 const std::string defaultChebyshevAlgorithm =
"first";
362 const bool defaultComputeMaxResNorm =
false;
363 const bool defaultComputeSpectralRadius =
true;
364 const bool defaultCkUseNativeSpMV = MV::node_type::is_gpu;
365 const bool defaultPreAllocateTempVector =
true;
366 const bool defaultDebug =
false;
372 RCP<const V> userInvDiagCopy;
373 ST lambdaMax = defaultLambdaMax;
374 ST lambdaMin = defaultLambdaMin;
375 ST eigRatio = defaultEigRatio;
376 MT boostFactor = defaultBoostFactor;
377 ST minDiagVal = defaultMinDiagVal;
378 int numIters = defaultNumIters;
379 int eigMaxIters = defaultEigMaxIters;
380 MT eigRelTolerance = defaultEigRelTolerance;
381 bool eigKeepVectors = defaultEigKeepVectors;
382 int eigNormalizationFreq = defaultEigNormalizationFreq;
383 bool zeroStartingSolution = defaultZeroStartingSolution;
384 bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
385 std::string chebyshevAlgorithm = defaultChebyshevAlgorithm;
386 bool computeMaxResNorm = defaultComputeMaxResNorm;
387 bool computeSpectralRadius = defaultComputeSpectralRadius;
388 bool ckUseNativeSpMV = defaultCkUseNativeSpMV;
389 bool preAllocateTempVector = defaultPreAllocateTempVector;
390 bool debug = defaultDebug;
397 if (plist.
isType<
bool>(
"debug")) {
398 debug = plist.
get<
bool>(
"debug");
399 }
else if (plist.
isType<
int>(
"debug")) {
400 const int debugInt = plist.
get<
bool>(
"debug");
401 debug = debugInt != 0;
414 const char opInvDiagLabel[] =
"chebyshev: operator inv diagonal";
417 RCP<const V> userInvDiag;
419 if (plist.
isType<
const V*>(opInvDiagLabel)) {
420 const V* rawUserInvDiag =
421 plist.
get<
const V*>(opInvDiagLabel);
423 userInvDiag =
rcp(rawUserInvDiag,
false);
424 }
else if (plist.
isType<
const V*>(opInvDiagLabel)) {
425 V* rawUserInvDiag = plist.
get<V*>(opInvDiagLabel);
427 userInvDiag =
rcp(const_cast<const V*>(rawUserInvDiag),
false);
428 }
else if (plist.
isType<RCP<const V>>(opInvDiagLabel)) {
429 userInvDiag = plist.
get<RCP<const V>>(opInvDiagLabel);
430 }
else if (plist.
isType<RCP<V>>(opInvDiagLabel)) {
431 RCP<V> userInvDiagNonConst =
432 plist.
get<RCP<V>>(opInvDiagLabel);
433 userInvDiag = rcp_const_cast<
const V>(userInvDiagNonConst);
434 }
else if (plist.
isType<
const V>(opInvDiagLabel)) {
435 const V& userInvDiagRef = plist.
get<
const V>(opInvDiagLabel);
437 userInvDiag = userInvDiagCopy;
438 }
else if (plist.
isType<V>(opInvDiagLabel)) {
439 V& userInvDiagNonConstRef = plist.
get<V>(opInvDiagLabel);
440 const V& userInvDiagRef =
const_cast<const V&
>(userInvDiagNonConstRef);
442 userInvDiag = userInvDiagCopy;
452 if (!userInvDiag.is_null() && userInvDiagCopy.is_null()) {
463 if (plist.
isParameter(
"chebyshev: use native spmv"))
464 ckUseNativeSpMV = plist.
get(
"chebyshev: use native spmv", ckUseNativeSpMV);
467 if (plist.
isParameter(
"chebyshev: pre-allocate temp vector"))
468 preAllocateTempVector = plist.
get(
"chebyshev: pre-allocate temp vector", preAllocateTempVector);
473 if (plist.
isParameter(
"chebyshev: max eigenvalue")) {
474 if (plist.
isType<
double>(
"chebyshev: max eigenvalue"))
475 lambdaMax = plist.
get<
double>(
"chebyshev: max eigenvalue");
477 lambdaMax = plist.
get<
ST>(
"chebyshev: max eigenvalue");
479 STS::isnaninf(lambdaMax), std::invalid_argument,
480 "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" "
481 "parameter is NaN or Inf. This parameter is optional, but if you "
482 "choose to supply it, it must have a finite value.");
484 if (plist.
isParameter(
"chebyshev: min eigenvalue")) {
485 if (plist.
isType<
double>(
"chebyshev: min eigenvalue"))
486 lambdaMin = plist.
get<
double>(
"chebyshev: min eigenvalue");
488 lambdaMin = plist.
get<
ST>(
"chebyshev: min eigenvalue");
490 STS::isnaninf(lambdaMin), std::invalid_argument,
491 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" "
492 "parameter is NaN or Inf. This parameter is optional, but if you "
493 "choose to supply it, it must have a finite value.");
497 if (plist.
isParameter(
"smoother: Chebyshev alpha")) {
498 if (plist.
isType<
double>(
"smoother: Chebyshev alpha"))
499 eigRatio = plist.
get<
double>(
"smoother: Chebyshev alpha");
501 eigRatio = plist.
get<
ST>(
"smoother: Chebyshev alpha");
504 eigRatio = plist.
get(
"chebyshev: ratio eigenvalue", eigRatio);
506 STS::isnaninf(eigRatio), std::invalid_argument,
507 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" "
508 "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. "
509 "This parameter is optional, but if you choose to supply it, it must have "
517 STS::real(eigRatio) < STS::real(STS::one()),
518 std::invalid_argument,
519 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\""
520 "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, "
521 "but you supplied the value "
527 const char paramName[] =
"chebyshev: boost factor";
531 boostFactor = plist.
get<
MT>(paramName);
532 }
else if (!std::is_same<double, MT>::value &&
533 plist.
isType<
double>(paramName)) {
534 const double dblBF = plist.
get<
double>(paramName);
535 boostFactor =
static_cast<MT>(dblBF);
538 "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\""
539 "parameter must have type magnitude_type (MT) or double.");
548 plist.
set(paramName, defaultBoostFactor);
551 "Ifpack2::Chebyshev::setParameters: \"" << paramName <<
"\" parameter "
552 "must be >= 1, but you supplied the value "
553 << boostFactor <<
".");
557 minDiagVal = plist.
get(
"chebyshev: min diagonal value", minDiagVal);
559 STS::isnaninf(minDiagVal), std::invalid_argument,
560 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" "
561 "parameter is NaN or Inf. This parameter is optional, but if you choose "
562 "to supply it, it must have a finite value.");
566 numIters = plist.
get<
int>(
"smoother: sweeps");
569 numIters = plist.
get<
int>(
"relaxation: sweeps");
571 numIters = plist.
get(
"chebyshev: degree", numIters);
573 numIters < 0, std::invalid_argument,
574 "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also "
575 "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a "
576 "nonnegative integer. You gave a value of "
580 if (plist.
isParameter(
"eigen-analysis: iterations")) {
581 eigMaxIters = plist.
get<
int>(
"eigen-analysis: iterations");
583 eigMaxIters = plist.
get(
"chebyshev: eigenvalue max iterations", eigMaxIters);
585 eigMaxIters < 0, std::invalid_argument,
586 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations"
587 "\" parameter (also called \"eigen-analysis: iterations\") must be a "
588 "nonnegative integer. You gave a value of "
589 << eigMaxIters <<
".");
591 if (plist.
isType<
double>(
"chebyshev: eigenvalue relative tolerance"))
592 eigRelTolerance = Teuchos::as<MT>(plist.
get<
double>(
"chebyshev: eigenvalue relative tolerance"));
593 else if (plist.
isType<
MT>(
"chebyshev: eigenvalue relative tolerance"))
594 eigRelTolerance = plist.
get<
MT>(
"chebyshev: eigenvalue relative tolerance");
595 else if (plist.
isType<
ST>(
"chebyshev: eigenvalue relative tolerance"))
598 eigKeepVectors = plist.
get(
"chebyshev: eigenvalue keep vectors", eigKeepVectors);
600 eigNormalizationFreq = plist.
get(
"chebyshev: eigenvalue normalization frequency", eigNormalizationFreq);
602 eigNormalizationFreq < 0, std::invalid_argument,
603 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue normalization frequency"
604 "\" parameter must be a "
605 "nonnegative integer. You gave a value of "
606 << eigNormalizationFreq <<
".")
608 zeroStartingSolution = plist.
get(
"chebyshev: zero starting solution",
609 zeroStartingSolution);
610 assumeMatrixUnchanged = plist.
get(
"chebyshev: assume matrix does not change",
611 assumeMatrixUnchanged);
616 chebyshevAlgorithm = plist.
get<std::string>(
"chebyshev: algorithm");
618 chebyshevAlgorithm !=
"first" &&
619 chebyshevAlgorithm !=
"textbook" &&
620 chebyshevAlgorithm !=
"fourth" &&
621 chebyshevAlgorithm !=
"opt_fourth",
622 std::invalid_argument,
623 "Ifpack2::Chebyshev: Ifpack2 only supports \"first\", \"textbook\", \"fourth\", and \"opt_fourth\", for \"chebyshev: algorithm\".");
626 #ifdef IFPACK2_ENABLE_DEPRECATED_CODE
630 if (plist.
isParameter(
"chebyshev: textbook algorithm")) {
631 const bool textbookAlgorithm = plist.
get<
bool>(
"chebyshev: textbook algorithm");
632 if (textbookAlgorithm) {
633 chebyshevAlgorithm =
"textbook";
635 chebyshevAlgorithm =
"first";
641 if (plist.
isParameter(
"chebyshev: compute max residual norm")) {
642 computeMaxResNorm = plist.
get<
bool>(
"chebyshev: compute max residual norm");
644 if (plist.
isParameter(
"chebyshev: compute spectral radius")) {
645 computeSpectralRadius = plist.
get<
bool>(
"chebyshev: compute spectral radius");
652 !plist.
get<
bool>(
"chebyshev: use block mode"),
653 std::invalid_argument,
654 "Ifpack2::Chebyshev requires that if you set \"chebyshev: use "
655 "block mode\" at all, you must set it to false. "
656 "Ifpack2::Chebyshev does not implement Ifpack's block mode.");
658 !plist.
get<
bool>(
"chebyshev: solve normal equations"),
659 std::invalid_argument,
660 "Ifpack2::Chebyshev does not and will never implement the Ifpack "
661 "parameter \"chebyshev: solve normal equations\". If you want to "
662 "solve the normal equations, construct a Tpetra::Operator that "
663 "implements A^* A, and use Chebyshev to solve A^* A x = A^* b.");
671 std::string eigenAnalysisType(
"power-method");
673 eigenAnalysisType = plist.
get<std::string>(
"eigen-analysis: type");
675 eigenAnalysisType !=
"power-method" &&
676 eigenAnalysisType !=
"power method" &&
677 eigenAnalysisType !=
"cg",
678 std::invalid_argument,
679 "Ifpack2::Chebyshev: Ifpack2 only supports \"power method\" and \"cg\" for \"eigen-analysis: type\".");
683 userInvDiag_ = userInvDiagCopy;
684 userLambdaMax_ = lambdaMax;
685 userLambdaMin_ = lambdaMin;
686 userEigRatio_ = eigRatio;
687 boostFactor_ =
static_cast<MT>(boostFactor);
688 minDiagVal_ = minDiagVal;
689 numIters_ = numIters;
690 eigMaxIters_ = eigMaxIters;
691 eigRelTolerance_ = eigRelTolerance;
692 eigKeepVectors_ = eigKeepVectors;
693 eigNormalizationFreq_ = eigNormalizationFreq;
694 eigenAnalysisType_ = eigenAnalysisType;
695 zeroStartingSolution_ = zeroStartingSolution;
696 assumeMatrixUnchanged_ = assumeMatrixUnchanged;
697 chebyshevAlgorithm_ = chebyshevAlgorithm;
698 computeMaxResNorm_ = computeMaxResNorm;
699 computeSpectralRadius_ = computeSpectralRadius;
700 ckUseNativeSpMV_ = ckUseNativeSpMV;
701 preAllocateTempVector_ = preAllocateTempVector;
712 myRank = A_->getComm()->getRank();
716 out_ = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
720 out_ = Teuchos::getFancyOStream(blackHole);
724 out_ = Teuchos::null;
728 template <
class ScalarType,
class MV>
732 diagOffsets_ = offsets_type();
733 savedDiagOffsets_ =
false;
735 computedLambdaMax_ = STS::nan();
736 computedLambdaMin_ = STS::nan();
737 eigVector_ = Teuchos::null;
738 eigVector2_ = Teuchos::null;
741 template <
class ScalarType,
class MV>
745 if (!assumeMatrixUnchanged_) {
762 myRank = A->getComm()->getRank();
766 out_ = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
769 out_ = Teuchos::getFancyOStream(blackHole);
773 out_ = Teuchos::null;
778 template <
class ScalarType,
class MV>
785 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
786 typename MV::local_ordinal_type,
787 typename MV::global_ordinal_type,
788 typename MV::node_type>
792 A_.
is_null(), std::runtime_error,
793 "Ifpack2::Chebyshev::compute: The input "
794 "matrix A is null. Please call setMatrix() with a nonnull input matrix "
795 "before calling this method.");
810 if (userInvDiag_.is_null()) {
812 Teuchos::rcp_dynamic_cast<
const crs_matrix_type>(A_);
814 if (!A_crsMat.
is_null() && A_crsMat->isFillComplete()) {
816 const size_t lclNumRows = A_crsMat->getLocalNumRows();
817 if (diagOffsets_.extent(0) < lclNumRows) {
818 diagOffsets_ = offsets_type();
819 diagOffsets_ = offsets_type(
"offsets", lclNumRows);
821 A_crsMat->getCrsGraph()->getLocalDiagOffsets(diagOffsets_);
822 savedDiagOffsets_ =
true;
823 D_ = makeInverseDiagonal(*A_,
true);
825 D_ = makeInverseDiagonal(*A_);
827 }
else if (!assumeMatrixUnchanged_) {
828 if (!A_crsMat.
is_null() && A_crsMat->isFillComplete()) {
831 if (!savedDiagOffsets_) {
832 const size_t lclNumRows = A_crsMat->getLocalNumRows();
833 if (diagOffsets_.extent(0) < lclNumRows) {
834 diagOffsets_ = offsets_type();
835 diagOffsets_ = offsets_type(
"offsets", lclNumRows);
837 A_crsMat->getCrsGraph()->getLocalDiagOffsets(diagOffsets_);
838 savedDiagOffsets_ =
true;
841 D_ = makeInverseDiagonal(*A_,
true);
843 D_ = makeInverseDiagonal(*A_);
847 D_ = makeRangeMapVectorConst(userInvDiag_);
851 const bool computedEigenvalueEstimates =
852 STS::isnaninf(computedLambdaMax_) || STS::isnaninf(computedLambdaMin_);
862 if (!assumeMatrixUnchanged_ ||
863 (!computedEigenvalueEstimates && STS::isnaninf(userLambdaMax_))) {
864 ST computedLambdaMax;
865 if ((eigenAnalysisType_ ==
"power method") || (eigenAnalysisType_ ==
"power-method")) {
867 if (eigVector_.is_null()) {
871 PowerMethod::computeInitialGuessForPowerMethod(*x,
false);
876 if (eigVector2_.is_null()) {
877 y =
rcp(
new V(A_->getRangeMap()));
884 computedLambdaMax = PowerMethod::powerMethodWithInitGuess(*A_, *D_, eigMaxIters_, x, y,
885 eigRelTolerance_, eigNormalizationFreq_, stream,
886 computeSpectralRadius_);
888 computedLambdaMax = cgMethod(*A_, *D_, eigMaxIters_);
891 STS::isnaninf(computedLambdaMax),
893 "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue "
894 "of D^{-1} A failed, by producing Inf or NaN. This probably means that "
895 "the matrix contains Inf or NaN values, or that it is badly scaled.");
897 STS::isnaninf(userEigRatio_),
899 "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
901 <<
"This should be impossible." << endl
902 <<
"Please report this bug to the Ifpack2 developers.");
908 const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
911 computedLambdaMax_ = computedLambdaMax;
912 computedLambdaMin_ = computedLambdaMin;
915 STS::isnaninf(userLambdaMax_) && STS::isnaninf(computedLambdaMax_),
917 "Ifpack2::Chebyshev::compute: " << endl
918 <<
"Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
920 <<
"This should be impossible." << endl
921 <<
"Please report this bug to the Ifpack2 developers.");
929 lambdaMaxForApply_ = STS::isnaninf(userLambdaMax_) ? computedLambdaMax_ : userLambdaMax_;
942 lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
943 eigRatioForApply_ = userEigRatio_;
945 if (chebyshevAlgorithm_ ==
"first") {
948 const ST one = Teuchos::as<ST>(1);
951 if (STS::magnitude(lambdaMaxForApply_ - one) < Teuchos::as<MT>(1.0e-6)) {
952 lambdaMinForApply_ = one;
953 lambdaMaxForApply_ = lambdaMinForApply_;
954 eigRatioForApply_ = one;
959 if (preAllocateTempVector_ && !D_.is_null()) {
960 makeTempMultiVector(*D_);
961 if (chebyshevAlgorithm_ ==
"fourth" || chebyshevAlgorithm_ ==
"opt_fourth") {
962 makeSecondTempMultiVector(*D_);
966 if (chebyshevAlgorithm_ ==
"textbook") {
972 if (ckUseNativeSpMV_) {
973 ck_->setAuxiliaryVectors(1);
978 template <
class ScalarType,
class MV>
982 return lambdaMaxForApply_;
985 template <
class ScalarType,
class MV>
988 const char prefix[] =
"Ifpack2::Chebyshev::apply: ";
991 *out_ <<
"apply: " << std::endl;
994 " Please call setMatrix() with a nonnull input matrix, and then call "
995 "compute(), before calling this method.");
997 prefix <<
"There is no estimate for the max eigenvalue."
999 << computeBeforeApplyReminder);
1000 TEUCHOS_TEST_FOR_EXCEPTION(STS::isnaninf(lambdaMinForApply_), std::runtime_error,
1001 prefix <<
"There is no estimate for the min eigenvalue."
1003 << computeBeforeApplyReminder);
1004 TEUCHOS_TEST_FOR_EXCEPTION(STS::isnaninf(eigRatioForApply_), std::runtime_error,
1005 prefix <<
"There is no estimate for the ratio of the max "
1006 "eigenvalue to the min eigenvalue."
1008 << computeBeforeApplyReminder);
1009 TEUCHOS_TEST_FOR_EXCEPTION(D_.is_null(), std::runtime_error, prefix <<
"The vector of inverse "
1010 "diagonal entries of the matrix has not yet been computed."
1012 << computeBeforeApplyReminder);
1014 if (chebyshevAlgorithm_ ==
"fourth" || chebyshevAlgorithm_ ==
"opt_fourth") {
1015 fourthKindApplyImpl(*A_, B, X, numIters_, lambdaMaxForApply_, *D_);
1016 }
else if (chebyshevAlgorithm_ ==
"textbook") {
1017 textbookApplyImpl(*A_, B, X, numIters_, lambdaMaxForApply_,
1018 lambdaMinForApply_, eigRatioForApply_, *D_);
1020 ifpackApplyImpl(*A_, B, X, numIters_, lambdaMaxForApply_,
1021 lambdaMinForApply_, eigRatioForApply_, *D_);
1024 if (computeMaxResNorm_ && B.getNumVectors() > 0) {
1025 MV R(B.getMap(), B.getNumVectors());
1026 computeResidual(R, B, *A_, X);
1029 return *std::max_element(norms.begin(), norms.end());
1035 template <
class ScalarType,
class MV>
1038 using Teuchos::rcpFromRef;
1039 this->describe(*(Teuchos::getFancyOStream(rcpFromRef(out))),
1043 template <
class ScalarType,
class MV>
1046 const ScalarType& alpha,
1050 solve(W, alpha, D_inv, B);
1051 Tpetra::deep_copy(X, W);
1054 template <
class ScalarType,
class MV>
1058 Tpetra::Details::residual(A, X, B, R);
1061 template <
class ScalarType,
class MV>
1062 void Chebyshev<ScalarType, MV>::
1063 solve(MV& Z,
const V& D_inv,
const MV& R) {
1064 Z.elementWiseMultiply(STS::one(), D_inv, R, STS::zero());
1067 template <
class ScalarType,
class MV>
1068 void Chebyshev<ScalarType, MV>::
1069 solve(MV& Z,
const ST alpha,
const V& D_inv,
const MV& R) {
1070 Z.elementWiseMultiply(alpha, D_inv, R, STS::zero());
1073 template <
class ScalarType,
class MV>
1075 Chebyshev<ScalarType, MV>::
1076 makeInverseDiagonal(
const row_matrix_type& A,
const bool useDiagOffsets)
const {
1078 using Teuchos::rcp_dynamic_cast;
1079 using Teuchos::rcpFromRef;
1082 if (!D_.is_null() &&
1083 D_->getMap()->isSameAs(*(A.getRowMap()))) {
1085 *out_ <<
"Reusing pre-existing vector for diagonal extraction" << std::endl;
1086 D_rowMap = Teuchos::rcp_const_cast<V>(D_);
1090 *out_ <<
"Allocated new vector for diagonal extraction" << std::endl;
1092 if (useDiagOffsets) {
1096 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
1097 typename MV::local_ordinal_type,
1098 typename MV::global_ordinal_type,
1099 typename MV::node_type>
1101 RCP<const crs_matrix_type> A_crsMat =
1102 rcp_dynamic_cast<
const crs_matrix_type>(rcpFromRef(A));
1103 if (!A_crsMat.is_null()) {
1105 !savedDiagOffsets_, std::logic_error,
1106 "Ifpack2::Details::Chebyshev::makeInverseDiagonal: "
1107 "It is not allowed to call this method with useDiagOffsets=true, "
1108 "if you have not previously saved offsets of diagonal entries. "
1109 "This situation should never arise if this class is used properly. "
1110 "Please report this bug to the Ifpack2 developers.");
1111 A_crsMat->getLocalDiagCopy(*D_rowMap, diagOffsets_);
1116 A.getLocalDiagCopy(*D_rowMap);
1118 RCP<V> D_rangeMap = makeRangeMapVector(D_rowMap);
1124 bool foundNonpositiveValue =
false;
1126 auto D_lcl = D_rangeMap->getLocalViewHost(Tpetra::Access::ReadOnly);
1127 auto D_lcl_1d = Kokkos::subview(D_lcl, Kokkos::ALL(), 0);
1129 typedef typename MV::impl_scalar_type IST;
1130 typedef typename MV::local_ordinal_type LO;
1131 typedef Kokkos::ArithTraits<IST> ATS;
1132 typedef Kokkos::ArithTraits<typename ATS::mag_type> STM;
1134 const LO lclNumRows =
static_cast<LO
>(D_rangeMap->getLocalLength());
1135 for (LO i = 0; i < lclNumRows; ++i) {
1136 if (STS::real(D_lcl_1d(i)) <= STM::zero()) {
1137 foundNonpositiveValue =
true;
1143 using Teuchos::outArg;
1145 using Teuchos::reduceAll;
1147 const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1148 int gblSuccess = lclSuccess;
1149 if (!D_rangeMap->getMap().is_null() && D_rangeMap->getMap()->getComm().is_null()) {
1151 reduceAll<int, int>(comm,
REDUCE_MIN, lclSuccess, outArg(gblSuccess));
1153 if (gblSuccess == 1) {
1154 *out_ <<
"makeInverseDiagonal: The matrix's diagonal entries all have "
1155 "positive real part (this is good for Chebyshev)."
1158 *out_ <<
"makeInverseDiagonal: The matrix's diagonal has at least one "
1159 "entry with nonpositive real part, on at least one process in the "
1160 "matrix's communicator. This is bad for Chebyshev."
1167 reciprocal_threshold(*D_rangeMap, minDiagVal_);
1168 return Teuchos::rcp_const_cast<
const V>(D_rangeMap);
1171 template <
class ScalarType,
class MV>
1173 Chebyshev<ScalarType, MV>::
1177 typedef Tpetra::Export<
typename MV::local_ordinal_type,
1178 typename MV::global_ordinal_type,
1179 typename MV::node_type>
1185 A_.
is_null(), std::logic_error,
1186 "Ifpack2::Details::Chebyshev::"
1187 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() "
1188 "with a nonnull input matrix before calling this method. This is probably "
1189 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1191 D.
is_null(), std::logic_error,
1192 "Ifpack2::Details::Chebyshev::"
1193 "makeRangeMapVector: The input Vector D is null. This is probably "
1194 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1196 RCP<const map_type> sourceMap = D->getMap();
1197 RCP<const map_type> rangeMap = A_->getRangeMap();
1198 RCP<const map_type> rowMap = A_->getRowMap();
1200 if (rangeMap->isSameAs(*sourceMap)) {
1205 RCP<const export_type> exporter;
1209 if (sourceMap->isSameAs(*rowMap)) {
1211 exporter = A_->getGraph()->getExporter();
1213 exporter =
rcp(
new export_type(sourceMap, rangeMap));
1216 if (exporter.is_null()) {
1220 D_out->doExport(*D, *exporter, Tpetra::ADD);
1221 return Teuchos::rcp_const_cast<
const V>(D_out);
1226 template <
class ScalarType,
class MV>
1228 Chebyshev<ScalarType, MV>::
1230 using Teuchos::rcp_const_cast;
1231 return rcp_const_cast<V>(makeRangeMapVectorConst(rcp_const_cast<V>(D)));
1234 template <
class ScalarType,
class MV>
1235 void Chebyshev<ScalarType, MV>::
1236 textbookApplyImpl(
const op_type& A,
1243 const V& D_inv)
const {
1245 const ST myLambdaMin = lambdaMax / eigRatio;
1247 const ST zero = Teuchos::as<ST>(0);
1248 const ST one = Teuchos::as<ST>(1);
1249 const ST two = Teuchos::as<ST>(2);
1250 const ST d = (lambdaMax + myLambdaMin) / two;
1251 const ST c = (lambdaMax - myLambdaMin) / two;
1253 if (zeroStartingSolution_ && numIters > 0) {
1257 MV R(B.getMap(), B.getNumVectors(),
false);
1258 MV P(B.getMap(), B.getNumVectors(),
false);
1259 MV Z(B.getMap(), B.getNumVectors(),
false);
1261 for (
int i = 0; i < numIters; ++i) {
1262 computeResidual(R, B, A, X);
1271 beta = alpha * (c / two) * (c / two);
1272 alpha = one / (d - beta);
1273 P.update(one, Z, beta);
1275 X.update(alpha, P, one);
1282 template <
class ScalarType,
class MV>
1283 void Chebyshev<ScalarType, MV>::
1284 fourthKindApplyImpl(
const op_type& A,
1291 std::vector<ScalarType> betas(numIters, 1.0);
1292 if (chebyshevAlgorithm_ ==
"opt_fourth") {
1293 betas = optimalWeightsImpl<ScalarType>(numIters);
1296 const ST invEig = MT(1) / (lambdaMax * boostFactor_);
1308 if (!zeroStartingSolution_) {
1310 Tpetra::deep_copy(X4, X);
1312 if (ck_.is_null()) {
1314 ck_ =
Teuchos::rcp(
new ChebyshevKernel<op_type>(A_op, ckUseNativeSpMV_));
1318 ck_->compute(Z, MT(4.0 / 3.0) * invEig, const_cast<V&>(D_inv),
1319 const_cast<MV&>(B), X4, STS::zero());
1322 X.update(betas[0], Z, STS::one());
1325 firstIterationWithZeroStartingSolution(Z, MT(4.0 / 3.0) * invEig, D_inv, B, X4);
1328 X.update(betas[0], Z, STS::zero());
1331 if (numIters > 1 && ck_.is_null()) {
1333 ck_ =
Teuchos::rcp(
new ChebyshevKernel<op_type>(A_op, ckUseNativeSpMV_));
1336 for (
int i = 1; i < numIters; ++i) {
1337 const ST zScale = (2.0 * i - 1.0) / (2.0 * i + 3.0);
1338 const ST rScale = MT((8.0 * i + 4.0) / (2.0 * i + 3.0)) * invEig;
1342 ck_->compute(Z, rScale, const_cast<V&>(D_inv),
1343 const_cast<MV&>(B), (X4), zScale);
1346 X.update(betas[i], Z, STS::one());
1350 template <
class ScalarType,
class MV>
1351 typename Chebyshev<ScalarType, MV>::MT
1352 Chebyshev<ScalarType, MV>::maxNormInf(
const MV& X) {
1355 return *std::max_element(norms.begin(), norms.end());
1358 template <
class ScalarType,
class MV>
1359 void Chebyshev<ScalarType, MV>::
1360 ifpackApplyImpl(
const op_type& A,
1369 #ifdef HAVE_IFPACK2_DEBUG
1370 const bool debug = debug_;
1372 const bool debug =
false;
1376 *out_ <<
" \\|B\\|_{\\infty} = " << maxNormInf(B) << endl;
1377 *out_ <<
" \\|X\\|_{\\infty} = " << maxNormInf(X) << endl;
1380 if (numIters <= 0) {
1383 const ST zero =
static_cast<ST
>(0.0);
1384 const ST one =
static_cast<ST
>(1.0);
1385 const ST two =
static_cast<ST
>(2.0);
1388 if (lambdaMin == one && lambdaMax == lambdaMin) {
1394 const ST alpha = lambdaMax / eigRatio;
1395 const ST beta = boostFactor_ * lambdaMax;
1396 const ST delta = two / (beta - alpha);
1397 const ST theta = (beta + alpha) / two;
1398 const ST s1 = theta * delta;
1401 *out_ <<
" alpha = " << alpha << endl
1402 <<
" beta = " << beta << endl
1403 <<
" delta = " << delta << endl
1404 <<
" theta = " << theta << endl
1405 <<
" s1 = " << s1 << endl;
1413 *out_ <<
" Iteration " << 1 <<
":" << endl
1414 <<
" - \\|D\\|_{\\infty} = " << D_->normInf() << endl;
1418 if (!zeroStartingSolution_) {
1421 if (ck_.is_null()) {
1423 ck_ =
Teuchos::rcp(
new ChebyshevKernel<op_type>(A_op, ckUseNativeSpMV_));
1427 ck_->compute(W, one / theta, const_cast<V&>(D_inv),
1428 const_cast<MV&>(B), X, zero);
1431 firstIterationWithZeroStartingSolution(W, one / theta, D_inv, B, X);
1435 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf(W) << endl
1436 <<
" - \\|X\\|_{\\infty} = " << maxNormInf(X) << endl;
1439 if (numIters > 1 && ck_.is_null()) {
1441 ck_ =
Teuchos::rcp(
new ChebyshevKernel<op_type>(A_op, ckUseNativeSpMV_));
1446 ST rhokp1, dtemp1, dtemp2;
1447 for (
int deg = 1; deg < numIters; ++deg) {
1449 *out_ <<
" Iteration " << deg + 1 <<
":" << endl
1450 <<
" - \\|D\\|_{\\infty} = " << D_->normInf() << endl
1451 <<
" - \\|B\\|_{\\infty} = " << maxNormInf(B) << endl
1452 <<
" - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm()
1454 <<
" - rhok = " << rhok << endl;
1457 rhokp1 = one / (two * s1 - rhok);
1458 dtemp1 = rhokp1 * rhok;
1459 dtemp2 = two * rhokp1 * delta;
1463 *out_ <<
" - dtemp1 = " << dtemp1 << endl
1464 <<
" - dtemp2 = " << dtemp2 << endl;
1469 ck_->compute(W, dtemp2, const_cast<V&>(D_inv),
1470 const_cast<MV&>(B), (X), dtemp1);
1473 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf(W) << endl
1474 <<
" - \\|X\\|_{\\infty} = " << maxNormInf(X) << endl;
1479 template <
class ScalarType,
class MV>
1480 typename Chebyshev<ScalarType, MV>::ST
1481 Chebyshev<ScalarType, MV>::
1482 cgMethodWithInitGuess(
const op_type& A,
1487 using MagnitudeType =
typename STS::magnitudeType;
1489 *out_ <<
" cgMethodWithInitGuess:" << endl;
1492 const ST one = STS::one();
1493 ST beta, betaOld = one, pAp, pApOld = one, alpha, rHz, rHzOld, rHzOld2 = one, lambdaMax;
1498 offdiag.
resize(numIters - 1);
1500 p =
rcp(
new V(A.getRangeMap()));
1501 z =
rcp(
new V(A.getRangeMap()));
1502 Ap =
rcp(
new V(A.getRangeMap()));
1505 solve(*p, D_inv, r);
1508 for (
int iter = 0; iter < numIters; ++iter) {
1510 *out_ <<
" Iteration " << iter << endl;
1515 r.update(-alpha, *Ap, one);
1517 solve(*z, D_inv, r);
1519 beta = rHz / rHzOld;
1520 p->update(one, *z, beta);
1522 diag[iter] = STS::real((betaOld * betaOld * pApOld + pAp) / rHzOld);
1523 offdiag[iter - 1] = -STS::real((betaOld * pApOld / (sqrt(rHzOld * rHzOld2))));
1525 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1526 *out_ <<
" offdiag[" << iter - 1 <<
"] = " << offdiag[iter - 1] << endl;
1527 *out_ <<
" rHz = " << rHz << endl;
1528 *out_ <<
" alpha = " << alpha << endl;
1529 *out_ <<
" beta = " << beta << endl;
1532 diag[iter] = STS::real(pAp / rHzOld);
1534 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1535 *out_ <<
" rHz = " << rHz << endl;
1536 *out_ <<
" alpha = " << alpha << endl;
1537 *out_ <<
" beta = " << beta << endl;
1545 lambdaMax = LapackHelper<ST>::tri_diag_spectral_radius(diag, offdiag);
1550 template <
class ScalarType,
class MV>
1551 typename Chebyshev<ScalarType, MV>::ST
1552 Chebyshev<ScalarType, MV>::
1553 cgMethod(
const op_type& A,
const V& D_inv,
const int numIters) {
1557 *out_ <<
"cgMethod:" << endl;
1561 if (eigVector_.is_null()) {
1562 r =
rcp(
new V(A.getDomainMap()));
1563 if (eigKeepVectors_)
1566 Details::computeInitialGuessForCG(D_inv, *r);
1570 ST lambdaMax = cgMethodWithInitGuess(A, D_inv, numIters, *r);
1575 template <
class ScalarType,
class MV>
1581 template <
class ScalarType,
class MV>
1588 template <
class ScalarType,
class MV>
1596 const size_t B_numVecs = B.getNumVectors();
1597 if (W_.is_null() || W_->getNumVectors() != B_numVecs) {
1598 W_ =
Teuchos::rcp(
new MV(B.getMap(), B_numVecs,
false));
1603 template <
class ScalarType,
class MV>
1605 Chebyshev<ScalarType, MV>::
1606 makeSecondTempMultiVector(
const MV& B) {
1611 const size_t B_numVecs = B.getNumVectors();
1612 if (W2_.is_null() || W2_->getNumVectors() != B_numVecs) {
1613 W2_ =
Teuchos::rcp(
new MV(B.getMap(), B_numVecs,
false));
1618 template <
class ScalarType,
class MV>
1622 std::ostringstream oss;
1625 oss <<
"\"Ifpack2::Details::Chebyshev\":"
1627 <<
"degree: " << numIters_
1628 <<
", lambdaMax: " << lambdaMaxForApply_
1629 <<
", alpha: " << eigRatioForApply_
1630 <<
", lambdaMin: " << lambdaMinForApply_
1631 <<
", boost factor: " << boostFactor_
1632 <<
", algorithm: " << chebyshevAlgorithm_;
1633 if (!userInvDiag_.is_null())
1634 oss <<
", diagonal: user-supplied";
1639 template <
class ScalarType,
class MV>
1667 myRank = A_->getComm()->getRank();
1672 out <<
"\"Ifpack2::Details::Chebyshev\":" << endl;
1678 out <<
"degree: " << numIters_ << endl
1679 <<
"lambdaMax: " << lambdaMaxForApply_ << endl
1680 <<
"alpha: " << eigRatioForApply_ << endl
1681 <<
"lambdaMin: " << lambdaMinForApply_ << endl
1682 <<
"boost factor: " << boostFactor_ << endl;
1690 out <<
"Template parameters:" << endl;
1693 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name() << endl
1694 <<
"MV: " << TypeNameTraits<MV>::name() << endl;
1700 out <<
"Computed parameters:" << endl;
1714 out <<
"unset" << endl;
1718 out <<
"set" << endl;
1726 D_->describe(out, vl);
1731 out <<
"W_: " << (W_.is_null() ?
"unset" :
"set") << endl
1732 <<
"computedLambdaMax_: " << computedLambdaMax_ << endl
1733 <<
"computedLambdaMin_: " << computedLambdaMin_ << endl
1734 <<
"lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1735 <<
"lambdaMinForApply_: " << lambdaMinForApply_ << endl
1736 <<
"eigRatioForApply_: " << eigRatioForApply_ << endl;
1741 out <<
"User parameters:" << endl;
1747 out <<
"userInvDiag_: ";
1748 if (userInvDiag_.is_null()) {
1749 out <<
"unset" << endl;
1751 out <<
"set" << endl;
1756 userInvDiag_->describe(out, vl);
1759 out <<
"userLambdaMax_: " << userLambdaMax_ << endl
1760 <<
"userLambdaMin_: " << userLambdaMin_ << endl
1761 <<
"userEigRatio_: " << userEigRatio_ << endl
1762 <<
"numIters_: " << numIters_ << endl
1763 <<
"eigMaxIters_: " << eigMaxIters_ << endl
1764 <<
"eigRelTolerance_: " << eigRelTolerance_ << endl
1765 <<
"eigNormalizationFreq_: " << eigNormalizationFreq_ << endl
1766 <<
"zeroStartingSolution_: " << zeroStartingSolution_ << endl
1767 <<
"assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1768 <<
"chebyshevAlgorithm_: " << chebyshevAlgorithm_ << endl
1769 <<
"computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1777 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S, LO, GO, N) \
1778 template class Ifpack2::Details::Chebyshev<S, Tpetra::MultiVector<S, LO, GO, N>>;
1780 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:328
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:987
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:1641
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:743
Compute scaled damped residual for Chebyshev.
Definition: Ifpack2_Details_ChebyshevKernel_decl.hpp:45
#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:1577
ParameterList & set(std::string const &name, T &&value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
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:263
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:165
void print(std::ostream &out)
Print instance data to the given output stream.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1037
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:85
ScalarType ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:81
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:75
Declaration of Chebyshev implementation.
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1621
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:1583
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:779
basic_oblackholestream< char, std::char_traits< char > > oblackholestream