44 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP 45 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP 56 #include "Ifpack2_Details_ChebyshevKernel.hpp" 57 #include "Kokkos_ArithTraits.hpp" 58 #include "Teuchos_FancyOStream.hpp" 59 #include "Teuchos_oblackholestream.hpp" 60 #include "Tpetra_Details_residual.hpp" 61 #include "Teuchos_LAPACK.hpp" 62 #include "Ifpack2_Details_LapackSupportsScalar.hpp" 72 const char computeBeforeApplyReminder[] =
73 "This means one of the following:\n" 74 " - you have not yet called compute() on this instance, or \n" 75 " - you didn't call compute() after calling setParameters().\n\n" 76 "After creating an Ifpack2::Chebyshev instance,\n" 77 "you must _always_ call compute() at least once before calling apply().\n" 78 "After calling compute() once, you do not need to call it again,\n" 79 "unless the matrix has changed or you have changed parameters\n" 80 "(by calling setParameters()).";
86 template<
class XV,
class SizeType =
typename XV::
size_type>
87 struct V_ReciprocalThresholdSelfFunctor
89 typedef typename XV::execution_space execution_space;
90 typedef typename XV::non_const_value_type value_type;
91 typedef SizeType size_type;
92 typedef Kokkos::Details::ArithTraits<value_type> KAT;
93 typedef typename KAT::mag_type mag_type;
96 const value_type minVal_;
97 const mag_type minValMag_;
99 V_ReciprocalThresholdSelfFunctor (
const XV& X,
100 const value_type& min_val) :
103 minValMag_ (KAT::abs (min_val))
106 KOKKOS_INLINE_FUNCTION
107 void operator() (
const size_type& i)
const 109 const mag_type X_i_abs = KAT::abs (X_[i]);
111 if (X_i_abs < minValMag_) {
115 X_[i] = KAT::one () / X_[i];
120 template<
class XV,
class SizeType =
typename XV::
size_type>
121 struct LocalReciprocalThreshold {
123 compute (
const XV& X,
124 const typename XV::non_const_value_type& minVal)
126 typedef typename XV::execution_space execution_space;
127 Kokkos::RangePolicy<execution_space, SizeType> policy (0, X.extent (0));
128 V_ReciprocalThresholdSelfFunctor<XV, SizeType> op (X, minVal);
129 Kokkos::parallel_for (policy, op);
133 template <
class TpetraVectorType,
134 const bool classic = TpetraVectorType::node_type::classic>
135 struct GlobalReciprocalThreshold {};
137 template <
class TpetraVectorType>
138 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
140 compute (TpetraVectorType& V,
141 const typename TpetraVectorType::scalar_type& min_val)
143 typedef typename TpetraVectorType::scalar_type scalar_type;
144 typedef typename TpetraVectorType::mag_type mag_type;
145 typedef Kokkos::Details::ArithTraits<scalar_type> STS;
147 const scalar_type ONE = STS::one ();
148 const mag_type min_val_abs = STS::abs (min_val);
150 Teuchos::ArrayRCP<scalar_type> V_0 = V.getDataNonConst (0);
151 const size_t lclNumRows = V.getLocalLength ();
153 for (
size_t i = 0; i < lclNumRows; ++i) {
154 const scalar_type V_0i = V_0[i];
155 if (STS::abs (V_0i) < min_val_abs) {
164 template <
class TpetraVectorType>
165 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
167 compute (TpetraVectorType& X,
168 const typename TpetraVectorType::scalar_type& minVal)
170 typedef typename TpetraVectorType::impl_scalar_type value_type;
172 const value_type minValS =
static_cast<value_type
> (minVal);
173 auto X_0 = Kokkos::subview (X.getLocalViewDevice (Tpetra::Access::ReadWrite),
175 LocalReciprocalThreshold<decltype (X_0) >::compute (X_0, minValS);
180 template <
typename S,
typename L,
typename G,
typename N>
182 reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V,
const S& minVal)
184 GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
192 template<
class OneDViewType,
193 class LocalOrdinal =
typename OneDViewType::size_type>
194 class PositivizeVector {
195 static_assert (Kokkos::Impl::is_view<OneDViewType>::value,
196 "OneDViewType must be a 1-D Kokkos::View.");
197 static_assert (static_cast<int> (OneDViewType::rank) == 1,
198 "This functor only works with a 1-D View.");
199 static_assert (std::is_integral<LocalOrdinal>::value,
200 "The second template parameter, LocalOrdinal, " 201 "must be an integer type.");
203 PositivizeVector (
const OneDViewType& x) : x_ (x) {}
205 KOKKOS_INLINE_FUNCTION
void 206 operator () (
const LocalOrdinal& i)
const 208 typedef typename OneDViewType::non_const_value_type IST;
209 typedef Kokkos::Details::ArithTraits<IST> STS;
210 typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
212 if (STS::real (x_(i)) < STM::zero ()) {
224 template<class ScalarType, const bool lapackSupportsScalarType = LapackSupportsScalar<ScalarType>::value>
228 tri_diag_spectral_radius(Teuchos::ArrayRCP<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
229 Teuchos::ArrayRCP<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> offdiag) {
230 throw std::runtime_error(
"LAPACK does not support the scalar type.");
235 template<
class ScalarType>
236 struct LapackHelper<ScalarType,true> {
239 tri_diag_spectral_radius(Teuchos::ArrayRCP<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
240 Teuchos::ArrayRCP<
typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> offdiag) {
241 using STS = Teuchos::ScalarTraits<ScalarType>;
242 using MagnitudeType =
typename STS::magnitudeType;
244 const int N = diag.size ();
245 ScalarType scalar_dummy;
246 std::vector<MagnitudeType> mag_dummy(4*N);
250 ScalarType lambdaMax = STS::one();
252 Teuchos::LAPACK<int,ScalarType> lapack;
253 lapack.PTEQR (char_N, N, diag.getRawPtr (), offdiag.getRawPtr (),
254 &scalar_dummy, 1, &mag_dummy[0], &info);
255 TEUCHOS_TEST_FOR_EXCEPTION
256 (info < 0, std::logic_error,
"Ifpack2::Details::LapackHelper::tri_diag_spectral_radius:" 257 "LAPACK's _PTEQR failed with info = " 258 << info <<
" < 0. This suggests there might be a bug in the way Ifpack2 " 259 "is calling LAPACK. Please report this to the Ifpack2 developers.");
261 lambdaMax = Teuchos::as<ScalarType> (diag[0]);
268 template<
class ScalarType,
class MV>
269 void Chebyshev<ScalarType, MV>::checkInputMatrix ()
const 271 TEUCHOS_TEST_FOR_EXCEPTION(
272 ! A_.is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
273 std::invalid_argument,
274 "Ifpack2::Chebyshev: The input matrix A must be square. " 275 "A has " << A_->getGlobalNumRows () <<
" rows and " 276 << A_->getGlobalNumCols () <<
" columns.");
280 if (debug_ && ! A_.is_null ()) {
281 Teuchos::RCP<const map_type> domainMap = A_->getDomainMap ();
282 Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap ();
287 TEUCHOS_TEST_FOR_EXCEPTION(
288 ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
289 "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be " 290 "the same (in the sense of isSameAs())." << std::endl <<
"We only check " 291 "for this in debug mode.");
296 template<
class ScalarType,
class MV>
298 Chebyshev<ScalarType, MV>::
299 checkConstructorInput ()
const 304 if (STS::isComplex) {
305 TEUCHOS_TEST_FOR_EXCEPTION
306 (
true, std::logic_error,
"Ifpack2::Chebyshev: This class' implementation " 307 "of Chebyshev iteration only works for real-valued, symmetric positive " 308 "definite matrices. However, you instantiated this class for ScalarType" 309 " = " << Teuchos::TypeNameTraits<ScalarType>::name () <<
", which is a " 310 "complex-valued type. While this may be algorithmically correct if all " 311 "of the complex numbers in the matrix have zero imaginary part, we " 312 "forbid using complex ScalarType altogether in order to remind you of " 313 "the limitations of our implementation (and of the algorithm itself).");
320 template<
class ScalarType,
class MV>
324 savedDiagOffsets_ (false),
325 computedLambdaMax_ (
STS::nan ()),
326 computedLambdaMin_ (
STS::nan ()),
327 lambdaMaxForApply_ (
STS::nan ()),
328 lambdaMinForApply_ (
STS::nan ()),
329 eigRatioForApply_ (
STS::nan ()),
330 userLambdaMax_ (
STS::nan ()),
331 userLambdaMin_ (
STS::nan ()),
333 minDiagVal_ (
STS::eps ()),
336 eigRelTolerance_(
Teuchos::ScalarTraits<
MT>::zero ()),
337 eigKeepVectors_(false),
338 eigenAnalysisType_(
"power method"),
339 eigNormalizationFreq_(1),
340 zeroStartingSolution_ (true),
341 assumeMatrixUnchanged_ (false),
342 textbookAlgorithm_ (false),
343 computeMaxResNorm_ (false),
346 checkConstructorInput ();
349 template<
class ScalarType,
class MV>
352 Teuchos::ParameterList& params) :
354 savedDiagOffsets_ (false),
355 computedLambdaMax_ (
STS::nan ()),
356 computedLambdaMin_ (
STS::nan ()),
357 lambdaMaxForApply_ (
STS::nan ()),
358 boostFactor_ (static_cast<
MT> (1.1)),
359 lambdaMinForApply_ (
STS::nan ()),
360 eigRatioForApply_ (
STS::nan ()),
361 userLambdaMax_ (
STS::nan ()),
362 userLambdaMin_ (
STS::nan ()),
364 minDiagVal_ (
STS::eps ()),
367 eigRelTolerance_(
Teuchos::ScalarTraits<
MT>::zero ()),
368 eigKeepVectors_(false),
369 eigenAnalysisType_(
"power method"),
370 eigNormalizationFreq_(1),
371 zeroStartingSolution_ (true),
372 assumeMatrixUnchanged_ (false),
373 textbookAlgorithm_ (false),
374 computeMaxResNorm_ (false),
377 checkConstructorInput ();
378 setParameters (params);
381 template<
class ScalarType,
class MV>
388 using Teuchos::rcp_const_cast;
400 const ST defaultLambdaMax = STS::nan ();
401 const ST defaultLambdaMin = STS::nan ();
408 const ST defaultEigRatio = Teuchos::as<ST> (30);
409 const MT defaultBoostFactor =
static_cast<MT> (1.1);
410 const ST defaultMinDiagVal = STS::eps ();
411 const int defaultNumIters = 1;
412 const int defaultEigMaxIters = 10;
413 const MT defaultEigRelTolerance = Teuchos::ScalarTraits<MT>::zero ();
414 const bool defaultEigKeepVectors =
false;
415 const int defaultEigNormalizationFreq = 1;
416 const bool defaultZeroStartingSolution =
true;
417 const bool defaultAssumeMatrixUnchanged =
false;
418 const bool defaultTextbookAlgorithm =
false;
419 const bool defaultComputeMaxResNorm =
false;
420 const bool defaultDebug =
false;
426 RCP<const V> userInvDiagCopy;
427 ST lambdaMax = defaultLambdaMax;
428 ST lambdaMin = defaultLambdaMin;
429 ST eigRatio = defaultEigRatio;
430 MT boostFactor = defaultBoostFactor;
431 ST minDiagVal = defaultMinDiagVal;
432 int numIters = defaultNumIters;
433 int eigMaxIters = defaultEigMaxIters;
434 MT eigRelTolerance = defaultEigRelTolerance;
435 bool eigKeepVectors = defaultEigKeepVectors;
436 int eigNormalizationFreq = defaultEigNormalizationFreq;
437 bool zeroStartingSolution = defaultZeroStartingSolution;
438 bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
439 bool textbookAlgorithm = defaultTextbookAlgorithm;
440 bool computeMaxResNorm = defaultComputeMaxResNorm;
441 bool debug = defaultDebug;
448 if (plist.isType<
bool> (
"debug")) {
449 debug = plist.get<
bool> (
"debug");
451 else if (plist.isType<
int> (
"debug")) {
452 const int debugInt = plist.get<
bool> (
"debug");
453 debug = debugInt != 0;
466 const char opInvDiagLabel[] =
"chebyshev: operator inv diagonal";
467 if (plist.isParameter (opInvDiagLabel)) {
469 RCP<const V> userInvDiag;
471 if (plist.isType<
const V*> (opInvDiagLabel)) {
472 const V* rawUserInvDiag =
473 plist.get<
const V*> (opInvDiagLabel);
475 userInvDiag = rcp (rawUserInvDiag,
false);
477 else if (plist.isType<
const V*> (opInvDiagLabel)) {
478 V* rawUserInvDiag = plist.get<
V*> (opInvDiagLabel);
480 userInvDiag = rcp (const_cast<const V*> (rawUserInvDiag),
false);
482 else if (plist.isType<RCP<const V>> (opInvDiagLabel)) {
483 userInvDiag = plist.get<RCP<const V> > (opInvDiagLabel);
485 else if (plist.isType<RCP<V>> (opInvDiagLabel)) {
486 RCP<V> userInvDiagNonConst =
487 plist.get<RCP<V> > (opInvDiagLabel);
488 userInvDiag = rcp_const_cast<
const V> (userInvDiagNonConst);
490 else if (plist.isType<
const V> (opInvDiagLabel)) {
491 const V& userInvDiagRef = plist.get<
const V> (opInvDiagLabel);
492 userInvDiagCopy = rcp (
new V (userInvDiagRef, Teuchos::Copy));
493 userInvDiag = userInvDiagCopy;
495 else if (plist.isType<
V> (opInvDiagLabel)) {
496 V& userInvDiagNonConstRef = plist.get<
V> (opInvDiagLabel);
497 const V& userInvDiagRef =
const_cast<const V&
> (userInvDiagNonConstRef);
498 userInvDiagCopy = rcp (
new V (userInvDiagRef, Teuchos::Copy));
499 userInvDiag = userInvDiagCopy;
509 if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
510 userInvDiagCopy = rcp (
new V (*userInvDiag, Teuchos::Copy));
522 if (plist.isParameter (
"chebyshev: max eigenvalue")) {
523 if (plist.isType<
double>(
"chebyshev: max eigenvalue"))
524 lambdaMax = plist.get<
double> (
"chebyshev: max eigenvalue");
526 lambdaMax = plist.get<
ST> (
"chebyshev: max eigenvalue");
527 TEUCHOS_TEST_FOR_EXCEPTION(
528 STS::isnaninf (lambdaMax), std::invalid_argument,
529 "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" " 530 "parameter is NaN or Inf. This parameter is optional, but if you " 531 "choose to supply it, it must have a finite value.");
533 if (plist.isParameter (
"chebyshev: min eigenvalue")) {
534 if (plist.isType<
double>(
"chebyshev: min eigenvalue"))
535 lambdaMin = plist.get<
double> (
"chebyshev: min eigenvalue");
537 lambdaMin = plist.get<
ST> (
"chebyshev: min eigenvalue");
538 TEUCHOS_TEST_FOR_EXCEPTION(
539 STS::isnaninf (lambdaMin), std::invalid_argument,
540 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" " 541 "parameter is NaN or Inf. This parameter is optional, but if you " 542 "choose to supply it, it must have a finite value.");
546 if (plist.isParameter (
"smoother: Chebyshev alpha")) {
547 if (plist.isType<
double>(
"smoother: Chebyshev alpha"))
548 eigRatio = plist.get<
double> (
"smoother: Chebyshev alpha");
550 eigRatio = plist.get<
ST> (
"smoother: Chebyshev alpha");
553 eigRatio = plist.get (
"chebyshev: ratio eigenvalue", eigRatio);
554 TEUCHOS_TEST_FOR_EXCEPTION(
555 STS::isnaninf (eigRatio), std::invalid_argument,
556 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" " 557 "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. " 558 "This parameter is optional, but if you choose to supply it, it must have " 565 TEUCHOS_TEST_FOR_EXCEPTION(
566 STS::real (eigRatio) < STS::real (STS::one ()),
567 std::invalid_argument,
568 "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\"" 569 "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, " 570 "but you supplied the value " << eigRatio <<
".");
575 const char paramName[] =
"chebyshev: boost factor";
577 if (plist.isParameter (paramName)) {
578 if (plist.isType<
MT> (paramName)) {
579 boostFactor = plist.get<
MT> (paramName);
581 else if (! std::is_same<double, MT>::value &&
582 plist.isType<
double> (paramName)) {
583 const double dblBF = plist.get<
double> (paramName);
584 boostFactor =
static_cast<MT> (dblBF);
587 TEUCHOS_TEST_FOR_EXCEPTION
588 (
true, std::invalid_argument,
589 "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\"" 590 "parameter must have type magnitude_type (MT) or double.");
600 plist.set (paramName, defaultBoostFactor);
602 TEUCHOS_TEST_FOR_EXCEPTION
603 (boostFactor < Teuchos::ScalarTraits<MT>::one (), std::invalid_argument,
604 "Ifpack2::Chebyshev::setParameters: \"" << paramName <<
"\" parameter " 605 "must be >= 1, but you supplied the value " << boostFactor <<
".");
609 minDiagVal = plist.get (
"chebyshev: min diagonal value", minDiagVal);
610 TEUCHOS_TEST_FOR_EXCEPTION(
611 STS::isnaninf (minDiagVal), std::invalid_argument,
612 "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" " 613 "parameter is NaN or Inf. This parameter is optional, but if you choose " 614 "to supply it, it must have a finite value.");
617 if (plist.isParameter (
"smoother: sweeps")) {
618 numIters = plist.get<
int> (
"smoother: sweeps");
620 if (plist.isParameter (
"relaxation: sweeps")) {
621 numIters = plist.get<
int> (
"relaxation: sweeps");
623 numIters = plist.get (
"chebyshev: degree", numIters);
624 TEUCHOS_TEST_FOR_EXCEPTION(
625 numIters < 0, std::invalid_argument,
626 "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also " 627 "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a " 628 "nonnegative integer. You gave a value of " << numIters <<
".");
631 if (plist.isParameter (
"eigen-analysis: iterations")) {
632 eigMaxIters = plist.get<
int> (
"eigen-analysis: iterations");
634 eigMaxIters = plist.get (
"chebyshev: eigenvalue max iterations", eigMaxIters);
635 TEUCHOS_TEST_FOR_EXCEPTION(
636 eigMaxIters < 0, std::invalid_argument,
637 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations" 638 "\" parameter (also called \"eigen-analysis: iterations\") must be a " 639 "nonnegative integer. You gave a value of " << eigMaxIters <<
".");
641 if (plist.isType<
double>(
"chebyshev: eigenvalue relative tolerance"))
642 eigRelTolerance = Teuchos::as<MT>(plist.get<
double> (
"chebyshev: eigenvalue relative tolerance"));
643 else if (plist.isType<
MT>(
"chebyshev: eigenvalue relative tolerance"))
644 eigRelTolerance = plist.get<
MT> (
"chebyshev: eigenvalue relative tolerance");
645 else if (plist.isType<
ST>(
"chebyshev: eigenvalue relative tolerance"))
646 eigRelTolerance = Teuchos::ScalarTraits<ST>::magnitude(plist.get<
ST> (
"chebyshev: eigenvalue relative tolerance"));
648 eigKeepVectors = plist.get (
"chebyshev: eigenvalue keep vectors", eigKeepVectors);
650 eigNormalizationFreq = plist.get (
"chebyshev: eigenvalue normalization frequency", eigNormalizationFreq);
651 TEUCHOS_TEST_FOR_EXCEPTION(
652 eigNormalizationFreq < 0, std::invalid_argument,
653 "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue normalization frequency" 654 "\" parameter must be a " 655 "nonnegative integer. You gave a value of " << eigNormalizationFreq <<
".")
657 zeroStartingSolution = plist.get (
"chebyshev: zero starting solution",
658 zeroStartingSolution);
659 assumeMatrixUnchanged = plist.get (
"chebyshev: assume matrix does not change",
660 assumeMatrixUnchanged);
664 if (plist.isParameter (
"chebyshev: textbook algorithm")) {
665 textbookAlgorithm = plist.get<
bool> (
"chebyshev: textbook algorithm");
667 if (plist.isParameter (
"chebyshev: compute max residual norm")) {
668 computeMaxResNorm = plist.get<
bool> (
"chebyshev: compute max residual norm");
674 TEUCHOS_TEST_FOR_EXCEPTION
675 (plist.isType<
bool> (
"chebyshev: use block mode") &&
676 ! plist.get<
bool> (
"chebyshev: use block mode"),
677 std::invalid_argument,
678 "Ifpack2::Chebyshev requires that if you set \"chebyshev: use " 679 "block mode\" at all, you must set it to false. " 680 "Ifpack2::Chebyshev does not implement Ifpack's block mode.");
681 TEUCHOS_TEST_FOR_EXCEPTION
682 (plist.isType<
bool> (
"chebyshev: solve normal equations") &&
683 ! plist.get<
bool> (
"chebyshev: solve normal equations"),
684 std::invalid_argument,
685 "Ifpack2::Chebyshev does not and will never implement the Ifpack " 686 "parameter \"chebyshev: solve normal equations\". If you want to " 687 "solve the normal equations, construct a Tpetra::Operator that " 688 "implements A^* A, and use Chebyshev to solve A^* A x = A^* b.");
696 std::string eigenAnalysisType (
"power-method");
697 if (plist.isParameter (
"eigen-analysis: type")) {
698 eigenAnalysisType = plist.get<std::string> (
"eigen-analysis: type");
699 TEUCHOS_TEST_FOR_EXCEPTION(
700 eigenAnalysisType !=
"power-method" &&
701 eigenAnalysisType !=
"power method" &&
702 eigenAnalysisType !=
"cg",
703 std::invalid_argument,
704 "Ifpack2::Chebyshev: Ifpack2 only supports \"power method\" and \"cg\" for \"eigen-analysis: type\".");
708 userInvDiag_ = userInvDiagCopy;
709 userLambdaMax_ = lambdaMax;
710 userLambdaMin_ = lambdaMin;
711 userEigRatio_ = eigRatio;
712 boostFactor_ =
static_cast<MT> (boostFactor);
713 minDiagVal_ = minDiagVal;
714 numIters_ = numIters;
715 eigMaxIters_ = eigMaxIters;
716 eigRelTolerance_ = eigRelTolerance;
717 eigKeepVectors_ = eigKeepVectors;
718 eigNormalizationFreq_ = eigNormalizationFreq;
719 eigenAnalysisType_ = eigenAnalysisType;
720 zeroStartingSolution_ = zeroStartingSolution;
721 assumeMatrixUnchanged_ = assumeMatrixUnchanged;
722 textbookAlgorithm_ = textbookAlgorithm;
723 computeMaxResNorm_ = computeMaxResNorm;
729 if (A_.is_null () || A_->getComm ().is_null ()) {
735 myRank = A_->getComm ()->getRank ();
739 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
742 using Teuchos::oblackholestream;
743 RCP<oblackholestream> blackHole (
new oblackholestream ());
744 out_ = Teuchos::getFancyOStream (blackHole);
749 out_ = Teuchos::null;
754 template<
class ScalarType,
class MV>
760 diagOffsets_ = offsets_type ();
761 savedDiagOffsets_ =
false;
763 computedLambdaMax_ = STS::nan ();
764 computedLambdaMin_ = STS::nan ();
765 eigVector_ = Teuchos::null;
766 eigVector2_ = Teuchos::null;
770 template<
class ScalarType,
class MV>
773 setMatrix (
const Teuchos::RCP<const row_matrix_type>& A)
775 if (A.getRawPtr () != A_.getRawPtr ()) {
776 if (! assumeMatrixUnchanged_) {
788 if (A.is_null () || A->getComm ().is_null ()) {
794 myRank = A->getComm ()->getRank ();
798 out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
801 Teuchos::RCP<Teuchos::oblackholestream> blackHole (
new Teuchos::oblackholestream ());
802 out_ = Teuchos::getFancyOStream (blackHole);
807 out_ = Teuchos::null;
812 template<
class ScalarType,
class MV>
821 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
822 typename MV::local_ordinal_type,
823 typename MV::global_ordinal_type,
824 typename MV::node_type> crs_matrix_type;
826 TEUCHOS_TEST_FOR_EXCEPTION(
827 A_.is_null (), std::runtime_error,
"Ifpack2::Chebyshev::compute: The input " 828 "matrix A is null. Please call setMatrix() with a nonnull input matrix " 829 "before calling this method.");
844 if (userInvDiag_.is_null ()) {
845 Teuchos::RCP<const crs_matrix_type> A_crsMat =
846 Teuchos::rcp_dynamic_cast<
const crs_matrix_type> (A_);
849 if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
851 const size_t lclNumRows = A_crsMat->getNodeNumRows ();
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;
858 D_ = makeInverseDiagonal (*A_,
true);
861 D_ = makeInverseDiagonal (*A_);
864 else if (! assumeMatrixUnchanged_) {
865 if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
868 if (! savedDiagOffsets_) {
869 const size_t lclNumRows = A_crsMat->getNodeNumRows ();
870 if (diagOffsets_.extent (0) < lclNumRows) {
871 diagOffsets_ = offsets_type ();
872 diagOffsets_ = offsets_type (
"offsets", lclNumRows);
874 A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
875 savedDiagOffsets_ =
true;
878 D_ = makeInverseDiagonal (*A_,
true);
881 D_ = makeInverseDiagonal (*A_);
886 D_ = makeRangeMapVectorConst (userInvDiag_);
890 const bool computedEigenvalueEstimates =
891 STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
900 if (! assumeMatrixUnchanged_ ||
901 (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
902 ST computedLambdaMax;
903 if ((eigenAnalysisType_ ==
"power method") || (eigenAnalysisType_ ==
"power-method"))
904 computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
906 computedLambdaMax = cgMethod (*A_, *D_, eigMaxIters_);
907 TEUCHOS_TEST_FOR_EXCEPTION(
908 STS::isnaninf (computedLambdaMax),
910 "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue " 911 "of D^{-1} A failed, by producing Inf or NaN. This probably means that " 912 "the matrix contains Inf or NaN values, or that it is badly scaled.");
913 TEUCHOS_TEST_FOR_EXCEPTION(
914 STS::isnaninf (userEigRatio_),
916 "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN." 917 << endl <<
"This should be impossible." << endl <<
918 "Please report this bug to the Ifpack2 developers.");
924 const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
927 computedLambdaMax_ = computedLambdaMax;
928 computedLambdaMin_ = computedLambdaMin;
930 TEUCHOS_TEST_FOR_EXCEPTION(
931 STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
933 "Ifpack2::Chebyshev::compute: " << endl <<
934 "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN." 935 << endl <<
"This should be impossible." << endl <<
936 "Please report this bug to the Ifpack2 developers.");
944 lambdaMaxForApply_ = STS::isnaninf (userLambdaMax_) ? computedLambdaMax_ : userLambdaMax_;
957 lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
958 eigRatioForApply_ = userEigRatio_;
960 if (! textbookAlgorithm_) {
963 const ST one = Teuchos::as<ST> (1);
966 if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
967 lambdaMinForApply_ = one;
968 lambdaMaxForApply_ = lambdaMinForApply_;
969 eigRatioForApply_ = one;
975 template<
class ScalarType,
class MV>
979 return lambdaMaxForApply_;
983 template<
class ScalarType,
class MV>
987 const char prefix[] =
"Ifpack2::Chebyshev::apply: ";
990 *out_ <<
"apply: " << std::endl;
992 TEUCHOS_TEST_FOR_EXCEPTION
993 (A_.is_null (), std::runtime_error, prefix <<
"The input matrix A is null. " 994 " Please call setMatrix() with a nonnull input matrix, and then call " 995 "compute(), before calling this method.");
996 TEUCHOS_TEST_FOR_EXCEPTION
997 (STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
998 prefix <<
"There is no estimate for the max eigenvalue." 999 << std::endl << computeBeforeApplyReminder);
1000 TEUCHOS_TEST_FOR_EXCEPTION
1001 (STS::isnaninf (lambdaMinForApply_), std::runtime_error,
1002 prefix <<
"There is no estimate for the min eigenvalue." 1003 << std::endl << computeBeforeApplyReminder);
1004 TEUCHOS_TEST_FOR_EXCEPTION
1005 (STS::isnaninf (eigRatioForApply_), std::runtime_error,
1006 prefix <<
"There is no estimate for the ratio of the max " 1007 "eigenvalue to the min eigenvalue." 1008 << std::endl << computeBeforeApplyReminder);
1009 TEUCHOS_TEST_FOR_EXCEPTION
1010 (D_.is_null (), std::runtime_error, prefix <<
"The vector of inverse " 1011 "diagonal entries of the matrix has not yet been computed." 1012 << std::endl << computeBeforeApplyReminder);
1014 if (textbookAlgorithm_) {
1015 textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1016 lambdaMinForApply_, eigRatioForApply_, *D_);
1019 ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1020 lambdaMinForApply_, eigRatioForApply_, *D_);
1023 if (computeMaxResNorm_ && B.getNumVectors () > 0) {
1024 MV R (B.getMap (), B.getNumVectors ());
1025 computeResidual (R, B, *A_, X);
1026 Teuchos::Array<MT> norms (B.getNumVectors ());
1028 return *std::max_element (norms.begin (), norms.end ());
1031 return Teuchos::ScalarTraits<MT>::zero ();
1035 template<
class ScalarType,
class MV>
1040 using Teuchos::rcpFromRef;
1041 this->describe (* (Teuchos::getFancyOStream (rcpFromRef (out))),
1042 Teuchos::VERB_MEDIUM);
1045 template<
class ScalarType,
class MV>
1049 const ScalarType& alpha,
1054 solve (W, alpha, D_inv, B);
1055 Tpetra::deep_copy (X, W);
1058 template<
class ScalarType,
class MV>
1062 const Teuchos::ETransp mode)
1064 Tpetra::Details::residual(A,X,B,R);
1067 template <
class ScalarType,
class MV>
1069 Chebyshev<ScalarType, MV>::
1070 solve (MV& Z,
const V& D_inv,
const MV& R) {
1071 Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
1074 template<
class ScalarType,
class MV>
1076 Chebyshev<ScalarType, MV>::
1077 solve (MV& Z,
const ST alpha,
const V& D_inv,
const MV& R) {
1078 Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
1081 template<
class ScalarType,
class MV>
1082 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
1083 Chebyshev<ScalarType, MV>::
1084 makeInverseDiagonal (
const row_matrix_type& A,
const bool useDiagOffsets)
const 1087 using Teuchos::rcpFromRef;
1088 using Teuchos::rcp_dynamic_cast;
1091 if (!D_.is_null() &&
1092 D_->getMap()->isSameAs(*(A.getGraph ()->getRowMap ()))) {
1094 *out_ <<
"Reusing pre-existing vector for diagonal extraction" << std::endl;
1095 D_rowMap = Teuchos::rcp_const_cast<V>(D_);
1097 D_rowMap = Teuchos::rcp(
new V (A.getGraph ()->getRowMap (),
false));
1099 *out_ <<
"Allocated new vector for diagonal extraction" << std::endl;
1101 if (useDiagOffsets) {
1105 typedef Tpetra::CrsMatrix<
typename MV::scalar_type,
1106 typename MV::local_ordinal_type,
1107 typename MV::global_ordinal_type,
1108 typename MV::node_type> crs_matrix_type;
1109 RCP<const crs_matrix_type> A_crsMat =
1110 rcp_dynamic_cast<
const crs_matrix_type> (rcpFromRef (A));
1111 if (! A_crsMat.is_null ()) {
1112 TEUCHOS_TEST_FOR_EXCEPTION(
1113 ! savedDiagOffsets_, std::logic_error,
1114 "Ifpack2::Details::Chebyshev::makeInverseDiagonal: " 1115 "It is not allowed to call this method with useDiagOffsets=true, " 1116 "if you have not previously saved offsets of diagonal entries. " 1117 "This situation should never arise if this class is used properly. " 1118 "Please report this bug to the Ifpack2 developers.");
1119 A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_);
1125 A.getLocalDiagCopy (*D_rowMap);
1127 RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
1133 bool foundNonpositiveValue =
false;
1135 auto D_lcl = D_rangeMap->getLocalViewHost (Tpetra::Access::ReadOnly);
1136 auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1138 typedef typename MV::impl_scalar_type IST;
1139 typedef typename MV::local_ordinal_type LO;
1140 typedef Kokkos::Details::ArithTraits<IST> STS;
1141 typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
1143 const LO lclNumRows =
static_cast<LO
> (D_rangeMap->getLocalLength ());
1144 for (LO i = 0; i < lclNumRows; ++i) {
1145 if (STS::real (D_lcl_1d(i)) <= STM::zero ()) {
1146 foundNonpositiveValue =
true;
1152 using Teuchos::outArg;
1153 using Teuchos::REDUCE_MIN;
1154 using Teuchos::reduceAll;
1156 const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1157 int gblSuccess = lclSuccess;
1158 if (! D_rangeMap->getMap ().is_null () && D_rangeMap->getMap ()->getComm ().is_null ()) {
1159 const Teuchos::Comm<int>& comm = * (D_rangeMap->getMap ()->getComm ());
1160 reduceAll<int, int> (comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
1162 if (gblSuccess == 1) {
1163 *out_ <<
"makeInverseDiagonal: The matrix's diagonal entries all have " 1164 "positive real part (this is good for Chebyshev)." << std::endl;
1167 *out_ <<
"makeInverseDiagonal: The matrix's diagonal has at least one " 1168 "entry with nonpositive real part, on at least one process in the " 1169 "matrix's communicator. This is bad for Chebyshev." << std::endl;
1175 reciprocal_threshold (*D_rangeMap, minDiagVal_);
1176 return Teuchos::rcp_const_cast<
const V> (D_rangeMap);
1180 template<
class ScalarType,
class MV>
1181 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
1182 Chebyshev<ScalarType, MV>::
1183 makeRangeMapVectorConst (
const Teuchos::RCP<const V>& D)
const 1187 typedef Tpetra::Export<
typename MV::local_ordinal_type,
1188 typename MV::global_ordinal_type,
1189 typename MV::node_type> export_type;
1193 TEUCHOS_TEST_FOR_EXCEPTION(
1194 A_.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::" 1195 "makeRangeMapVector: The input matrix A is null. Please call setMatrix() " 1196 "with a nonnull input matrix before calling this method. This is probably " 1197 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1198 TEUCHOS_TEST_FOR_EXCEPTION(
1199 D.is_null (), std::logic_error,
"Ifpack2::Details::Chebyshev::" 1200 "makeRangeMapVector: The input Vector D is null. This is probably " 1201 "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1203 RCP<const map_type> sourceMap = D->getMap ();
1204 RCP<const map_type> rangeMap = A_->getRangeMap ();
1205 RCP<const map_type> rowMap = A_->getRowMap ();
1207 if (rangeMap->isSameAs (*sourceMap)) {
1213 RCP<const export_type> exporter;
1217 if (sourceMap->isSameAs (*rowMap)) {
1219 exporter = A_->getGraph ()->getExporter ();
1222 exporter = rcp (
new export_type (sourceMap, rangeMap));
1225 if (exporter.is_null ()) {
1229 RCP<V> D_out = rcp (
new V (*D, Teuchos::Copy));
1230 D_out->doExport (*D, *exporter, Tpetra::ADD);
1231 return Teuchos::rcp_const_cast<
const V> (D_out);
1237 template<
class ScalarType,
class MV>
1238 Teuchos::RCP<typename Chebyshev<ScalarType, MV>::V>
1239 Chebyshev<ScalarType, MV>::
1240 makeRangeMapVector (
const Teuchos::RCP<V>& D)
const 1242 using Teuchos::rcp_const_cast;
1243 return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1247 template<
class ScalarType,
class MV>
1249 Chebyshev<ScalarType, MV>::
1250 textbookApplyImpl (
const op_type& A,
1257 const V& D_inv)
const 1260 const ST myLambdaMin = lambdaMax / eigRatio;
1262 const ST zero = Teuchos::as<ST> (0);
1263 const ST one = Teuchos::as<ST> (1);
1264 const ST two = Teuchos::as<ST> (2);
1265 const ST d = (lambdaMax + myLambdaMin) / two;
1266 const ST c = (lambdaMax - myLambdaMin) / two;
1268 if (zeroStartingSolution_ && numIters > 0) {
1272 MV R (B.getMap (), B.getNumVectors (),
false);
1273 MV P (B.getMap (), B.getNumVectors (),
false);
1274 MV Z (B.getMap (), B.getNumVectors (),
false);
1276 for (
int i = 0; i < numIters; ++i) {
1277 computeResidual (R, B, A, X);
1278 solve (Z, D_inv, R);
1286 beta = alpha * (c/two) * (c/two);
1287 alpha = one / (d - beta);
1288 P.update (one, Z, beta);
1290 X.update (alpha, P, one);
1297 template<
class ScalarType,
class MV>
1298 typename Chebyshev<ScalarType, MV>::MT
1299 Chebyshev<ScalarType, MV>::maxNormInf (
const MV& X) {
1300 Teuchos::Array<MT> norms (X.getNumVectors ());
1301 X.normInf (norms());
1302 return *std::max_element (norms.begin (), norms.end ());
1305 template<
class ScalarType,
class MV>
1307 Chebyshev<ScalarType, MV>::
1308 ifpackApplyImpl (
const op_type& A,
1318 #ifdef HAVE_IFPACK2_DEBUG 1319 const bool debug = debug_;
1321 const bool debug =
false;
1325 *out_ <<
" \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1326 *out_ <<
" \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1329 if (numIters <= 0) {
1332 const ST zero =
static_cast<ST
> (0.0);
1333 const ST one =
static_cast<ST
> (1.0);
1334 const ST two =
static_cast<ST
> (2.0);
1337 if (lambdaMin == one && lambdaMax == lambdaMin) {
1338 solve (X, D_inv, B);
1343 const ST alpha = lambdaMax / eigRatio;
1344 const ST beta = boostFactor_ * lambdaMax;
1345 const ST delta = two / (beta - alpha);
1346 const ST theta = (beta + alpha) / two;
1347 const ST s1 = theta * delta;
1350 *out_ <<
" alpha = " << alpha << endl
1351 <<
" beta = " << beta << endl
1352 <<
" delta = " << delta << endl
1353 <<
" theta = " << theta << endl
1354 <<
" s1 = " << s1 << endl;
1358 Teuchos::RCP<MV> W_ptr = makeTempMultiVector (B);
1362 *out_ <<
" Iteration " << 1 <<
":" << endl
1363 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1367 if (! zeroStartingSolution_) {
1370 if (ck_.is_null ()) {
1371 Teuchos::RCP<const op_type> A_op = A_;
1372 ck_ = Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op));
1376 ck_->compute (W, one/theta, const_cast<V&> (D_inv),
1377 const_cast<MV&> (B), X, zero);
1381 firstIterationWithZeroStartingSolution (W, one/theta, D_inv, B, X);
1385 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1386 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1389 if (numIters > 1 && ck_.is_null ()) {
1390 Teuchos::RCP<const op_type> A_op = A_;
1391 ck_ = Teuchos::rcp (
new ChebyshevKernel<op_type> (A_op));
1396 ST rhokp1, dtemp1, dtemp2;
1397 for (
int deg = 1; deg < numIters; ++deg) {
1399 *out_ <<
" Iteration " << deg+1 <<
":" << endl
1400 <<
" - \\|D\\|_{\\infty} = " << D_->normInf () << endl
1401 <<
" - \\|B\\|_{\\infty} = " << maxNormInf (B) << endl
1402 <<
" - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm ()
1403 << endl <<
" - rhok = " << rhok << endl;
1406 rhokp1 = one / (two * s1 - rhok);
1407 dtemp1 = rhokp1 * rhok;
1408 dtemp2 = two * rhokp1 * delta;
1412 *out_ <<
" - dtemp1 = " << dtemp1 << endl
1413 <<
" - dtemp2 = " << dtemp2 << endl;
1418 ck_->compute (W, dtemp2, const_cast<V&> (D_inv),
1419 const_cast<MV&> (B), (X), dtemp1);
1422 *out_ <<
" - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1423 <<
" - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1428 template<
class ScalarType,
class MV>
1429 typename Chebyshev<ScalarType, MV>::ST
1430 Chebyshev<ScalarType, MV>::
1431 powerMethodWithInitGuess (
const op_type& A,
1438 *out_ <<
" powerMethodWithInitGuess:" << endl;
1441 const ST zero =
static_cast<ST
> (0.0);
1442 const ST one =
static_cast<ST
> (1.0);
1443 ST lambdaMax = zero;
1444 ST lambdaMaxOld = one;
1448 if (eigVector2_.is_null()) {
1449 y = rcp(
new V(A.getRangeMap ()));
1450 if (eigKeepVectors_)
1455 TEUCHOS_TEST_FOR_EXCEPTION
1456 (norm == zero, std::runtime_error,
1457 "Ifpack2::Chebyshev::powerMethodWithInitGuess: The initial guess " 1458 "has zero norm. This could be either because Tpetra::Vector::" 1459 "randomize filled the vector with zeros (if that was used to " 1460 "compute the initial guess), or because the norm2 method has a " 1461 "bug. The first is not impossible, but unlikely.");
1464 *out_ <<
" Original norm1(x): " << x.norm1 ()
1465 <<
", norm2(x): " << norm << endl;
1468 x.scale (one / norm);
1471 *out_ <<
" norm1(x.scale(one/norm)): " << x.norm1 () << endl;
1474 for (
int iter = 0; iter < numIters-1; ++iter) {
1476 *out_ <<
" Iteration " << iter << endl;
1479 solve (x, D_inv, *y);
1481 if (((iter+1) % eigNormalizationFreq_ == 0) && (iter < numIters-2)) {
1485 *out_ <<
" norm is zero; returning zero" << endl;
1486 *out_ <<
" Power method terminated after "<< iter <<
" iterations." << endl;
1490 lambdaMaxOld = lambdaMax;
1491 lambdaMax = pow(norm, Teuchos::ScalarTraits<MT>::one() / eigNormalizationFreq_);
1492 if (Teuchos::ScalarTraits<ST>::magnitude(lambdaMax-lambdaMaxOld) < eigRelTolerance_ * Teuchos::ScalarTraits<ST>::magnitude(lambdaMax)) {
1494 *out_ <<
" lambdaMax: " << lambdaMax << endl;
1495 *out_ <<
" Power method terminated after "<< iter <<
" iterations." << endl;
1498 }
else if (debug_) {
1499 *out_ <<
" lambdaMaxOld: " << lambdaMaxOld << endl;
1500 *out_ <<
" lambdaMax: " << lambdaMax << endl;
1501 *out_ <<
" |lambdaMax-lambdaMaxOld|/|lambdaMax|: " << Teuchos::ScalarTraits<ST>::magnitude(lambdaMax-lambdaMaxOld)/Teuchos::ScalarTraits<ST>::magnitude(lambdaMax) << endl;
1504 x.scale (one / norm);
1508 *out_ <<
" lambdaMax: " << lambdaMax << endl;
1514 *out_ <<
" norm is zero; returning zero" << endl;
1515 *out_ <<
" Power method terminated after "<< numIters <<
" iterations." << endl;
1519 x.scale (one / norm);
1521 solve (*y, D_inv, *y);
1522 lambdaMax = y->dot (x);
1524 *out_ <<
" lambdaMax: " << lambdaMax << endl;
1525 *out_ <<
" Power method terminated after "<< numIters <<
" iterations." << endl;
1531 template<
class ScalarType,
class MV>
1533 Chebyshev<ScalarType, MV>::
1534 computeInitialGuessForPowerMethod (V& x,
const bool nonnegativeRealParts)
const 1536 typedef typename MV::device_type::execution_space dev_execution_space;
1537 typedef typename MV::local_ordinal_type LO;
1541 if (nonnegativeRealParts) {
1542 auto x_lcl = x.getLocalViewDevice (Tpetra::Access::ReadWrite);
1543 auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
1545 const LO lclNumRows =
static_cast<LO
> (x.getLocalLength ());
1546 Kokkos::RangePolicy<dev_execution_space, LO> range (0, lclNumRows);
1547 PositivizeVector<decltype (x_lcl_1d), LO> functor (x_lcl_1d);
1548 Kokkos::parallel_for (range, functor);
1552 template<
class ScalarType,
class MV>
1553 typename Chebyshev<ScalarType, MV>::ST
1554 Chebyshev<ScalarType, MV>::
1555 powerMethod (
const op_type& A,
const V& D_inv,
const int numIters)
1559 *out_ <<
"powerMethod:" << endl;
1562 const ST zero =
static_cast<ST
> (0.0);
1564 if (eigVector_.is_null()) {
1565 x = rcp(
new V(A.getDomainMap ()));
1566 if (eigKeepVectors_)
1571 computeInitialGuessForPowerMethod (*x,
false);
1575 ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, *x);
1584 if (STS::real (lambdaMax) < STS::real (zero)) {
1586 *out_ <<
"real(lambdaMax) = " << STS::real (lambdaMax) <<
" < 0; " 1587 "try again with a different random initial guess" << endl;
1596 computeInitialGuessForPowerMethod (*x,
true);
1597 lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, *x);
1603 template<
class ScalarType,
class MV>
1604 typename Chebyshev<ScalarType, MV>::ST
1605 Chebyshev<ScalarType, MV>::
1606 cgMethodWithInitGuess (
const op_type& A,
1612 using STS = Teuchos::ScalarTraits<ST>;
1613 using MagnitudeType =
typename STS::magnitudeType;
1615 *out_ <<
" cgMethodWithInitGuess:" << endl;
1618 const ST one = STS::one();
1619 ST beta, betaOld = one, pAp, pApOld = one, alpha, rHz, rHzOld, rHzOld2 = one, lambdaMax;
1621 Teuchos::ArrayRCP<MagnitudeType> diag, offdiag;
1622 Teuchos::RCP<V> p, z, Ap;
1623 diag.resize(numIters);
1624 offdiag.resize(numIters-1);
1626 p = rcp(
new V(A.getRangeMap ()));
1627 z = rcp(
new V(A.getRangeMap ()));
1628 Ap = rcp(
new V(A.getRangeMap ()));
1631 solve (*p, D_inv, r);
1634 for (
int iter = 0; iter < numIters; ++iter) {
1636 *out_ <<
" Iteration " << iter << endl;
1641 r.update (-alpha, *Ap, one);
1643 solve (*z, D_inv, r);
1645 beta = rHz / rHzOld;
1646 p->update (one, *z, beta);
1648 diag[iter] = STS::real((betaOld*betaOld * pApOld + pAp) / rHzOld);
1649 offdiag[iter-1] = -STS::real((betaOld * pApOld / (sqrt(rHzOld * rHzOld2))));
1651 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1652 *out_ <<
" offdiag["<< iter-1 <<
"] = " << offdiag[iter-1] << endl;
1655 diag[iter] = STS::real(pAp/rHzOld);
1657 *out_ <<
" diag[" << iter <<
"] = " << diag[iter] << endl;
1665 lambdaMax = LapackHelper<ST>::tri_diag_spectral_radius(diag, offdiag);
1671 template<
class ScalarType,
class MV>
1672 typename Chebyshev<ScalarType, MV>::ST
1673 Chebyshev<ScalarType, MV>::
1674 cgMethod (
const op_type& A,
const V& D_inv,
const int numIters)
1678 *out_ <<
"cgMethod:" << endl;
1682 if (eigVector_.is_null()) {
1683 r = rcp(
new V(A.getDomainMap ()));
1684 if (eigKeepVectors_)
1689 computeInitialGuessForPowerMethod (*r,
false);
1693 ST lambdaMax = cgMethodWithInitGuess (A, D_inv, numIters, *r);
1698 template<
class ScalarType,
class MV>
1699 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::row_matrix_type>
1704 template<
class ScalarType,
class MV>
1712 template<
class ScalarType,
class MV>
1721 const size_t B_numVecs = B.getNumVectors ();
1722 if (W_.is_null () || W_->getNumVectors () != B_numVecs) {
1723 W_ = Teuchos::rcp (
new MV (B.getMap (), B_numVecs,
false));
1728 template<
class ScalarType,
class MV>
1730 Chebyshev<ScalarType, MV>::
1731 description ()
const {
1732 std::ostringstream oss;
1735 oss <<
"\"Ifpack2::Details::Chebyshev\":" 1737 <<
"degree: " << numIters_
1738 <<
", lambdaMax: " << lambdaMaxForApply_
1739 <<
", alpha: " << eigRatioForApply_
1740 <<
", lambdaMin: " << lambdaMinForApply_
1741 <<
", boost factor: " << boostFactor_;
1742 if (!userInvDiag_.is_null())
1743 oss <<
", diagonal: user-supplied";
1748 template<
class ScalarType,
class MV>
1752 const Teuchos::EVerbosityLevel verbLevel)
const 1754 using Teuchos::TypeNameTraits;
1757 const Teuchos::EVerbosityLevel vl =
1758 (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1759 if (vl == Teuchos::VERB_NONE) {
1769 Teuchos::OSTab tab0 (out);
1775 if (A_.is_null () || A_->getComm ().is_null ()) {
1779 myRank = A_->getComm ()->getRank ();
1784 out <<
"\"Ifpack2::Details::Chebyshev\":" << endl;
1786 Teuchos::OSTab tab1 (out);
1788 if (vl == Teuchos::VERB_LOW) {
1790 out <<
"degree: " << numIters_ << endl
1791 <<
"lambdaMax: " << lambdaMaxForApply_ << endl
1792 <<
"alpha: " << eigRatioForApply_ << endl
1793 <<
"lambdaMin: " << lambdaMinForApply_ << endl
1794 <<
"boost factor: " << boostFactor_ << endl;
1802 out <<
"Template parameters:" << endl;
1804 Teuchos::OSTab tab2 (out);
1805 out <<
"ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1806 <<
"MV: " << TypeNameTraits<MV>::name () << endl;
1812 out <<
"Computed parameters:" << endl;
1818 Teuchos::OSTab tab2 (out);
1824 if (D_.is_null ()) {
1826 out <<
"unset" << endl;
1829 else if (vl <= Teuchos::VERB_HIGH) {
1831 out <<
"set" << endl;
1840 D_->describe (out, vl);
1845 out <<
"W_: " << (W_.is_null () ?
"unset" :
"set") << endl
1846 <<
"computedLambdaMax_: " << computedLambdaMax_ << endl
1847 <<
"computedLambdaMin_: " << computedLambdaMin_ << endl
1848 <<
"lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1849 <<
"lambdaMinForApply_: " << lambdaMinForApply_ << endl
1850 <<
"eigRatioForApply_: " << eigRatioForApply_ << endl;
1855 out <<
"User parameters:" << endl;
1860 Teuchos::OSTab tab2 (out);
1861 out <<
"userInvDiag_: ";
1862 if (userInvDiag_.is_null ()) {
1863 out <<
"unset" << endl;
1865 else if (vl <= Teuchos::VERB_HIGH) {
1866 out <<
"set" << endl;
1872 userInvDiag_->describe (out, vl);
1875 out <<
"userLambdaMax_: " << userLambdaMax_ << endl
1876 <<
"userLambdaMin_: " << userLambdaMin_ << endl
1877 <<
"userEigRatio_: " << userEigRatio_ << endl
1878 <<
"numIters_: " << numIters_ << endl
1879 <<
"eigMaxIters_: " << eigMaxIters_ << endl
1880 <<
"eigRelTolerance_: " << eigRelTolerance_ << endl
1881 <<
"eigNormalizationFreq_: " << eigNormalizationFreq_ << endl
1882 <<
"zeroStartingSolution_: " << zeroStartingSolution_ << endl
1883 <<
"assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1884 <<
"textbookAlgorithm_: " << textbookAlgorithm_ << endl
1885 <<
"computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1893 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \ 1894 template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >; 1896 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:322
Ifpack2 implementation details.
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:199
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:116
scalar_type ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:112
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:106
Declaration of Chebyshev implementation.
Definition: Ifpack2_Container_decl.hpp:576
Teuchos::ScalarTraits< scalar_type > STS
Traits class for ST.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:114
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
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