Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_Details_Chebyshev_def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
4 //
5 // Copyright 2009 NTESS and the Ifpack2 contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
11 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
12 
19 
20 #include "Ifpack2_PowerMethod.hpp"
23 // #include "Ifpack2_Details_ScaledDampedResidual.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"
29 #include "Teuchos_LAPACK.hpp"
30 #include "Ifpack2_Details_LapackSupportsScalar.hpp"
31 #include <cmath>
32 #include <iostream>
33 
34 namespace Ifpack2 {
35 namespace Details {
36 
37 namespace { // (anonymous)
38 
39 // We use this text a lot in error messages.
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()).";
49 
50 } // namespace
51 
52 // ReciprocalThreshold stuff below needs to be in a namspace visible outside
53 // of this file
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;
61 
62  XV X_;
63  const value_type minVal_;
64  const mag_type minValMag_;
65 
66  V_ReciprocalThresholdSelfFunctor(const XV& X,
67  const value_type& min_val)
68  : X_(X)
69  , minVal_(min_val)
70  , minValMag_(KAT::abs(min_val)) {}
71 
72  KOKKOS_INLINE_FUNCTION
73  void operator()(const size_type& i) const {
74  const mag_type X_i_abs = KAT::abs(X_[i]);
75 
76  if (X_i_abs < minValMag_) {
77  X_[i] = minVal_;
78  } else {
79  X_[i] = KAT::one() / X_[i];
80  }
81  }
82 };
83 
84 template <class XV, class SizeType = typename XV::size_type>
85 struct LocalReciprocalThreshold {
86  static void
87  compute(const XV& X,
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);
93  }
94 };
95 
96 template <class TpetraVectorType,
97  const bool classic = TpetraVectorType::node_type::classic>
98 struct GlobalReciprocalThreshold {};
99 
100 template <class TpetraVectorType>
101 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
102  static void
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;
108 
109  const scalar_type ONE = STS::one();
110  const mag_type min_val_abs = STS::abs(min_val);
111 
112  Teuchos::ArrayRCP<scalar_type> V_0 = V.getDataNonConst(0);
113  const size_t lclNumRows = V.getLocalLength();
114 
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) {
118  V_0[i] = min_val;
119  } else {
120  V_0[i] = ONE / V_0i;
121  }
122  }
123  }
124 };
125 
126 template <class TpetraVectorType>
127 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
128  static void
129  compute(TpetraVectorType& X,
130  const typename TpetraVectorType::scalar_type& minVal) {
131  typedef typename TpetraVectorType::impl_scalar_type value_type;
132 
133  const value_type minValS = static_cast<value_type>(minVal);
134  auto X_0 = Kokkos::subview(X.getLocalViewDevice(Tpetra::Access::ReadWrite),
135  Kokkos::ALL(), 0);
136  LocalReciprocalThreshold<decltype(X_0)>::compute(X_0, minValS);
137  }
138 };
139 
140 // Utility function for inverting diagonal with threshold.
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);
144 }
145 
146 template <class ScalarType, const bool lapackSupportsScalarType = LapackSupportsScalar<ScalarType>::value>
147 struct LapackHelper {
148  static ScalarType
149  tri_diag_spectral_radius(Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
151  throw std::runtime_error("LAPACK does not support the scalar type.");
152  }
153 };
154 
155 template <class V>
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>;
159 
160  // Initial randomization of the vector
161  x.randomize();
162 
163  // Zero the stuff that where the diagonal is equal to one. These are assumed to
164  // correspond to OAZ rows in the matrix.
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);
168 
171 
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)
175  x_view(i, 0) = ZERO;
176  });
177 }
178 
179 template <class ScalarType>
180 struct LapackHelper<ScalarType, true> {
181  static ScalarType
182  tri_diag_spectral_radius(Teuchos::ArrayRCP<typename Teuchos::ScalarTraits<ScalarType>::magnitudeType> diag,
185  using MagnitudeType = typename STS::magnitudeType;
186  int info = 0;
187  const int N = diag.size();
188  ScalarType scalar_dummy;
189  std::vector<MagnitudeType> mag_dummy(4 * N);
190  char char_N = 'N';
191 
192  // lambdaMin = one;
193  ScalarType lambdaMax = STS::one();
194  if (N > 2) {
196  lapack.PTEQR(char_N, N, diag.getRawPtr(), offdiag.getRawPtr(),
197  &scalar_dummy, 1, &mag_dummy[0], &info);
198  TEUCHOS_TEST_FOR_EXCEPTION(info < 0, std::logic_error,
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.");
203  // lambdaMin = Teuchos::as<ScalarType> (diag[N-1]);
204  lambdaMax = Teuchos::as<ScalarType>(diag[0]);
205  }
206  return lambdaMax;
207  }
208 };
209 
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. "
216  "A has "
217  << A_->getGlobalNumRows() << " rows and "
218  << A_->getGlobalNumCols() << " columns.");
219 
220  // In debug mode, test that the domain and range Maps of the matrix
221  // are the same.
222  if (debug_ && !A_.is_null()) {
223  Teuchos::RCP<const map_type> domainMap = A_->getDomainMap();
224  Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap();
225 
226  // isSameAs is a collective, but if the two pointers are the same,
227  // isSameAs will assume that they are the same on all processes, and
228  // return true without an all-reduce.
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())."
233  << std::endl
234  << "We only check "
235  "for this in debug mode.");
236  }
237 }
238 
239 template <class ScalarType, class MV>
240 void Chebyshev<ScalarType, MV>::
241  checkConstructorInput() const {
242  // mfh 12 Aug 2016: The if statement avoids an "unreachable
243  // statement" warning for the checkInputMatrix() call, when
244  // STS::isComplex is false.
245  if (STS::isComplex) {
246  TEUCHOS_TEST_FOR_EXCEPTION(true, std::logic_error,
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"
250  " = "
251  << Teuchos::TypeNameTraits<ScalarType>::name() << ", which is a "
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).");
256  } else {
257  checkInputMatrix();
258  }
259 }
260 
261 template <class ScalarType, class MV>
264  : A_(A)
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())
275  , numIters_(1)
276  , eigMaxIters_(10)
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)
288  , debug_(false) {
289  checkConstructorInput();
290 }
291 
292 template <class ScalarType, class MV>
295  Teuchos::ParameterList& params)
296  : A_(A)
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())
308  , numIters_(1)
309  , eigMaxIters_(10)
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)
321  , debug_(false) {
322  checkConstructorInput();
323  setParameters(params);
324 }
325 
326 template <class ScalarType, class MV>
329  using Teuchos::RCP;
330  using Teuchos::rcp;
331  using Teuchos::rcp_const_cast;
332 
333  // Note to developers: The logic for this method is complicated,
334  // because we want to accept Ifpack and ML parameters whenever
335  // possible, but we don't want to add their default values to the
336  // user's ParameterList. That's why we do all the isParameter()
337  // checks, instead of using the two-argument version of get()
338  // everywhere. The min and max eigenvalue parameters are also a
339  // special case, because we decide whether or not to do eigenvalue
340  // analysis based on whether the user supplied the max eigenvalue.
341 
342  // Default values of all the parameters.
343  const ST defaultLambdaMax = STS::nan();
344  const ST defaultLambdaMin = STS::nan();
345  // 30 is Ifpack::Chebyshev's default. ML has a different default
346  // eigRatio for smoothers and the coarse-grid solve (if using
347  // Chebyshev for that). The former uses 20; the latter uses 30.
348  // We're testing smoothers here, so use 20. (However, if you give
349  // ML an Epetra matrix, it will use Ifpack for Chebyshev, in which
350  // case it would defer to Ifpack's default settings.)
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;
356  const MT defaultEigRelTolerance = Teuchos::ScalarTraits<MT>::zero();
357  const bool defaultEigKeepVectors = false;
358  const int defaultEigNormalizationFreq = 1;
359  const bool defaultZeroStartingSolution = true; // Ifpack::Chebyshev default
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;
367 
368  // We'll set the instance data transactionally, after all reads
369  // from the ParameterList. That way, if any of the ParameterList
370  // reads fail (e.g., due to the wrong parameter type), we will not
371  // have left the instance data in a half-changed state.
372  RCP<const V> userInvDiagCopy; // if nonnull: deep copy of user's Vector
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;
391 
392  // Fetch the parameters from the ParameterList. Defer all
393  // externally visible side effects until we have finished all
394  // ParameterList interaction. This makes the method satisfy the
395  // strong exception guarantee.
396 
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;
402  }
403 
404  // Get the user-supplied inverse diagonal.
405  //
406  // Check for a raw pointer (const V* or V*), for Ifpack
407  // compatibility, as well as for RCP<const V>, RCP<V>, const V, or
408  // V. We'll make a deep copy of the vector at the end of this
409  // method anyway, so its const-ness doesn't matter. We handle the
410  // latter two cases ("const V" or "V") specially (copy them into
411  // userInvDiagCopy first, which is otherwise null at the end of the
412  // long if-then chain) to avoid an extra copy.
413 
414  const char opInvDiagLabel[] = "chebyshev: operator inv diagonal";
415  if (plist.isParameter(opInvDiagLabel)) {
416  // Pointer to the user's Vector, if provided.
417  RCP<const V> userInvDiag;
418 
419  if (plist.isType<const V*>(opInvDiagLabel)) {
420  const V* rawUserInvDiag =
421  plist.get<const V*>(opInvDiagLabel);
422  // Nonowning reference (we'll make a deep copy below)
423  userInvDiag = rcp(rawUserInvDiag, false);
424  } else if (plist.isType<const V*>(opInvDiagLabel)) {
425  V* rawUserInvDiag = plist.get<V*>(opInvDiagLabel);
426  // Nonowning reference (we'll make a deep copy below)
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);
436  userInvDiagCopy = rcp(new V(userInvDiagRef, Teuchos::Copy));
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);
441  userInvDiagCopy = rcp(new V(userInvDiagRef, Teuchos::Copy));
442  userInvDiag = userInvDiagCopy;
443  }
444 
445  // NOTE: If the user's parameter has some strange type that we
446  // didn't test above, userInvDiag might still be null. You may
447  // want to add an error test for this condition. Currently, we
448  // just say in this case that the user didn't give us a Vector.
449 
450  // If we have userInvDiag but don't have a deep copy yet, make a
451  // deep copy now.
452  if (!userInvDiag.is_null() && userInvDiagCopy.is_null()) {
453  userInvDiagCopy = rcp(new V(*userInvDiag, Teuchos::Copy));
454  }
455 
456  // NOTE: userInvDiag, if provided, is a row Map version of the
457  // Vector. We don't necessarily have a range Map yet. compute()
458  // would be the proper place to compute the range Map version of
459  // userInvDiag.
460  }
461 
462  // Load the kernel fuse override from the parameter list
463  if (plist.isParameter("chebyshev: use native spmv"))
464  ckUseNativeSpMV = plist.get("chebyshev: use native spmv", ckUseNativeSpMV);
465 
466  // Load the pre-allocate overrride from the parameter list
467  if (plist.isParameter("chebyshev: pre-allocate temp vector"))
468  preAllocateTempVector = plist.get("chebyshev: pre-allocate temp vector", preAllocateTempVector);
469 
470  // Don't fill in defaults for the max or min eigenvalue, because
471  // this class uses the existence of those parameters to determine
472  // whether it should do eigenanalysis.
473  if (plist.isParameter("chebyshev: max eigenvalue")) {
474  if (plist.isType<double>("chebyshev: max eigenvalue"))
475  lambdaMax = plist.get<double>("chebyshev: max eigenvalue");
476  else
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.");
483  }
484  if (plist.isParameter("chebyshev: min eigenvalue")) {
485  if (plist.isType<double>("chebyshev: min eigenvalue"))
486  lambdaMin = plist.get<double>("chebyshev: min eigenvalue");
487  else
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.");
494  }
495 
496  // Only fill in Ifpack2's name for the default parameter, not ML's.
497  if (plist.isParameter("smoother: Chebyshev alpha")) { // ML compatibility
498  if (plist.isType<double>("smoother: Chebyshev alpha"))
499  eigRatio = plist.get<double>("smoother: Chebyshev alpha");
500  else
501  eigRatio = plist.get<ST>("smoother: Chebyshev alpha");
502  }
503  // Ifpack2's name overrides ML's name.
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 "
510  "a finite value.");
511  // mfh 11 Feb 2013: This class is currently only correct for real
512  // Scalar types, but we still want it to build for complex Scalar
513  // type so that users of Ifpack2::Factory can build their
514  // executables for real or complex Scalar type. Thus, we take the
515  // real parts here, which are always less-than comparable.
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 "
522  << eigRatio << ".");
523 
524  // See Github Issue #234. This parameter may be either MT
525  // (preferred) or double. We check both.
526  {
527  const char paramName[] = "chebyshev: boost factor";
528 
529  if (plist.isParameter(paramName)) {
530  if (plist.isType<MT>(paramName)) { // MT preferred
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);
536  } else {
537  TEUCHOS_TEST_FOR_EXCEPTION(true, std::invalid_argument,
538  "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\""
539  "parameter must have type magnitude_type (MT) or double.");
540  }
541  } else { // parameter not in the list
542  // mfh 12 Aug 2016: To preserve current behavior (that fills in
543  // any parameters not in the input ParameterList with their
544  // default values), we call set() here. I don't actually like
545  // this behavior; I prefer the Belos model, where the input
546  // ParameterList is a delta from current behavior. However, I
547  // don't want to break things.
548  plist.set(paramName, defaultBoostFactor);
549  }
550  TEUCHOS_TEST_FOR_EXCEPTION(boostFactor < Teuchos::ScalarTraits<MT>::one(), std::invalid_argument,
551  "Ifpack2::Chebyshev::setParameters: \"" << paramName << "\" parameter "
552  "must be >= 1, but you supplied the value "
553  << boostFactor << ".");
554  }
555 
556  // Same name in Ifpack2 and Ifpack.
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.");
563 
564  // Only fill in Ifpack2's name, not ML's or Ifpack's.
565  if (plist.isParameter("smoother: sweeps")) { // ML compatibility
566  numIters = plist.get<int>("smoother: sweeps");
567  } // Ifpack's name overrides ML's name.
568  if (plist.isParameter("relaxation: sweeps")) { // Ifpack compatibility
569  numIters = plist.get<int>("relaxation: sweeps");
570  } // Ifpack2's name overrides Ifpack's name.
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 "
577  << numIters << ".");
578 
579  // The last parameter name overrides the first.
580  if (plist.isParameter("eigen-analysis: iterations")) { // ML compatibility
581  eigMaxIters = plist.get<int>("eigen-analysis: iterations");
582  } // Ifpack2's name overrides ML's name.
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 << ".");
590 
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"))
596  eigRelTolerance = Teuchos::ScalarTraits<ST>::magnitude(plist.get<ST>("chebyshev: eigenvalue relative tolerance"));
597 
598  eigKeepVectors = plist.get("chebyshev: eigenvalue keep vectors", eigKeepVectors);
599 
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 << ".")
607 
608  zeroStartingSolution = plist.get("chebyshev: zero starting solution",
609  zeroStartingSolution);
610  assumeMatrixUnchanged = plist.get("chebyshev: assume matrix does not change",
611  assumeMatrixUnchanged);
612 
613  // We don't want to fill these parameters in, because they shouldn't
614  // be visible to Ifpack2::Chebyshev users.
615  if (plist.isParameter("chebyshev: algorithm")) {
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\".");
624  }
625 
626 #ifdef IFPACK2_ENABLE_DEPRECATED_CODE
627  // to preserve behavior with previous input decks, only read "chebyshev:textbook algorithm" setting
628  // if a user has not specified "chebyshev: algorithm"
629  if (!plist.isParameter("chebyshev: algorithm")) {
630  if (plist.isParameter("chebyshev: textbook algorithm")) {
631  const bool textbookAlgorithm = plist.get<bool>("chebyshev: textbook algorithm");
632  if (textbookAlgorithm) {
633  chebyshevAlgorithm = "textbook";
634  } else {
635  chebyshevAlgorithm = "first";
636  }
637  }
638  }
639 #endif
640 
641  if (plist.isParameter("chebyshev: compute max residual norm")) {
642  computeMaxResNorm = plist.get<bool>("chebyshev: compute max residual norm");
643  }
644  if (plist.isParameter("chebyshev: compute spectral radius")) {
645  computeSpectralRadius = plist.get<bool>("chebyshev: compute spectral radius");
646  }
647 
648  // Test for Ifpack parameters that we won't ever implement here.
649  // Be careful to use the one-argument version of get(), since the
650  // two-argment version adds the parameter if it's not there.
651  TEUCHOS_TEST_FOR_EXCEPTION(plist.isType<bool>("chebyshev: use block mode") &&
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.");
657  TEUCHOS_TEST_FOR_EXCEPTION(plist.isType<bool>("chebyshev: solve normal equations") &&
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.");
664 
665  // Test for Ifpack parameters that we haven't implemented yet.
666  //
667  // For now, we only check that this ML parameter, if provided, has
668  // the one value that we support. We consider other values "invalid
669  // arguments" rather than "logic errors," because Ifpack does not
670  // implement eigenanalyses other than the power method.
671  std::string eigenAnalysisType("power-method");
672  if (plist.isParameter("eigen-analysis: type")) {
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\".");
680  }
681 
682  // We've validated all the parameters, so it's safe now to "commit" them.
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;
702  debug_ = debug;
703 
704  if (debug_) {
705  // Only print if myRank == 0.
706  int myRank = -1;
707  if (A_.is_null() || A_->getComm().is_null()) {
708  // We don't have a communicator (yet), so we assume that
709  // everybody can print. Revise this expectation in setMatrix().
710  myRank = 0;
711  } else {
712  myRank = A_->getComm()->getRank();
713  }
714 
715  if (myRank == 0) {
716  out_ = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
717  } else {
718  using Teuchos::oblackholestream; // prints nothing
719  RCP<oblackholestream> blackHole(new oblackholestream());
720  out_ = Teuchos::getFancyOStream(blackHole);
721  }
722  } else { // NOT debug
723  // free the "old" output stream, if there was one
724  out_ = Teuchos::null;
725  }
726 }
727 
728 template <class ScalarType, class MV>
730  ck_ = Teuchos::null;
731  D_ = Teuchos::null;
732  diagOffsets_ = offsets_type();
733  savedDiagOffsets_ = false;
734  W_ = Teuchos::null;
735  computedLambdaMax_ = STS::nan();
736  computedLambdaMin_ = STS::nan();
737  eigVector_ = Teuchos::null;
738  eigVector2_ = Teuchos::null;
739 }
740 
741 template <class ScalarType, class MV>
744  if (A.getRawPtr() != A_.getRawPtr()) {
745  if (!assumeMatrixUnchanged_) {
746  reset();
747  }
748  A_ = A;
749  ck_ = Teuchos::null; // constructed on demand
750 
751  // The communicator may have changed, or we may not have had a
752  // communicator before. Thus, we may have to reset the debug
753  // output stream.
754  if (debug_) {
755  // Only print if myRank == 0.
756  int myRank = -1;
757  if (A.is_null() || A->getComm().is_null()) {
758  // We don't have a communicator (yet), so we assume that
759  // everybody can print. Revise this expectation in setMatrix().
760  myRank = 0;
761  } else {
762  myRank = A->getComm()->getRank();
763  }
764 
765  if (myRank == 0) {
766  out_ = Teuchos::getFancyOStream(Teuchos::rcpFromRef(std::cerr));
767  } else {
769  out_ = Teuchos::getFancyOStream(blackHole); // print nothing on other processes
770  }
771  } else { // NOT debug
772  // free the "old" output stream, if there was one
773  out_ = Teuchos::null;
774  }
775  }
776 }
777 
778 template <class ScalarType, class MV>
780  using std::endl;
781  // Some of the optimizations below only work if A_ is a
782  // Tpetra::CrsMatrix. We'll make our best guess about its type
783  // here, since we have no way to get back the original fifth
784  // template parameter.
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>
789  crs_matrix_type;
790 
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.");
796 
797  // If A_ is a CrsMatrix and its graph is constant, we presume that
798  // the user plans to reuse the structure of A_, but possibly change
799  // A_'s values before each compute() call. This is the intended use
800  // case for caching the offsets of the diagonal entries of A_, to
801  // speed up extraction of diagonal entries on subsequent compute()
802  // calls.
803 
804  // FIXME (mfh 22 Jan 2013, 10 Feb 2013) In all cases when we use
805  // isnaninf() in this method, we really only want to check if the
806  // number is NaN. Inf means something different. However,
807  // Teuchos::ScalarTraits doesn't distinguish the two cases.
808 
809  // makeInverseDiagonal() returns a range Map Vector.
810  if (userInvDiag_.is_null()) {
812  Teuchos::rcp_dynamic_cast<const crs_matrix_type>(A_);
813  if (D_.is_null()) { // We haven't computed D_ before
814  if (!A_crsMat.is_null() && A_crsMat->isFillComplete()) {
815  // It's a CrsMatrix with a const graph; cache diagonal offsets.
816  const size_t lclNumRows = A_crsMat->getLocalNumRows();
817  if (diagOffsets_.extent(0) < lclNumRows) {
818  diagOffsets_ = offsets_type(); // clear first to save memory
819  diagOffsets_ = offsets_type("offsets", lclNumRows);
820  }
821  A_crsMat->getCrsGraph()->getLocalDiagOffsets(diagOffsets_);
822  savedDiagOffsets_ = true;
823  D_ = makeInverseDiagonal(*A_, true);
824  } else { // either A_ is not a CrsMatrix, or its graph is nonconst
825  D_ = makeInverseDiagonal(*A_);
826  }
827  } else if (!assumeMatrixUnchanged_) { // D_ exists but A_ may have changed
828  if (!A_crsMat.is_null() && A_crsMat->isFillComplete()) {
829  // It's a CrsMatrix with a const graph; cache diagonal offsets
830  // if we haven't already.
831  if (!savedDiagOffsets_) {
832  const size_t lclNumRows = A_crsMat->getLocalNumRows();
833  if (diagOffsets_.extent(0) < lclNumRows) {
834  diagOffsets_ = offsets_type(); // clear first to save memory
835  diagOffsets_ = offsets_type("offsets", lclNumRows);
836  }
837  A_crsMat->getCrsGraph()->getLocalDiagOffsets(diagOffsets_);
838  savedDiagOffsets_ = true;
839  }
840  // Now we're guaranteed to have cached diagonal offsets.
841  D_ = makeInverseDiagonal(*A_, true);
842  } else { // either A_ is not a CrsMatrix, or its graph is nonconst
843  D_ = makeInverseDiagonal(*A_);
844  }
845  }
846  } else { // the user provided an inverse diagonal
847  D_ = makeRangeMapVectorConst(userInvDiag_);
848  }
849 
850  // Have we estimated eigenvalues before?
851  const bool computedEigenvalueEstimates =
852  STS::isnaninf(computedLambdaMax_) || STS::isnaninf(computedLambdaMin_);
853 
854  // Only recompute the eigenvalue estimates if
855  // - we are supposed to assume that the matrix may have changed, or
856  // - they haven't been computed before, and the user hasn't given
857  // us at least an estimate of the max eigenvalue.
858  //
859  // We at least need an estimate of the max eigenvalue. This is the
860  // most important one if using Chebyshev as a smoother.
861 
862  if (!assumeMatrixUnchanged_ ||
863  (!computedEigenvalueEstimates && STS::isnaninf(userLambdaMax_))) {
864  ST computedLambdaMax;
865  if ((eigenAnalysisType_ == "power method") || (eigenAnalysisType_ == "power-method")) {
866  Teuchos::RCP<V> x;
867  if (eigVector_.is_null()) {
868  x = Teuchos::rcp(new V(A_->getDomainMap()));
869  if (eigKeepVectors_)
870  eigVector_ = x;
871  PowerMethod::computeInitialGuessForPowerMethod(*x, false);
872  } else
873  x = eigVector_;
874 
875  Teuchos::RCP<V> y;
876  if (eigVector2_.is_null()) {
877  y = rcp(new V(A_->getRangeMap()));
878  if (eigKeepVectors_)
879  eigVector2_ = y;
880  } else
881  y = eigVector2_;
882 
883  Teuchos::RCP<Teuchos::FancyOStream> stream = (debug_ ? out_ : Teuchos::null);
884  computedLambdaMax = PowerMethod::powerMethodWithInitGuess(*A_, *D_, eigMaxIters_, x, y,
885  eigRelTolerance_, eigNormalizationFreq_, stream,
886  computeSpectralRadius_);
887  } else {
888  computedLambdaMax = cgMethod(*A_, *D_, eigMaxIters_);
889  }
891  STS::isnaninf(computedLambdaMax),
892  std::runtime_error,
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_),
898  std::logic_error,
899  "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
900  << endl
901  << "This should be impossible." << endl
902  << "Please report this bug to the Ifpack2 developers.");
903 
904  // The power method doesn't estimate the min eigenvalue, so we
905  // do our best to provide an estimate. userEigRatio_ has a
906  // reasonable default value, and if the user provided it, we
907  // have already checked that its value is finite and >= 1.
908  const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
909 
910  // Defer "committing" results until all computations succeeded.
911  computedLambdaMax_ = computedLambdaMax;
912  computedLambdaMin_ = computedLambdaMin;
913  } else {
915  STS::isnaninf(userLambdaMax_) && STS::isnaninf(computedLambdaMax_),
916  std::logic_error,
917  "Ifpack2::Chebyshev::compute: " << endl
918  << "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
919  << endl
920  << "This should be impossible." << endl
921  << "Please report this bug to the Ifpack2 developers.");
922  }
923 
925  // Figure out the eigenvalue estimates that apply() will use.
927 
928  // Always favor the user's max eigenvalue estimate, if provided.
929  lambdaMaxForApply_ = STS::isnaninf(userLambdaMax_) ? computedLambdaMax_ : userLambdaMax_;
930 
931  // mfh 11 Feb 2013: For now, we imitate Ifpack by ignoring the
932  // user's min eigenvalue estimate, and using the given eigenvalue
933  // ratio to estimate the min eigenvalue. We could instead do this:
934  // favor the user's eigenvalue ratio estimate, but if it's not
935  // provided, use lambdaMax / lambdaMin. However, we want Chebyshev
936  // to have sensible smoother behavior if the user did not provide
937  // eigenvalue estimates. Ifpack's behavior attempts to push down
938  // the error terms associated with the largest eigenvalues, while
939  // expecting that users will only want a small number of iterations,
940  // so that error terms associated with the smallest eigenvalues
941  // won't grow too much. This is sensible behavior for a smoother.
942  lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
943  eigRatioForApply_ = userEigRatio_;
944 
945  if (chebyshevAlgorithm_ == "first") {
946  // Ifpack has a special-case modification of the eigenvalue bounds
947  // for the case where the max eigenvalue estimate is close to one.
948  const ST one = Teuchos::as<ST>(1);
949  // FIXME (mfh 20 Nov 2013) Should scale this 1.0e-6 term
950  // appropriately for MT's machine precision.
951  if (STS::magnitude(lambdaMaxForApply_ - one) < Teuchos::as<MT>(1.0e-6)) {
952  lambdaMinForApply_ = one;
953  lambdaMaxForApply_ = lambdaMinForApply_;
954  eigRatioForApply_ = one; // Ifpack doesn't include this line.
955  }
956  }
957 
958  // Allocate temporary vector
959  if (preAllocateTempVector_ && !D_.is_null()) {
960  makeTempMultiVector(*D_);
961  if (chebyshevAlgorithm_ == "fourth" || chebyshevAlgorithm_ == "opt_fourth") {
962  makeSecondTempMultiVector(*D_);
963  }
964  }
965 
966  if (chebyshevAlgorithm_ == "textbook") {
967  // no-op
968  } else {
969  if (ck_.is_null()) {
970  ck_ = Teuchos::rcp(new ChebyshevKernel<op_type>(A_, ckUseNativeSpMV_));
971  }
972  if (ckUseNativeSpMV_) {
973  ck_->setAuxiliaryVectors(1);
974  }
975  }
976 }
977 
978 template <class ScalarType, class MV>
979 ScalarType
981  getLambdaMaxForApply() const {
982  return lambdaMaxForApply_;
983 }
984 
985 template <class ScalarType, class MV>
987 Chebyshev<ScalarType, MV>::apply(const MV& B, MV& X) {
988  const char prefix[] = "Ifpack2::Chebyshev::apply: ";
989 
990  if (debug_) {
991  *out_ << "apply: " << std::endl;
992  }
993  TEUCHOS_TEST_FOR_EXCEPTION(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(STS::isnaninf(lambdaMaxForApply_), std::runtime_error,
997  prefix << "There is no estimate for the max eigenvalue."
998  << std::endl
999  << computeBeforeApplyReminder);
1000  TEUCHOS_TEST_FOR_EXCEPTION(STS::isnaninf(lambdaMinForApply_), std::runtime_error,
1001  prefix << "There is no estimate for the min eigenvalue."
1002  << std::endl
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."
1007  << std::endl
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."
1011  << std::endl
1012  << computeBeforeApplyReminder);
1013 
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_);
1019  } else {
1020  ifpackApplyImpl(*A_, B, X, numIters_, lambdaMaxForApply_,
1021  lambdaMinForApply_, eigRatioForApply_, *D_);
1022  }
1023 
1024  if (computeMaxResNorm_ && B.getNumVectors() > 0) {
1025  MV R(B.getMap(), B.getNumVectors());
1026  computeResidual(R, B, *A_, X);
1027  Teuchos::Array<MT> norms(B.getNumVectors());
1028  R.norm2(norms());
1029  return *std::max_element(norms.begin(), norms.end());
1030  } else {
1032  }
1033 }
1034 
1035 template <class ScalarType, class MV>
1037  print(std::ostream& out) {
1038  using Teuchos::rcpFromRef;
1039  this->describe(*(Teuchos::getFancyOStream(rcpFromRef(out))),
1041 }
1042 
1043 template <class ScalarType, class MV>
1046  const ScalarType& alpha,
1047  const V& D_inv,
1048  const MV& B,
1049  MV& X) {
1050  solve(W, alpha, D_inv, B); // W = alpha*D_inv*B
1051  Tpetra::deep_copy(X, W); // X = 0 + W
1052 }
1053 
1054 template <class ScalarType, class MV>
1056  computeResidual(MV& R, const MV& B, const op_type& A, const MV& X,
1057  const Teuchos::ETransp mode) {
1058  Tpetra::Details::residual(A, X, B, R);
1059 }
1060 
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());
1065 }
1066 
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());
1071 }
1072 
1073 template <class ScalarType, class MV>
1075 Chebyshev<ScalarType, MV>::
1076  makeInverseDiagonal(const row_matrix_type& A, const bool useDiagOffsets) const {
1077  using Teuchos::RCP;
1078  using Teuchos::rcp_dynamic_cast;
1079  using Teuchos::rcpFromRef;
1080 
1081  RCP<V> D_rowMap;
1082  if (!D_.is_null() &&
1083  D_->getMap()->isSameAs(*(A.getRowMap()))) {
1084  if (debug_)
1085  *out_ << "Reusing pre-existing vector for diagonal extraction" << std::endl;
1086  D_rowMap = Teuchos::rcp_const_cast<V>(D_);
1087  } else {
1088  D_rowMap = Teuchos::rcp(new V(A.getRowMap(), /*zeroOut=*/false));
1089  if (debug_)
1090  *out_ << "Allocated new vector for diagonal extraction" << std::endl;
1091  }
1092  if (useDiagOffsets) {
1093  // The optimizations below only work if A_ is a Tpetra::CrsMatrix.
1094  // We'll make our best guess about its type here, since we have no
1095  // way to get back the original fifth template parameter.
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>
1100  crs_matrix_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_);
1112  }
1113  } else {
1114  // This always works for a Tpetra::RowMatrix, even if it is not a
1115  // Tpetra::CrsMatrix. We just don't have offsets in this case.
1116  A.getLocalDiagCopy(*D_rowMap);
1117  }
1118  RCP<V> D_rangeMap = makeRangeMapVector(D_rowMap);
1119 
1120  if (debug_) {
1121  // In debug mode, make sure that all diagonal entries are
1122  // positive, on all processes. Note that *out_ only prints on
1123  // Process 0 of the matrix's communicator.
1124  bool foundNonpositiveValue = false;
1125  {
1126  auto D_lcl = D_rangeMap->getLocalViewHost(Tpetra::Access::ReadOnly);
1127  auto D_lcl_1d = Kokkos::subview(D_lcl, Kokkos::ALL(), 0);
1128 
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;
1133 
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;
1138  break;
1139  }
1140  }
1141  }
1142 
1143  using Teuchos::outArg;
1144  using Teuchos::REDUCE_MIN;
1145  using Teuchos::reduceAll;
1146 
1147  const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1148  int gblSuccess = lclSuccess; // to be overwritten
1149  if (!D_rangeMap->getMap().is_null() && D_rangeMap->getMap()->getComm().is_null()) {
1150  const Teuchos::Comm<int>& comm = *(D_rangeMap->getMap()->getComm());
1151  reduceAll<int, int>(comm, REDUCE_MIN, lclSuccess, outArg(gblSuccess));
1152  }
1153  if (gblSuccess == 1) {
1154  *out_ << "makeInverseDiagonal: The matrix's diagonal entries all have "
1155  "positive real part (this is good for Chebyshev)."
1156  << std::endl;
1157  } else {
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."
1161  << std::endl;
1162  }
1163  }
1164 
1165  // Invert the diagonal entries, replacing entries less (in
1166  // magnitude) than the user-specified value with that value.
1167  reciprocal_threshold(*D_rangeMap, minDiagVal_);
1168  return Teuchos::rcp_const_cast<const V>(D_rangeMap);
1169 }
1170 
1171 template <class ScalarType, class MV>
1173 Chebyshev<ScalarType, MV>::
1174  makeRangeMapVectorConst(const Teuchos::RCP<const V>& D) const {
1175  using Teuchos::RCP;
1176  using Teuchos::rcp;
1177  typedef Tpetra::Export<typename MV::local_ordinal_type,
1178  typename MV::global_ordinal_type,
1179  typename MV::node_type>
1180  export_type;
1181  // This throws logic_error instead of runtime_error, because the
1182  // methods that call makeRangeMapVector should all have checked
1183  // whether A_ is null before calling this method.
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.");
1195 
1196  RCP<const map_type> sourceMap = D->getMap();
1197  RCP<const map_type> rangeMap = A_->getRangeMap();
1198  RCP<const map_type> rowMap = A_->getRowMap();
1199 
1200  if (rangeMap->isSameAs(*sourceMap)) {
1201  // The given vector's Map is the same as the matrix's range Map.
1202  // That means we don't need to Export. This is the fast path.
1203  return D;
1204  } else { // We need to Export.
1205  RCP<const export_type> exporter;
1206  // Making an Export object from scratch is expensive enough that
1207  // it's worth the O(1) global reductions to call isSameAs(), to
1208  // see if we can avoid that cost.
1209  if (sourceMap->isSameAs(*rowMap)) {
1210  // We can reuse the matrix's Export object, if there is one.
1211  exporter = A_->getGraph()->getExporter();
1212  } else { // We have to make a new Export object.
1213  exporter = rcp(new export_type(sourceMap, rangeMap));
1214  }
1215 
1216  if (exporter.is_null()) {
1217  return D; // Row Map and range Map are the same; no need to Export.
1218  } else { // Row Map and range Map are _not_ the same; must Export.
1219  RCP<V> D_out = rcp(new V(*D, Teuchos::Copy));
1220  D_out->doExport(*D, *exporter, Tpetra::ADD);
1221  return Teuchos::rcp_const_cast<const V>(D_out);
1222  }
1223  }
1224 }
1225 
1226 template <class ScalarType, class MV>
1228 Chebyshev<ScalarType, MV>::
1229  makeRangeMapVector(const Teuchos::RCP<V>& D) const {
1230  using Teuchos::rcp_const_cast;
1231  return rcp_const_cast<V>(makeRangeMapVectorConst(rcp_const_cast<V>(D)));
1232 }
1233 
1234 template <class ScalarType, class MV>
1235 void Chebyshev<ScalarType, MV>::
1236  textbookApplyImpl(const op_type& A,
1237  const MV& B,
1238  MV& X,
1239  const int numIters,
1240  const ST lambdaMax,
1241  const ST lambdaMin,
1242  const ST eigRatio,
1243  const V& D_inv) const {
1244  (void)lambdaMin; // Forestall compiler warning.
1245  const ST myLambdaMin = lambdaMax / eigRatio;
1246 
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; // Ifpack2 calls this theta
1251  const ST c = (lambdaMax - myLambdaMin) / two; // Ifpack2 calls this 1/delta
1252 
1253  if (zeroStartingSolution_ && numIters > 0) {
1254  // If zero iterations, then input X is output X.
1255  X.putScalar(zero);
1256  }
1257  MV R(B.getMap(), B.getNumVectors(), false);
1258  MV P(B.getMap(), B.getNumVectors(), false);
1259  MV Z(B.getMap(), B.getNumVectors(), false);
1260  ST alpha, beta;
1261  for (int i = 0; i < numIters; ++i) {
1262  computeResidual(R, B, A, X); // R = B - A*X
1263  solve(Z, D_inv, R); // z = D_inv * R, that is, D \ R.
1264  if (i == 0) {
1265  P = Z;
1266  alpha = two / d;
1267  } else {
1268  // beta = (c * alpha / two)^2;
1269  // const ST sqrtBeta = c * alpha / two;
1270  // beta = sqrtBeta * sqrtBeta;
1271  beta = alpha * (c / two) * (c / two);
1272  alpha = one / (d - beta);
1273  P.update(one, Z, beta); // P = Z + beta*P
1274  }
1275  X.update(alpha, P, one); // X = X + alpha*P
1276  // If we compute the residual here, we could either do R = B -
1277  // A*X, or R = R - alpha*A*P. Since we choose the former, we
1278  // can move the computeResidual call to the top of the loop.
1279  }
1280 }
1281 
1282 template <class ScalarType, class MV>
1283 void Chebyshev<ScalarType, MV>::
1284  fourthKindApplyImpl(const op_type& A,
1285  const MV& B,
1286  MV& X,
1287  const int numIters,
1288  const ST lambdaMax,
1289  const V& D_inv) {
1290  // standard 4th kind Chebyshev smoother has \beta_i := 1
1291  std::vector<ScalarType> betas(numIters, 1.0);
1292  if (chebyshevAlgorithm_ == "opt_fourth") {
1293  betas = optimalWeightsImpl<ScalarType>(numIters);
1294  }
1295 
1296  const ST invEig = MT(1) / (lambdaMax * boostFactor_);
1297 
1298  // Fetch cached temporary (multi)vector.
1299  Teuchos::RCP<MV> Z_ptr = makeTempMultiVector(B);
1300  MV& Z = *Z_ptr;
1301 
1302  // Store 4th-kind result (needed as temporary for bootstrapping opt. 4th-kind Chebyshev)
1303  // Fetch the second cached temporary (multi)vector.
1304  Teuchos::RCP<MV> X4_ptr = makeSecondTempMultiVector(B);
1305  MV& X4 = *X4_ptr;
1306 
1307  // Special case for the first iteration.
1308  if (!zeroStartingSolution_) {
1309  // X4 = X
1310  Tpetra::deep_copy(X4, X);
1311 
1312  if (ck_.is_null()) {
1313  Teuchos::RCP<const op_type> A_op = A_;
1314  ck_ = Teuchos::rcp(new ChebyshevKernel<op_type>(A_op, ckUseNativeSpMV_));
1315  }
1316  // Z := (4/3 * invEig)*D_inv*(B-A*X4)
1317  // X4 := X4 + Z
1318  ck_->compute(Z, MT(4.0 / 3.0) * invEig, const_cast<V&>(D_inv),
1319  const_cast<MV&>(B), X4, STS::zero());
1320 
1321  // X := X + beta[0] * Z
1322  X.update(betas[0], Z, STS::one());
1323  } else {
1324  // Z := (4/3 * invEig)*D_inv*B and X := 0 + Z.
1325  firstIterationWithZeroStartingSolution(Z, MT(4.0 / 3.0) * invEig, D_inv, B, X4);
1326 
1327  // X := 0 + beta * Z
1328  X.update(betas[0], Z, STS::zero());
1329  }
1330 
1331  if (numIters > 1 && ck_.is_null()) {
1332  Teuchos::RCP<const op_type> A_op = A_;
1333  ck_ = Teuchos::rcp(new ChebyshevKernel<op_type>(A_op, ckUseNativeSpMV_));
1334  }
1335 
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;
1339 
1340  // Z := rScale*D_inv*(B - A*X4) + zScale*Z.
1341  // X4 := X4 + Z
1342  ck_->compute(Z, rScale, const_cast<V&>(D_inv),
1343  const_cast<MV&>(B), (X4), zScale);
1344 
1345  // X := X + beta[i] * Z
1346  X.update(betas[i], Z, STS::one());
1347  }
1348 }
1349 
1350 template <class ScalarType, class MV>
1351 typename Chebyshev<ScalarType, MV>::MT
1352 Chebyshev<ScalarType, MV>::maxNormInf(const MV& X) {
1353  Teuchos::Array<MT> norms(X.getNumVectors());
1354  X.normInf(norms());
1355  return *std::max_element(norms.begin(), norms.end());
1356 }
1357 
1358 template <class ScalarType, class MV>
1359 void Chebyshev<ScalarType, MV>::
1360  ifpackApplyImpl(const op_type& A,
1361  const MV& B,
1362  MV& X,
1363  const int numIters,
1364  const ST lambdaMax,
1365  const ST lambdaMin,
1366  const ST eigRatio,
1367  const V& D_inv) {
1368  using std::endl;
1369 #ifdef HAVE_IFPACK2_DEBUG
1370  const bool debug = debug_;
1371 #else
1372  const bool debug = false;
1373 #endif
1374 
1375  if (debug) {
1376  *out_ << " \\|B\\|_{\\infty} = " << maxNormInf(B) << endl;
1377  *out_ << " \\|X\\|_{\\infty} = " << maxNormInf(X) << endl;
1378  }
1379 
1380  if (numIters <= 0) {
1381  return;
1382  }
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);
1386 
1387  // Quick solve when the matrix A is the identity.
1388  if (lambdaMin == one && lambdaMax == lambdaMin) {
1389  solve(X, D_inv, B);
1390  return;
1391  }
1392 
1393  // Initialize coefficients
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;
1399 
1400  if (debug) {
1401  *out_ << " alpha = " << alpha << endl
1402  << " beta = " << beta << endl
1403  << " delta = " << delta << endl
1404  << " theta = " << theta << endl
1405  << " s1 = " << s1 << endl;
1406  }
1407 
1408  // Fetch cached temporary (multi)vector.
1409  Teuchos::RCP<MV> W_ptr = makeTempMultiVector(B);
1410  MV& W = *W_ptr;
1411 
1412  if (debug) {
1413  *out_ << " Iteration " << 1 << ":" << endl
1414  << " - \\|D\\|_{\\infty} = " << D_->normInf() << endl;
1415  }
1416 
1417  // Special case for the first iteration.
1418  if (!zeroStartingSolution_) {
1419  // mfh 22 May 2019: Tests don't actually exercise this path.
1420 
1421  if (ck_.is_null()) {
1422  Teuchos::RCP<const op_type> A_op = A_;
1423  ck_ = Teuchos::rcp(new ChebyshevKernel<op_type>(A_op, ckUseNativeSpMV_));
1424  }
1425  // W := (1/theta)*D_inv*(B-A*X) and X := X + W.
1426  // X := X + W
1427  ck_->compute(W, one / theta, const_cast<V&>(D_inv),
1428  const_cast<MV&>(B), X, zero);
1429  } else {
1430  // W := (1/theta)*D_inv*B and X := 0 + W.
1431  firstIterationWithZeroStartingSolution(W, one / theta, D_inv, B, X);
1432  }
1433 
1434  if (debug) {
1435  *out_ << " - \\|W\\|_{\\infty} = " << maxNormInf(W) << endl
1436  << " - \\|X\\|_{\\infty} = " << maxNormInf(X) << endl;
1437  }
1438 
1439  if (numIters > 1 && ck_.is_null()) {
1440  Teuchos::RCP<const op_type> A_op = A_;
1441  ck_ = Teuchos::rcp(new ChebyshevKernel<op_type>(A_op, ckUseNativeSpMV_));
1442  }
1443 
1444  // The rest of the iterations.
1445  ST rhok = one / s1;
1446  ST rhokp1, dtemp1, dtemp2;
1447  for (int deg = 1; deg < numIters; ++deg) {
1448  if (debug) {
1449  *out_ << " Iteration " << deg + 1 << ":" << endl
1450  << " - \\|D\\|_{\\infty} = " << D_->normInf() << endl
1451  << " - \\|B\\|_{\\infty} = " << maxNormInf(B) << endl
1452  << " - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm()
1453  << endl
1454  << " - rhok = " << rhok << endl;
1455  }
1456 
1457  rhokp1 = one / (two * s1 - rhok);
1458  dtemp1 = rhokp1 * rhok;
1459  dtemp2 = two * rhokp1 * delta;
1460  rhok = rhokp1;
1461 
1462  if (debug) {
1463  *out_ << " - dtemp1 = " << dtemp1 << endl
1464  << " - dtemp2 = " << dtemp2 << endl;
1465  }
1466 
1467  // W := dtemp2*D_inv*(B - A*X) + dtemp1*W.
1468  // X := X + W
1469  ck_->compute(W, dtemp2, const_cast<V&>(D_inv),
1470  const_cast<MV&>(B), (X), dtemp1);
1471 
1472  if (debug) {
1473  *out_ << " - \\|W\\|_{\\infty} = " << maxNormInf(W) << endl
1474  << " - \\|X\\|_{\\infty} = " << maxNormInf(X) << endl;
1475  }
1476  }
1477 }
1478 
1479 template <class ScalarType, class MV>
1480 typename Chebyshev<ScalarType, MV>::ST
1481 Chebyshev<ScalarType, MV>::
1482  cgMethodWithInitGuess(const op_type& A,
1483  const V& D_inv,
1484  const int numIters,
1485  V& r) {
1486  using std::endl;
1487  using MagnitudeType = typename STS::magnitudeType;
1488  if (debug_) {
1489  *out_ << " cgMethodWithInitGuess:" << endl;
1490  }
1491 
1492  const ST one = STS::one();
1493  ST beta, betaOld = one, pAp, pApOld = one, alpha, rHz, rHzOld, rHzOld2 = one, lambdaMax;
1494  // ST lambdaMin;
1495  Teuchos::ArrayRCP<MagnitudeType> diag, offdiag;
1496  Teuchos::RCP<V> p, z, Ap;
1497  diag.resize(numIters);
1498  offdiag.resize(numIters - 1);
1499 
1500  p = rcp(new V(A.getRangeMap()));
1501  z = rcp(new V(A.getRangeMap()));
1502  Ap = rcp(new V(A.getRangeMap()));
1503 
1504  // Tpetra::Details::residual (A, x, *b, *r);
1505  solve(*p, D_inv, r);
1506  rHz = r.dot(*p);
1507 
1508  for (int iter = 0; iter < numIters; ++iter) {
1509  if (debug_) {
1510  *out_ << " Iteration " << iter << endl;
1511  }
1512  A.apply(*p, *Ap);
1513  pAp = p->dot(*Ap);
1514  alpha = rHz / pAp;
1515  r.update(-alpha, *Ap, one);
1516  rHzOld = rHz;
1517  solve(*z, D_inv, r);
1518  rHz = r.dot(*z);
1519  beta = rHz / rHzOld;
1520  p->update(one, *z, beta);
1521  if (iter > 0) {
1522  diag[iter] = STS::real((betaOld * betaOld * pApOld + pAp) / rHzOld);
1523  offdiag[iter - 1] = -STS::real((betaOld * pApOld / (sqrt(rHzOld * rHzOld2))));
1524  if (debug_) {
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;
1530  }
1531  } else {
1532  diag[iter] = STS::real(pAp / rHzOld);
1533  if (debug_) {
1534  *out_ << " diag[" << iter << "] = " << diag[iter] << endl;
1535  *out_ << " rHz = " << rHz << endl;
1536  *out_ << " alpha = " << alpha << endl;
1537  *out_ << " beta = " << beta << endl;
1538  }
1539  }
1540  rHzOld2 = rHzOld;
1541  betaOld = beta;
1542  pApOld = pAp;
1543  }
1544 
1545  lambdaMax = LapackHelper<ST>::tri_diag_spectral_radius(diag, offdiag);
1546 
1547  return lambdaMax;
1548 }
1549 
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) {
1554  using std::endl;
1555 
1556  if (debug_) {
1557  *out_ << "cgMethod:" << endl;
1558  }
1559 
1560  Teuchos::RCP<V> r;
1561  if (eigVector_.is_null()) {
1562  r = rcp(new V(A.getDomainMap()));
1563  if (eigKeepVectors_)
1564  eigVector_ = r;
1565  // For CG, we need to get the BCs right and we'll use D_inv to get that
1566  Details::computeInitialGuessForCG(D_inv, *r);
1567  } else
1568  r = eigVector_;
1569 
1570  ST lambdaMax = cgMethodWithInitGuess(A, D_inv, numIters, *r);
1571 
1572  return lambdaMax;
1573 }
1574 
1575 template <class ScalarType, class MV>
1578  return A_;
1579 }
1580 
1581 template <class ScalarType, class MV>
1584  // Technically, this is true, because the matrix must be symmetric.
1585  return true;
1586 }
1587 
1588 template <class ScalarType, class MV>
1591  makeTempMultiVector(const MV& B) {
1592  // ETP 02/08/17: We must check not only if the temporary vectors are
1593  // null, but also if the number of columns match, since some multi-RHS
1594  // solvers (e.g., Belos) may call apply() with different numbers of columns.
1595 
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));
1599  }
1600  return W_;
1601 }
1602 
1603 template <class ScalarType, class MV>
1605 Chebyshev<ScalarType, MV>::
1606  makeSecondTempMultiVector(const MV& B) {
1607  // ETP 02/08/17: We must check not only if the temporary vectors are
1608  // null, but also if the number of columns match, since some multi-RHS
1609  // solvers (e.g., Belos) may call apply() with different numbers of columns.
1610 
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));
1614  }
1615  return W2_;
1616 }
1617 
1618 template <class ScalarType, class MV>
1619 std::string
1621  description() const {
1622  std::ostringstream oss;
1623  // YAML requires quoting the key in this case, to distinguish
1624  // key's colons from the colon that separates key from value.
1625  oss << "\"Ifpack2::Details::Chebyshev\":"
1626  << "{"
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";
1635  oss << "}";
1636  return oss.str();
1637 }
1638 
1639 template <class ScalarType, class MV>
1642  const Teuchos::EVerbosityLevel verbLevel) const {
1643  using std::endl;
1645 
1646  const Teuchos::EVerbosityLevel vl =
1647  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1648  if (vl == Teuchos::VERB_NONE) {
1649  return; // print NOTHING
1650  }
1651 
1652  // By convention, describe() starts with a tab.
1653  //
1654  // This does affect all processes on which it's valid to print to
1655  // 'out'. However, it does not actually print spaces to 'out'
1656  // unless operator<< gets called, so it's safe to use on all
1657  // processes.
1658  Teuchos::OSTab tab0(out);
1659 
1660  // We only print on Process 0 of the matrix's communicator. If
1661  // the matrix isn't set, we don't have a communicator, so we have
1662  // to assume that every process can print.
1663  int myRank = -1;
1664  if (A_.is_null() || A_->getComm().is_null()) {
1665  myRank = 0;
1666  } else {
1667  myRank = A_->getComm()->getRank();
1668  }
1669  if (myRank == 0) {
1670  // YAML requires quoting the key in this case, to distinguish
1671  // key's colons from the colon that separates key from value.
1672  out << "\"Ifpack2::Details::Chebyshev\":" << endl;
1673  }
1674  Teuchos::OSTab tab1(out);
1675 
1676  if (vl == Teuchos::VERB_LOW) {
1677  if (myRank == 0) {
1678  out << "degree: " << numIters_ << endl
1679  << "lambdaMax: " << lambdaMaxForApply_ << endl
1680  << "alpha: " << eigRatioForApply_ << endl
1681  << "lambdaMin: " << lambdaMinForApply_ << endl
1682  << "boost factor: " << boostFactor_ << endl;
1683  }
1684  return;
1685  }
1686 
1687  // vl > Teuchos::VERB_LOW
1688 
1689  if (myRank == 0) {
1690  out << "Template parameters:" << endl;
1691  {
1692  Teuchos::OSTab tab2(out);
1693  out << "ScalarType: " << TypeNameTraits<ScalarType>::name() << endl
1694  << "MV: " << TypeNameTraits<MV>::name() << endl;
1695  }
1696 
1697  // "Computed parameters" literally means "parameters whose
1698  // values were computed by compute()."
1699  if (myRank == 0) {
1700  out << "Computed parameters:" << endl;
1701  }
1702  }
1703 
1704  // Print computed parameters
1705  {
1706  Teuchos::OSTab tab2(out);
1707  // Users might want to see the values in the computed inverse
1708  // diagonal, so we print them out at the highest verbosity.
1709  if (myRank == 0) {
1710  out << "D_: ";
1711  }
1712  if (D_.is_null()) {
1713  if (myRank == 0) {
1714  out << "unset" << endl;
1715  }
1716  } else if (vl <= Teuchos::VERB_HIGH) {
1717  if (myRank == 0) {
1718  out << "set" << endl;
1719  }
1720  } else { // D_ not null and vl > Teuchos::VERB_HIGH
1721  if (myRank == 0) {
1722  out << endl;
1723  }
1724  // By convention, describe() first indents, then prints.
1725  // We can rely on other describe() implementations to do that.
1726  D_->describe(out, vl);
1727  }
1728  if (myRank == 0) {
1729  // W_ is scratch space; its values are irrelevant.
1730  // All that matters is whether or not they have been set.
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;
1737  }
1738  } // print computed parameters
1739 
1740  if (myRank == 0) {
1741  out << "User parameters:" << endl;
1742  }
1743 
1744  // Print user parameters
1745  {
1746  Teuchos::OSTab tab2(out);
1747  out << "userInvDiag_: ";
1748  if (userInvDiag_.is_null()) {
1749  out << "unset" << endl;
1750  } else if (vl <= Teuchos::VERB_HIGH) {
1751  out << "set" << endl;
1752  } else { // userInvDiag_ not null and vl > Teuchos::VERB_HIGH
1753  if (myRank == 0) {
1754  out << endl;
1755  }
1756  userInvDiag_->describe(out, vl);
1757  }
1758  if (myRank == 0) {
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;
1770  }
1771  } // print user parameters
1772 }
1773 
1774 } // namespace Details
1775 } // namespace Ifpack2
1776 
1777 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S, LO, GO, N) \
1778  template class Ifpack2::Details::Chebyshev<S, Tpetra::MultiVector<S, LO, GO, N>>;
1779 
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
T * getRawPtr() 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&#39;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
bool is_null() const