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 "Kokkos_ArithTraits.hpp"
56 #include "Teuchos_FancyOStream.hpp"
57 #include "Teuchos_oblackholestream.hpp"
58 #include <cmath>
59 #include <iostream>
60 
61 namespace Ifpack2 {
62 namespace Details {
63 
64 namespace { // (anonymous)
65 
66 // We use this text a lot in error messages.
67 const char computeBeforeApplyReminder[] =
68  "This means one of the following:\n"
69  " - you have not yet called compute() on this instance, or \n"
70  " - you didn't call compute() after calling setParameters().\n\n"
71  "After creating an Ifpack2::Chebyshev instance,\n"
72  "you must _always_ call compute() at least once before calling apply().\n"
73  "After calling compute() once, you do not need to call it again,\n"
74  "unless the matrix has changed or you have changed parameters\n"
75  "(by calling setParameters()).";
76 
77 } // namespace (anonymous)
78 
79 // ReciprocalThreshold stuff below needs to be in a namspace visible outside
80 // of this file
81 template<class XV, class SizeType = typename XV::size_type>
82 struct V_ReciprocalThresholdSelfFunctor
83 {
84  typedef typename XV::execution_space execution_space;
85  typedef typename XV::non_const_value_type value_type;
86  typedef SizeType size_type;
87  typedef Kokkos::Details::ArithTraits<value_type> KAT;
88  typedef typename KAT::mag_type mag_type;
89 
90  XV X_;
91  const value_type minVal_;
92  const mag_type minValMag_;
93 
94  V_ReciprocalThresholdSelfFunctor (const XV& X,
95  const value_type& min_val) :
96  X_ (X),
97  minVal_ (min_val),
98  minValMag_ (KAT::abs (min_val))
99  {}
100 
101  KOKKOS_INLINE_FUNCTION
102  void operator() (const size_type& i) const
103  {
104  const mag_type X_i_abs = KAT::abs (X_[i]);
105 
106  if (X_i_abs < minValMag_) {
107  X_[i] = minVal_;
108  }
109  else {
110  X_[i] = KAT::one () / X_[i];
111  }
112  }
113 };
114 
115 template<class XV, class SizeType = typename XV::size_type>
116 struct LocalReciprocalThreshold {
117  static void
118  compute (const XV& X,
119  const typename XV::non_const_value_type& minVal)
120  {
121  typedef typename XV::execution_space execution_space;
122  Kokkos::RangePolicy<execution_space, SizeType> policy (0, X.extent (0));
123  V_ReciprocalThresholdSelfFunctor<XV, SizeType> op (X, minVal);
124  Kokkos::parallel_for (policy, op);
125  }
126 };
127 
128 template <class TpetraVectorType,
129  const bool classic = TpetraVectorType::node_type::classic>
130 struct GlobalReciprocalThreshold {};
131 
132 template <class TpetraVectorType>
133 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
134  static void
135  compute (TpetraVectorType& V,
136  const typename TpetraVectorType::scalar_type& min_val)
137  {
138  typedef typename TpetraVectorType::scalar_type scalar_type;
139  typedef typename TpetraVectorType::mag_type mag_type;
140  typedef Kokkos::Details::ArithTraits<scalar_type> STS;
141 
142  const scalar_type ONE = STS::one ();
143  const mag_type min_val_abs = STS::abs (min_val);
144 
145  Teuchos::ArrayRCP<scalar_type> V_0 = V.getDataNonConst (0);
146  const size_t lclNumRows = V.getLocalLength ();
147 
148  for (size_t i = 0; i < lclNumRows; ++i) {
149  const scalar_type V_0i = V_0[i];
150  if (STS::abs (V_0i) < min_val_abs) {
151  V_0[i] = min_val;
152  } else {
153  V_0[i] = ONE / V_0i;
154  }
155  }
156  }
157 };
158 
159 template <class TpetraVectorType>
160 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
161  static void
162  compute (TpetraVectorType& X,
163  const typename TpetraVectorType::scalar_type& minVal)
164  {
165  typedef typename TpetraVectorType::impl_scalar_type value_type;
166  typedef typename TpetraVectorType::device_type::memory_space memory_space;
167 
168  X.template sync<memory_space> ();
169  X.template modify<memory_space> ();
170 
171  const value_type minValS = static_cast<value_type> (minVal);
172  auto X_0 = Kokkos::subview (X.template getLocalView<memory_space> (),
173  Kokkos::ALL (), 0);
174  LocalReciprocalThreshold<decltype (X_0) >::compute (X_0, minValS);
175  }
176 };
177 
178 // Utility function for inverting diagonal with threshold.
179 template <typename S, typename L, typename G, typename N>
180 void
181 reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V, const S& minVal)
182 {
183  GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
184 }
185 
186 namespace { // (anonymous)
187 
188 // Functor for making sure the real parts of all entries of a vector
189 // are nonnegative. We use this in computeInitialGuessForPowerMethod
190 // below.
191 template<class OneDViewType,
192  class LocalOrdinal = typename OneDViewType::size_type>
193 class PositivizeVector {
194  static_assert (Kokkos::Impl::is_view<OneDViewType>::value,
195  "OneDViewType must be a 1-D Kokkos::View.");
196  static_assert (static_cast<int> (OneDViewType::rank) == 1,
197  "This functor only works with a 1-D View.");
198  static_assert (std::is_integral<LocalOrdinal>::value,
199  "The second template parameter, LocalOrdinal, "
200  "must be an integer type.");
201 public:
202  PositivizeVector (const OneDViewType& x) : x_ (x) {}
203 
204  KOKKOS_INLINE_FUNCTION void
205  operator () (const LocalOrdinal& i) const
206  {
207  typedef typename OneDViewType::non_const_value_type IST;
208  typedef Kokkos::Details::ArithTraits<IST> STS;
209  typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
210 
211  if (STS::real (x_(i)) < STM::zero ()) {
212  x_(i) = -x_(i);
213  }
214  }
215 
216 private:
217  OneDViewType x_;
218 };
219 
220 } // namespace (anonymous)
221 
222 template<class ScalarType, class MV>
223 void Chebyshev<ScalarType, MV>::checkInputMatrix () const
224 {
226  ! A_.is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
227  std::invalid_argument,
228  "Ifpack2::Chebyshev: The input matrix A must be square. "
229  "A has " << A_->getGlobalNumRows () << " rows and "
230  << A_->getGlobalNumCols () << " columns.");
231 
232  // In debug mode, test that the domain and range Maps of the matrix
233  // are the same.
234  if (debug_ && ! A_.is_null ()) {
235  Teuchos::RCP<const map_type> domainMap = A_->getDomainMap ();
236  Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap ();
237 
238  // isSameAs is a collective, but if the two pointers are the same,
239  // isSameAs will assume that they are the same on all processes, and
240  // return true without an all-reduce.
242  ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
243  "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be "
244  "the same (in the sense of isSameAs())." << std::endl << "We only check "
245  "for this in debug mode.");
246  }
247 }
248 
249 
250 template<class ScalarType, class MV>
251 void
252 Chebyshev<ScalarType, MV>::
253 checkConstructorInput () const
254 {
255  // mfh 12 Aug 2016: The if statement avoids an "unreachable
256  // statement" warning for the checkInputMatrix() call, when
257  // STS::isComplex is false.
258  if (STS::isComplex) {
260  (true, std::logic_error, "Ifpack2::Chebyshev: This class' implementation "
261  "of Chebyshev iteration only works for real-valued, symmetric positive "
262  "definite matrices. However, you instantiated this class for ScalarType"
263  " = " << Teuchos::TypeNameTraits<ScalarType>::name () << ", which is a "
264  "complex-valued type. While this may be algorithmically correct if all "
265  "of the complex numbers in the matrix have zero imaginary part, we "
266  "forbid using complex ScalarType altogether in order to remind you of "
267  "the limitations of our implementation (and of the algorithm itself).");
268  }
269  else {
270  checkInputMatrix ();
271  }
272 }
273 
274 template<class ScalarType, class MV>
277  A_ (A),
278  savedDiagOffsets_ (false),
279  computedLambdaMax_ (STS::nan ()),
280  computedLambdaMin_ (STS::nan ()),
281  lambdaMaxForApply_ (STS::nan ()),
282  lambdaMinForApply_ (STS::nan ()),
283  eigRatioForApply_ (STS::nan ()),
284  userLambdaMax_ (STS::nan ()),
285  userLambdaMin_ (STS::nan ()),
286  userEigRatio_ (Teuchos::as<ST> (30)),
287  minDiagVal_ (STS::eps ()),
288  numIters_ (1),
289  eigMaxIters_ (10),
290  zeroStartingSolution_ (true),
291  assumeMatrixUnchanged_ (false),
292  textbookAlgorithm_ (false),
293  computeMaxResNorm_ (false),
294  debug_ (false)
295 {
296  checkConstructorInput ();
297 }
298 
299 template<class ScalarType, class MV>
302  A_ (A),
303  savedDiagOffsets_ (false),
304  computedLambdaMax_ (STS::nan ()),
305  computedLambdaMin_ (STS::nan ()),
306  lambdaMaxForApply_ (STS::nan ()),
307  boostFactor_ (static_cast<MT> (1.1)),
308  lambdaMinForApply_ (STS::nan ()),
309  eigRatioForApply_ (STS::nan ()),
310  userLambdaMax_ (STS::nan ()),
311  userLambdaMin_ (STS::nan ()),
312  userEigRatio_ (Teuchos::as<ST> (30)),
313  minDiagVal_ (STS::eps ()),
314  numIters_ (1),
315  eigMaxIters_ (10),
316  zeroStartingSolution_ (true),
317  assumeMatrixUnchanged_ (false),
318  textbookAlgorithm_ (false),
319  computeMaxResNorm_ (false),
320  debug_ (false)
321 {
322  checkConstructorInput ();
323  setParameters (params);
324 }
325 
326 template<class ScalarType, class MV>
327 void
330 {
331  using Teuchos::RCP;
332  using Teuchos::rcp;
333  using Teuchos::rcp_const_cast;
334 
335  // Note to developers: The logic for this method is complicated,
336  // because we want to accept Ifpack and ML parameters whenever
337  // possible, but we don't want to add their default values to the
338  // user's ParameterList. That's why we do all the isParameter()
339  // checks, instead of using the two-argument version of get()
340  // everywhere. The min and max eigenvalue parameters are also a
341  // special case, because we decide whether or not to do eigenvalue
342  // analysis based on whether the user supplied the max eigenvalue.
343 
344  // Default values of all the parameters.
345  const ST defaultLambdaMax = STS::nan ();
346  const ST defaultLambdaMin = STS::nan ();
347  // 30 is Ifpack::Chebyshev's default. ML has a different default
348  // eigRatio for smoothers and the coarse-grid solve (if using
349  // Chebyshev for that). The former uses 20; the latter uses 30.
350  // We're testing smoothers here, so use 20. (However, if you give
351  // ML an Epetra matrix, it will use Ifpack for Chebyshev, in which
352  // case it would defer to Ifpack's default settings.)
353  const ST defaultEigRatio = Teuchos::as<ST> (30);
354  const MT defaultBoostFactor = static_cast<MT> (1.1);
355  const ST defaultMinDiagVal = STS::eps ();
356  const int defaultNumIters = 1;
357  const int defaultEigMaxIters = 10;
358  const bool defaultZeroStartingSolution = true; // Ifpack::Chebyshev default
359  const bool defaultAssumeMatrixUnchanged = false;
360  const bool defaultTextbookAlgorithm = false;
361  const bool defaultComputeMaxResNorm = false;
362  const bool defaultDebug = false;
363 
364  // We'll set the instance data transactionally, after all reads
365  // from the ParameterList. That way, if any of the ParameterList
366  // reads fail (e.g., due to the wrong parameter type), we will not
367  // have left the instance data in a half-changed state.
368  RCP<const V> userInvDiagCopy; // if nonnull: deep copy of user's Vector
369  ST lambdaMax = defaultLambdaMax;
370  ST lambdaMin = defaultLambdaMin;
371  ST eigRatio = defaultEigRatio;
372  MT boostFactor = defaultBoostFactor;
373  ST minDiagVal = defaultMinDiagVal;
374  int numIters = defaultNumIters;
375  int eigMaxIters = defaultEigMaxIters;
376  bool zeroStartingSolution = defaultZeroStartingSolution;
377  bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
378  bool textbookAlgorithm = defaultTextbookAlgorithm;
379  bool computeMaxResNorm = defaultComputeMaxResNorm;
380  bool debug = defaultDebug;
381 
382  // Fetch the parameters from the ParameterList. Defer all
383  // externally visible side effects until we have finished all
384  // ParameterList interaction. This makes the method satisfy the
385  // strong exception guarantee.
386 
387  if (plist.isParameter ("debug")) {
388  try { // a bool
389  debug = plist.get<bool> ("debug");
390  }
392 
393  try { // an int instead of a bool
394  int debugInt = plist.get<bool> ("debug");
395  debug = debugInt != 0;
396  }
398  }
399 
400  // Get the user-supplied inverse diagonal.
401  //
402  // Check for a raw pointer (const V* or V*), for Ifpack
403  // compatibility, as well as for RCP<const V>, RCP<V>, const V, or
404  // V. We'll make a deep copy of the vector at the end of this
405  // method anyway, so its const-ness doesn't matter. We handle the
406  // latter two cases ("const V" or "V") specially (copy them into
407  // userInvDiagCopy first, which is otherwise null at the end of the
408  // long if-then chain) to avoid an extra copy.
409  if (plist.isParameter ("chebyshev: operator inv diagonal")) {
410  // Pointer to the user's Vector, if provided.
411  RCP<const V> userInvDiag;
412 
413  try { // Could the type be const V*?
414  const V* rawUserInvDiag =
415  plist.get<const V*> ("chebyshev: operator inv diagonal");
416  // Nonowning reference (we'll make a deep copy below)
417  userInvDiag = rcp (rawUserInvDiag, false);
419  }
420  if (userInvDiag.is_null ()) {
421  try { // Could the type be V*?
422  V* rawUserInvDiag = plist.get<V*> ("chebyshev: operator inv diagonal");
423  // Nonowning reference (we'll make a deep copy below)
424  userInvDiag = rcp (const_cast<const V*> (rawUserInvDiag), false);
426  }
427  }
428  if (userInvDiag.is_null ()) {
429  try { // Could the type be RCP<const V>?
430  userInvDiag =
431  plist.get<RCP<const V> > ("chebyshev: operator inv diagonal");
433  }
434  }
435  if (userInvDiag.is_null ()) {
436  try { // Could the type be RCP<V>?
437  RCP<V> userInvDiagNonConst =
438  plist.get<RCP<V> > ("chebyshev: operator inv diagonal");
439  userInvDiag = rcp_const_cast<const V> (userInvDiagNonConst);
441  }
442  }
443  if (userInvDiag.is_null ()) {
444 #ifndef _MSC_VER
445  try { // Could the type be const V?
446  // ParameterList::get() returns by reference. Thus, we don't
447  // have to invoke Vector's copy constructor here. It's good
448  // practice not to make an RCP to this reference, even though
449  // it should be valid as long as the ParameterList that holds
450  // it is valid. Thus, we make our deep copy here, rather than
451  // waiting to do it below.
452  const V& userInvDiagRef =
453  plist.get<const V> ("chebyshev: operator inv diagonal");
454  userInvDiagCopy = rcp (new V (userInvDiagRef, Teuchos::Copy));
455  // Tell the if-chain below not to keep trying.
456  userInvDiag = userInvDiagCopy;
458  }
459 #else
461  true, std::runtime_error,
462  "Ifpack2::Chebyshev::setParameters: \"chebyshev: operator inv diagonal\" "
463  "plist.get<const V> does not compile around return held == other_held "
464  "in Teuchos::any in Visual Studio. Can't fix it now, so throwing "
465  "in case someone builds there.");
466 #endif
467  }
468  if (userInvDiag.is_null ()) {
469 #ifndef _MSC_VER
470  try { // Could the type be V?
471  // ParameterList::get() returns by reference. Thus, we don't
472  // have to invoke Vector's copy constructor here. It's good
473  // practice not to make an RCP to this reference, even though
474  // it should be valid as long as the ParameterList that holds
475  // it is valid. Thus, we make our deep copy here, rather than
476  // waiting to do it below.
477  V& userInvDiagNonConstRef =
478  plist.get<V> ("chebyshev: operator inv diagonal");
479  const V& userInvDiagRef = const_cast<const V&> (userInvDiagNonConstRef);
480  userInvDiagCopy = rcp (new V (userInvDiagRef, Teuchos::Copy));
481  // Tell the if-chain below not to keep trying.
482  userInvDiag = userInvDiagCopy;
484  }
485 #else
487  true, std::runtime_error,
488  "Ifpack2::Chebyshev::setParameters: \"chebyshev: operator inv diagonal\" "
489  "plist.get<V> does not compile around return held == other_held "
490  "in Teuchos::any in Visual Studio. Can't fix it now, so throwing "
491  "in case someone builds there.");
492 #endif
493  }
494 
495  // NOTE: If the user's parameter has some strange type that we
496  // didn't test above, userInvDiag might still be null. You may
497  // want to add an error test for this condition. Currently, we
498  // just say in this case that the user didn't give us a Vector.
499 
500  // If we have userInvDiag but don't have a deep copy yet, make a
501  // deep copy now.
502  if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
503  userInvDiagCopy = rcp (new V (*userInvDiag, Teuchos::Copy));
504  }
505 
506  // NOTE: userInvDiag, if provided, is a row Map version of the
507  // Vector. We don't necessarily have a range Map yet. compute()
508  // would be the proper place to compute the range Map version of
509  // userInvDiag.
510  }
511 
512  // Don't fill in defaults for the max or min eigenvalue, because
513  // this class uses the existence of those parameters to determine
514  // whether it should do eigenanalysis.
515  if (plist.isParameter ("chebyshev: max eigenvalue")) {
516  if (plist.isType<double>("chebyshev: max eigenvalue"))
517  lambdaMax = plist.get<double> ("chebyshev: max eigenvalue");
518  else
519  lambdaMax = plist.get<ST> ("chebyshev: max eigenvalue");
521  STS::isnaninf (lambdaMax), std::invalid_argument,
522  "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" "
523  "parameter is NaN or Inf. This parameter is optional, but if you "
524  "choose to supply it, it must have a finite value.");
525  }
526  if (plist.isParameter ("chebyshev: min eigenvalue")) {
527  if (plist.isType<double>("chebyshev: min eigenvalue"))
528  lambdaMin = plist.get<double> ("chebyshev: min eigenvalue");
529  else
530  lambdaMin = plist.get<ST> ("chebyshev: min eigenvalue");
532  STS::isnaninf (lambdaMin), std::invalid_argument,
533  "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" "
534  "parameter is NaN or Inf. This parameter is optional, but if you "
535  "choose to supply it, it must have a finite value.");
536  }
537 
538  // Only fill in Ifpack2's name for the default parameter, not ML's.
539  if (plist.isParameter ("smoother: Chebyshev alpha")) { // ML compatibility
540  if (plist.isType<double>("smoother: Chebyshev alpha"))
541  eigRatio = plist.get<double> ("smoother: Chebyshev alpha");
542  else
543  eigRatio = plist.get<ST> ("smoother: Chebyshev alpha");
544  }
545  // Ifpack2's name overrides ML's name.
546  eigRatio = plist.get ("chebyshev: ratio eigenvalue", eigRatio);
548  STS::isnaninf (eigRatio), std::invalid_argument,
549  "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" "
550  "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. "
551  "This parameter is optional, but if you choose to supply it, it must have "
552  "a finite value.");
553  // mfh 11 Feb 2013: This class is currently only correct for real
554  // Scalar types, but we still want it to build for complex Scalar
555  // type so that users of Ifpack2::Factory can build their
556  // executables for real or complex Scalar type. Thus, we take the
557  // real parts here, which are always less-than comparable.
559  STS::real (eigRatio) < STS::real (STS::one ()),
560  std::invalid_argument,
561  "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\""
562  "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, "
563  "but you supplied the value " << eigRatio << ".");
564 
565  // See Github Issue #234. This parameter may be either MT
566  // (preferred) or double. We check both.
567  {
568  const char paramName[] = "chebyshev: boost factor";
569 
570  if (plist.isParameter (paramName)) {
571  if (plist.isType<MT> (paramName)) { // MT preferred
572  boostFactor = plist.get<MT> (paramName);
573  }
574  else if (! std::is_same<double, MT>::value &&
575  plist.isType<double> (paramName)) {
576  const double dblBF = plist.get<double> (paramName);
577  boostFactor = static_cast<MT> (dblBF);
578  }
579  else {
581  (true, std::invalid_argument,
582  "Ifpack2::Chebyshev::setParameters: \"chebyshev: boost factor\""
583  "parameter must have type magnitude_type (MT) or double.");
584  }
585  }
586  else { // parameter not in the list
587  // mfh 12 Aug 2016: To preserve current behavior (that fills in
588  // any parameters not in the input ParameterList with their
589  // default values), we call set() here. I don't actually like
590  // this behavior; I prefer the Belos model, where the input
591  // ParameterList is a delta from current behavior. However, I
592  // don't want to break things.
593  plist.set (paramName, defaultBoostFactor);
594  }
596  (boostFactor < Teuchos::ScalarTraits<MT>::one (), std::invalid_argument,
597  "Ifpack2::Chebyshev::setParameters: \"" << paramName << "\" parameter "
598  "must be >= 1, but you supplied the value " << boostFactor << ".");
599  }
600 
601  // Same name in Ifpack2 and Ifpack.
602  minDiagVal = plist.get ("chebyshev: min diagonal value", minDiagVal);
604  STS::isnaninf (minDiagVal), std::invalid_argument,
605  "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" "
606  "parameter is NaN or Inf. This parameter is optional, but if you choose "
607  "to supply it, it must have a finite value.");
608 
609  // Only fill in Ifpack2's name, not ML's or Ifpack's.
610  if (plist.isParameter ("smoother: sweeps")) { // ML compatibility
611  numIters = plist.get<int> ("smoother: sweeps");
612  } // Ifpack's name overrides ML's name.
613  if (plist.isParameter ("relaxation: sweeps")) { // Ifpack compatibility
614  numIters = plist.get<int> ("relaxation: sweeps");
615  } // Ifpack2's name overrides Ifpack's name.
616  numIters = plist.get ("chebyshev: degree", numIters);
618  numIters < 0, std::invalid_argument,
619  "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also "
620  "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a "
621  "nonnegative integer. You gave a value of " << numIters << ".");
622 
623  // The last parameter name overrides the first.
624  if (plist.isParameter ("eigen-analysis: iterations")) { // ML compatibility
625  eigMaxIters = plist.get<int> ("eigen-analysis: iterations");
626  } // Ifpack2's name overrides ML's name.
627  eigMaxIters = plist.get ("chebyshev: eigenvalue max iterations", eigMaxIters);
629  eigMaxIters < 0, std::invalid_argument,
630  "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations"
631  "\" parameter (also called \"eigen-analysis: iterations\") must be a "
632  "nonnegative integer. You gave a value of " << eigMaxIters << ".");
633 
634  zeroStartingSolution = plist.get ("chebyshev: zero starting solution",
635  zeroStartingSolution);
636  assumeMatrixUnchanged = plist.get ("chebyshev: assume matrix does not change",
637  assumeMatrixUnchanged);
638 
639  // We don't want to fill these parameters in, because they shouldn't
640  // be visible to Ifpack2::Chebyshev users.
641  if (plist.isParameter ("chebyshev: textbook algorithm")) {
642  textbookAlgorithm = plist.get<bool> ("chebyshev: textbook algorithm");
643  }
644  if (plist.isParameter ("chebyshev: compute max residual norm")) {
645  computeMaxResNorm = plist.get<bool> ("chebyshev: compute max residual norm");
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.
652  plist.isParameter ("chebyshev: use block mode") &&
653  ! plist.get<bool> ("chebyshev: use block mode"),
654  std::invalid_argument,
655  "Ifpack2::Chebyshev requires that if you supply the Ifpack parameter "
656  "\"chebyshev: use block mode\", it must be set to false. Ifpack2's "
657  "Chebyshev does not implement Ifpack's block mode.");
659  plist.isParameter ("chebyshev: solve normal equations") &&
660  ! plist.get<bool> ("chebyshev: solve normal equations"),
661  std::invalid_argument,
662  "Ifpack2::Chebyshev does not and will never implement the Ifpack "
663  "parameter \"chebyshev: solve normal equations\". If you want to solve "
664  "the normal equations, construct a Tpetra::Operator that implements "
665  "A^* A, and use Chebyshev to solve A^* A x = A^* b.");
666 
667  // Test for Ifpack parameters that we haven't implemented yet.
668  //
669  // For now, we only check that this ML parameter, if provided, has
670  // the one value that we support. We consider other values "invalid
671  // arguments" rather than "logic errors," because Ifpack does not
672  // implement eigenanalyses other than the power method.
673  std::string eigenAnalysisType ("power-method");
674  if (plist.isParameter ("eigen-analysis: type")) {
675  eigenAnalysisType = plist.get<std::string> ("eigen-analysis: type");
677  eigenAnalysisType != "power-method" &&
678  eigenAnalysisType != "power method",
679  std::invalid_argument,
680  "Ifpack2::Chebyshev: This class supports the ML parameter \"eigen-"
681  "analysis: type\" for backwards compatibility. However, Ifpack2 "
682  "currently only supports the \"power-method\" option for this "
683  "parameter. This imitates Ifpack, which only implements the power "
684  "method for eigenanalysis.");
685  }
686 
687  // We've validated all the parameters, so it's safe now to "commit" them.
688  userInvDiag_ = userInvDiagCopy;
689  userLambdaMax_ = lambdaMax;
690  userLambdaMin_ = lambdaMin;
691  userEigRatio_ = eigRatio;
692  boostFactor_ = static_cast<MT> (boostFactor);
693  minDiagVal_ = minDiagVal;
694  numIters_ = numIters;
695  eigMaxIters_ = eigMaxIters;
696  zeroStartingSolution_ = zeroStartingSolution;
697  assumeMatrixUnchanged_ = assumeMatrixUnchanged;
698  textbookAlgorithm_ = textbookAlgorithm;
699  computeMaxResNorm_ = computeMaxResNorm;
700  debug_ = debug;
701 
702  if (debug_) {
703  // Only print if myRank == 0.
704  int myRank = -1;
705  if (A_.is_null () || A_->getComm ().is_null ()) {
706  // We don't have a communicator (yet), so we assume that
707  // everybody can print. Revise this expectation in setMatrix().
708  myRank = 0;
709  }
710  else {
711  myRank = A_->getComm ()->getRank ();
712  }
713 
714  if (myRank == 0) {
715  out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
716  }
717  else {
719  out_ = Teuchos::getFancyOStream (blackHole); // print nothing on other processes
720  }
721  }
722  else { // NOT debug
723  // free the "old" output stream, if there was one
724  out_ = Teuchos::null;
725  }
726 }
727 
728 
729 template<class ScalarType, class MV>
730 void
732 {
733  D_ = Teuchos::null;
734  diagOffsets_ = offsets_type ();
735  savedDiagOffsets_ = false;
736  V_ = Teuchos::null;
737  W_ = Teuchos::null;
738  computedLambdaMax_ = STS::nan ();
739  computedLambdaMin_ = STS::nan ();
740 }
741 
742 
743 template<class ScalarType, class MV>
744 void
746 {
747  if (A.getRawPtr () != A_.getRawPtr ()) {
748  if (! assumeMatrixUnchanged_) {
749  reset ();
750  }
751  A_ = A;
752 
753  // The communicator may have changed, or we may not have had a
754  // communicator before. Thus, we may have to reset the debug
755  // output stream.
756  if (debug_) {
757  // Only print if myRank == 0.
758  int myRank = -1;
759  if (A.is_null () || A->getComm ().is_null ()) {
760  // We don't have a communicator (yet), so we assume that
761  // everybody can print. Revise this expectation in setMatrix().
762  myRank = 0;
763  }
764  else {
765  myRank = A->getComm ()->getRank ();
766  }
767 
768  if (myRank == 0) {
769  out_ = Teuchos::getFancyOStream (Teuchos::rcpFromRef (std::cerr));
770  }
771  else {
773  out_ = Teuchos::getFancyOStream (blackHole); // print nothing on other processes
774  }
775  }
776  else { // NOT debug
777  // free the "old" output stream, if there was one
778  out_ = Teuchos::null;
779  }
780  }
781 }
782 
783 
784 template<class ScalarType, class MV>
785 void
787 {
788  using std::endl;
789  // Some of the optimizations below only work if A_ is a
790  // Tpetra::CrsMatrix. We'll make our best guess about its type
791  // here, since we have no way to get back the original fifth
792  // template parameter.
793  typedef Tpetra::CrsMatrix<typename MV::scalar_type,
794  typename MV::local_ordinal_type,
795  typename MV::global_ordinal_type,
796  typename MV::node_type> crs_matrix_type;
797 
799  A_.is_null (), std::runtime_error, "Ifpack2::Chebyshev::compute: The input "
800  "matrix A is null. Please call setMatrix() with a nonnull input matrix "
801  "before calling this method.");
802 
803  // If A_ is a CrsMatrix and its graph is constant, we presume that
804  // the user plans to reuse the structure of A_, but possibly change
805  // A_'s values before each compute() call. This is the intended use
806  // case for caching the offsets of the diagonal entries of A_, to
807  // speed up extraction of diagonal entries on subsequent compute()
808  // calls.
809 
810  // FIXME (mfh 22 Jan 2013, 10 Feb 2013) In all cases when we use
811  // isnaninf() in this method, we really only want to check if the
812  // number is NaN. Inf means something different. However,
813  // Teuchos::ScalarTraits doesn't distinguish the two cases.
814 
815  // makeInverseDiagonal() returns a range Map Vector.
816  if (userInvDiag_.is_null ()) {
818  Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
819 
820  if (D_.is_null ()) { // We haven't computed D_ before
821  if (! A_crsMat.is_null () && A_crsMat->isStaticGraph ()) {
822  // It's a CrsMatrix with a const graph; cache diagonal offsets.
823  const size_t lclNumRows = A_crsMat->getNodeNumRows ();
824  if (diagOffsets_.extent (0) < lclNumRows) {
825  diagOffsets_ = offsets_type (); // clear first to save memory
826  diagOffsets_ = offsets_type ("offsets", lclNumRows);
827  }
828  A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
829  savedDiagOffsets_ = true;
830  D_ = makeInverseDiagonal (*A_, true);
831  }
832  else { // either A_ is not a CrsMatrix, or its graph is nonconst
833  D_ = makeInverseDiagonal (*A_);
834  }
835  }
836  else if (! assumeMatrixUnchanged_) { // D_ exists but A_ may have changed
837  if (! A_crsMat.is_null () && A_crsMat->isStaticGraph ()) {
838  // It's a CrsMatrix with a const graph; cache diagonal offsets
839  // if we haven't already.
840  if (! savedDiagOffsets_) {
841  const size_t lclNumRows = A_crsMat->getNodeNumRows ();
842  if (diagOffsets_.extent (0) < lclNumRows) {
843  diagOffsets_ = offsets_type (); // clear first to save memory
844  diagOffsets_ = offsets_type ("offsets", lclNumRows);
845  }
846  A_crsMat->getCrsGraph ()->getLocalDiagOffsets (diagOffsets_);
847  savedDiagOffsets_ = true;
848  }
849  // Now we're guaranteed to have cached diagonal offsets.
850  D_ = makeInverseDiagonal (*A_, true);
851  }
852  else { // either A_ is not a CrsMatrix, or its graph is nonconst
853  D_ = makeInverseDiagonal (*A_);
854  }
855  }
856  }
857  else { // the user provided an inverse diagonal
858  D_ = makeRangeMapVectorConst (userInvDiag_);
859  }
860 
861  // Have we estimated eigenvalues before?
862  const bool computedEigenvalueEstimates =
863  STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
864 
865  // Only recompute the eigenvalue estimates if
866  // - we are supposed to assume that the matrix may have changed, or
867  // - they haven't been computed before, and the user hasn't given
868  // us at least an estimate of the max eigenvalue.
869  //
870  // We at least need an estimate of the max eigenvalue. This is the
871  // most important one if using Chebyshev as a smoother.
872  if (! assumeMatrixUnchanged_ ||
873  (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
874  const ST computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
876  STS::isnaninf (computedLambdaMax),
877  std::runtime_error,
878  "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue "
879  "of D^{-1} A failed, by producing Inf or NaN. This probably means that "
880  "the matrix contains Inf or NaN values, or that it is badly scaled.");
882  STS::isnaninf (userEigRatio_),
883  std::logic_error,
884  "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
885  << endl << "This should be impossible." << endl <<
886  "Please report this bug to the Ifpack2 developers.");
887 
888  // The power method doesn't estimate the min eigenvalue, so we
889  // do our best to provide an estimate. userEigRatio_ has a
890  // reasonable default value, and if the user provided it, we
891  // have already checked that its value is finite and >= 1.
892  const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
893 
894  // Defer "committing" results until all computations succeeded.
895  computedLambdaMax_ = computedLambdaMax;
896  computedLambdaMin_ = computedLambdaMin;
897  } else {
899  STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
900  std::logic_error,
901  "Ifpack2::Chebyshev::compute: " << endl <<
902  "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
903  << endl << "This should be impossible." << endl <<
904  "Please report this bug to the Ifpack2 developers.");
905  }
906 
908  // Figure out the eigenvalue estimates that apply() will use.
910 
911  // Always favor the user's max eigenvalue estimate, if provided.
912  if (STS::isnaninf (userLambdaMax_)) {
913  lambdaMaxForApply_ = computedLambdaMax_;
914  } else {
915  lambdaMaxForApply_ = userLambdaMax_;
916  }
917  // mfh 11 Feb 2013: For now, we imitate Ifpack by ignoring the
918  // user's min eigenvalue estimate, and using the given eigenvalue
919  // ratio to estimate the min eigenvalue. We could instead do this:
920  // favor the user's eigenvalue ratio estimate, but if it's not
921  // provided, use lambdaMax / lambdaMin. However, we want Chebyshev
922  // to have sensible smoother behavior if the user did not provide
923  // eigenvalue estimates. Ifpack's behavior attempts to push down
924  // the error terms associated with the largest eigenvalues, while
925  // expecting that users will only want a small number of iterations,
926  // so that error terms associated with the smallest eigenvalues
927  // won't grow too much. This is sensible behavior for a smoother.
928  lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
929  eigRatioForApply_ = userEigRatio_;
930 
931  if (! textbookAlgorithm_) {
932  // Ifpack has a special-case modification of the eigenvalue bounds
933  // for the case where the max eigenvalue estimate is close to one.
934  const ST one = Teuchos::as<ST> (1);
935  // FIXME (mfh 20 Nov 2013) Should scale this 1.0e-6 term
936  // appropriately for MT's machine precision.
937  if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
938  lambdaMinForApply_ = one;
939  lambdaMaxForApply_ = lambdaMinForApply_;
940  eigRatioForApply_ = one; // Ifpack doesn't include this line.
941  }
942  }
943 }
944 
945 
946 template<class ScalarType, class MV>
947 ScalarType
949 getLambdaMaxForApply() const {
950  return lambdaMaxForApply_;
951 }
952 
953 
954 template<class ScalarType, class MV>
957 {
958  const char prefix[] = "Ifpack2::Chebyshev::apply: ";
959 
960  if (debug_) {
961  *out_ << "apply: " << std::endl;
962  }
964  (A_.is_null (), std::runtime_error, prefix << "The input matrix A is null. "
965  " Please call setMatrix() with a nonnull input matrix, and then call "
966  "compute(), before calling this method.");
968  (STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
969  prefix << "There is no estimate for the max eigenvalue."
970  << std::endl << computeBeforeApplyReminder);
972  (STS::isnaninf (lambdaMinForApply_), std::runtime_error,
973  prefix << "There is no estimate for the min eigenvalue."
974  << std::endl << computeBeforeApplyReminder);
976  (STS::isnaninf (eigRatioForApply_), std::runtime_error,
977  prefix <<"There is no estimate for the ratio of the max "
978  "eigenvalue to the min eigenvalue."
979  << std::endl << computeBeforeApplyReminder);
981  (D_.is_null (), std::runtime_error, prefix << "The vector of inverse "
982  "diagonal entries of the matrix has not yet been computed."
983  << std::endl << computeBeforeApplyReminder);
984 
985  if (textbookAlgorithm_) {
986  textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
987  lambdaMinForApply_, eigRatioForApply_, *D_);
988  }
989  else {
990  ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
991  lambdaMinForApply_, eigRatioForApply_, *D_);
992  }
993 
994  if (computeMaxResNorm_ && B.getNumVectors () > 0) {
995  MV R (B.getMap (), B.getNumVectors ());
996  computeResidual (R, B, *A_, X);
997  Teuchos::Array<MT> norms (B.getNumVectors ());
998  R.norm2 (norms ());
999  return *std::max_element (norms.begin (), norms.end ());
1000  }
1001  else {
1003  }
1004 }
1005 
1006 template<class ScalarType, class MV>
1007 void
1009 print (std::ostream& out) {
1010  this->describe (* (Teuchos::getFancyOStream (Teuchos::rcpFromRef (out))),
1012 }
1013 
1014 template<class ScalarType, class MV>
1015 void
1017 computeResidual (MV& R, const MV& B, const op_type& A, const MV& X,
1018  const Teuchos::ETransp mode)
1019 {
1020  Tpetra::deep_copy(R, B);
1021  A.apply (X, R, mode, -STS::one(), STS::one());
1022 }
1023 
1024 template<class ScalarType, class MV>
1025 void
1027 solve (MV& Z, const V& D_inv, const MV& R) {
1028  Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
1029 }
1030 
1031 template<class ScalarType, class MV>
1032 void
1033 Chebyshev<ScalarType, MV>::
1034 solve (MV& Z, const ST alpha, const V& D_inv, const MV& R) {
1035  Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
1036 }
1037 
1038 template<class ScalarType, class MV>
1040 Chebyshev<ScalarType, MV>::
1041 makeInverseDiagonal (const row_matrix_type& A, const bool useDiagOffsets) const
1042 {
1043  using Teuchos::RCP;
1044  using Teuchos::rcpFromRef;
1045  using Teuchos::rcp_dynamic_cast;
1046 
1047  RCP<V> D_rowMap (new V (A.getGraph ()->getRowMap ()));
1048  if (useDiagOffsets) {
1049  // The optimizations below only work if A_ is a Tpetra::CrsMatrix.
1050  // We'll make our best guess about its type here, since we have no
1051  // way to get back the original fifth template parameter.
1052  typedef Tpetra::CrsMatrix<typename MV::scalar_type,
1053  typename MV::local_ordinal_type,
1054  typename MV::global_ordinal_type,
1055  typename MV::node_type> crs_matrix_type;
1056  RCP<const crs_matrix_type> A_crsMat =
1057  rcp_dynamic_cast<const crs_matrix_type> (rcpFromRef (A));
1058  if (! A_crsMat.is_null ()) {
1060  ! savedDiagOffsets_, std::logic_error,
1061  "Ifpack2::Details::Chebyshev::makeInverseDiagonal: "
1062  "It is not allowed to call this method with useDiagOffsets=true, "
1063  "if you have not previously saved offsets of diagonal entries. "
1064  "This situation should never arise if this class is used properly. "
1065  "Please report this bug to the Ifpack2 developers.");
1066  A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_);
1067  }
1068  }
1069  else {
1070  // This always works for a Tpetra::RowMatrix, even if it is not a
1071  // Tpetra::CrsMatrix. We just don't have offsets in this case.
1072  A.getLocalDiagCopy (*D_rowMap);
1073  }
1074  RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
1075 
1076  if (debug_) {
1077  // In debug mode, make sure that all diagonal entries are
1078  // positive, on all processes. Note that *out_ only prints on
1079  // Process 0 of the matrix's communicator.
1080  D_rangeMap->template sync<Kokkos::HostSpace> ();
1081  auto D_lcl = D_rangeMap->template getLocalView<Kokkos::HostSpace> ();
1082  auto D_lcl_1d = Kokkos::subview (D_lcl, Kokkos::ALL (), 0);
1083 
1084  typedef typename MV::impl_scalar_type IST;
1085  typedef typename MV::local_ordinal_type LO;
1086  typedef Kokkos::Details::ArithTraits<IST> STS;
1087  typedef Kokkos::Details::ArithTraits<typename STS::mag_type> STM;
1088 
1089  const LO lclNumRows = static_cast<LO> (D_rangeMap->getLocalLength ());
1090  bool foundNonpositiveValue = false;
1091  for (LO i = 0; i < lclNumRows; ++i) {
1092  if (STS::real (D_lcl_1d(i)) <= STM::zero ()) {
1093  foundNonpositiveValue = true;
1094  break;
1095  }
1096  }
1097 
1098  using Teuchos::outArg;
1099  using Teuchos::REDUCE_MIN;
1100  using Teuchos::reduceAll;
1101 
1102  const int lclSuccess = foundNonpositiveValue ? 0 : 1;
1103  int gblSuccess = lclSuccess; // to be overwritten
1104  if (! D_rangeMap->getMap ().is_null () && D_rangeMap->getMap ()->getComm ().is_null ()) {
1105  const Teuchos::Comm<int>& comm = * (D_rangeMap->getMap ()->getComm ());
1106  reduceAll<int, int> (comm, REDUCE_MIN, lclSuccess, outArg (gblSuccess));
1107  }
1108  if (gblSuccess == 1) {
1109  *out_ << "makeInverseDiagonal: The matrix's diagonal entries all have "
1110  "positive real part (this is good for Chebyshev)." << std::endl;
1111  }
1112  else {
1113  *out_ << "makeInverseDiagonal: The matrix's diagonal has at least one "
1114  "entry with nonpositive real part, on at least one process in the "
1115  "matrix's communicator. This is bad for Chebyshev." << std::endl;
1116  }
1117  }
1118 
1119  // Invert the diagonal entries, replacing entries less (in
1120  // magnitude) than the user-specified value with that value.
1121  reciprocal_threshold (*D_rangeMap, minDiagVal_);
1122  return Teuchos::rcp_const_cast<const V> (D_rangeMap);
1123 }
1124 
1125 
1126 template<class ScalarType, class MV>
1128 Chebyshev<ScalarType, MV>::
1129 makeRangeMapVectorConst (const Teuchos::RCP<const V>& D) const
1130 {
1131  using Teuchos::RCP;
1132  using Teuchos::rcp;
1133  typedef Tpetra::Export<typename MV::local_ordinal_type,
1134  typename MV::global_ordinal_type,
1135  typename MV::node_type> export_type;
1136  // This throws logic_error instead of runtime_error, because the
1137  // methods that call makeRangeMapVector should all have checked
1138  // whether A_ is null before calling this method.
1140  A_.is_null (), std::logic_error, "Ifpack2::Details::Chebyshev::"
1141  "makeRangeMapVector: The input matrix A is null. Please call setMatrix() "
1142  "with a nonnull input matrix before calling this method. This is probably "
1143  "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1145  D.is_null (), std::logic_error, "Ifpack2::Details::Chebyshev::"
1146  "makeRangeMapVector: The input Vector D is null. This is probably "
1147  "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
1148 
1149  RCP<const map_type> sourceMap = D->getMap ();
1150  RCP<const map_type> rangeMap = A_->getRangeMap ();
1151  RCP<const map_type> rowMap = A_->getRowMap ();
1152 
1153  if (rangeMap->isSameAs (*sourceMap)) {
1154  // The given vector's Map is the same as the matrix's range Map.
1155  // That means we don't need to Export. This is the fast path.
1156  return D;
1157  }
1158  else { // We need to Export.
1159  RCP<const export_type> exporter;
1160  // Making an Export object from scratch is expensive enough that
1161  // it's worth the O(1) global reductions to call isSameAs(), to
1162  // see if we can avoid that cost.
1163  if (sourceMap->isSameAs (*rowMap)) {
1164  // We can reuse the matrix's Export object, if there is one.
1165  exporter = A_->getGraph ()->getExporter ();
1166  }
1167  else { // We have to make a new Export object.
1168  exporter = rcp (new export_type (sourceMap, rangeMap));
1169  }
1170 
1171  if (exporter.is_null ()) {
1172  return D; // Row Map and range Map are the same; no need to Export.
1173  }
1174  else { // Row Map and range Map are _not_ the same; must Export.
1175  RCP<V> D_out = rcp (new V (*D, Teuchos::Copy));
1176  D_out->doExport (*D, *exporter, Tpetra::ADD);
1177  return Teuchos::rcp_const_cast<const V> (D_out);
1178  }
1179  }
1180 }
1181 
1182 
1183 template<class ScalarType, class MV>
1185 Chebyshev<ScalarType, MV>::
1186 makeRangeMapVector (const Teuchos::RCP<V>& D) const
1187 {
1188  using Teuchos::rcp_const_cast;
1189  return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1190 }
1191 
1192 
1193 template<class ScalarType, class MV>
1194 void
1195 Chebyshev<ScalarType, MV>::
1196 textbookApplyImpl (const op_type& A,
1197  const MV& B,
1198  MV& X,
1199  const int numIters,
1200  const ST lambdaMax,
1201  const ST lambdaMin,
1202  const ST eigRatio,
1203  const V& D_inv) const
1204 {
1205  (void) lambdaMin; // Forestall compiler warning.
1206  const ST myLambdaMin = lambdaMax / eigRatio;
1207 
1208  const ST zero = Teuchos::as<ST> (0);
1209  const ST one = Teuchos::as<ST> (1);
1210  const ST two = Teuchos::as<ST> (2);
1211  const ST d = (lambdaMax + myLambdaMin) / two; // Ifpack2 calls this theta
1212  const ST c = (lambdaMax - myLambdaMin) / two; // Ifpack2 calls this 1/delta
1213 
1214  if (zeroStartingSolution_ && numIters > 0) {
1215  // If zero iterations, then input X is output X.
1216  X.putScalar (zero);
1217  }
1218  MV R (B.getMap (), B.getNumVectors (), false);
1219  MV P (B.getMap (), B.getNumVectors (), false);
1220  MV Z (B.getMap (), B.getNumVectors (), false);
1221  ST alpha, beta;
1222  for (int i = 0; i < numIters; ++i) {
1223  computeResidual (R, B, A, X); // R = B - A*X
1224  solve (Z, D_inv, R); // z = D_inv * R, that is, D \ R.
1225  if (i == 0) {
1226  P = Z;
1227  alpha = two / d;
1228  } else {
1229  //beta = (c * alpha / two)^2;
1230  //const ST sqrtBeta = c * alpha / two;
1231  //beta = sqrtBeta * sqrtBeta;
1232  beta = alpha * (c/two) * (c/two);
1233  alpha = one / (d - beta);
1234  P.update (one, Z, beta); // P = Z + beta*P
1235  }
1236  X.update (alpha, P, one); // X = X + alpha*P
1237  // If we compute the residual here, we could either do R = B -
1238  // A*X, or R = R - alpha*A*P. Since we choose the former, we
1239  // can move the computeResidual call to the top of the loop.
1240  }
1241 }
1242 
1243 template<class ScalarType, class MV>
1244 typename Chebyshev<ScalarType, MV>::MT
1245 Chebyshev<ScalarType, MV>::maxNormInf (const MV& X) {
1246  Teuchos::Array<MT> norms (X.getNumVectors ());
1247  X.normInf (norms());
1248  return *std::max_element (norms.begin (), norms.end ());
1249 }
1250 
1251 template<class ScalarType, class MV>
1252 void
1253 Chebyshev<ScalarType, MV>::
1254 ifpackApplyImpl (const op_type& A,
1255  const MV& B,
1256  MV& X,
1257  const int numIters,
1258  const ST lambdaMax,
1259  const ST lambdaMin,
1260  const ST eigRatio,
1261  const V& D_inv)
1262 {
1263  using std::endl;
1264 #ifdef HAVE_IFPACK2_DEBUG
1265  const bool debug = debug_;
1266 #else
1267  const bool debug = false;
1268 #endif
1269 
1270  if (debug) {
1271  *out_ << " \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1272  *out_ << " \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1273  }
1274 
1275  if (numIters <= 0) {
1276  return;
1277  }
1278  const ST one = static_cast<ST> (1.0);
1279  const ST two = static_cast<ST> (2.0);
1280 
1281  // Quick solve when the matrix A is the identity.
1282  if (lambdaMin == one && lambdaMax == lambdaMin) {
1283  solve (X, D_inv, B);
1284  return;
1285  }
1286 
1287  // Initialize coefficients
1288  const ST alpha = lambdaMax / eigRatio;
1289  const ST beta = boostFactor_ * lambdaMax;
1290  const ST delta = two / (beta - alpha);
1291  const ST theta = (beta + alpha) / two;
1292  const ST s1 = theta * delta;
1293 
1294  if (debug) {
1295  *out_ << " alpha = " << alpha << endl
1296  << " beta = " << beta << endl
1297  << " delta = " << delta << endl
1298  << " theta = " << theta << endl
1299  << " s1 = " << s1 << endl;
1300  }
1301 
1302  // Fetch cached temporary vectors.
1303  Teuchos::RCP<MV> V_ptr, W_ptr;
1304  makeTempMultiVectors (V_ptr, W_ptr, B);
1305 
1306  // mfh 28 Jan 2013: We write V1 instead of V, so as not to confuse
1307  // the multivector V with the typedef V (for Tpetra::Vector).
1308  //MV V1 (B.getMap (), B.getNumVectors (), false);
1309  //MV W (B.getMap (), B.getNumVectors (), false);
1310  MV& V1 = *V_ptr;
1311  MV& W = *W_ptr;
1312 
1313  if (debug) {
1314  *out_ << " Iteration " << 1 << ":" << endl
1315  << " - \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1316  }
1317 
1318  // Special case for the first iteration.
1319  if (! zeroStartingSolution_) {
1320  computeResidual (V1, B, A, X); // V1 = B - A*X
1321  if (debug) {
1322  *out_ << " - \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1323  }
1324 
1325  solve (W, one/theta, D_inv, V1); // W = (1/theta)*D_inv*(B-A*X)
1326  if (debug) {
1327  *out_ << " - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1328  }
1329  X.update (one, W, one); // X = X + W
1330  }
1331  else {
1332  solve (W, one/theta, D_inv, B); // W = (1/theta)*D_inv*B
1333  if (debug) {
1334  *out_ << " - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1335  }
1336  Tpetra::deep_copy (X, W); // X = 0 + W
1337  }
1338  if (debug) {
1339  *out_ << " - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1340  }
1341 
1342  // The rest of the iterations.
1343  ST rhok = one / s1;
1344  ST rhokp1, dtemp1, dtemp2;
1345  for (int deg = 1; deg < numIters; ++deg) {
1346  if (debug) {
1347  *out_ << " Iteration " << deg+1 << ":" << endl
1348  << " - \\|D\\|_{\\infty} = " << D_->normInf () << endl
1349  << " - \\|B\\|_{\\infty} = " << maxNormInf (B) << endl
1350  << " - \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm () << endl
1351  << " - rhok = " << rhok << endl;
1352  V1.putScalar (STS::zero ()); // ???????
1353  }
1354 
1355  computeResidual (V1, B, A, X); // V1 = B - A*X
1356  if (debug) {
1357  *out_ << " - \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1358  }
1359 
1360  rhokp1 = one / (two * s1 - rhok);
1361  dtemp1 = rhokp1 * rhok;
1362  dtemp2 = two * rhokp1 * delta;
1363  rhok = rhokp1;
1364 
1365  if (debug) {
1366  *out_ << " - dtemp1 = " << dtemp1 << endl
1367  << " - dtemp2 = " << dtemp2 << endl;
1368  }
1369 
1370 
1371  W.elementWiseMultiply (dtemp2, D_inv, V1, dtemp1);
1372  X.update (one, W, one);
1373 
1374  if (debug) {
1375  *out_ << " - \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1376  *out_ << " - \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1377  }
1378  }
1379 }
1380 
1381 template<class ScalarType, class MV>
1382 typename Chebyshev<ScalarType, MV>::ST
1383 Chebyshev<ScalarType, MV>::
1384 powerMethodWithInitGuess (const op_type& A,
1385  const V& D_inv,
1386  const int numIters,
1387  V& x)
1388 {
1389  using std::endl;
1390  if (debug_) {
1391  *out_ << " powerMethodWithInitGuess:" << endl;
1392  }
1393 
1394  const ST zero = static_cast<ST> (0.0);
1395  const ST one = static_cast<ST> (1.0);
1396  ST lambdaMax = zero;
1397  ST RQ_top, RQ_bottom, norm;
1398 
1399  V y (A.getRangeMap ());
1400  norm = x.norm2 ();
1402  (norm == zero, std::runtime_error,
1403  "Ifpack2::Chebyshev::powerMethodWithInitGuess: "
1404  "The initial guess has zero norm. This could be either because Tpetra::"
1405  "Vector::randomize() filled the vector with zeros (if that was used to "
1406  "compute the initial guess), or because the norm2 method has a bug. The "
1407  "first is not impossible, but unlikely.");
1408 
1409  if (debug_) {
1410  *out_ << " Original norm1(x): " << x.norm1 () << ", norm2(x): " << norm << endl;
1411  }
1412 
1413  x.scale (one / norm);
1414 
1415  if (debug_) {
1416  *out_ << " norm1(x.scale(one/norm)): " << x.norm1 () << endl;
1417  }
1418 
1419  for (int iter = 0; iter < numIters; ++iter) {
1420  if (debug_) {
1421  *out_ << " Iteration " << iter << endl;
1422  }
1423  A.apply (x, y);
1424  solve (y, D_inv, y);
1425  RQ_top = y.dot (x);
1426  RQ_bottom = x.dot (x);
1427  if (debug_) {
1428  *out_ << " RQ_top: " << RQ_top << ", RQ_bottom: " << RQ_bottom << endl;
1429  }
1430  lambdaMax = RQ_top / RQ_bottom;
1431  norm = y.norm2 ();
1432  if (norm == zero) { // Return something reasonable.
1433  if (debug_) {
1434  *out_ << " norm is zero; returning zero" << endl;
1435  }
1436  return zero;
1437  }
1438  x.update (one / norm, y, zero);
1439  }
1440  if (debug_) {
1441  *out_ << " lambdaMax: " << lambdaMax << endl;
1442  }
1443  return lambdaMax;
1444 }
1445 
1446 template<class ScalarType, class MV>
1447 void
1448 Chebyshev<ScalarType, MV>::
1449 computeInitialGuessForPowerMethod (V& x, const bool nonnegativeRealParts) const
1450 {
1451  typedef typename MV::device_type::execution_space dev_execution_space;
1452  typedef typename MV::device_type::memory_space dev_memory_space;
1453  typedef typename MV::local_ordinal_type LO;
1454 
1455  x.randomize ();
1456 
1457  if (nonnegativeRealParts) {
1458  x.template modify<dev_memory_space> ();
1459  auto x_lcl = x.template getLocalView<dev_memory_space> ();
1460  auto x_lcl_1d = Kokkos::subview (x_lcl, Kokkos::ALL (), 0);
1461 
1462  const LO lclNumRows = static_cast<LO> (x.getLocalLength ());
1463  Kokkos::RangePolicy<dev_execution_space, LO> range (0, lclNumRows);
1464  PositivizeVector<decltype (x_lcl_1d), LO> functor (x_lcl_1d);
1465  Kokkos::parallel_for (range, functor);
1466  }
1467 }
1468 
1469 template<class ScalarType, class MV>
1470 typename Chebyshev<ScalarType, MV>::ST
1471 Chebyshev<ScalarType, MV>::
1472 powerMethod (const op_type& A, const V& D_inv, const int numIters)
1473 {
1474  using std::endl;
1475  if (debug_) {
1476  *out_ << "powerMethod:" << endl;
1477  }
1478 
1479  const ST zero = static_cast<ST> (0.0);
1480  V x (A.getDomainMap ());
1481  // For the first pass, just let the pseudorandom number generator
1482  // fill x with whatever values it wants; don't try to make its
1483  // entries nonnegative.
1484  computeInitialGuessForPowerMethod (x, false);
1485 
1486  ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1487 
1488  // mfh 07 Jan 2015: Taking the real part here is only a concession
1489  // to the compiler, so that this class can build with ScalarType =
1490  // std::complex<T>. Our Chebyshev implementation only works with
1491  // real, symmetric positive definite matrices. The right thing to
1492  // do would be what Belos does, which is provide a partial
1493  // specialization for ScalarType = std::complex<T> with a stub
1494  // implementation (that builds, but whose constructor throws).
1495  if (STS::real (lambdaMax) < STS::real (zero)) {
1496  if (debug_) {
1497  *out_ << "real(lambdaMax) = " << STS::real (lambdaMax) << " < 0; try "
1498  "again with a different random initial guess" << endl;
1499  }
1500  // Max eigenvalue estimate was negative. Perhaps we got unlucky
1501  // with the random initial guess. Try again with a different (but
1502  // still random) initial guess. Only try again once, so that the
1503  // run time is bounded.
1504 
1505  // For the second pass, make all the entries of the initial guess
1506  // vector have nonnegative real parts.
1507  computeInitialGuessForPowerMethod (x, true);
1508  lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1509  }
1510  return lambdaMax;
1511 }
1512 
1513 template<class ScalarType, class MV>
1516  return A_;
1517 }
1518 
1519 template<class ScalarType, class MV>
1520 bool
1523  // Technically, this is true, because the matrix must be symmetric.
1524  return true;
1525 }
1526 
1527 template<class ScalarType, class MV>
1528 void
1531  Teuchos::RCP<MV>& W,
1532  const MV& X)
1533 {
1534  // ETP 02/08/17: We must check not only if the temporary vectors are
1535  // null, but also if the number of columns match, since some multi-RHS
1536  // solvers (e.g., Belos) may call apply() with different numbers of columns.
1537  if (V_.is_null () || V_->getNumVectors () != X.getNumVectors ()) {
1538  V_ = Teuchos::rcp (new MV (X.getMap (), X.getNumVectors (), false));
1539  }
1540  //W must be initialized to zero when it is used as a multigrid smoother.
1541  if (W_.is_null () || W_->getNumVectors () != X.getNumVectors ()) {
1542  W_ = Teuchos::rcp (new MV (X.getMap (), X.getNumVectors (), true));
1543  }
1544  V1 = V_;
1545  W = W_;
1546 }
1547 
1548 template<class ScalarType, class MV>
1549 std::string
1550 Chebyshev<ScalarType, MV>::
1551 description () const {
1552  std::ostringstream oss;
1553  // YAML requires quoting the key in this case, to distinguish
1554  // key's colons from the colon that separates key from value.
1555  oss << "\"Ifpack2::Details::Chebyshev\":"
1556  << "{"
1557  << "degree: " << numIters_
1558  << ", lambdaMax: " << lambdaMaxForApply_
1559  << ", alpha: " << eigRatioForApply_
1560  << ", lambdaMin: " << lambdaMinForApply_
1561  << ", boost factor: " << boostFactor_
1562  << "}";
1563  return oss.str();
1564 }
1565 
1566 template<class ScalarType, class MV>
1567 void
1570  const Teuchos::EVerbosityLevel verbLevel) const
1571 {
1573  using std::endl;
1574 
1575  const Teuchos::EVerbosityLevel vl =
1576  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1577  if (vl == Teuchos::VERB_NONE) {
1578  return; // print NOTHING
1579  }
1580 
1581  // By convention, describe() starts with a tab.
1582  //
1583  // This does affect all processes on which it's valid to print to
1584  // 'out'. However, it does not actually print spaces to 'out'
1585  // unless operator<< gets called, so it's safe to use on all
1586  // processes.
1587  Teuchos::OSTab tab0 (out);
1588 
1589  // We only print on Process 0 of the matrix's communicator. If
1590  // the matrix isn't set, we don't have a communicator, so we have
1591  // to assume that every process can print.
1592  int myRank = -1;
1593  if (A_.is_null () || A_->getComm ().is_null ()) {
1594  myRank = 0;
1595  }
1596  else {
1597  myRank = A_->getComm ()->getRank ();
1598  }
1599  if (myRank == 0) {
1600  // YAML requires quoting the key in this case, to distinguish
1601  // key's colons from the colon that separates key from value.
1602  out << "\"Ifpack2::Details::Chebyshev\":" << endl;
1603  }
1604  Teuchos::OSTab tab1 (out);
1605 
1606  if (vl == Teuchos::VERB_LOW) {
1607  if (myRank == 0) {
1608  out << "degree: " << numIters_ << endl
1609  << "lambdaMax: " << lambdaMaxForApply_ << endl
1610  << "alpha: " << eigRatioForApply_ << endl
1611  << "lambdaMin: " << lambdaMinForApply_ << endl
1612  << "boost factor: " << boostFactor_ << endl;
1613  }
1614  return;
1615  }
1616 
1617  // vl > Teuchos::VERB_LOW
1618 
1619  if (myRank == 0) {
1620  out << "Template parameters:" << endl;
1621  {
1622  Teuchos::OSTab tab2 (out);
1623  out << "ScalarType: " << TypeNameTraits<ScalarType>::name () << endl
1624  << "MV: " << TypeNameTraits<MV>::name () << endl;
1625  }
1626 
1627  // "Computed parameters" literally means "parameters whose
1628  // values were computed by compute()."
1629  if (myRank == 0) {
1630  out << "Computed parameters:" << endl;
1631  }
1632  }
1633 
1634  // Print computed parameters
1635  {
1636  Teuchos::OSTab tab2 (out);
1637  // Users might want to see the values in the computed inverse
1638  // diagonal, so we print them out at the highest verbosity.
1639  if (myRank == 0) {
1640  out << "D_: ";
1641  }
1642  if (D_.is_null ()) {
1643  if (myRank == 0) {
1644  out << "unset" << endl;
1645  }
1646  }
1647  else if (vl <= Teuchos::VERB_HIGH) {
1648  if (myRank == 0) {
1649  out << "set" << endl;
1650  }
1651  }
1652  else { // D_ not null and vl > Teuchos::VERB_HIGH
1653  if (myRank == 0) {
1654  out << endl;
1655  }
1656  // By convention, describe() first indents, then prints.
1657  // We can rely on other describe() implementations to do that.
1658  D_->describe (out, vl);
1659  }
1660  if (myRank == 0) {
1661  // V_ and W_ are scratch space; their values are irrelevant.
1662  // All that matters is whether or not they have been set.
1663  out << "V_: " << (V_.is_null () ? "unset" : "set") << endl
1664  << "W_: " << (W_.is_null () ? "unset" : "set") << endl
1665  << "computedLambdaMax_: " << computedLambdaMax_ << endl
1666  << "computedLambdaMin_: " << computedLambdaMin_ << endl
1667  << "lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1668  << "lambdaMinForApply_: " << lambdaMinForApply_ << endl
1669  << "eigRatioForApply_: " << eigRatioForApply_ << endl;
1670  }
1671  } // print computed parameters
1672 
1673  if (myRank == 0) {
1674  out << "User parameters:" << endl;
1675  }
1676 
1677  // Print user parameters
1678  {
1679  Teuchos::OSTab tab2 (out);
1680  out << "userInvDiag_: ";
1681  if (userInvDiag_.is_null ()) {
1682  out << "unset" << endl;
1683  }
1684  else if (vl <= Teuchos::VERB_HIGH) {
1685  out << "set" << endl;
1686  }
1687  else { // userInvDiag_ not null and vl > Teuchos::VERB_HIGH
1688  if (myRank == 0) {
1689  out << endl;
1690  }
1691  userInvDiag_->describe (out, vl);
1692  }
1693  if (myRank == 0) {
1694  out << "userLambdaMax_: " << userLambdaMax_ << endl
1695  << "userLambdaMin_: " << userLambdaMin_ << endl
1696  << "userEigRatio_: " << userEigRatio_ << endl
1697  << "numIters_: " << numIters_ << endl
1698  << "eigMaxIters_: " << eigMaxIters_ << endl
1699  << "zeroStartingSolution_: " << zeroStartingSolution_ << endl
1700  << "assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1701  << "textbookAlgorithm_: " << textbookAlgorithm_ << endl
1702  << "computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1703  }
1704  } // print user parameters
1705 }
1706 
1707 } // namespace Details
1708 } // namespace Ifpack2
1709 
1710 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \
1711  template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >;
1712 
1713 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
VERB_DEFAULT
VERB_LOW
VERB_NONE
T & get(const std::string &name, T def_value)
ParameterList & set(std::string const &name, T const &value, std::string const &docString="", RCP< const ParameterEntryValidator > const &validator=null)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:276
bool is_null() const
VERB_MEDIUM
Ifpack2 implementation details.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
VERB_HIGH
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:199
LinearOp zero(const VectorSpace &vs)
bool isType(const std::string &name) const
TEUCHOS_DEPRECATED void reduceAll(const Comm< Ordinal > &comm, const EReductionType reductType, const Packet &send, Packet *globalReduct)
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:111
scalar_type ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:107
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:101
Declaration of Chebyshev implementation.
TypeTo as(const TypeFrom &t)
T * getRawPtr() const
bool isParameter(const std::string &name) const
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
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:126