44#ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
45#define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
57#include "Ifpack2_Details_ChebyshevKernel.hpp"
58#include "Kokkos_ArithTraits.hpp"
59#include "Teuchos_FancyOStream.hpp"
60#include "Teuchos_oblackholestream.hpp"
61#include "Tpetra_Details_residual.hpp"
62#include "Teuchos_LAPACK.hpp"
63#include "Ifpack2_Details_LapackSupportsScalar.hpp"
73const char computeBeforeApplyReminder[] =
74 "This means one of the following:\n"
75 " - you have not yet called compute() on this instance, or \n"
76 " - you didn't call compute() after calling setParameters().\n\n"
77 "After creating an Ifpack2::Chebyshev instance,\n"
78 "you must _always_ call compute() at least once before calling apply().\n"
79 "After calling compute() once, you do not need to call it again,\n"
80 "unless the matrix has changed or you have changed parameters\n"
81 "(by calling setParameters()).";
87template<
class XV,
class SizeType =
typename XV::
size_type>
88struct V_ReciprocalThresholdSelfFunctor
90 typedef typename XV::execution_space execution_space;
91 typedef typename XV::non_const_value_type value_type;
92 typedef SizeType size_type;
93 typedef Kokkos::Details::ArithTraits<value_type> KAT;
94 typedef typename KAT::mag_type mag_type;
97 const value_type minVal_;
98 const mag_type minValMag_;
100 V_ReciprocalThresholdSelfFunctor (
const XV& X,
101 const value_type& min_val) :
104 minValMag_ (KAT::abs (min_val))
107 KOKKOS_INLINE_FUNCTION
108 void operator() (
const size_type& i)
const
110 const mag_type X_i_abs = KAT::abs (X_[i]);
112 if (X_i_abs < minValMag_) {
116 X_[i] = KAT::one () / X_[i];
121template<
class XV,
class SizeType =
typename XV::
size_type>
122struct LocalReciprocalThreshold {
124 compute (
const XV& X,
125 const typename XV::non_const_value_type& minVal)
127 typedef typename XV::execution_space execution_space;
128 Kokkos::RangePolicy<execution_space, SizeType> policy (0, X.extent (0));
129 V_ReciprocalThresholdSelfFunctor<XV, SizeType> op (X, minVal);
130 Kokkos::parallel_for (policy, op);
134template <
class TpetraVectorType,
135 const bool classic = TpetraVectorType::node_type::classic>
136struct GlobalReciprocalThreshold {};
138template <
class TpetraVectorType>
139struct GlobalReciprocalThreshold<TpetraVectorType, true> {
141 compute (TpetraVectorType& V,
142 const typename TpetraVectorType::scalar_type& min_val)
144 typedef typename TpetraVectorType::scalar_type scalar_type;
145 typedef typename TpetraVectorType::mag_type mag_type;
146 typedef Kokkos::Details::ArithTraits<scalar_type>
STS;
148 const scalar_type ONE = STS::one ();
149 const mag_type min_val_abs = STS::abs (min_val);
151 Teuchos::ArrayRCP<scalar_type> V_0 =
V.getDataNonConst (0);
152 const size_t lclNumRows =
V.getLocalLength ();
154 for (
size_t i = 0; i < lclNumRows; ++i) {
155 const scalar_type V_0i = V_0[i];
156 if (STS::abs (V_0i) < min_val_abs) {
165template <
class TpetraVectorType>
166struct GlobalReciprocalThreshold<TpetraVectorType, false> {
168 compute (TpetraVectorType& X,
169 const typename TpetraVectorType::scalar_type& minVal)
171 typedef typename TpetraVectorType::impl_scalar_type value_type;
173 const value_type minValS =
static_cast<value_type
> (minVal);
174 auto X_0 = Kokkos::subview (X.getLocalViewDevice (Tpetra::Access::ReadWrite),
176 LocalReciprocalThreshold<
decltype (X_0) >::compute (X_0, minValS);
181template <
typename S,
typename L,
typename G,
typename N>
183reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V,
const S& minVal)
185 GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
189template<class ScalarType, const bool lapackSupportsScalarType = LapackSupportsScalar<ScalarType>::value>
193 tri_diag_spectral_radius(Teuchos::ArrayRCP<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
194 Teuchos::ArrayRCP<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> offdiag) {
195 throw std::runtime_error(
"LAPACK does not support the scalar type.");
200template<
class ScalarType>
201struct LapackHelper<ScalarType,true> {
204 tri_diag_spectral_radius(Teuchos::ArrayRCP<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
205 Teuchos::ArrayRCP<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> offdiag) {
206 using STS = Teuchos::ScalarTraits<ScalarType>;
207 using MagnitudeType =
typename STS::magnitudeType;
209 const int N = diag.size ();
210 ScalarType scalar_dummy;
211 std::vector<MagnitudeType> mag_dummy(4*N);
215 ScalarType lambdaMax = STS::one();
217 Teuchos::LAPACK<int,ScalarType> lapack;
218 lapack.PTEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
219 &scalar_dummy, 1, &mag_dummy[0], &info);
220 TEUCHOS_TEST_FOR_EXCEPTION
221 (info < 0, std::logic_error,
"Ifpack2::Details::LapackHelper::tri_diag_spectral_radius:"
222 "LAPACK's _PTEQR failed with info = "
223 << info <<
" < 0. This suggests there might be a bug in the way Ifpack2 "
224 "is calling LAPACK. Please report this to the Ifpack2 developers.");
226 lambdaMax = Teuchos::as<ScalarType> (diag[0]);
233template<
class ScalarType,
class MV>
234void Chebyshev<ScalarType, MV>::checkInputMatrix ()
const
236 TEUCHOS_TEST_FOR_EXCEPTION(
237 ! A_.is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
238 std::invalid_argument,
239 "Ifpack2::Chebyshev: The input matrix A must be square. "
240 "A has " << A_->getGlobalNumRows () <<
" rows and "
241 << A_->getGlobalNumCols () <<
" columns.");
245 if (debug_ && ! A_.is_null ()) {
246 Teuchos::RCP<const map_type> domainMap = A_->getDomainMap ();
247 Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap ();
252 TEUCHOS_TEST_FOR_EXCEPTION(
253 ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
254 "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be "
255 "the same (in the sense of isSameAs())." << std::endl <<
"We only check "
256 "for this in debug mode.");
261template<
class ScalarType,
class MV>
263Chebyshev<ScalarType, MV>::
264checkConstructorInput ()
const
269 if (STS::isComplex) {
270 TEUCHOS_TEST_FOR_EXCEPTION
271 (
true, std::logic_error,
"Ifpack2::Chebyshev: This class' implementation "
272 "of Chebyshev iteration only works for real-valued, symmetric positive "
273 "definite matrices. However, you instantiated this class for ScalarType"
274 " = " << Teuchos::TypeNameTraits<ScalarType>::name () <<
", which is a "
275 "complex-valued type. While this may be algorithmically correct if all "
276 "of the complex numbers in the matrix have zero imaginary part, we "
277 "forbid using complex ScalarType altogether in order to remind you of "
278 "the limitations of our implementation (and of the algorithm itself).");
285template<
class ScalarType,
class MV>
287Chebyshev (Teuchos::RCP<const row_matrix_type> A) :
289 savedDiagOffsets_ (false),
290 computedLambdaMax_ (
STS::nan ()),
291 computedLambdaMin_ (
STS::nan ()),
292 lambdaMaxForApply_ (
STS::nan ()),
293 lambdaMinForApply_ (
STS::nan ()),
294 eigRatioForApply_ (
STS::nan ()),
295 userLambdaMax_ (
STS::nan ()),
296 userLambdaMin_ (
STS::nan ()),
297 userEigRatio_ (Teuchos::as<
ST> (30)),
298 minDiagVal_ (
STS::eps ()),
301 eigRelTolerance_(Teuchos::ScalarTraits<
MT>::zero ()),
302 eigKeepVectors_(false),
303 eigenAnalysisType_(
"power method"),
304 eigNormalizationFreq_(1),
305 zeroStartingSolution_ (true),
306 assumeMatrixUnchanged_ (false),
307 textbookAlgorithm_ (false),
308 computeMaxResNorm_ (false),
309 computeSpectralRadius_(true),
310 ckUseNativeSpMV_(false),
313 checkConstructorInput ();
316template<
class ScalarType,
class MV>
318Chebyshev (Teuchos::RCP<const row_matrix_type> A,
319 Teuchos::ParameterList& params) :
321 savedDiagOffsets_ (false),
322 computedLambdaMax_ (
STS::nan ()),
323 computedLambdaMin_ (
STS::nan ()),
324 lambdaMaxForApply_ (
STS::nan ()),
325 boostFactor_ (static_cast<
MT> (1.1)),
326 lambdaMinForApply_ (
STS::nan ()),
327 eigRatioForApply_ (
STS::nan ()),
328 userLambdaMax_ (
STS::nan ()),
329 userLambdaMin_ (
STS::nan ()),
330 userEigRatio_ (Teuchos::as<
ST> (30)),
331 minDiagVal_ (
STS::eps ()),
334 eigRelTolerance_(Teuchos::ScalarTraits<
MT>::zero ()),
335 eigKeepVectors_(false),
336 eigenAnalysisType_(
"power method"),
337 eigNormalizationFreq_(1),
338 zeroStartingSolution_ (true),
339 assumeMatrixUnchanged_ (false),
340 textbookAlgorithm_ (false),
341 computeMaxResNorm_ (false),
342 computeSpectralRadius_(true),
343 ckUseNativeSpMV_(false),
346 checkConstructorInput ();
347 setParameters (params);
350template<
class ScalarType,
class MV>
357 using Teuchos::rcp_const_cast;
369 const ST defaultLambdaMax = STS::nan ();
370 const ST defaultLambdaMin = STS::nan ();
377 const ST defaultEigRatio = Teuchos::as<ST> (30);
378 const MT defaultBoostFactor =
static_cast<MT> (1.1);
379 const ST defaultMinDiagVal = STS::eps ();
380 const int defaultNumIters = 1;
381 const int defaultEigMaxIters = 10;
382 const MT defaultEigRelTolerance = Teuchos::ScalarTraits<MT>::zero ();
383 const bool defaultEigKeepVectors =
false;
384 const int defaultEigNormalizationFreq = 1;
385 const bool defaultZeroStartingSolution =
true;
386 const bool defaultAssumeMatrixUnchanged =
false;
387 const bool defaultTextbookAlgorithm =
false;
388 const bool defaultComputeMaxResNorm =
false;
389 const bool defaultComputeSpectralRadius =
true;
390 const bool defaultCkUseNativeSpMV =
false;
391 const bool defaultDebug =
false;
397 RCP<const V> userInvDiagCopy;
398 ST lambdaMax = defaultLambdaMax;
399 ST lambdaMin = defaultLambdaMin;
400 ST eigRatio = defaultEigRatio;
401 MT boostFactor = defaultBoostFactor;
402 ST minDiagVal = defaultMinDiagVal;
403 int numIters = defaultNumIters;
404 int eigMaxIters = defaultEigMaxIters;
405 MT eigRelTolerance = defaultEigRelTolerance;
406 bool eigKeepVectors = defaultEigKeepVectors;
407 int eigNormalizationFreq = defaultEigNormalizationFreq;
408 bool zeroStartingSolution = defaultZeroStartingSolution;
409 bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
410 bool textbookAlgorithm = defaultTextbookAlgorithm;
411 bool computeMaxResNorm = defaultComputeMaxResNorm;
412 bool computeSpectralRadius = defaultComputeSpectralRadius;
413 bool ckUseNativeSpMV = defaultCkUseNativeSpMV;
414 bool debug = defaultDebug;
421 if (plist.isType<
bool> (
"debug")) {
422 debug = plist.get<
bool> (
"debug");
424 else if (plist.isType<
int> (
"debug")) {
425 const int debugInt = plist.get<
bool> (
"debug");
426 debug = debugInt != 0;
439 const char opInvDiagLabel[] =
"chebyshev: operator inv diagonal";
440 if (plist.isParameter (opInvDiagLabel)) {
442 RCP<const V> userInvDiag;
444 if (plist.isType<
const V*> (opInvDiagLabel)) {
445 const V* rawUserInvDiag =
446 plist.get<
const V*> (opInvDiagLabel);
448 userInvDiag = rcp (rawUserInvDiag,
false);
450 else if (plist.isType<
const V*> (opInvDiagLabel)) {
451 V* rawUserInvDiag = plist.get<
V*> (opInvDiagLabel);
453 userInvDiag = rcp (
const_cast<const V*
> (rawUserInvDiag),
false);
455 else if (plist.isType<RCP<const V>> (opInvDiagLabel)) {
456 userInvDiag = plist.get<RCP<const V> > (opInvDiagLabel);
458 else if (plist.isType<RCP<V>> (opInvDiagLabel)) {
459 RCP<V> userInvDiagNonConst =
460 plist.get<RCP<V> > (opInvDiagLabel);
461 userInvDiag = rcp_const_cast<const V> (userInvDiagNonConst);
463 else if (plist.isType<
const V> (opInvDiagLabel)) {
464 const V& userInvDiagRef = plist.get<
const V> (opInvDiagLabel);
465 userInvDiagCopy = rcp (
new V (userInvDiagRef, Teuchos::Copy));
466 userInvDiag = userInvDiagCopy;
468 else if (plist.isType<
V> (opInvDiagLabel)) {
469 V& userInvDiagNonConstRef = plist.get<
V> (opInvDiagLabel);
470 const V& userInvDiagRef =
const_cast<const V&
> (userInvDiagNonConstRef);
471 userInvDiagCopy = rcp (
new V (userInvDiagRef, Teuchos::Copy));
472 userInvDiag = userInvDiagCopy;
482 if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
483 userInvDiagCopy = rcp (
new V (*userInvDiag, Teuchos::Copy));
493 if (plist.isParameter (
"chebyshev: use native spmv"))
494 ckUseNativeSpMV = plist.get(
"chebyshev: use native spmv", ckUseNativeSpMV);
499 if (plist.isParameter (
"chebyshev: max eigenvalue")) {
500 if (plist.isType<
double>(
"chebyshev: max eigenvalue"))
501 lambdaMax = plist.get<
double> (
"chebyshev: max eigenvalue");
503 lambdaMax = plist.get<
ST> (
"chebyshev: max eigenvalue");
504 TEUCHOS_TEST_FOR_EXCEPTION(
505 STS::isnaninf (lambdaMax), std::invalid_argument,
506 "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" "
507 "parameter is NaN or Inf. This parameter is optional, but if you "
508 "choose to supply it, it must have a finite value.");
510 if (plist.isParameter (
"chebyshev: min eigenvalue")) {
511 if (plist.isType<
double>(
"chebyshev: min eigenvalue"))
512 lambdaMin = plist.get<
double> (
"chebyshev: min eigenvalue");
514 lambdaMin = plist.get<
ST> (
"chebyshev: min eigenvalue");
515 TEUCHOS_TEST_FOR_EXCEPTION(
516 STS::isnaninf (lambdaMin), std::invalid_argument,
517 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" "
518 "parameter is NaN or Inf. This parameter is optional, but if you "
519 "choose to supply it, it must have a finite value.");
523 if (plist.isParameter (
"smoother: Chebyshev alpha")) {
524 if (plist.isType<
double>(
"smoother: Chebyshev alpha"))
525 eigRatio = plist.get<
double> (
"smoother: Chebyshev alpha");
527 eigRatio = plist.get<
ST> (
"smoother: Chebyshev alpha");
530 eigRatio = plist.get (
"chebyshev: ratio eigenvalue", eigRatio);
531 TEUCHOS_TEST_FOR_EXCEPTION(
532 STS::isnaninf (eigRatio), std::invalid_argument,
533 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" "
534 "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. "
535 "This parameter is optional, but if you choose to supply it, it must have "
542 TEUCHOS_TEST_FOR_EXCEPTION(
543 STS::real (eigRatio) < STS::real (STS::one ()),
544 std::invalid_argument,
545 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\""
546 "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, "
547 "but you supplied the value " << eigRatio <<
".");
552 const char paramName[] =
"chebyshev: boost factor";
554 if (plist.isParameter (paramName)) {
555 if (plist.isType<
MT> (paramName)) {
556 boostFactor = plist.get<
MT> (paramName);
558 else if (! std::is_same<double, MT>::value &&
559 plist.isType<
double> (paramName)) {
560 const double dblBF = plist.get<
double> (paramName);
561 boostFactor =
static_cast<MT> (dblBF);
564 TEUCHOS_TEST_FOR_EXCEPTION
565 (
true, std::invalid_argument,
566 "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\""
567 "parameter must have type magnitude_type (MT) or double.");
577 plist.set (paramName, defaultBoostFactor);
579 TEUCHOS_TEST_FOR_EXCEPTION
580 (boostFactor < Teuchos::ScalarTraits<MT>::one (), std::invalid_argument,
581 "Ifpack2::Chebyshev::setParameters: \"" << paramName <<
"\" parameter "
582 "must be >= 1, but you supplied the value " << boostFactor <<
".");
586 minDiagVal = plist.get (
"chebyshev: min diagonal value", minDiagVal);
587 TEUCHOS_TEST_FOR_EXCEPTION(
588 STS::isnaninf (minDiagVal), std::invalid_argument,
589 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" "
590 "parameter is NaN or Inf. This parameter is optional, but if you choose "
591 "to supply it, it must have a finite value.");
594 if (plist.isParameter (
"smoother: sweeps")) {
595 numIters = plist.get<
int> (
"smoother: sweeps");
597 if (plist.isParameter (
"relaxation: sweeps")) {
598 numIters = plist.get<
int> (
"relaxation: sweeps");
600 numIters = plist.get (
"chebyshev: degree", numIters);
601 TEUCHOS_TEST_FOR_EXCEPTION(
602 numIters < 0, std::invalid_argument,
603 "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also "
604 "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a "
605 "nonnegative integer. You gave a value of " << numIters <<
".");
608 if (plist.isParameter (
"eigen-analysis: iterations")) {
609 eigMaxIters = plist.get<
int> (
"eigen-analysis: iterations");
611 eigMaxIters = plist.get (
"chebyshev: eigenvalue max iterations", eigMaxIters);
612 TEUCHOS_TEST_FOR_EXCEPTION(
613 eigMaxIters < 0, std::invalid_argument,
614 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations"
615 "\" parameter (also called \"eigen-analysis: iterations\") must be a "
616 "nonnegative integer. You gave a value of " << eigMaxIters <<
".");
618 if (plist.isType<
double>(
"chebyshev: eigenvalue relative tolerance"))
619 eigRelTolerance = Teuchos::as<MT>(plist.get<
double> (
"chebyshev: eigenvalue relative tolerance"));
620 else if (plist.isType<
MT>(
"chebyshev: eigenvalue relative tolerance"))
621 eigRelTolerance = plist.get<
MT> (
"chebyshev: eigenvalue relative tolerance");
622 else if (plist.isType<
ST>(
"chebyshev: eigenvalue relative tolerance"))
623 eigRelTolerance = Teuchos::ScalarTraits<ST>::magnitude(plist.get<
ST> (
"chebyshev: eigenvalue relative tolerance"));
625 eigKeepVectors = plist.get (
"chebyshev: eigenvalue keep vectors", eigKeepVectors);
627 eigNormalizationFreq = plist.get (
"chebyshev: eigenvalue normalization frequency", eigNormalizationFreq);
628 TEUCHOS_TEST_FOR_EXCEPTION(
629 eigNormalizationFreq < 0, std::invalid_argument,
630 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue normalization frequency"
631 "\" parameter must be a "
632 "nonnegative integer. You gave a value of " << eigNormalizationFreq <<
".")
634 zeroStartingSolution = plist.get (
"chebyshev: zero starting solution",
635 zeroStartingSolution);
636 assumeMatrixUnchanged = plist.get (
"chebyshev: assume matrix does not change",
637 assumeMatrixUnchanged);
641 if (plist.isParameter (
"chebyshev: textbook algorithm")) {
642 textbookAlgorithm = plist.get<
bool> (
"chebyshev: textbook algorithm");
644 if (plist.isParameter (
"chebyshev: compute max residual norm")) {
645 computeMaxResNorm = plist.get<
bool> (
"chebyshev: compute max residual norm");
647 if (plist.isParameter (
"chebyshev: compute spectral radius")) {
648 computeSpectralRadius = plist.get<
bool> (
"chebyshev: compute spectral radius");
654 TEUCHOS_TEST_FOR_EXCEPTION
655 (plist.isType<
bool> (
"chebyshev: use block mode") &&
656 ! plist.get<
bool> (
"chebyshev: use block mode"),
657 std::invalid_argument,
658 "Ifpack2::Chebyshev requires that if you set \"chebyshev: use "
659 "block mode\" at all, you must set it to false. "
660 "Ifpack2::Chebyshev does not implement Ifpack's block mode.");
661 TEUCHOS_TEST_FOR_EXCEPTION
662 (plist.isType<
bool> (
"chebyshev: solve normal equations") &&
663 ! plist.get<
bool> (
"chebyshev: solve normal equations"),
664 std::invalid_argument,
665 "Ifpack2::Chebyshev does not and will never implement the Ifpack "
666 "parameter \"chebyshev: solve normal equations\". If you want to "
667 "solve the normal equations, construct a Tpetra::Operator that "
668 "implements A^* A, and use Chebyshev to solve A^* A x = A^* b.");
676 std::string eigenAnalysisType (
"power-method");
677 if (plist.isParameter (
"eigen-analysis: type")) {
678 eigenAnalysisType = plist.get<std::string> (
"eigen-analysis: type");
679 TEUCHOS_TEST_FOR_EXCEPTION(
680 eigenAnalysisType !=
"power-method" &&
681 eigenAnalysisType !=
"power method" &&
682 eigenAnalysisType !=
"cg",
683 std::invalid_argument,
684 "Ifpack2::Chebyshev: Ifpack2 only supports \"power method\" and \"cg\" for \"eigen-analysis: type\".");
688 userInvDiag_ = userInvDiagCopy;
689 userLambdaMax_ = lambdaMax;
690 userLambdaMin_ = lambdaMin;
691 userEigRatio_ = eigRatio;
692 boostFactor_ =
static_cast<MT> (boostFactor);
693 minDiagVal_ = minDiagVal;
694 numIters_ = numIters;
695 eigMaxIters_ = eigMaxIters;
696 eigRelTolerance_ = eigRelTolerance;
697 eigKeepVectors_ = eigKeepVectors;
698 eigNormalizationFreq_ = eigNormalizationFreq;
699 eigenAnalysisType_ = eigenAnalysisType;
700 zeroStartingSolution_ = zeroStartingSolution;
701 assumeMatrixUnchanged_ = assumeMatrixUnchanged;
702 textbookAlgorithm_ = textbookAlgorithm;
703 computeMaxResNorm_ = computeMaxResNorm;
704 computeSpectralRadius_ = computeSpectralRadius;
705 ckUseNativeSpMV_ = ckUseNativeSpMV;
711 if (A_.is_null () || A_->getComm ().is_null ()) {
717 myRank = A_->getComm ()->getRank ();
721 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
724 using Teuchos::oblackholestream;
725 RCP<oblackholestream> blackHole (
new oblackholestream ());
726 out_ = Teuchos::getFancyOStream (blackHole);
731 out_ = Teuchos::null;
736template<
class ScalarType,
class MV>
742 diagOffsets_ = offsets_type ();
743 savedDiagOffsets_ =
false;
745 computedLambdaMax_ = STS::nan ();
746 computedLambdaMin_ = STS::nan ();
747 eigVector_ = Teuchos::null;
748 eigVector2_ = Teuchos::null;
752template<
class ScalarType,
class MV>
755setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
757 if (A.getRawPtr () != A_.getRawPtr ()) {
758 if (! assumeMatrixUnchanged_) {
770 if (A.is_null () || A->getComm ().is_null ()) {
776 myRank = A->getComm ()->getRank ();
780 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
783 Teuchos::RCP<Teuchos::oblackholestream> blackHole (
new Teuchos::oblackholestream ());
784 out_ = Teuchos::getFancyOStream (blackHole);
789 out_ = Teuchos::null;
794template<
class ScalarType,
class MV>
803 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
804 typename MV::local_ordinal_type,
805 typename MV::global_ordinal_type,
806 typename MV::node_type> crs_matrix_type;
808 TEUCHOS_TEST_FOR_EXCEPTION(
809 A_.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::compute: The input "
810 "matrix A is null. Please call setMatrix() with a nonnull input matrix "
811 "before calling this method.");
826 if (userInvDiag_.is_null ()) {
827 Teuchos::RCP<const crs_matrix_type> A_crsMat =
828 Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
831 if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
833 const size_t lclNumRows = A_crsMat->getLocalNumRows ();
834 if (diagOffsets_.extent (0) < lclNumRows) {
835 diagOffsets_ = offsets_type ();
836 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
838 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
839 savedDiagOffsets_ =
true;
840 D_ = makeInverseDiagonal (*A_,
true);
843 D_ = makeInverseDiagonal (*A_);
846 else if (! assumeMatrixUnchanged_) {
847 if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
850 if (! savedDiagOffsets_) {
851 const size_t lclNumRows = A_crsMat->getLocalNumRows ();
852 if (diagOffsets_.extent (0) < lclNumRows) {
853 diagOffsets_ = offsets_type ();
854 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
856 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
857 savedDiagOffsets_ =
true;
860 D_ = makeInverseDiagonal (*A_,
true);
863 D_ = makeInverseDiagonal (*A_);
868 D_ = makeRangeMapVectorConst (userInvDiag_);
872 const bool computedEigenvalueEstimates =
873 STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
883 if (! assumeMatrixUnchanged_ ||
884 (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
885 ST computedLambdaMax;
886 if ((eigenAnalysisType_ ==
"power method") || (eigenAnalysisType_ ==
"power-method")) {
888 if (eigVector_.is_null()) {
889 x = Teuchos::rcp(
new V(A_->getDomainMap ()));
892 PowerMethod::computeInitialGuessForPowerMethod (*x,
false);
897 if (eigVector2_.is_null()) {
898 y = rcp(
new V(A_->getRangeMap ()));
904 Teuchos::RCP<Teuchos::FancyOStream> stream = (debug_ ? out_ : Teuchos::null);
905 computedLambdaMax = PowerMethod::powerMethodWithInitGuess (*A_, *D_, eigMaxIters_, x, y,
906 eigRelTolerance_, eigNormalizationFreq_, stream,
907 computeSpectralRadius_);
910 computedLambdaMax = cgMethod (*A_, *D_, eigMaxIters_);
911 TEUCHOS_TEST_FOR_EXCEPTION(
912 STS::isnaninf (computedLambdaMax),
914 "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue "
915 "of D^{-1} A failed, by producing Inf or NaN. This probably means that "
916 "the matrix contains Inf or NaN values, or that it is badly scaled.");
917 TEUCHOS_TEST_FOR_EXCEPTION(
918 STS::isnaninf (userEigRatio_),
920 "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
921 << endl <<
"This should be impossible." << endl <<
922 "Please report this bug to the Ifpack2 developers.");
928 const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
931 computedLambdaMax_ = computedLambdaMax;
932 computedLambdaMin_ = computedLambdaMin;
934 TEUCHOS_TEST_FOR_EXCEPTION(
935 STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
937 "Ifpack2::Chebyshev::compute: " << endl <<
938 "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
939 << endl <<
"This should be impossible." << endl <<
940 "Please report this bug to the Ifpack2 developers.");
948 lambdaMaxForApply_ = STS::isnaninf (userLambdaMax_) ? computedLambdaMax_ : userLambdaMax_;
961 lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
962 eigRatioForApply_ = userEigRatio_;
964 if (! textbookAlgorithm_) {
967 const ST one = Teuchos::as<ST> (1);
970 if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
971 lambdaMinForApply_ = one;
972 lambdaMaxForApply_ = lambdaMinForApply_;
973 eigRatioForApply_ = one;
979template<
class ScalarType,
class MV>
983 return lambdaMaxForApply_;
987template<
class ScalarType,
class MV>
991 const char prefix[] =
"Ifpack2::Chebyshev::apply: ";
994 *out_ <<
"apply: " << std::endl;
996 TEUCHOS_TEST_FOR_EXCEPTION
997 (A_.is_null (), std::runtime_error, prefix <<
"The input matrix A is null. "
998 " Please call setMatrix() with a nonnull input matrix, and then call "
999 "compute(), before calling this method.");
1000 TEUCHOS_TEST_FOR_EXCEPTION
1001 (STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
1002 prefix <<
"There is no estimate for the max eigenvalue."
1003 << std::endl << computeBeforeApplyReminder);
1004 TEUCHOS_TEST_FOR_EXCEPTION
1005 (STS::isnaninf (lambdaMinForApply_), std::runtime_error,
1006 prefix <<
"There is no estimate for the min eigenvalue."
1007 << std::endl << computeBeforeApplyReminder);
1008 TEUCHOS_TEST_FOR_EXCEPTION
1009 (STS::isnaninf (eigRatioForApply_), std::runtime_error,
1010 prefix <<
"There is no estimate for the ratio of the max "
1011 "eigenvalue to the min eigenvalue."
1012 << std::endl << computeBeforeApplyReminder);
1013 TEUCHOS_TEST_FOR_EXCEPTION
1014 (D_.is_null (), std::runtime_error, prefix <<
"The vector of inverse "
1015 "diagonal entries of the matrix has not yet been computed."
1016 << std::endl << computeBeforeApplyReminder);
1018 if (textbookAlgorithm_) {
1019 textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1020 lambdaMinForApply_, eigRatioForApply_, *D_);
1023 ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1024 lambdaMinForApply_, eigRatioForApply_, *D_);
1027 if (computeMaxResNorm_ && B.getNumVectors () > 0) {
1028 MV R (B.getMap (), B.getNumVectors ());
1029 computeResidual (R, B, *A_, X);
1030 Teuchos::Array<MT> norms (B.getNumVectors ());
1032 return *std::max_element (norms.begin (), norms.end ());
1035 return Teuchos::ScalarTraits<MT>::zero ();
1039template<
class ScalarType,
class MV>
1042print (std::ostream& out)
1044 using Teuchos::rcpFromRef;
1045 this->describe (* (Teuchos::getFancyOStream (rcpFromRef (out))),
1046 Teuchos::VERB_MEDIUM);
1049template<
class ScalarType,
class MV>
1053 const ScalarType& alpha,
1058 solve (W, alpha, D_inv, B);
1059 Tpetra::deep_copy (X, W);
1062template<
class ScalarType,
class MV>
1066 const Teuchos::ETransp mode)
1068 Tpetra::Details::residual(A,X,B,R);
1071template <
class ScalarType,
class MV>
1073Chebyshev<ScalarType, MV>::
1074solve (MV& Z,
const V& D_inv,
const MV& R) {
1075 Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
1078template<
class ScalarType,
class MV>
1080Chebyshev<ScalarType, MV>::
1081solve (MV& Z,
const ST alpha,
const V& D_inv,
const MV& R) {
1082 Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
1085template<
class ScalarType,
class MV>
1086Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
1087Chebyshev<ScalarType, MV>::
1088makeInverseDiagonal (
const row_matrix_type& A,
const bool useDiagOffsets)
const
1091 using Teuchos::rcpFromRef;
1092 using Teuchos::rcp_dynamic_cast;
1095 if (!D_.is_null() &&
1096 D_->getMap()->isSameAs(*(A.getRowMap ()))) {
1098 *out_ <<
"Reusing pre-existing vector for diagonal extraction" << std::endl;
1099 D_rowMap = Teuchos::rcp_const_cast<V>(D_);
1101 D_rowMap = Teuchos::rcp(
new V (A.getRowMap (),
false));
1103 *out_ <<
"Allocated new vector for diagonal extraction" << std::endl;
1105 if (useDiagOffsets) {
1109 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
1110 typename MV::local_ordinal_type,
1111 typename MV::global_ordinal_type,
1112 typename MV::node_type> crs_matrix_type;
1113 RCP<const crs_matrix_type> A_crsMat =
1114 rcp_dynamic_cast<const crs_matrix_type> (rcpFromRef (A));
1115 if (! A_crsMat.is_null ()) {
1116 TEUCHOS_TEST_FOR_EXCEPTION(
1117 ! savedDiagOffsets_, std::logic_error,
1118 "Ifpack2::Details::Chebyshev::makeInverseDiagonal: "
1119 "It is not allowed to call this method with useDiagOffsets=true, "
1120 "if you have not previously saved offsets of diagonal entries. "
1121 "This situation should never arise if this class is used properly. "
1122 "Please report this bug to the Ifpack2 developers.");
1123 A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_);
1129 A.getLocalDiagCopy (*D_rowMap);
1131 RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
1137 bool foundNonpositiveValue =
false;
1139 auto D_lcl = D_rangeMap->getLocalViewHost (Tpetra::Access::ReadOnly);
1140 auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1142 typedef typename MV::impl_scalar_type IST;
1143 typedef typename MV::local_ordinal_type LO;
1144 typedef Kokkos::Details::ArithTraits<IST> ATS;
1145 typedef Kokkos::Details::ArithTraits<typename ATS::mag_type> STM;
1147 const LO lclNumRows =
static_cast<LO
> (D_rangeMap->getLocalLength ());
1148 for (LO i = 0; i < lclNumRows; ++i) {
1149 if (STS::real (D_lcl_1d(i)) <= STM::zero ()) {
1150 foundNonpositiveValue =
true;
1156 using Teuchos::outArg;
1157 using Teuchos::REDUCE_MIN;
1158 using Teuchos::reduceAll;
1160 const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1161 int gblSuccess = lclSuccess;
1162 if (! D_rangeMap->getMap ().is_null () && D_rangeMap->getMap ()->getComm ().is_null ()) {
1163 const Teuchos::Comm<int>& comm = * (D_rangeMap->getMap ()->getComm ());
1164 reduceAll<int, int> (comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
1166 if (gblSuccess == 1) {
1167 *out_ <<
"makeInverseDiagonal: The matrix's diagonal entries all have "
1168 "positive real part (this is good for Chebyshev)." << std::endl;
1171 *out_ <<
"makeInverseDiagonal: The matrix's diagonal has at least one "
1172 "entry with nonpositive real part, on at least one process in the "
1173 "matrix's communicator. This is bad for Chebyshev." << std::endl;
1179 reciprocal_threshold (*D_rangeMap, minDiagVal_);
1180 return Teuchos::rcp_const_cast<const V> (D_rangeMap);
1184template<
class ScalarType,
class MV>
1185Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
1186Chebyshev<ScalarType, MV>::
1187makeRangeMapVectorConst (
const Teuchos::RCP<const V>& D)
const
1191 typedef Tpetra::Export<
typename MV::local_ordinal_type,
1192 typename MV::global_ordinal_type,
1193 typename MV::node_type> export_type;
1197 TEUCHOS_TEST_FOR_EXCEPTION(
1198 A_.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1199 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() "
1200 "with a nonnull input matrix before calling this method. This is probably "
1201 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1202 TEUCHOS_TEST_FOR_EXCEPTION(
1203 D.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::"
1204 "makeRangeMapVector: The input Vector D is null. This is probably "
1205 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1207 RCP<const map_type> sourceMap = D->getMap ();
1208 RCP<const map_type> rangeMap = A_->getRangeMap ();
1209 RCP<const map_type> rowMap = A_->getRowMap ();
1211 if (rangeMap->isSameAs (*sourceMap)) {
1217 RCP<const export_type> exporter;
1221 if (sourceMap->isSameAs (*rowMap)) {
1223 exporter = A_->getGraph ()->getExporter ();
1226 exporter = rcp (
new export_type (sourceMap, rangeMap));
1229 if (exporter.is_null ()) {
1233 RCP<V> D_out = rcp (
new V (*D, Teuchos::Copy));
1234 D_out->doExport (*D, *exporter, Tpetra::ADD);
1235 return Teuchos::rcp_const_cast<const V> (D_out);
1241template<
class ScalarType,
class MV>
1242Teuchos::RCP<typename Chebyshev<ScalarType, MV>::V>
1243Chebyshev<ScalarType, MV>::
1244makeRangeMapVector (
const Teuchos::RCP<V>& D)
const
1246 using Teuchos::rcp_const_cast;
1247 return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1251template<
class ScalarType,
class MV>
1253Chebyshev<ScalarType, MV>::
1254textbookApplyImpl (
const op_type& A,
1261 const V& D_inv)
const
1264 const ST myLambdaMin = lambdaMax / eigRatio;
1266 const ST zero = Teuchos::as<ST> (0);
1267 const ST one = Teuchos::as<ST> (1);
1268 const ST two = Teuchos::as<ST> (2);
1269 const ST d = (lambdaMax + myLambdaMin) / two;
1270 const ST c = (lambdaMax - myLambdaMin) / two;
1272 if (zeroStartingSolution_ && numIters > 0) {
1276 MV R (B.getMap (), B.getNumVectors (),
false);
1277 MV P (B.getMap (), B.getNumVectors (),
false);
1278 MV Z (B.getMap (), B.getNumVectors (),
false);
1280 for (
int i = 0; i < numIters; ++i) {
1281 computeResidual (R, B, A, X);
1282 solve (Z, D_inv, R);
1290 beta = alpha * (c/two) * (c/two);
1291 alpha = one / (d - beta);
1292 P.update (one, Z, beta);
1294 X.update (alpha, P, one);
1301template<
class ScalarType,
class MV>
1303Chebyshev<ScalarType, MV>::maxNormInf (
const MV& X) {
1304 Teuchos::Array<MT> norms (X.getNumVectors ());
1305 X.normInf (norms());
1306 return *std::max_element (norms.begin (), norms.end ());
1309template<
class ScalarType,
class MV>
1311Chebyshev<ScalarType, MV>::
1312ifpackApplyImpl (
const op_type& A,
1322#ifdef HAVE_IFPACK2_DEBUG
1323 const bool debug = debug_;
1325 const bool debug =
false;
1329 *out_ <<
" \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1330 *out_ <<
" \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1333 if (numIters <= 0) {
1336 const ST zero =
static_cast<ST
> (0.0);
1337 const ST one =
static_cast<ST
> (1.0);
1338 const ST two =
static_cast<ST
> (2.0);
1341 if (lambdaMin == one && lambdaMax == lambdaMin) {
1342 solve (X, D_inv, B);
1347 const ST alpha = lambdaMax / eigRatio;
1348 const ST beta = boostFactor_ * lambdaMax;
1349 const ST delta = two / (beta - alpha);
1350 const ST theta = (beta + alpha) / two;
1351 const ST s1 = theta * delta;
1354 *out_ <<
" alpha = " << alpha << endl
1355 <<
" beta = " << beta << endl
1356 <<
" delta = " << delta << endl
1357 <<
" theta = " << theta << endl
1358 <<
" s1 = " << s1 << endl;
1362 Teuchos::RCP<MV> W_ptr = makeTempMultiVector (B);
1366 *out_ <<
" Iteration " << 1 <<
":" << endl
1367 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1371 if (! zeroStartingSolution_) {
1374 if (ck_.is_null ()) {
1375 Teuchos::RCP<const op_type> A_op = A_;
1376 ck_ = Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1380 ck_->compute (W, one/theta,
const_cast<V&
> (D_inv),
1381 const_cast<MV&
> (B), X, zero);
1385 firstIterationWithZeroStartingSolution (W, one/theta, D_inv, B, X);
1389 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1390 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1393 if (numIters > 1 && ck_.is_null ()) {
1394 Teuchos::RCP<const op_type> A_op = A_;
1395 ck_ = Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op, ckUseNativeSpMV_));
1400 ST rhokp1, dtemp1, dtemp2;
1401 for (
int deg = 1; deg < numIters; ++deg) {
1403 *out_ <<
" Iteration " << deg+1 <<
":" << endl
1404 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl
1405 <<
" - \\|B\\|_{\\infty} = " << maxNormInf (B) << endl
1406 <<
" - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm ()
1407 << endl <<
" - rhok = " << rhok << endl;
1410 rhokp1 = one / (two * s1 - rhok);
1411 dtemp1 = rhokp1 * rhok;
1412 dtemp2 = two * rhokp1 * delta;
1416 *out_ <<
" - dtemp1 = " << dtemp1 << endl
1417 <<
" - dtemp2 = " << dtemp2 << endl;
1422 ck_->compute (W, dtemp2,
const_cast<V&
> (D_inv),
1423 const_cast<MV&
> (B), (X), dtemp1);
1426 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1427 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1433template<
class ScalarType,
class MV>
1435Chebyshev<ScalarType, MV>::
1436cgMethodWithInitGuess (
const op_type& A,
1442 using MagnitudeType =
typename STS::magnitudeType;
1444 *out_ <<
" cgMethodWithInitGuess:" << endl;
1447 const ST one = STS::one();
1448 ST beta, betaOld = one, pAp, pApOld = one, alpha, rHz, rHzOld, rHzOld2 = one, lambdaMax;
1450 Teuchos::ArrayRCP<MagnitudeType> diag, offdiag;
1451 Teuchos::RCP<V> p, z, Ap;
1452 diag.resize(numIters);
1453 offdiag.resize(numIters-1);
1455 p = rcp(
new V(A.getRangeMap ()));
1456 z = rcp(
new V(A.getRangeMap ()));
1457 Ap = rcp(
new V(A.getRangeMap ()));
1460 solve (*p, D_inv, r);
1463 for (
int iter = 0; iter < numIters; ++iter) {
1465 *out_ <<
" Iteration " << iter << endl;
1470 r.update (-alpha, *Ap, one);
1472 solve (*z, D_inv, r);
1474 beta = rHz / rHzOld;
1475 p->update (one, *z, beta);
1477 diag[iter] = STS::real((betaOld*betaOld * pApOld + pAp) / rHzOld);
1478 offdiag[iter-1] = -STS::real((betaOld * pApOld / (sqrt(rHzOld * rHzOld2))));
1480 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1481 *out_ <<
" offdiag["<< iter-1 <<
"] = " << offdiag[iter-1] << endl;
1484 diag[iter] = STS::real(pAp/rHzOld);
1486 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1494 lambdaMax = LapackHelper<ST>::tri_diag_spectral_radius(diag, offdiag);
1500template<
class ScalarType,
class MV>
1502Chebyshev<ScalarType, MV>::
1503cgMethod (
const op_type& A,
const V& D_inv,
const int numIters)
1507 *out_ <<
"cgMethod:" << endl;
1511 if (eigVector_.is_null()) {
1512 r = rcp(
new V(A.getDomainMap ()));
1513 if (eigKeepVectors_)
1518 PowerMethod::computeInitialGuessForPowerMethod (*r,
false);
1522 ST lambdaMax = cgMethodWithInitGuess (A, D_inv, numIters, *r);
1527template<
class ScalarType,
class MV>
1528Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::row_matrix_type>
1533template<
class ScalarType,
class MV>
1541template<
class ScalarType,
class MV>
1550 const size_t B_numVecs = B.getNumVectors ();
1551 if (W_.is_null () || W_->getNumVectors () != B_numVecs) {
1552 W_ = Teuchos::rcp (
new MV (B.getMap (), B_numVecs,
false));
1557template<
class ScalarType,
class MV>
1561 std::ostringstream oss;
1564 oss <<
"\"Ifpack2::Details::Chebyshev\":"
1566 <<
"degree: " << numIters_
1567 <<
", lambdaMax: " << lambdaMaxForApply_
1568 <<
", alpha: " << eigRatioForApply_
1569 <<
", lambdaMin: " << lambdaMinForApply_
1570 <<
", boost factor: " << boostFactor_;
1571 if (!userInvDiag_.is_null())
1572 oss <<
", diagonal: user-supplied";
1577template<
class ScalarType,
class MV>
1580describe (Teuchos::FancyOStream& out,
1581 const Teuchos::EVerbosityLevel verbLevel)
const
1583 using Teuchos::TypeNameTraits;
1586 const Teuchos::EVerbosityLevel vl =
1587 (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1588 if (vl == Teuchos::VERB_NONE) {
1598 Teuchos::OSTab tab0 (out);
1604 if (A_.is_null () || A_->getComm ().is_null ()) {
1608 myRank = A_->getComm ()->getRank ();
1613 out <<
"\"Ifpack2::Details::Chebyshev\":" << endl;
1615 Teuchos::OSTab tab1 (out);
1617 if (vl == Teuchos::VERB_LOW) {
1619 out <<
"degree: " << numIters_ << endl
1620 <<
"lambdaMax: " << lambdaMaxForApply_ << endl
1621 <<
"alpha: " << eigRatioForApply_ << endl
1622 <<
"lambdaMin: " << lambdaMinForApply_ << endl
1623 <<
"boost factor: " << boostFactor_ << endl;
1631 out <<
"Template parameters:" << endl;
1633 Teuchos::OSTab tab2 (out);
1634 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1635 <<
"MV: " << TypeNameTraits<MV>::name () << endl;
1641 out <<
"Computed parameters:" << endl;
1647 Teuchos::OSTab tab2 (out);
1653 if (D_.is_null ()) {
1655 out <<
"unset" << endl;
1658 else if (vl <= Teuchos::VERB_HIGH) {
1660 out <<
"set" << endl;
1669 D_->describe (out, vl);
1674 out <<
"W_: " << (W_.is_null () ?
"unset" :
"set") << endl
1675 <<
"computedLambdaMax_: " << computedLambdaMax_ << endl
1676 <<
"computedLambdaMin_: " << computedLambdaMin_ << endl
1677 <<
"lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1678 <<
"lambdaMinForApply_: " << lambdaMinForApply_ << endl
1679 <<
"eigRatioForApply_: " << eigRatioForApply_ << endl;
1684 out <<
"User parameters:" << endl;
1689 Teuchos::OSTab tab2 (out);
1690 out <<
"userInvDiag_: ";
1691 if (userInvDiag_.is_null ()) {
1692 out <<
"unset" << endl;
1694 else if (vl <= Teuchos::VERB_HIGH) {
1695 out <<
"set" << endl;
1701 userInvDiag_->describe (out, vl);
1704 out <<
"userLambdaMax_: " << userLambdaMax_ << endl
1705 <<
"userLambdaMin_: " << userLambdaMin_ << endl
1706 <<
"userEigRatio_: " << userEigRatio_ << endl
1707 <<
"numIters_: " << numIters_ << endl
1708 <<
"eigMaxIters_: " << eigMaxIters_ << endl
1709 <<
"eigRelTolerance_: " << eigRelTolerance_ << endl
1710 <<
"eigNormalizationFreq_: " << eigNormalizationFreq_ << endl
1711 <<
"zeroStartingSolution_: " << zeroStartingSolution_ << endl
1712 <<
"assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1713 <<
"textbookAlgorithm_: " << textbookAlgorithm_ << endl
1714 <<
"computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1722#define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \
1723 template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >;
Declaration of Chebyshev implementation.
Definition of power methods.
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:208
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:106
Tpetra::Vector< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > V
Specialization of Tpetra::Vector.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:131
Teuchos::ScalarTraits< scalar_type > STS
Traits class for ST.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:114
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:796
void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:353
void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set the matrix.
Definition: Ifpack2_Details_Chebyshev_def.hpp:755
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1560
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:287
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:1580
void print(std::ostream &out)
Print instance data to the given output stream.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1042
ScalarType ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:112
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:116
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:989
bool hasTransposeApply() const
Whether it's possible to apply the transpose of this operator.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1536
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the matrix given to the constructor.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1529
Ifpack2 implementation details.
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74