Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_Details_Chebyshev_def.hpp
Go to the documentation of this file.
1 /*
2 //@HEADER
3 // ***********************************************************************
4 //
5 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
6 // Copyright (2009) Sandia Corporation
7 //
8 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
9 // license for use of this work by or on behalf of the U.S. Government.
10 //
11 // Redistribution and use in source and binary forms, with or without
12 // modification, are permitted provided that the following conditions are
13 // met:
14 //
15 // 1. Redistributions of source code must retain the above copyright
16 // notice, this list of conditions and the following disclaimer.
17 //
18 // 2. Redistributions in binary form must reproduce the above copyright
19 // notice, this list of conditions and the following disclaimer in the
20 // documentation and/or other materials provided with the distribution.
21 //
22 // 3. Neither the name of the Corporation nor the names of the
23 // contributors may be used to endorse or promote products derived from
24 // this software without specific prior written permission.
25 //
26 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
27 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
28 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
29 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
30 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
31 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
32 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
33 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
34 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
35 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
36 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
37 //
38 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
39 //
40 // ***********************************************************************
41 //@HEADER
42 */
43 
44 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
45 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
46 
53 
55 // #include "Ifpack2_Details_ScaledDampedResidual.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"
63 #include <cmath>
64 #include <iostream>
65 
66 namespace Ifpack2 {
67 namespace Details {
68 
69 namespace { // (anonymous)
70 
71 // We use this text a lot in error messages.
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()).";
81 
82 } // namespace (anonymous)
83 
84 // ReciprocalThreshold stuff below needs to be in a namspace visible outside
85 // of this file
86 template<class XV, class SizeType = typename XV::size_type>
87 struct V_ReciprocalThresholdSelfFunctor
88 {
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;
94 
95  XV X_;
96  const value_type minVal_;
97  const mag_type minValMag_;
98 
99  V_ReciprocalThresholdSelfFunctor (const XV& X,
100  const value_type& min_val) :
101  X_ (X),
102  minVal_ (min_val),
103  minValMag_ (KAT::abs (min_val))
104  {}
105 
106  KOKKOS_INLINE_FUNCTION
107  void operator() (const size_type& i) const
108  {
109  const mag_type X_i_abs = KAT::abs (X_[i]);
110 
111  if (X_i_abs < minValMag_) {
112  X_[i] = minVal_;
113  }
114  else {
115  X_[i] = KAT::one () / X_[i];
116  }
117  }
118 };
119 
120 template<class XV, class SizeType = typename XV::size_type>
121 struct LocalReciprocalThreshold {
122  static void
123  compute (const XV& X,
124  const typename XV::non_const_value_type& minVal)
125  {
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);
130  }
131 };
132 
133 template <class TpetraVectorType,
134  const bool classic = TpetraVectorType::node_type::classic>
135 struct GlobalReciprocalThreshold {};
136 
137 template <class TpetraVectorType>
138 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
139  static void
140  compute (TpetraVectorType& V,
141  const typename TpetraVectorType::scalar_type& min_val)
142  {
143  typedef typename TpetraVectorType::scalar_type scalar_type;
144  typedef typename TpetraVectorType::mag_type mag_type;
145  typedef Kokkos::Details::ArithTraits<scalar_type> STS;
146 
147  const scalar_type ONE = STS::one ();
148  const mag_type min_val_abs = STS::abs (min_val);
149 
150  Teuchos::ArrayRCP<scalar_type> V_0 = V.getDataNonConst (0);
151  const size_t lclNumRows = V.getLocalLength ();
152 
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) {
156  V_0[i] = min_val;
157  } else {
158  V_0[i] = ONE / V_0i;
159  }
160  }
161  }
162 };
163 
164 template <class TpetraVectorType>
165 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
166  static void
167  compute (TpetraVectorType& X,
168  const typename TpetraVectorType::scalar_type& minVal)
169  {
170  typedef typename TpetraVectorType::impl_scalar_type value_type;
171 
172  const value_type minValS = static_cast<value_type> (minVal);
173  auto X_0 = Kokkos::subview (X.getLocalViewDevice (Tpetra::Access::ReadWrite),
174  Kokkos::ALL (), 0);
175  LocalReciprocalThreshold<decltype (X_0) >::compute (X_0, minValS);
176  }
177 };
178 
179 // Utility function for inverting diagonal with threshold.
180 template <typename S, typename L, typename G, typename N>
181 void
182 reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V, const S& minVal)
183 {
184  GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
185 }
186 
187 namespace { // (anonymous)
188 
189 // Functor for making sure the real parts of all entries of a vector
190 // are nonnegative. We use this in computeInitialGuessForPowerMethod
191 // below.
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.");
202 public:
203  PositivizeVector (const OneDViewType& x) : x_ (x) {}
204 
205  KOKKOS_INLINE_FUNCTION void
206  operator () (const LocalOrdinal& i) const
207  {
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;
211 
212  if (STS::real (x_(i)) < STM::zero ()) {
213  x_(i) = -x_(i);
214  }
215  }
216 
217 private:
218  OneDViewType x_;
219 };
220 
221 } // namespace (anonymous)
222 
223 
224 template<class ScalarType, const bool lapackSupportsScalarType = LapackSupportsScalar<ScalarType>::value>
225 struct LapackHelper{
226  static
227  ScalarType
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.");
231  }
232 };
233 
234 
235 template<class ScalarType>
236 struct LapackHelper<ScalarType,true> {
237  static
238  ScalarType
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;
243  int info = 0;
244  const int N = diag.size ();
245  ScalarType scalar_dummy;
246  std::vector<MagnitudeType> mag_dummy(4*N);
247  char char_N = 'N';
248 
249  // lambdaMin = one;
250  ScalarType lambdaMax = STS::one();
251  if( N > 2 ) {
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.");
260  // lambdaMin = Teuchos::as<ScalarType> (diag[N-1]);
261  lambdaMax = Teuchos::as<ScalarType> (diag[0]);
262  }
263  return lambdaMax;
264  }
265 };
266 
267 
268 template<class ScalarType, class MV>
269 void Chebyshev<ScalarType, MV>::checkInputMatrix () const
270 {
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.");
277 
278  // In debug mode, test that the domain and range Maps of the matrix
279  // are the same.
280  if (debug_ && ! A_.is_null ()) {
281  Teuchos::RCP<const map_type> domainMap = A_->getDomainMap ();
282  Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap ();
283 
284  // isSameAs is a collective, but if the two pointers are the same,
285  // isSameAs will assume that they are the same on all processes, and
286  // return true without an all-reduce.
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.");
292  }
293 }
294 
295 
296 template<class ScalarType, class MV>
297 void
298 Chebyshev<ScalarType, MV>::
299 checkConstructorInput () const
300 {
301  // mfh 12 Aug 2016: The if statement avoids an "unreachable
302  // statement" warning for the checkInputMatrix() call, when
303  // STS::isComplex is false.
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).");
314  }
315  else {
316  checkInputMatrix ();
317  }
318 }
319 
320 template<class ScalarType, class MV>
322 Chebyshev (Teuchos::RCP<const row_matrix_type> A) :
323  A_ (A),
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 ()),
332  userEigRatio_ (Teuchos::as<ST> (30)),
333  minDiagVal_ (STS::eps ()),
334  numIters_ (1),
335  eigMaxIters_ (10),
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),
344  debug_ (false)
345 {
346  checkConstructorInput ();
347 }
348 
349 template<class ScalarType, class MV>
351 Chebyshev (Teuchos::RCP<const row_matrix_type> A,
352  Teuchos::ParameterList& params) :
353  A_ (A),
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 ()),
363  userEigRatio_ (Teuchos::as<ST> (30)),
364  minDiagVal_ (STS::eps ()),
365  numIters_ (1),
366  eigMaxIters_ (10),
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),
375  debug_ (false)
376 {
377  checkConstructorInput ();
378  setParameters (params);
379 }
380 
381 template<class ScalarType, class MV>
382 void
384 setParameters (Teuchos::ParameterList& plist)
385 {
386  using Teuchos::RCP;
387  using Teuchos::rcp;
388  using Teuchos::rcp_const_cast;
389 
390  // Note to developers: The logic for this method is complicated,
391  // because we want to accept Ifpack and ML parameters whenever
392  // possible, but we don't want to add their default values to the
393  // user's ParameterList. That's why we do all the isParameter()
394  // checks, instead of using the two-argument version of get()
395  // everywhere. The min and max eigenvalue parameters are also a
396  // special case, because we decide whether or not to do eigenvalue
397  // analysis based on whether the user supplied the max eigenvalue.
398 
399  // Default values of all the parameters.
400  const ST defaultLambdaMax = STS::nan ();
401  const ST defaultLambdaMin = STS::nan ();
402  // 30 is Ifpack::Chebyshev's default. ML has a different default
403  // eigRatio for smoothers and the coarse-grid solve (if using
404  // Chebyshev for that). The former uses 20; the latter uses 30.
405  // We're testing smoothers here, so use 20. (However, if you give
406  // ML an Epetra matrix, it will use Ifpack for Chebyshev, in which
407  // case it would defer to Ifpack's default settings.)
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; // Ifpack::Chebyshev default
417  const bool defaultAssumeMatrixUnchanged = false;
418  const bool defaultTextbookAlgorithm = false;
419  const bool defaultComputeMaxResNorm = false;
420  const bool defaultDebug = false;
421 
422  // We'll set the instance data transactionally, after all reads
423  // from the ParameterList. That way, if any of the ParameterList
424  // reads fail (e.g., due to the wrong parameter type), we will not
425  // have left the instance data in a half-changed state.
426  RCP<const V> userInvDiagCopy; // if nonnull: deep copy of user's Vector
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;
442 
443  // Fetch the parameters from the ParameterList. Defer all
444  // externally visible side effects until we have finished all
445  // ParameterList interaction. This makes the method satisfy the
446  // strong exception guarantee.
447 
448  if (plist.isType<bool> ("debug")) {
449  debug = plist.get<bool> ("debug");
450  }
451  else if (plist.isType<int> ("debug")) {
452  const int debugInt = plist.get<bool> ("debug");
453  debug = debugInt != 0;
454  }
455 
456  // Get the user-supplied inverse diagonal.
457  //
458  // Check for a raw pointer (const V* or V*), for Ifpack
459  // compatibility, as well as for RCP<const V>, RCP<V>, const V, or
460  // V. We'll make a deep copy of the vector at the end of this
461  // method anyway, so its const-ness doesn't matter. We handle the
462  // latter two cases ("const V" or "V") specially (copy them into
463  // userInvDiagCopy first, which is otherwise null at the end of the
464  // long if-then chain) to avoid an extra copy.
465 
466  const char opInvDiagLabel[] = "chebyshev: operator inv diagonal";
467  if (plist.isParameter (opInvDiagLabel)) {
468  // Pointer to the user's Vector, if provided.
469  RCP<const V> userInvDiag;
470 
471  if (plist.isType<const V*> (opInvDiagLabel)) {
472  const V* rawUserInvDiag =
473  plist.get<const V*> (opInvDiagLabel);
474  // Nonowning reference (we'll make a deep copy below)
475  userInvDiag = rcp (rawUserInvDiag, false);
476  }
477  else if (plist.isType<const V*> (opInvDiagLabel)) {
478  V* rawUserInvDiag = plist.get<V*> (opInvDiagLabel);
479  // Nonowning reference (we'll make a deep copy below)
480  userInvDiag = rcp (const_cast<const V*> (rawUserInvDiag), false);
481  }
482  else if (plist.isType<RCP<const V>> (opInvDiagLabel)) {
483  userInvDiag = plist.get<RCP<const V> > (opInvDiagLabel);
484  }
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);
489  }
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;
494  }
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;
500  }
501 
502  // NOTE: If the user's parameter has some strange type that we
503  // didn't test above, userInvDiag might still be null. You may
504  // want to add an error test for this condition. Currently, we
505  // just say in this case that the user didn't give us a Vector.
506 
507  // If we have userInvDiag but don't have a deep copy yet, make a
508  // deep copy now.
509  if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
510  userInvDiagCopy = rcp (new V (*userInvDiag, Teuchos::Copy));
511  }
512 
513  // NOTE: userInvDiag, if provided, is a row Map version of the
514  // Vector. We don't necessarily have a range Map yet. compute()
515  // would be the proper place to compute the range Map version of
516  // userInvDiag.
517  }
518 
519  // Don't fill in defaults for the max or min eigenvalue, because
520  // this class uses the existence of those parameters to determine
521  // whether it should do eigenanalysis.
522  if (plist.isParameter ("chebyshev: max eigenvalue")) {
523  if (plist.isType<double>("chebyshev: max eigenvalue"))
524  lambdaMax = plist.get<double> ("chebyshev: max eigenvalue");
525  else
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.");
532  }
533  if (plist.isParameter ("chebyshev: min eigenvalue")) {
534  if (plist.isType<double>("chebyshev: min eigenvalue"))
535  lambdaMin = plist.get<double> ("chebyshev: min eigenvalue");
536  else
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.");
543  }
544 
545  // Only fill in Ifpack2's name for the default parameter, not ML's.
546  if (plist.isParameter ("smoother: Chebyshev alpha")) { // ML compatibility
547  if (plist.isType<double>("smoother: Chebyshev alpha"))
548  eigRatio = plist.get<double> ("smoother: Chebyshev alpha");
549  else
550  eigRatio = plist.get<ST> ("smoother: Chebyshev alpha");
551  }
552  // Ifpack2's name overrides ML's name.
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 "
559  "a finite value.");
560  // mfh 11 Feb 2013: This class is currently only correct for real
561  // Scalar types, but we still want it to build for complex Scalar
562  // type so that users of Ifpack2::Factory can build their
563  // executables for real or complex Scalar type. Thus, we take the
564  // real parts here, which are always less-than comparable.
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 << ".");
571 
572  // See Github Issue #234. This parameter may be either MT
573  // (preferred) or double. We check both.
574  {
575  const char paramName[] = "chebyshev: boost factor";
576 
577  if (plist.isParameter (paramName)) {
578  if (plist.isType<MT> (paramName)) { // MT preferred
579  boostFactor = plist.get<MT> (paramName);
580  }
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);
585  }
586  else {
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.");
591  }
592  }
593  else { // parameter not in the list
594  // mfh 12 Aug 2016: To preserve current behavior (that fills in
595  // any parameters not in the input ParameterList with their
596  // default values), we call set() here. I don't actually like
597  // this behavior; I prefer the Belos model, where the input
598  // ParameterList is a delta from current behavior. However, I
599  // don't want to break things.
600  plist.set (paramName, defaultBoostFactor);
601  }
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 << ".");
606  }
607 
608  // Same name in Ifpack2 and Ifpack.
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.");
615 
616  // Only fill in Ifpack2's name, not ML's or Ifpack's.
617  if (plist.isParameter ("smoother: sweeps")) { // ML compatibility
618  numIters = plist.get<int> ("smoother: sweeps");
619  } // Ifpack's name overrides ML's name.
620  if (plist.isParameter ("relaxation: sweeps")) { // Ifpack compatibility
621  numIters = plist.get<int> ("relaxation: sweeps");
622  } // Ifpack2's name overrides Ifpack's name.
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 << ".");
629 
630  // The last parameter name overrides the first.
631  if (plist.isParameter ("eigen-analysis: iterations")) { // ML compatibility
632  eigMaxIters = plist.get<int> ("eigen-analysis: iterations");
633  } // Ifpack2's name overrides ML's name.
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 << ".");
640 
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"));
647 
648  eigKeepVectors = plist.get ("chebyshev: eigenvalue keep vectors", eigKeepVectors);
649 
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 << ".")
656 
657  zeroStartingSolution = plist.get ("chebyshev: zero starting solution",
658  zeroStartingSolution);
659  assumeMatrixUnchanged = plist.get ("chebyshev: assume matrix does not change",
660  assumeMatrixUnchanged);
661 
662  // We don't want to fill these parameters in, because they shouldn't
663  // be visible to Ifpack2::Chebyshev users.
664  if (plist.isParameter ("chebyshev: textbook algorithm")) {
665  textbookAlgorithm = plist.get<bool> ("chebyshev: textbook algorithm");
666  }
667  if (plist.isParameter ("chebyshev: compute max residual norm")) {
668  computeMaxResNorm = plist.get<bool> ("chebyshev: compute max residual norm");
669  }
670 
671  // Test for Ifpack parameters that we won't ever implement here.
672  // Be careful to use the one-argument version of get(), since the
673  // two-argment version adds the parameter if it's not there.
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.");
689 
690  // Test for Ifpack parameters that we haven't implemented yet.
691  //
692  // For now, we only check that this ML parameter, if provided, has
693  // the one value that we support. We consider other values "invalid
694  // arguments" rather than "logic errors," because Ifpack does not
695  // implement eigenanalyses other than the power method.
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\".");
705  }
706 
707  // We've validated all the parameters, so it's safe now to "commit" them.
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;
724  debug_ = debug;
725 
726  if (debug_) {
727  // Only print if myRank == 0.
728  int myRank = -1;
729  if (A_.is_null () || A_->getComm ().is_null ()) {
730  // We don't have a communicator (yet), so we assume that
731  // everybody can print. Revise this expectation in setMatrix().
732  myRank = 0;
733  }
734  else {
735  myRank = A_->getComm ()->getRank ();
736  }
737 
738  if (myRank == 0) {
739  out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
740  }
741  else {
742  using Teuchos::oblackholestream; // prints nothing
743  RCP<oblackholestream> blackHole (new oblackholestream ());
744  out_ = Teuchos::getFancyOStream (blackHole);
745  }
746  }
747  else { // NOT debug
748  // free the "old" output stream, if there was one
749  out_ = Teuchos::null;
750  }
751 }
752 
753 
754 template<class ScalarType, class MV>
755 void
757 {
758  ck_ = Teuchos::null;
759  D_ = Teuchos::null;
760  diagOffsets_ = offsets_type ();
761  savedDiagOffsets_ = false;
762  W_ = Teuchos::null;
763  computedLambdaMax_ = STS::nan ();
764  computedLambdaMin_ = STS::nan ();
765  eigVector_ = Teuchos::null;
766  eigVector2_ = Teuchos::null;
767 }
768 
769 
770 template<class ScalarType, class MV>
771 void
773 setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
774 {
775  if (A.getRawPtr () != A_.getRawPtr ()) {
776  if (! assumeMatrixUnchanged_) {
777  reset ();
778  }
779  A_ = A;
780  ck_ = Teuchos::null; // constructed on demand
781 
782  // The communicator may have changed, or we may not have had a
783  // communicator before. Thus, we may have to reset the debug
784  // output stream.
785  if (debug_) {
786  // Only print if myRank == 0.
787  int myRank = -1;
788  if (A.is_null () || A->getComm ().is_null ()) {
789  // We don't have a communicator (yet), so we assume that
790  // everybody can print. Revise this expectation in setMatrix().
791  myRank = 0;
792  }
793  else {
794  myRank = A->getComm ()->getRank ();
795  }
796 
797  if (myRank == 0) {
798  out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
799  }
800  else {
801  Teuchos::RCP<Teuchos::oblackholestream> blackHole (new Teuchos::oblackholestream ());
802  out_ = Teuchos::getFancyOStream (blackHole); // print nothing on other processes
803  }
804  }
805  else { // NOT debug
806  // free the "old" output stream, if there was one
807  out_ = Teuchos::null;
808  }
809  }
810 }
811 
812 template<class ScalarType, class MV>
813 void
815 {
816  using std::endl;
817  // Some of the optimizations below only work if A_ is a
818  // Tpetra::CrsMatrix. We'll make our best guess about its type
819  // here, since we have no way to get back the original fifth
820  // template parameter.
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;
825 
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.");
830 
831  // If A_ is a CrsMatrix and its graph is constant, we presume that
832  // the user plans to reuse the structure of A_, but possibly change
833  // A_'s values before each compute() call. This is the intended use
834  // case for caching the offsets of the diagonal entries of A_, to
835  // speed up extraction of diagonal entries on subsequent compute()
836  // calls.
837 
838  // FIXME (mfh 22 Jan 2013, 10 Feb 2013) In all cases when we use
839  // isnaninf() in this method, we really only want to check if the
840  // number is NaN. Inf means something different. However,
841  // Teuchos::ScalarTraits doesn't distinguish the two cases.
842 
843  // makeInverseDiagonal() returns a range Map Vector.
844  if (userInvDiag_.is_null ()) {
845  Teuchos::RCP<const crs_matrix_type> A_crsMat =
846  Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
847 
848  if (D_.is_null ()) { // We haven't computed D_ before
849  if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
850  // It's a CrsMatrix with a const graph; cache diagonal offsets.
851  const size_t lclNumRows = A_crsMat->getNodeNumRows ();
852  if (diagOffsets_.extent (0) < lclNumRows) {
853  diagOffsets_ = offsets_type (); // clear first to save memory
854  diagOffsets_ = offsets_type ("offsets", lclNumRows);
855  }
856  A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
857  savedDiagOffsets_ = true;
858  D_ = makeInverseDiagonal (*A_, true);
859  }
860  else { // either A_ is not a CrsMatrix, or its graph is nonconst
861  D_ = makeInverseDiagonal (*A_);
862  }
863  }
864  else if (! assumeMatrixUnchanged_) { // D_ exists but A_ may have changed
865  if (! A_crsMat.is_null () && A_crsMat->isFillComplete ()) {
866  // It's a CrsMatrix with a const graph; cache diagonal offsets
867  // if we haven't already.
868  if (! savedDiagOffsets_) {
869  const size_t lclNumRows = A_crsMat->getNodeNumRows ();
870  if (diagOffsets_.extent (0) < lclNumRows) {
871  diagOffsets_ = offsets_type (); // clear first to save memory
872  diagOffsets_ = offsets_type ("offsets", lclNumRows);
873  }
874  A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
875  savedDiagOffsets_ = true;
876  }
877  // Now we're guaranteed to have cached diagonal offsets.
878  D_ = makeInverseDiagonal (*A_, true);
879  }
880  else { // either A_ is not a CrsMatrix, or its graph is nonconst
881  D_ = makeInverseDiagonal (*A_);
882  }
883  }
884  }
885  else { // the user provided an inverse diagonal
886  D_ = makeRangeMapVectorConst (userInvDiag_);
887  }
888 
889  // Have we estimated eigenvalues before?
890  const bool computedEigenvalueEstimates =
891  STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
892 
893  // Only recompute the eigenvalue estimates if
894  // - we are supposed to assume that the matrix may have changed, or
895  // - they haven't been computed before, and the user hasn't given
896  // us at least an estimate of the max eigenvalue.
897  //
898  // We at least need an estimate of the max eigenvalue. This is the
899  // most important one if using Chebyshev as a smoother.
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_);
905  else
906  computedLambdaMax = cgMethod (*A_, *D_, eigMaxIters_);
907  TEUCHOS_TEST_FOR_EXCEPTION(
908  STS::isnaninf (computedLambdaMax),
909  std::runtime_error,
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_),
915  std::logic_error,
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.");
919 
920  // The power method doesn't estimate the min eigenvalue, so we
921  // do our best to provide an estimate. userEigRatio_ has a
922  // reasonable default value, and if the user provided it, we
923  // have already checked that its value is finite and >= 1.
924  const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
925 
926  // Defer "committing" results until all computations succeeded.
927  computedLambdaMax_ = computedLambdaMax;
928  computedLambdaMin_ = computedLambdaMin;
929  } else {
930  TEUCHOS_TEST_FOR_EXCEPTION(
931  STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
932  std::logic_error,
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.");
937  }
938 
940  // Figure out the eigenvalue estimates that apply() will use.
942 
943  // Always favor the user's max eigenvalue estimate, if provided.
944  lambdaMaxForApply_ = STS::isnaninf (userLambdaMax_) ? computedLambdaMax_ : userLambdaMax_;
945 
946  // mfh 11 Feb 2013: For now, we imitate Ifpack by ignoring the
947  // user's min eigenvalue estimate, and using the given eigenvalue
948  // ratio to estimate the min eigenvalue. We could instead do this:
949  // favor the user's eigenvalue ratio estimate, but if it's not
950  // provided, use lambdaMax / lambdaMin. However, we want Chebyshev
951  // to have sensible smoother behavior if the user did not provide
952  // eigenvalue estimates. Ifpack's behavior attempts to push down
953  // the error terms associated with the largest eigenvalues, while
954  // expecting that users will only want a small number of iterations,
955  // so that error terms associated with the smallest eigenvalues
956  // won't grow too much. This is sensible behavior for a smoother.
957  lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
958  eigRatioForApply_ = userEigRatio_;
959 
960  if (! textbookAlgorithm_) {
961  // Ifpack has a special-case modification of the eigenvalue bounds
962  // for the case where the max eigenvalue estimate is close to one.
963  const ST one = Teuchos::as<ST> (1);
964  // FIXME (mfh 20 Nov 2013) Should scale this 1.0e-6 term
965  // appropriately for MT's machine precision.
966  if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
967  lambdaMinForApply_ = one;
968  lambdaMaxForApply_ = lambdaMinForApply_;
969  eigRatioForApply_ = one; // Ifpack doesn't include this line.
970  }
971  }
972 }
973 
974 
975 template<class ScalarType, class MV>
976 ScalarType
978 getLambdaMaxForApply() const {
979  return lambdaMaxForApply_;
980 }
981 
982 
983 template<class ScalarType, class MV>
986 {
987  const char prefix[] = "Ifpack2::Chebyshev::apply: ";
988 
989  if (debug_) {
990  *out_ << "apply: " << std::endl;
991  }
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);
1013 
1014  if (textbookAlgorithm_) {
1015  textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1016  lambdaMinForApply_, eigRatioForApply_, *D_);
1017  }
1018  else {
1019  ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
1020  lambdaMinForApply_, eigRatioForApply_, *D_);
1021  }
1022 
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 ());
1027  R.norm2 (norms ());
1028  return *std::max_element (norms.begin (), norms.end ());
1029  }
1030  else {
1031  return Teuchos::ScalarTraits<MT>::zero ();
1032  }
1033 }
1034 
1035 template<class ScalarType, class MV>
1036 void
1038 print (std::ostream& out)
1039 {
1040  using Teuchos::rcpFromRef;
1041  this->describe (* (Teuchos::getFancyOStream (rcpFromRef (out))),
1042  Teuchos::VERB_MEDIUM);
1043 }
1044 
1045 template<class ScalarType, class MV>
1046 void
1049  const ScalarType& alpha,
1050  const V& D_inv,
1051  const MV& B,
1052  MV& X)
1053 {
1054  solve (W, alpha, D_inv, B); // W = alpha*D_inv*B
1055  Tpetra::deep_copy (X, W); // X = 0 + W
1056 }
1057 
1058 template<class ScalarType, class MV>
1059 void
1061 computeResidual (MV& R, const MV& B, const op_type& A, const MV& X,
1062  const Teuchos::ETransp mode)
1063 {
1064  Tpetra::Details::residual(A,X,B,R);
1065 }
1066 
1067 template <class ScalarType, class MV>
1068 void
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());
1072 }
1073 
1074 template<class ScalarType, class MV>
1075 void
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());
1079 }
1080 
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
1085 {
1086  using Teuchos::RCP;
1087  using Teuchos::rcpFromRef;
1088  using Teuchos::rcp_dynamic_cast;
1089 
1090  RCP<V> D_rowMap;
1091  if (!D_.is_null() &&
1092  D_->getMap()->isSameAs(*(A.getGraph ()->getRowMap ()))) {
1093  if (debug_)
1094  *out_ << "Reusing pre-existing vector for diagonal extraction" << std::endl;
1095  D_rowMap = Teuchos::rcp_const_cast<V>(D_);
1096  } else {
1097  D_rowMap = Teuchos::rcp(new V (A.getGraph ()->getRowMap (), /*zeroOut=*/false));
1098  if (debug_)
1099  *out_ << "Allocated new vector for diagonal extraction" << std::endl;
1100  }
1101  if (useDiagOffsets) {
1102  // The optimizations below only work if A_ is a Tpetra::CrsMatrix.
1103  // We'll make our best guess about its type here, since we have no
1104  // way to get back the original fifth template parameter.
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_);
1120  }
1121  }
1122  else {
1123  // This always works for a Tpetra::RowMatrix, even if it is not a
1124  // Tpetra::CrsMatrix. We just don't have offsets in this case.
1125  A.getLocalDiagCopy (*D_rowMap);
1126  }
1127  RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
1128 
1129  if (debug_) {
1130  // In debug mode, make sure that all diagonal entries are
1131  // positive, on all processes. Note that *out_ only prints on
1132  // Process 0 of the matrix's communicator.
1133  bool foundNonpositiveValue = false;
1134  {
1135  auto D_lcl = D_rangeMap->getLocalViewHost (Tpetra::Access::ReadOnly);
1136  auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1137 
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;
1142 
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;
1147  break;
1148  }
1149  }
1150  }
1151 
1152  using Teuchos::outArg;
1153  using Teuchos::REDUCE_MIN;
1154  using Teuchos::reduceAll;
1155 
1156  const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1157  int gblSuccess = lclSuccess; // to be overwritten
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));
1161  }
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;
1165  }
1166  else {
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;
1170  }
1171  }
1172 
1173  // Invert the diagonal entries, replacing entries less (in
1174  // magnitude) than the user-specified value with that value.
1175  reciprocal_threshold (*D_rangeMap, minDiagVal_);
1176  return Teuchos::rcp_const_cast<const V> (D_rangeMap);
1177 }
1178 
1179 
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
1184 {
1185  using Teuchos::RCP;
1186  using Teuchos::rcp;
1187  typedef Tpetra::Export<typename MV::local_ordinal_type,
1188  typename MV::global_ordinal_type,
1189  typename MV::node_type> export_type;
1190  // This throws logic_error instead of runtime_error, because the
1191  // methods that call makeRangeMapVector should all have checked
1192  // whether A_ is null before calling this method.
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.");
1202 
1203  RCP<const map_type> sourceMap = D->getMap ();
1204  RCP<const map_type> rangeMap = A_->getRangeMap ();
1205  RCP<const map_type> rowMap = A_->getRowMap ();
1206 
1207  if (rangeMap->isSameAs (*sourceMap)) {
1208  // The given vector's Map is the same as the matrix's range Map.
1209  // That means we don't need to Export. This is the fast path.
1210  return D;
1211  }
1212  else { // We need to Export.
1213  RCP<const export_type> exporter;
1214  // Making an Export object from scratch is expensive enough that
1215  // it's worth the O(1) global reductions to call isSameAs(), to
1216  // see if we can avoid that cost.
1217  if (sourceMap->isSameAs (*rowMap)) {
1218  // We can reuse the matrix's Export object, if there is one.
1219  exporter = A_->getGraph ()->getExporter ();
1220  }
1221  else { // We have to make a new Export object.
1222  exporter = rcp (new export_type (sourceMap, rangeMap));
1223  }
1224 
1225  if (exporter.is_null ()) {
1226  return D; // Row Map and range Map are the same; no need to Export.
1227  }
1228  else { // Row Map and range Map are _not_ the same; must Export.
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);
1232  }
1233  }
1234 }
1235 
1236 
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
1241 {
1242  using Teuchos::rcp_const_cast;
1243  return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1244 }
1245 
1246 
1247 template<class ScalarType, class MV>
1248 void
1249 Chebyshev<ScalarType, MV>::
1250 textbookApplyImpl (const op_type& A,
1251  const MV& B,
1252  MV& X,
1253  const int numIters,
1254  const ST lambdaMax,
1255  const ST lambdaMin,
1256  const ST eigRatio,
1257  const V& D_inv) const
1258 {
1259  (void) lambdaMin; // Forestall compiler warning.
1260  const ST myLambdaMin = lambdaMax / eigRatio;
1261 
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; // Ifpack2 calls this theta
1266  const ST c = (lambdaMax - myLambdaMin) / two; // Ifpack2 calls this 1/delta
1267 
1268  if (zeroStartingSolution_ && numIters > 0) {
1269  // If zero iterations, then input X is output X.
1270  X.putScalar (zero);
1271  }
1272  MV R (B.getMap (), B.getNumVectors (), false);
1273  MV P (B.getMap (), B.getNumVectors (), false);
1274  MV Z (B.getMap (), B.getNumVectors (), false);
1275  ST alpha, beta;
1276  for (int i = 0; i < numIters; ++i) {
1277  computeResidual (R, B, A, X); // R = B - A*X
1278  solve (Z, D_inv, R); // z = D_inv * R, that is, D \ R.
1279  if (i == 0) {
1280  P = Z;
1281  alpha = two / d;
1282  } else {
1283  //beta = (c * alpha / two)^2;
1284  //const ST sqrtBeta = c * alpha / two;
1285  //beta = sqrtBeta * sqrtBeta;
1286  beta = alpha * (c/two) * (c/two);
1287  alpha = one / (d - beta);
1288  P.update (one, Z, beta); // P = Z + beta*P
1289  }
1290  X.update (alpha, P, one); // X = X + alpha*P
1291  // If we compute the residual here, we could either do R = B -
1292  // A*X, or R = R - alpha*A*P. Since we choose the former, we
1293  // can move the computeResidual call to the top of the loop.
1294  }
1295 }
1296 
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 ());
1303 }
1304 
1305 template<class ScalarType, class MV>
1306 void
1307 Chebyshev<ScalarType, MV>::
1308 ifpackApplyImpl (const op_type& A,
1309  const MV& B,
1310  MV& X,
1311  const int numIters,
1312  const ST lambdaMax,
1313  const ST lambdaMin,
1314  const ST eigRatio,
1315  const V& D_inv)
1316 {
1317  using std::endl;
1318 #ifdef HAVE_IFPACK2_DEBUG
1319  const bool debug = debug_;
1320 #else
1321  const bool debug = false;
1322 #endif
1323 
1324  if (debug) {
1325  *out_ << " \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1326  *out_ << " \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1327  }
1328 
1329  if (numIters <= 0) {
1330  return;
1331  }
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);
1335 
1336  // Quick solve when the matrix A is the identity.
1337  if (lambdaMin == one && lambdaMax == lambdaMin) {
1338  solve (X, D_inv, B);
1339  return;
1340  }
1341 
1342  // Initialize coefficients
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;
1348 
1349  if (debug) {
1350  *out_ << " alpha = " << alpha << endl
1351  << " beta = " << beta << endl
1352  << " delta = " << delta << endl
1353  << " theta = " << theta << endl
1354  << " s1 = " << s1 << endl;
1355  }
1356 
1357  // Fetch cached temporary (multi)vector.
1358  Teuchos::RCP<MV> W_ptr = makeTempMultiVector (B);
1359  MV& W = *W_ptr;
1360 
1361  if (debug) {
1362  *out_ << " Iteration " << 1 << ":" << endl
1363  << " - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1364  }
1365 
1366  // Special case for the first iteration.
1367  if (! zeroStartingSolution_) {
1368  // mfh 22 May 2019: Tests don't actually exercise this path.
1369 
1370  if (ck_.is_null ()) {
1371  Teuchos::RCP<const op_type> A_op = A_;
1372  ck_ = Teuchos::rcp (new ChebyshevKernel<op_type> (A_op));
1373  }
1374  // W := (1/theta)*D_inv*(B-A*X) and X := X + W.
1375  // X := X + W
1376  ck_->compute (W, one/theta, const_cast<V&> (D_inv),
1377  const_cast<MV&> (B), X, zero);
1378  }
1379  else {
1380  // W := (1/theta)*D_inv*B and X := 0 + W.
1381  firstIterationWithZeroStartingSolution (W, one/theta, D_inv, B, X);
1382  }
1383 
1384  if (debug) {
1385  *out_ << " - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1386  << " - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1387  }
1388 
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));
1392  }
1393 
1394  // The rest of the iterations.
1395  ST rhok = one / s1;
1396  ST rhokp1, dtemp1, dtemp2;
1397  for (int deg = 1; deg < numIters; ++deg) {
1398  if (debug) {
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;
1404  }
1405 
1406  rhokp1 = one / (two * s1 - rhok);
1407  dtemp1 = rhokp1 * rhok;
1408  dtemp2 = two * rhokp1 * delta;
1409  rhok = rhokp1;
1410 
1411  if (debug) {
1412  *out_ << " - dtemp1 = " << dtemp1 << endl
1413  << " - dtemp2 = " << dtemp2 << endl;
1414  }
1415 
1416  // W := dtemp2*D_inv*(B - A*X) + dtemp1*W.
1417  // X := X + W
1418  ck_->compute (W, dtemp2, const_cast<V&> (D_inv),
1419  const_cast<MV&> (B), (X), dtemp1);
1420 
1421  if (debug) {
1422  *out_ << " - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl
1423  << " - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1424  }
1425  }
1426 }
1427 
1428 template<class ScalarType, class MV>
1429 typename Chebyshev<ScalarType, MV>::ST
1430 Chebyshev<ScalarType, MV>::
1431 powerMethodWithInitGuess (const op_type& A,
1432  const V& D_inv,
1433  const int numIters,
1434  V& x)
1435 {
1436  using std::endl;
1437  if (debug_) {
1438  *out_ << " powerMethodWithInitGuess:" << endl;
1439  }
1440 
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;
1445  ST norm;
1446 
1447  Teuchos::RCP<V> y;
1448  if (eigVector2_.is_null()) {
1449  y = rcp(new V(A.getRangeMap ()));
1450  if (eigKeepVectors_)
1451  eigVector2_ = y;
1452  } else
1453  y = eigVector2_;
1454  norm = x.norm2 ();
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.");
1462 
1463  if (debug_) {
1464  *out_ << " Original norm1(x): " << x.norm1 ()
1465  << ", norm2(x): " << norm << endl;
1466  }
1467 
1468  x.scale (one / norm);
1469 
1470  if (debug_) {
1471  *out_ << " norm1(x.scale(one/norm)): " << x.norm1 () << endl;
1472  }
1473 
1474  for (int iter = 0; iter < numIters-1; ++iter) {
1475  if (debug_) {
1476  *out_ << " Iteration " << iter << endl;
1477  }
1478  A.apply (x, *y);
1479  solve (x, D_inv, *y);
1480 
1481  if (((iter+1) % eigNormalizationFreq_ == 0) && (iter < numIters-2)) {
1482  norm = x.norm2 ();
1483  if (norm == zero) { // Return something reasonable.
1484  if (debug_) {
1485  *out_ << " norm is zero; returning zero" << endl;
1486  *out_ << " Power method terminated after "<< iter << " iterations." << endl;
1487  }
1488  return zero;
1489  } else {
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)) {
1493  if (debug_) {
1494  *out_ << " lambdaMax: " << lambdaMax << endl;
1495  *out_ << " Power method terminated after "<< iter << " iterations." << endl;
1496  }
1497  return lambdaMax;
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;
1502  }
1503  }
1504  x.scale (one / norm);
1505  }
1506  }
1507  if (debug_) {
1508  *out_ << " lambdaMax: " << lambdaMax << endl;
1509  }
1510 
1511  norm = x.norm2 ();
1512  if (norm == zero) { // Return something reasonable.
1513  if (debug_) {
1514  *out_ << " norm is zero; returning zero" << endl;
1515  *out_ << " Power method terminated after "<< numIters << " iterations." << endl;
1516  }
1517  return zero;
1518  }
1519  x.scale (one / norm);
1520  A.apply (x, *y);
1521  solve (*y, D_inv, *y);
1522  lambdaMax = y->dot (x);
1523  if (debug_) {
1524  *out_ << " lambdaMax: " << lambdaMax << endl;
1525  *out_ << " Power method terminated after "<< numIters << " iterations." << endl;
1526  }
1527 
1528  return lambdaMax;
1529 }
1530 
1531 template<class ScalarType, class MV>
1532 void
1533 Chebyshev<ScalarType, MV>::
1534 computeInitialGuessForPowerMethod (V& x, const bool nonnegativeRealParts) const
1535 {
1536  typedef typename MV::device_type::execution_space dev_execution_space;
1537  typedef typename MV::local_ordinal_type LO;
1538 
1539  x.randomize ();
1540 
1541  if (nonnegativeRealParts) {
1542  auto x_lcl = x.getLocalViewDevice (Tpetra::Access::ReadWrite);
1543  auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
1544 
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);
1549  }
1550 }
1551 
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)
1556 {
1557  using std::endl;
1558  if (debug_) {
1559  *out_ << "powerMethod:" << endl;
1560  }
1561 
1562  const ST zero = static_cast<ST> (0.0);
1563  Teuchos::RCP<V> x;
1564  if (eigVector_.is_null()) {
1565  x = rcp(new V(A.getDomainMap ()));
1566  if (eigKeepVectors_)
1567  eigVector_ = x;
1568  // For the first pass, just let the pseudorandom number generator
1569  // fill x with whatever values it wants; don't try to make its
1570  // entries nonnegative.
1571  computeInitialGuessForPowerMethod (*x, false);
1572  } else
1573  x = eigVector_;
1574 
1575  ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, *x);
1576 
1577  // mfh 07 Jan 2015: Taking the real part here is only a concession
1578  // to the compiler, so that this class can build with ScalarType =
1579  // std::complex<T>. Our Chebyshev implementation only works with
1580  // real, symmetric positive definite matrices. The right thing to
1581  // do would be what Belos does, which is provide a partial
1582  // specialization for ScalarType = std::complex<T> with a stub
1583  // implementation (that builds, but whose constructor throws).
1584  if (STS::real (lambdaMax) < STS::real (zero)) {
1585  if (debug_) {
1586  *out_ << "real(lambdaMax) = " << STS::real (lambdaMax) << " < 0; "
1587  "try again with a different random initial guess" << endl;
1588  }
1589  // Max eigenvalue estimate was negative. Perhaps we got unlucky
1590  // with the random initial guess. Try again with a different (but
1591  // still random) initial guess. Only try again once, so that the
1592  // run time is bounded.
1593 
1594  // For the second pass, make all the entries of the initial guess
1595  // vector have nonnegative real parts.
1596  computeInitialGuessForPowerMethod (*x, true);
1597  lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, *x);
1598  }
1599  return lambdaMax;
1600 }
1601 
1602 
1603 template<class ScalarType, class MV>
1604 typename Chebyshev<ScalarType, MV>::ST
1605 Chebyshev<ScalarType, MV>::
1606 cgMethodWithInitGuess (const op_type& A,
1607  const V& D_inv,
1608  const int numIters,
1609  V& r)
1610 {
1611  using std::endl;
1612  using STS = Teuchos::ScalarTraits<ST>;
1613  using MagnitudeType = typename STS::magnitudeType;
1614  if (debug_) {
1615  *out_ << " cgMethodWithInitGuess:" << endl;
1616  }
1617 
1618  const ST one = STS::one();
1619  ST beta, betaOld = one, pAp, pApOld = one, alpha, rHz, rHzOld, rHzOld2 = one, lambdaMax;
1620  // ST lambdaMin;
1621  Teuchos::ArrayRCP<MagnitudeType> diag, offdiag;
1622  Teuchos::RCP<V> p, z, Ap;
1623  diag.resize(numIters);
1624  offdiag.resize(numIters-1);
1625 
1626  p = rcp(new V(A.getRangeMap ()));
1627  z = rcp(new V(A.getRangeMap ()));
1628  Ap = rcp(new V(A.getRangeMap ()));
1629 
1630  // Tpetra::Details::residual (A, x, *b, *r);
1631  solve (*p, D_inv, r);
1632  rHz = r.dot (*p);
1633 
1634  for (int iter = 0; iter < numIters; ++iter) {
1635  if (debug_) {
1636  *out_ << " Iteration " << iter << endl;
1637  }
1638  A.apply (*p, *Ap);
1639  pAp = p->dot (*Ap);
1640  alpha = rHz/pAp;
1641  r.update (-alpha, *Ap, one);
1642  rHzOld = rHz;
1643  solve (*z, D_inv, r);
1644  rHz = r.dot (*z);
1645  beta = rHz / rHzOld;
1646  p->update (one, *z, beta);
1647  if (iter>0) {
1648  diag[iter] = STS::real((betaOld*betaOld * pApOld + pAp) / rHzOld);
1649  offdiag[iter-1] = -STS::real((betaOld * pApOld / (sqrt(rHzOld * rHzOld2))));
1650  if (debug_) {
1651  *out_ << " diag[" << iter << "] = " << diag[iter] << endl;
1652  *out_ << " offdiag["<< iter-1 << "] = " << offdiag[iter-1] << endl;
1653  }
1654  } else {
1655  diag[iter] = STS::real(pAp/rHzOld);
1656  if (debug_) {
1657  *out_ << " diag[" << iter << "] = " << diag[iter] << endl;
1658  }
1659  }
1660  rHzOld2 = rHzOld;
1661  betaOld = beta;
1662  pApOld = pAp;
1663  }
1664 
1665  lambdaMax = LapackHelper<ST>::tri_diag_spectral_radius(diag, offdiag);
1666 
1667  return lambdaMax;
1668 }
1669 
1670 
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)
1675 {
1676  using std::endl;
1677  if (debug_) {
1678  *out_ << "cgMethod:" << endl;
1679  }
1680 
1681  Teuchos::RCP<V> r;
1682  if (eigVector_.is_null()) {
1683  r = rcp(new V(A.getDomainMap ()));
1684  if (eigKeepVectors_)
1685  eigVector_ = r;
1686  // For the first pass, just let the pseudorandom number generator
1687  // fill x with whatever values it wants; don't try to make its
1688  // entries nonnegative.
1689  computeInitialGuessForPowerMethod (*r, false);
1690  } else
1691  r = eigVector_;
1692 
1693  ST lambdaMax = cgMethodWithInitGuess (A, D_inv, numIters, *r);
1694 
1695  return lambdaMax;
1696 }
1697 
1698 template<class ScalarType, class MV>
1699 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::row_matrix_type>
1701  return A_;
1702 }
1703 
1704 template<class ScalarType, class MV>
1705 bool
1708  // Technically, this is true, because the matrix must be symmetric.
1709  return true;
1710 }
1711 
1712 template<class ScalarType, class MV>
1713 Teuchos::RCP<MV>
1715 makeTempMultiVector (const MV& B)
1716 {
1717  // ETP 02/08/17: We must check not only if the temporary vectors are
1718  // null, but also if the number of columns match, since some multi-RHS
1719  // solvers (e.g., Belos) may call apply() with different numbers of columns.
1720 
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));
1724  }
1725  return W_;
1726 }
1727 
1728 template<class ScalarType, class MV>
1729 std::string
1730 Chebyshev<ScalarType, MV>::
1731 description () const {
1732  std::ostringstream oss;
1733  // YAML requires quoting the key in this case, to distinguish
1734  // key's colons from the colon that separates key from value.
1735  oss << "\"Ifpack2::Details::Chebyshev\":"
1736  << "{"
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";
1744  oss << "}";
1745  return oss.str();
1746 }
1747 
1748 template<class ScalarType, class MV>
1749 void
1751 describe (Teuchos::FancyOStream& out,
1752  const Teuchos::EVerbosityLevel verbLevel) const
1753 {
1754  using Teuchos::TypeNameTraits;
1755  using std::endl;
1756 
1757  const Teuchos::EVerbosityLevel vl =
1758  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1759  if (vl == Teuchos::VERB_NONE) {
1760  return; // print NOTHING
1761  }
1762 
1763  // By convention, describe() starts with a tab.
1764  //
1765  // This does affect all processes on which it's valid to print to
1766  // 'out'. However, it does not actually print spaces to 'out'
1767  // unless operator<< gets called, so it's safe to use on all
1768  // processes.
1769  Teuchos::OSTab tab0 (out);
1770 
1771  // We only print on Process 0 of the matrix's communicator. If
1772  // the matrix isn't set, we don't have a communicator, so we have
1773  // to assume that every process can print.
1774  int myRank = -1;
1775  if (A_.is_null () || A_->getComm ().is_null ()) {
1776  myRank = 0;
1777  }
1778  else {
1779  myRank = A_->getComm ()->getRank ();
1780  }
1781  if (myRank == 0) {
1782  // YAML requires quoting the key in this case, to distinguish
1783  // key's colons from the colon that separates key from value.
1784  out << "\"Ifpack2::Details::Chebyshev\":" << endl;
1785  }
1786  Teuchos::OSTab tab1 (out);
1787 
1788  if (vl == Teuchos::VERB_LOW) {
1789  if (myRank == 0) {
1790  out << "degree: " << numIters_ << endl
1791  << "lambdaMax: " << lambdaMaxForApply_ << endl
1792  << "alpha: " << eigRatioForApply_ << endl
1793  << "lambdaMin: " << lambdaMinForApply_ << endl
1794  << "boost factor: " << boostFactor_ << endl;
1795  }
1796  return;
1797  }
1798 
1799  // vl > Teuchos::VERB_LOW
1800 
1801  if (myRank == 0) {
1802  out << "Template parameters:" << endl;
1803  {
1804  Teuchos::OSTab tab2 (out);
1805  out << "ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1806  << "MV: " << TypeNameTraits<MV>::name () << endl;
1807  }
1808 
1809  // "Computed parameters" literally means "parameters whose
1810  // values were computed by compute()."
1811  if (myRank == 0) {
1812  out << "Computed parameters:" << endl;
1813  }
1814  }
1815 
1816  // Print computed parameters
1817  {
1818  Teuchos::OSTab tab2 (out);
1819  // Users might want to see the values in the computed inverse
1820  // diagonal, so we print them out at the highest verbosity.
1821  if (myRank == 0) {
1822  out << "D_: ";
1823  }
1824  if (D_.is_null ()) {
1825  if (myRank == 0) {
1826  out << "unset" << endl;
1827  }
1828  }
1829  else if (vl <= Teuchos::VERB_HIGH) {
1830  if (myRank == 0) {
1831  out << "set" << endl;
1832  }
1833  }
1834  else { // D_ not null and vl > Teuchos::VERB_HIGH
1835  if (myRank == 0) {
1836  out << endl;
1837  }
1838  // By convention, describe() first indents, then prints.
1839  // We can rely on other describe() implementations to do that.
1840  D_->describe (out, vl);
1841  }
1842  if (myRank == 0) {
1843  // W_ is scratch space; its values are irrelevant.
1844  // All that matters is whether or not they have been set.
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;
1851  }
1852  } // print computed parameters
1853 
1854  if (myRank == 0) {
1855  out << "User parameters:" << endl;
1856  }
1857 
1858  // Print user parameters
1859  {
1860  Teuchos::OSTab tab2 (out);
1861  out << "userInvDiag_: ";
1862  if (userInvDiag_.is_null ()) {
1863  out << "unset" << endl;
1864  }
1865  else if (vl <= Teuchos::VERB_HIGH) {
1866  out << "set" << endl;
1867  }
1868  else { // userInvDiag_ not null and vl > Teuchos::VERB_HIGH
1869  if (myRank == 0) {
1870  out << endl;
1871  }
1872  userInvDiag_->describe (out, vl);
1873  }
1874  if (myRank == 0) {
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;
1886  }
1887  } // print user parameters
1888 }
1889 
1890 } // namespace Details
1891 } // namespace Ifpack2
1892 
1893 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \
1894  template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >;
1895 
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