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: Tempated 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 // ***********************************************************************
45 //
46 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
47 // Copyright (2004) Sandia Corporation
48 //
49 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
50 // license for use of this work by or on behalf of the U.S. Government.
51 //
52 // This library is free software; you can redistribute it and/or modify
53 // it under the terms of the GNU Lesser General Public License as
54 // published by the Free Software Foundation; either version 2.1 of the
55 // License, or (at your option) any later version.
56 //
57 // This library is distributed in the hope that it will be useful, but
58 // WITHOUT ANY WARRANTY; without even the implied warranty of
59 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
60 // Lesser General Public License for more details.
61 //
62 // You should have received a copy of the GNU Lesser General Public
63 // License along with this library; if not, write to the Free Software
64 // Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301
65 // USA
66 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
67 //
68 // ***********************************************************************
69 
70 #ifndef IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
71 #define IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
72 
79 
81 #include <Kokkos_ArithTraits.hpp>
82 #include <cmath>
83 
84 // Uncommit the #define line below if you want Chebyshev to do extra
85 // debug checking and generate a lot of debugging output to stderr (on
86 // all processes in the communicator). Even if you uncomment this
87 // line, the debugging code will only be enabled if the CMake option
88 // Teuchos_ENABLE_DEBUG was set to ON when configuring Trilinos.
89 //#define IFPACK_DETAILS_CHEBYSHEV_DEBUG 1
90 
91 namespace Ifpack2 {
92 namespace Details {
93 
94 namespace { // (anonymous)
95 
96 // We use this text a lot in error messages.
97 const char computeBeforeApplyReminder[] =
98  "This means one of the following:\n"
99  " - you have not yet called compute() on this instance, or \n"
100  " - you didn't call compute() after calling setParameters().\n\n"
101  "After creating an Ifpack2::Chebyshev instance,\n"
102  "you must _always_ call compute() at least once before calling apply().\n"
103  "After calling compute() once, you do not need to call it again,\n"
104  "unless the matrix has changed or you have changed parameters\n"
105  "(by calling setParameters()).";
106 
107 } // namespace (anonymous)
108 
109 // ReciprocalThreshold stuff below needs to be in a namspace visible outside
110 // of this file
111 template<class XV, class SizeType = typename XV::size_type>
112 struct V_ReciprocalThresholdSelfFunctor
113 {
114  typedef typename XV::execution_space execution_space;
115  typedef typename XV::non_const_value_type value_type;
116  typedef SizeType size_type;
117  typedef Kokkos::Details::ArithTraits<value_type> KAT;
118  typedef typename KAT::mag_type mag_type;
119 
120  XV X_;
121  const value_type m_min_val;
122  const mag_type m_min_val_mag;
123 
124  V_ReciprocalThresholdSelfFunctor (const XV& X,
125  const value_type& min_val) :
126  X_ (X),
127  m_min_val (min_val),
128  m_min_val_mag (KAT::abs (min_val))
129  {}
130 
131  KOKKOS_INLINE_FUNCTION
132  void operator() (const size_type& i) const
133  {
134  if (KAT::abs (X_(i)) < m_min_val_mag) {
135  X_[i] = m_min_val;
136  }
137  else {
138  X_[i] = KAT::one () / X_[i];
139  }
140  }
141 };
142 
143 template<class XV, class SizeType = typename XV::size_type>
144 struct LocalReciprocalThreshold {
145  static void
146  compute (const XV& X,
147  const typename XV::non_const_value_type& minVal)
148  {
149  typedef typename XV::execution_space execution_space;
150  Kokkos::RangePolicy<execution_space, SizeType> policy (0, X.dimension_0 ());
151  V_ReciprocalThresholdSelfFunctor<XV, SizeType> op (X, minVal);
152  Kokkos::parallel_for (policy, op);
153  }
154 };
155 
156 template <class TpetraVectorType,
157  const bool classic = TpetraVectorType::node_type::classic>
158 struct GlobalReciprocalThreshold {};
159 
160 template <class TpetraVectorType>
161 struct GlobalReciprocalThreshold<TpetraVectorType, true> {
162  static void
163  compute (TpetraVectorType& V,
164  const typename TpetraVectorType::scalar_type& min_val)
165  {
166  typedef typename TpetraVectorType::scalar_type scalar_type;
167  typedef typename TpetraVectorType::mag_type mag_type;
168  typedef Kokkos::Details::ArithTraits<scalar_type> STS;
169 
170  const scalar_type ONE = STS::one ();
171  const mag_type min_val_abs = STS::abs (min_val);
172 
173  Teuchos::ArrayRCP<scalar_type> V_0 = V.getDataNonConst (0);
174  const size_t lclNumRows = V.getLocalLength ();
175 
176  for (size_t i = 0; i < lclNumRows; ++i) {
177  const scalar_type V_0i = V_0[i];
178  if (STS::abs (V_0i) < min_val_abs) {
179  V_0[i] = min_val;
180  } else {
181  V_0[i] = ONE / V_0i;
182  }
183  }
184  }
185 };
186 
187 template <class TpetraVectorType>
188 struct GlobalReciprocalThreshold<TpetraVectorType, false> {
189  static void
190  compute (TpetraVectorType& X,
191  const typename TpetraVectorType::scalar_type& minVal)
192  {
193  typedef typename TpetraVectorType::dual_view_type::non_const_value_type value_type;
194  typedef typename TpetraVectorType::execution_space::memory_space memory_space;
195 
196  auto X_lcl = X.getDualView ();
197  X_lcl.template sync<memory_space> ();
198  X_lcl.template modify<memory_space> ();
199 
200  const value_type minValS = static_cast<value_type> (minVal);
201 
202  auto X_0 = Kokkos::subview (X_lcl.d_view, Kokkos::ALL (), 0);
203  LocalReciprocalThreshold<decltype (X_0) >::compute (X_0, minValS);
204  }
205 };
206 
207 // Utility function for inverting diagonal with threshold.
208 template <typename S, typename L, typename G, typename N>
209 void
210 reciprocal_threshold (Tpetra::Vector<S,L,G,N>& V, const S& minVal)
211 {
212  GlobalReciprocalThreshold<Tpetra::Vector<S,L,G,N> >::compute (V, minVal);
213 }
214 
215 template<class ScalarType, class MV>
216 void Chebyshev<ScalarType, MV>::checkInputMatrix () const
217 {
218  TEUCHOS_TEST_FOR_EXCEPTION(
219  ! A_.is_null () && A_->getGlobalNumRows () != A_->getGlobalNumCols (),
220  std::invalid_argument,
221  "Ifpack2::Chebyshev: The input matrix A must be square. "
222  "A has " << A_->getGlobalNumRows () << " rows and "
223  << A_->getGlobalNumCols () << " columns.");
224 
225  // In a debug build, test that the domain and range Maps of the
226  // matrix are the same.
227 #ifdef HAVE_TEUCHOS_DEBUG
228  if (! A_.is_null ()) {
229  Teuchos::RCP<const map_type> domainMap = A_->getDomainMap ();
230  Teuchos::RCP<const map_type> rangeMap = A_->getRangeMap ();
231 
232  // isSameAs is a collective, but if the two pointers are the same,
233  // isSameAs will assume that they are the same on all processes, and
234  // return true without an all-reduce.
235  TEUCHOS_TEST_FOR_EXCEPTION(
236  ! domainMap->isSameAs (*rangeMap), std::invalid_argument,
237  "Ifpack2::Chebyshev: The domain Map and range Map of the matrix must be "
238  "the same (in the sense of isSameAs())." << std::endl << "We only check "
239  "for this if Trilinos was built with the CMake configuration option "
240  "Teuchos_ENABLE_DEBUG set to ON.");
241  }
242 #endif // HAVE_TEUCHOS_DEBUG
243 }
244 
245 
246 template<class ScalarType, class MV>
247 void
248 Chebyshev<ScalarType, MV>::
249 checkConstructorInput () const
250 {
251  TEUCHOS_TEST_FOR_EXCEPTION(STS::isComplex, std::logic_error,
252  "Ifpack2::Chebyshev: This class' implementation of Chebyshev iteration "
253  "only works for real-valued, symmetric positive definite matrices. "
254  "However, you instantiated this class for ScalarType = "
255  << Teuchos::TypeNameTraits<ScalarType>::name () << ", which is a complex-"
256  "valued type. While this may be algorithmically correct if all of the "
257  "complex numbers in the matrix have zero imaginary part, we forbid using "
258  "complex ScalarType altogether in order to remind you of the limitations "
259  "of our implementation (and of the algorithm itself).");
260 
261  checkInputMatrix ();
262 }
263 
264 template<class ScalarType, class MV>
266 Chebyshev (Teuchos::RCP<const row_matrix_type> A) :
267  A_ (A),
268  savedDiagOffsets_ (false),
269  computedLambdaMax_ (STS::nan ()),
270  computedLambdaMin_ (STS::nan ()),
271  lambdaMaxForApply_ (STS::nan ()),
272  lambdaMinForApply_ (STS::nan ()),
273  eigRatioForApply_ (STS::nan ()),
274  userLambdaMax_ (STS::nan ()),
275  userLambdaMin_ (STS::nan ()),
276  userEigRatio_ (Teuchos::as<ST> (30)),
277  minDiagVal_ (STS::eps ()),
278  numIters_ (1),
279  eigMaxIters_ (10),
280  zeroStartingSolution_ (true),
281  assumeMatrixUnchanged_ (false),
282  textbookAlgorithm_ (false),
283  computeMaxResNorm_ (false)
284 {
285  checkConstructorInput ();
286 }
287 
288 template<class ScalarType, class MV>
290 Chebyshev (Teuchos::RCP<const row_matrix_type> A, Teuchos::ParameterList& params) :
291  A_ (A),
292  savedDiagOffsets_ (false),
293  computedLambdaMax_ (STS::nan ()),
294  computedLambdaMin_ (STS::nan ()),
295  lambdaMaxForApply_ (STS::nan ()),
296  lambdaMinForApply_ (STS::nan ()),
297  eigRatioForApply_ (STS::nan ()),
298  userLambdaMax_ (STS::nan ()),
299  userLambdaMin_ (STS::nan ()),
300  userEigRatio_ (Teuchos::as<ST> (30)),
301  minDiagVal_ (STS::eps ()),
302  numIters_ (1),
303  eigMaxIters_ (10),
304  zeroStartingSolution_ (true),
305  assumeMatrixUnchanged_ (false),
306  textbookAlgorithm_ (false),
307  computeMaxResNorm_ (false)
308 {
309  checkConstructorInput ();
310  setParameters (params);
311 }
312 
313 template<class ScalarType, class MV>
314 void
316 setParameters (Teuchos::ParameterList& plist)
317 {
318  using Teuchos::RCP;
319  using Teuchos::rcp;
320  using Teuchos::rcp_const_cast;
321 
322  // Note to developers: The logic for this method is complicated,
323  // because we want to accept Ifpack and ML parameters whenever
324  // possible, but we don't want to add their default values to the
325  // user's ParameterList. That's why we do all the isParameter()
326  // checks, instead of using the two-argument version of get()
327  // everywhere. The min and max eigenvalue parameters are also a
328  // special case, because we decide whether or not to do eigenvalue
329  // analysis based on whether the user supplied the max eigenvalue.
330 
331  // Default values of all the parameters.
332  const ST defaultLambdaMax = STS::nan ();
333  const ST defaultLambdaMin = STS::nan ();
334  // 30 is Ifpack::Chebyshev's default. ML has a different default
335  // eigRatio for smoothers and the coarse-grid solve (if using
336  // Chebyshev for that). The former uses 20; the latter uses 30.
337  // We're testing smoothers here, so use 20. (However, if you give
338  // ML an Epetra matrix, it will use Ifpack for Chebyshev, in which
339  // case it would defer to Ifpack's default settings.)
340  const ST defaultEigRatio = Teuchos::as<ST> (30);
341  const ST defaultMinDiagVal = STS::eps ();
342  const int defaultNumIters = 1;
343  const int defaultEigMaxIters = 10;
344  const bool defaultZeroStartingSolution = true; // Ifpack::Chebyshev default
345  const bool defaultAssumeMatrixUnchanged = false;
346  const bool defaultTextbookAlgorithm = false;
347  const bool defaultComputeMaxResNorm = false;
348 
349  // We'll set the instance data transactionally, after all reads
350  // from the ParameterList. That way, if any of the ParameterList
351  // reads fail (e.g., due to the wrong parameter type), we will not
352  // have left the instance data in a half-changed state.
353  RCP<const V> userInvDiagCopy; // if nonnull: deep copy of user's Vector
354  ST lambdaMax = defaultLambdaMax;
355  ST lambdaMin = defaultLambdaMin;
356  ST eigRatio = defaultEigRatio;
357  ST minDiagVal = defaultMinDiagVal;
358  int numIters = defaultNumIters;
359  int eigMaxIters = defaultEigMaxIters;
360  bool zeroStartingSolution = defaultZeroStartingSolution;
361  bool assumeMatrixUnchanged = defaultAssumeMatrixUnchanged;
362  bool textbookAlgorithm = defaultTextbookAlgorithm;
363  bool computeMaxResNorm = defaultComputeMaxResNorm;
364 
365  // Fetch the parameters from the ParameterList. Defer all
366  // externally visible side effects until we have finished all
367  // ParameterList interaction. This makes the method satisfy the
368  // strong exception guarantee.
369 
370  // Get the user-supplied inverse diagonal.
371  //
372  // Check for a raw pointer (const V* or V*), for Ifpack
373  // compatibility, as well as for RCP<const V>, RCP<V>, const V, or
374  // V. We'll make a deep copy of the vector at the end of this
375  // method anyway, so its const-ness doesn't matter. We handle the
376  // latter two cases ("const V" or "V") specially (copy them into
377  // userInvDiagCopy first, which is otherwise null at the end of the
378  // long if-then chain) to avoid an extra copy.
379  if (plist.isParameter ("chebyshev: operator inv diagonal")) {
380  // Pointer to the user's Vector, if provided.
381  RCP<const V> userInvDiag;
382 
383  try { // Could the type be const V*?
384  const V* rawUserInvDiag =
385  plist.get<const V*> ("chebyshev: operator inv diagonal");
386  // Nonowning reference (we'll make a deep copy below)
387  userInvDiag = rcp (rawUserInvDiag, false);
388  } catch (Teuchos::Exceptions::InvalidParameterType&) {
389  }
390  if (userInvDiag.is_null ()) {
391  try { // Could the type be V*?
392  V* rawUserInvDiag = plist.get<V*> ("chebyshev: operator inv diagonal");
393  // Nonowning reference (we'll make a deep copy below)
394  userInvDiag = rcp (const_cast<const V*> (rawUserInvDiag), false);
395  } catch (Teuchos::Exceptions::InvalidParameterType&) {
396  }
397  }
398  if (userInvDiag.is_null ()) {
399  try { // Could the type be RCP<const V>?
400  userInvDiag =
401  plist.get<RCP<const V> > ("chebyshev: operator inv diagonal");
402  } catch (Teuchos::Exceptions::InvalidParameterType&) {
403  }
404  }
405  if (userInvDiag.is_null ()) {
406  try { // Could the type be RCP<V>?
407  RCP<V> userInvDiagNonConst =
408  plist.get<RCP<V> > ("chebyshev: operator inv diagonal");
409  userInvDiag = rcp_const_cast<const V> (userInvDiagNonConst);
410  } catch (Teuchos::Exceptions::InvalidParameterType&) {
411  }
412  }
413  if (userInvDiag.is_null ()) {
414 #ifndef _MSC_VER
415  try { // Could the type be const V?
416  // ParameterList::get() returns by reference. Thus, we don't
417  // have to invoke Vector's copy constructor here. It's good
418  // practice not to make an RCP to this reference, even though
419  // it should be valid as long as the ParameterList that holds
420  // it is valid. Thus, we make our deep copy here, rather than
421  // waiting to do it below.
422  const V& userInvDiagRef =
423  plist.get<const V> ("chebyshev: operator inv diagonal");
424  userInvDiagCopy = rcp (new V (userInvDiagRef, Teuchos::Copy));
425  // Tell the if-chain below not to keep trying.
426  userInvDiag = userInvDiagCopy;
427  } catch (Teuchos::Exceptions::InvalidParameterType&) {
428  }
429 #else
430  TEUCHOS_TEST_FOR_EXCEPTION(
431  true, std::runtime_error,
432  "Ifpack2::Chebyshev::setParameters: \"chebyshev: operator inv diagonal\" "
433  "plist.get<const V> does not compile around return held == other_held "
434  "in Teuchos::any in Visual Studio. Can't fix it now, so throwing "
435  "in case someone builds there.");
436 #endif
437  }
438  if (userInvDiag.is_null ()) {
439 #ifndef _MSC_VER
440  try { // Could the type be V?
441  // ParameterList::get() returns by reference. Thus, we don't
442  // have to invoke Vector's copy constructor here. It's good
443  // practice not to make an RCP to this reference, even though
444  // it should be valid as long as the ParameterList that holds
445  // it is valid. Thus, we make our deep copy here, rather than
446  // waiting to do it below.
447  V& userInvDiagNonConstRef =
448  plist.get<V> ("chebyshev: operator inv diagonal");
449  const V& userInvDiagRef = const_cast<const V&> (userInvDiagNonConstRef);
450  userInvDiagCopy = rcp (new V (userInvDiagRef, Teuchos::Copy));
451  // Tell the if-chain below not to keep trying.
452  userInvDiag = userInvDiagCopy;
453  } catch (Teuchos::Exceptions::InvalidParameterType&) {
454  }
455 #else
456  TEUCHOS_TEST_FOR_EXCEPTION(
457  true, std::runtime_error,
458  "Ifpack2::Chebyshev::setParameters: \"chebyshev: operator inv diagonal\" "
459  "plist.get<V> does not compile around return held == other_held "
460  "in Teuchos::any in Visual Studio. Can't fix it now, so throwing "
461  "in case someone builds there.");
462 #endif
463  }
464 
465  // NOTE: If the user's parameter has some strange type that we
466  // didn't test above, userInvDiag might still be null. You may
467  // want to add an error test for this condition. Currently, we
468  // just say in this case that the user didn't give us a Vector.
469 
470  // If we have userInvDiag but don't have a deep copy yet, make a
471  // deep copy now.
472  if (! userInvDiag.is_null () && userInvDiagCopy.is_null ()) {
473  userInvDiagCopy = rcp (new V (*userInvDiag, Teuchos::Copy));
474  }
475 
476  // NOTE: userInvDiag, if provided, is a row Map version of the
477  // Vector. We don't necessarily have a range Map yet. compute()
478  // would be the proper place to compute the range Map version of
479  // userInvDiag.
480  }
481 
482  // Don't fill in defaults for the max or min eigenvalue, because
483  // this class uses the existence of those parameters to determine
484  // whether it should do eigenanalysis.
485  if (plist.isParameter ("chebyshev: max eigenvalue")) {
486  if (plist.isType<double>("chebyshev: max eigenvalue"))
487  lambdaMax = plist.get<double> ("chebyshev: max eigenvalue");
488  else
489  lambdaMax = plist.get<ST> ("chebyshev: max eigenvalue");
490  TEUCHOS_TEST_FOR_EXCEPTION(
491  STS::isnaninf (lambdaMax), std::invalid_argument,
492  "Ifpack2::Chebyshev::setParameters: \"chebyshev: max eigenvalue\" "
493  "parameter is NaN or Inf. This parameter is optional, but if you "
494  "choose to supply it, it must have a finite value.");
495  }
496  if (plist.isParameter ("chebyshev: min eigenvalue")) {
497  if (plist.isType<double>("chebyshev: min eigenvalue"))
498  lambdaMin = plist.get<double> ("chebyshev: min eigenvalue");
499  else
500  lambdaMin = plist.get<ST> ("chebyshev: min eigenvalue");
501  TEUCHOS_TEST_FOR_EXCEPTION(
502  STS::isnaninf (lambdaMin), std::invalid_argument,
503  "Ifpack2::Chebyshev::setParameters: \"chebyshev: min eigenvalue\" "
504  "parameter is NaN or Inf. This parameter is optional, but if you "
505  "choose to supply it, it must have a finite value.");
506  }
507 
508  // Only fill in Ifpack2's name for the default parameter, not ML's.
509  if (plist.isParameter ("smoother: Chebyshev alpha")) { // ML compatibility
510  if (plist.isType<double>("smoother: Chebyshev alpha"))
511  eigRatio = plist.get<double> ("smoother: Chebyshev alpha");
512  else
513  eigRatio = plist.get<ST> ("smoother: Chebyshev alpha");
514  }
515  // Ifpack2's name overrides ML's name.
516  eigRatio = plist.get ("chebyshev: ratio eigenvalue", eigRatio);
517  TEUCHOS_TEST_FOR_EXCEPTION(
518  STS::isnaninf (eigRatio), std::invalid_argument,
519  "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\" "
520  "parameter (also called \"smoother: Chebyshev alpha\") is NaN or Inf. "
521  "This parameter is optional, but if you choose to supply it, it must have "
522  "a finite value.");
523  // mfh 11 Feb 2013: This class is currently only correct for real
524  // Scalar types, but we still want it to build for complex Scalar
525  // type so that users of Ifpack2::Factory can build their
526  // executables for real or complex Scalar type. Thus, we take the
527  // real parts here, which are always less-than comparable.
528  TEUCHOS_TEST_FOR_EXCEPTION(
529  STS::real (eigRatio) < STS::real (STS::one ()),
530  std::invalid_argument,
531  "Ifpack2::Chebyshev::setParameters: \"chebyshev: ratio eigenvalue\""
532  "parameter (also called \"smoother: Chebyshev alpha\") must be >= 1, "
533  "but you supplied the value " << eigRatio << ".");
534 
535  // Same name in Ifpack2 and Ifpack.
536  minDiagVal = plist.get ("chebyshev: min diagonal value", minDiagVal);
537  TEUCHOS_TEST_FOR_EXCEPTION(
538  STS::isnaninf (minDiagVal), std::invalid_argument,
539  "Ifpack2::Chebyshev::setParameters: \"chebyshev: min diagonal value\" "
540  "parameter is NaN or Inf. This parameter is optional, but if you choose "
541  "to supply it, it must have a finite value.");
542 
543  // Only fill in Ifpack2's name, not ML's or Ifpack's.
544  if (plist.isParameter ("smoother: sweeps")) { // ML compatibility
545  numIters = plist.get<int> ("smoother: sweeps");
546  } // Ifpack's name overrides ML's name.
547  if (plist.isParameter ("relaxation: sweeps")) { // Ifpack compatibility
548  numIters = plist.get<int> ("relaxation: sweeps");
549  } // Ifpack2's name overrides Ifpack's name.
550  numIters = plist.get ("chebyshev: degree", numIters);
551  TEUCHOS_TEST_FOR_EXCEPTION(
552  numIters < 0, std::invalid_argument,
553  "Ifpack2::Chebyshev::setParameters: \"chebyshev: degree\" parameter (also "
554  "called \"smoother: sweeps\" or \"relaxation: sweeps\") must be a "
555  "nonnegative integer. You gave a value of " << numIters << ".");
556 
557  // The last parameter name overrides the first.
558  if (plist.isParameter ("eigen-analysis: iterations")) { // ML compatibility
559  eigMaxIters = plist.get<int> ("eigen-analysis: iterations");
560  } // Ifpack2's name overrides ML's name.
561  eigMaxIters = plist.get ("chebyshev: eigenvalue max iterations", eigMaxIters);
562  TEUCHOS_TEST_FOR_EXCEPTION(
563  eigMaxIters < 0, std::invalid_argument,
564  "Ifpack2::Chebyshev::setParameters: \"chebyshev: eigenvalue max iterations"
565  "\" parameter (also called \"eigen-analysis: iterations\") must be a "
566  "nonnegative integer. You gave a value of " << eigMaxIters << ".");
567 
568  zeroStartingSolution = plist.get ("chebyshev: zero starting solution",
569  zeroStartingSolution);
570  assumeMatrixUnchanged = plist.get ("chebyshev: assume matrix does not change",
571  assumeMatrixUnchanged);
572 
573  // We don't want to fill these parameters in, because they shouldn't
574  // be visible to Ifpack2::Chebyshev users.
575  if (plist.isParameter ("chebyshev: textbook algorithm")) {
576  textbookAlgorithm = plist.get<bool> ("chebyshev: textbook algorithm");
577  }
578  if (plist.isParameter ("chebyshev: compute max residual norm")) {
579  computeMaxResNorm = plist.get<bool> ("chebyshev: compute max residual norm");
580  }
581 
582  // Test for Ifpack parameters that we won't ever implement here.
583  // Be careful to use the one-argument version of get(), since the
584  // two-argment version adds the parameter if it's not there.
585  TEUCHOS_TEST_FOR_EXCEPTION(
586  plist.isParameter ("chebyshev: use block mode") &&
587  ! plist.get<bool> ("chebyshev: use block mode"),
588  std::invalid_argument,
589  "Ifpack2::Chebyshev requires that if you supply the Ifpack parameter "
590  "\"chebyshev: use block mode\", it must be set to false. Ifpack2's "
591  "Chebyshev does not implement Ifpack's block mode.");
592  TEUCHOS_TEST_FOR_EXCEPTION(
593  plist.isParameter ("chebyshev: solve normal equations") &&
594  ! plist.get<bool> ("chebyshev: solve normal equations"),
595  std::invalid_argument,
596  "Ifpack2::Chebyshev does not and will never implement the Ifpack "
597  "parameter \"chebyshev: solve normal equations\". If you want to solve "
598  "the normal equations, construct a Tpetra::Operator that implements "
599  "A^* A, and use Chebyshev to solve A^* A x = A^* b.");
600 
601  // Test for Ifpack parameters that we haven't implemented yet.
602  //
603  // For now, we only check that this ML parameter, if provided, has
604  // the one value that we support. We consider other values "invalid
605  // arguments" rather than "logic errors," because Ifpack does not
606  // implement eigenanalyses other than the power method.
607  std::string eigenAnalysisType ("power-method");
608  if (plist.isParameter ("eigen-analysis: type")) {
609  eigenAnalysisType = plist.get<std::string> ("eigen-analysis: type");
610  TEUCHOS_TEST_FOR_EXCEPTION(
611  eigenAnalysisType != "power-method" &&
612  eigenAnalysisType != "power method",
613  std::invalid_argument,
614  "Ifpack2::Chebyshev: This class supports the ML parameter \"eigen-"
615  "analysis: type\" for backwards compatibility. However, Ifpack2 "
616  "currently only supports the \"power-method\" option for this "
617  "parameter. This imitates Ifpack, which only implements the power "
618  "method for eigenanalysis.");
619  }
620 
621  // We've validated all the parameters, so it's safe now to "commit" them.
622  userInvDiag_ = userInvDiagCopy;
623  userLambdaMax_ = lambdaMax;
624  userLambdaMin_ = lambdaMin;
625  userEigRatio_ = eigRatio;
626  minDiagVal_ = minDiagVal;
627  numIters_ = numIters;
628  eigMaxIters_ = eigMaxIters;
629  zeroStartingSolution_ = zeroStartingSolution;
630  assumeMatrixUnchanged_ = assumeMatrixUnchanged;
631  textbookAlgorithm_ = textbookAlgorithm;
632  computeMaxResNorm_ = computeMaxResNorm;
633 }
634 
635 
636 template<class ScalarType, class MV>
637 void
639 {
640  D_ = Teuchos::null;
641  diagOffsets_ = Teuchos::null;
642  savedDiagOffsets_ = false;
643  V_ = Teuchos::null;
644  W_ = Teuchos::null;
645  computedLambdaMax_ = STS::nan ();
646  computedLambdaMin_ = STS::nan ();
647 }
648 
649 
650 template<class ScalarType, class MV>
651 void
652 Chebyshev<ScalarType, MV>::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
653 {
654  if (A.getRawPtr () != A_.getRawPtr ()) {
655  if (! assumeMatrixUnchanged_) {
656  reset ();
657  }
658  A_ = A;
659  }
660 }
661 
662 
663 template<class ScalarType, class MV>
664 void
666 {
667  using std::endl;
668  // Some of the optimizations below only work if A_ is a
669  // Tpetra::CrsMatrix. We'll make our best guess about its type
670  // here, since we have no way to get back the original fifth
671  // template parameter.
672  typedef Tpetra::CrsMatrix<typename MV::scalar_type,
673  typename MV::local_ordinal_type,
674  typename MV::global_ordinal_type,
675  typename MV::node_type> crs_matrix_type;
676 
677  TEUCHOS_TEST_FOR_EXCEPTION(
678  A_.is_null (), std::runtime_error, "Ifpack2::Chebyshev::compute: The input "
679  "matrix A is null. Please call setMatrix() with a nonnull input matrix "
680  "before calling this method.");
681 
682  // If A_ is a CrsMatrix and its graph is constant, we presume that
683  // the user plans to reuse the structure of A_, but possibly change
684  // A_'s values before each compute() call. This is the intended use
685  // case for caching the offsets of the diagonal entries of A_, to
686  // speed up extraction of diagonal entries on subsequent compute()
687  // calls.
688 
689  // FIXME (mfh 22 Jan 2013, 10 Feb 2013) In all cases when we use
690  // isnaninf() in this method, we really only want to check if the
691  // number is NaN. Inf means something different. However,
692  // Teuchos::ScalarTraits doesn't distinguish the two cases.
693 
694  // makeInverseDiagonal() returns a range Map Vector.
695  if (userInvDiag_.is_null ()) {
696  Teuchos::RCP<const crs_matrix_type> A_crsMat =
697  Teuchos::rcp_dynamic_cast<const crs_matrix_type> (A_);
698 
699  if (D_.is_null ()) { // We haven't computed D_ before
700  if (! A_crsMat.is_null () && A_crsMat->isStaticGraph ()) {
701  // It's a CrsMatrix with a const graph; cache diagonal offsets.
702  A_crsMat->getLocalDiagOffsets (diagOffsets_);
703  savedDiagOffsets_ = true;
704  D_ = makeInverseDiagonal (*A_, true);
705  }
706  else { // either A_ is not a CrsMatrix, or its graph is nonconst
707  D_ = makeInverseDiagonal (*A_);
708  }
709  }
710  else if (! assumeMatrixUnchanged_) { // D_ exists but A_ may have changed
711  if (! A_crsMat.is_null () && A_crsMat->isStaticGraph ()) {
712  // It's a CrsMatrix with a const graph; cache diagonal offsets
713  // if we haven't already.
714  if (! savedDiagOffsets_) {
715  A_crsMat->getLocalDiagOffsets (diagOffsets_);
716  savedDiagOffsets_ = true;
717  }
718  // Now we're guaranteed to have cached diagonal offsets.
719  D_ = makeInverseDiagonal (*A_, true);
720  }
721  else { // either A_ is not a CrsMatrix, or its graph is nonconst
722  D_ = makeInverseDiagonal (*A_);
723  }
724  }
725  }
726  else { // the user provided an inverse diagonal
727  D_ = makeRangeMapVectorConst (userInvDiag_);
728  }
729 
730  // Have we estimated eigenvalues before?
731  const bool computedEigenvalueEstimates =
732  STS::isnaninf (computedLambdaMax_) || STS::isnaninf (computedLambdaMin_);
733 
734  // Only recompute the eigenvalue estimates if
735  // - we are supposed to assume that the matrix may have changed, or
736  // - they haven't been computed before, and the user hasn't given
737  // us at least an estimate of the max eigenvalue.
738  //
739  // We at least need an estimate of the max eigenvalue. This is the
740  // most important one if using Chebyshev as a smoother.
741  if (! assumeMatrixUnchanged_ ||
742  (! computedEigenvalueEstimates && STS::isnaninf (userLambdaMax_))) {
743  const ST computedLambdaMax = powerMethod (*A_, *D_, eigMaxIters_);
744  TEUCHOS_TEST_FOR_EXCEPTION(
745  STS::isnaninf (computedLambdaMax),
746  std::runtime_error,
747  "Ifpack2::Chebyshev::compute: Estimation of the max eigenvalue "
748  "of D^{-1} A failed, by producing Inf or NaN. This probably means that "
749  "the matrix contains Inf or NaN values, or that it is badly scaled.");
750  TEUCHOS_TEST_FOR_EXCEPTION(
751  STS::isnaninf (userEigRatio_),
752  std::logic_error,
753  "Ifpack2::Chebyshev::compute: userEigRatio_ is Inf or NaN."
754  << endl << "This should be impossible." << endl <<
755  "Please report this bug to the Ifpack2 developers.");
756 
757  // The power method doesn't estimate the min eigenvalue, so we
758  // do our best to provide an estimate. userEigRatio_ has a
759  // reasonable default value, and if the user provided it, we
760  // have already checked that its value is finite and >= 1.
761  const ST computedLambdaMin = computedLambdaMax / userEigRatio_;
762 
763  // Defer "committing" results until all computations succeeded.
764  computedLambdaMax_ = computedLambdaMax;
765  computedLambdaMin_ = computedLambdaMin;
766  } else {
767  TEUCHOS_TEST_FOR_EXCEPTION(
768  STS::isnaninf (userLambdaMax_) && STS::isnaninf (computedLambdaMax_),
769  std::logic_error,
770  "Ifpack2::Chebyshev::compute: " << endl <<
771  "Both userLambdaMax_ and computedLambdaMax_ are Inf or NaN."
772  << endl << "This should be impossible." << endl <<
773  "Please report this bug to the Ifpack2 developers.");
774  }
775 
777  // Figure out the eigenvalue estimates that apply() will use.
779 
780  // Always favor the user's max eigenvalue estimate, if provided.
781  if (STS::isnaninf (userLambdaMax_)) {
782  lambdaMaxForApply_ = computedLambdaMax_;
783  } else {
784  lambdaMaxForApply_ = userLambdaMax_;
785  }
786  // mfh 11 Feb 2013: For now, we imitate Ifpack by ignoring the
787  // user's min eigenvalue estimate, and using the given eigenvalue
788  // ratio to estimate the min eigenvalue. We could instead do this:
789  // favor the user's eigenvalue ratio estimate, but if it's not
790  // provided, use lambdaMax / lambdaMin. However, we want Chebyshev
791  // to have sensible smoother behavior if the user did not provide
792  // eigenvalue estimates. Ifpack's behavior attempts to push down
793  // the error terms associated with the largest eigenvalues, while
794  // expecting that users will only want a small number of iterations,
795  // so that error terms associated with the smallest eigenvalues
796  // won't grow too much. This is sensible behavior for a smoother.
797  lambdaMinForApply_ = lambdaMaxForApply_ / userEigRatio_;
798  eigRatioForApply_ = userEigRatio_;
799 
800  if (! textbookAlgorithm_) {
801  // Ifpack has a special-case modification of the eigenvalue bounds
802  // for the case where the max eigenvalue estimate is close to one.
803  const ST one = Teuchos::as<ST> (1);
804  // FIXME (mfh 20 Nov 2013) Should scale this 1.0e-6 term
805  // appropriately for MT's machine precision.
806  if (STS::magnitude (lambdaMaxForApply_ - one) < Teuchos::as<MT> (1.0e-6)) {
807  lambdaMinForApply_ = one;
808  lambdaMaxForApply_ = lambdaMinForApply_;
809  eigRatioForApply_ = one; // Ifpack doesn't include this line.
810  }
811  }
812 }
813 
814 
815 template<class ScalarType, class MV>
816 ScalarType
818 getLambdaMaxForApply() const {
819  return lambdaMaxForApply_;
820 }
821 
822 
823 template<class ScalarType, class MV>
826 {
827  TEUCHOS_TEST_FOR_EXCEPTION(
828  A_.is_null (), std::runtime_error, "Ifpack2::Chebyshev::apply: The input "
829  "matrix A is null. Please call setMatrix() with a nonnull input matrix, "
830  "and then call compute(), before calling this method.");
831  TEUCHOS_TEST_FOR_EXCEPTION(
832  STS::isnaninf (lambdaMaxForApply_), std::runtime_error,
833  "Ifpack2::Chebyshev::apply: There is no estimate for the max eigenvalue."
834  << std::endl << computeBeforeApplyReminder);
835  TEUCHOS_TEST_FOR_EXCEPTION(
836  STS::isnaninf (lambdaMinForApply_), std::runtime_error,
837  "Ifpack2::Chebyshev::apply: There is no estimate for the min eigenvalue."
838  << std::endl << computeBeforeApplyReminder);
839  TEUCHOS_TEST_FOR_EXCEPTION(
840  STS::isnaninf (eigRatioForApply_), std::runtime_error,
841  "Ifpack2::Chebyshev::apply: There is no estimate for the ratio of the max "
842  "eigenvalue to the min eigenvalue."
843  << std::endl << computeBeforeApplyReminder);
844  TEUCHOS_TEST_FOR_EXCEPTION(
845  D_.is_null (), std::runtime_error,
846  "Ifpack2::Chebyshev::apply: "
847  "The vector of inverse diagonal entries of the matrix has not yet been "
848  "computed." << std::endl << computeBeforeApplyReminder);
849 
850  if (textbookAlgorithm_) {
851  textbookApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
852  lambdaMinForApply_, eigRatioForApply_, *D_);
853  } else {
854  ifpackApplyImpl (*A_, B, X, numIters_, lambdaMaxForApply_,
855  lambdaMinForApply_, eigRatioForApply_, *D_);
856  }
857 
858  if (computeMaxResNorm_ && B.getNumVectors () > 0) {
859  MV R (B.getMap (), B.getNumVectors ());
860  computeResidual (R, B, *A_, X);
861  Teuchos::Array<MT> norms (B.getNumVectors ());
862  R.norm2 (norms ());
863  return *std::max_element (norms.begin (), norms.end ());
864  } else {
865  return Teuchos::ScalarTraits<MT>::zero ();
866  }
867 }
868 
869 template<class ScalarType, class MV>
870 void
872 print (std::ostream& out) {
873  this->describe (* (Teuchos::getFancyOStream (Teuchos::rcpFromRef (out))),
874  Teuchos::VERB_MEDIUM);
875 }
876 
877 template<class ScalarType, class MV>
878 void
880 computeResidual (MV& R, const MV& B, const op_type& A, const MV& X,
881  const Teuchos::ETransp mode)
882 {
883  Tpetra::deep_copy(R, B);
884  A.apply (X, R, mode, -STS::one(), STS::one());
885 }
886 
887 template<class ScalarType, class MV>
888 void
890 solve (MV& Z, const V& D_inv, const MV& R) {
891  Z.elementWiseMultiply (STS::one(), D_inv, R, STS::zero());
892 }
893 
894 template<class ScalarType, class MV>
895 void
897 solve (MV& Z, const ST alpha, const V& D_inv, const MV& R) {
898  Z.elementWiseMultiply (alpha, D_inv, R, STS::zero());
899 }
900 
901 template<class ScalarType, class MV>
902 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
904 makeInverseDiagonal (const row_matrix_type& A, const bool useDiagOffsets) const
905 {
906  using Teuchos::RCP;
907  using Teuchos::rcpFromRef;
908  using Teuchos::rcp_dynamic_cast;
909 
910  RCP<V> D_rowMap (new V (A.getGraph ()->getRowMap ()));
911  if (useDiagOffsets) {
912  // The optimizations below only work if A_ is a Tpetra::CrsMatrix.
913  // We'll make our best guess about its type here, since we have no
914  // way to get back the original fifth template parameter.
915  typedef Tpetra::CrsMatrix<typename MV::scalar_type,
916  typename MV::local_ordinal_type,
917  typename MV::global_ordinal_type,
918  typename MV::node_type> crs_matrix_type;
919  RCP<const crs_matrix_type> A_crsMat =
920  rcp_dynamic_cast<const crs_matrix_type> (rcpFromRef (A));
921  if (! A_crsMat.is_null ()) {
922  TEUCHOS_TEST_FOR_EXCEPTION(
923  ! savedDiagOffsets_, std::logic_error,
924  "Ifpack2::Details::Chebyshev::makeInverseDiagonal: "
925  "It is not allowed to call this method with useDiagOffsets=true, "
926  "if you have not previously saved offsets of diagonal entries. "
927  "This situation should never arise if this class is used properly. "
928  "Please report this bug to the Ifpack2 developers.");
929  A_crsMat->getLocalDiagCopy (*D_rowMap, diagOffsets_ ());
930  }
931  }
932  else {
933  // This always works for a Tpetra::RowMatrix, even if it is not a
934  // Tpetra::CrsMatrix. We just don't have offsets in this case.
935  A.getLocalDiagCopy (*D_rowMap);
936  }
937  RCP<V> D_rangeMap = makeRangeMapVector (D_rowMap);
938 
939  // Invert the diagonal entries, replacing entries less (in
940  // magnitude) than the user-specified value with that value.
941  reciprocal_threshold (*D_rangeMap, minDiagVal_);
942  return Teuchos::rcp_const_cast<const V> (D_rangeMap);
943 }
944 
945 
946 template<class ScalarType, class MV>
947 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::V>
949 makeRangeMapVectorConst (const Teuchos::RCP<const V>& D) const
950 {
951  using Teuchos::RCP;
952  using Teuchos::rcp;
953  typedef Tpetra::Export<typename MV::local_ordinal_type,
954  typename MV::global_ordinal_type,
955  typename MV::node_type> export_type;
956  // This throws logic_error instead of runtime_error, because the
957  // methods that call makeRangeMapVector should all have checked
958  // whether A_ is null before calling this method.
959  TEUCHOS_TEST_FOR_EXCEPTION(
960  A_.is_null (), std::logic_error, "Ifpack2::Details::Chebyshev::"
961  "makeRangeMapVector: The input matrix A is null. Please call setMatrix() "
962  "with a nonnull input matrix before calling this method. This is probably "
963  "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
964  TEUCHOS_TEST_FOR_EXCEPTION(
965  D.is_null (), std::logic_error, "Ifpack2::Details::Chebyshev::"
966  "makeRangeMapVector: The input Vector D is null. This is probably "
967  "a bug in Ifpack2; please report this bug to the Ifpack2 developers.");
968 
969  RCP<const map_type> sourceMap = D->getMap ();
970  RCP<const map_type> rangeMap = A_->getRangeMap ();
971  RCP<const map_type> rowMap = A_->getRowMap ();
972 
973  if (rangeMap->isSameAs (*sourceMap)) {
974  // The given vector's Map is the same as the matrix's range Map.
975  // That means we don't need to Export. This is the fast path.
976  return D;
977  }
978  else { // We need to Export.
979  RCP<const export_type> exporter;
980  // Making an Export object from scratch is expensive enough that
981  // it's worth the O(1) global reductions to call isSameAs(), to
982  // see if we can avoid that cost.
983  if (sourceMap->isSameAs (*rowMap)) {
984  // We can reuse the matrix's Export object, if there is one.
985  exporter = A_->getGraph ()->getExporter ();
986  }
987  else { // We have to make a new Export object.
988  exporter = rcp (new export_type (sourceMap, rangeMap));
989  }
990 
991  if (exporter.is_null ()) {
992  return D; // Row Map and range Map are the same; no need to Export.
993  }
994  else { // Row Map and range Map are _not_ the same; must Export.
995  RCP<V> D_out = rcp (new V (*D, Teuchos::Copy));
996  D_out->doExport (*D, *exporter, Tpetra::ADD);
997  return Teuchos::rcp_const_cast<const V> (D_out);
998  }
999  }
1000 }
1001 
1002 
1003 template<class ScalarType, class MV>
1004 Teuchos::RCP<typename Chebyshev<ScalarType, MV>::V>
1006 makeRangeMapVector (const Teuchos::RCP<V>& D) const
1007 {
1008  using Teuchos::rcp_const_cast;
1009  return rcp_const_cast<V> (makeRangeMapVectorConst (rcp_const_cast<V> (D)));
1010 }
1011 
1012 
1013 template<class ScalarType, class MV>
1014 void
1016 textbookApplyImpl (const op_type& A,
1017  const MV& B,
1018  MV& X,
1019  const int numIters,
1020  const ST lambdaMax,
1021  const ST lambdaMin,
1022  const ST eigRatio,
1023  const V& D_inv) const
1024 {
1025  (void) lambdaMin; // Forestall compiler warning.
1026  const ST myLambdaMin = lambdaMax / eigRatio;
1027 
1028  const ST zero = Teuchos::as<ST> (0);
1029  const ST one = Teuchos::as<ST> (1);
1030  const ST two = Teuchos::as<ST> (2);
1031  const ST d = (lambdaMax + myLambdaMin) / two; // Ifpack2 calls this theta
1032  const ST c = (lambdaMax - myLambdaMin) / two; // Ifpack2 calls this 1/delta
1033 
1034  if (zeroStartingSolution_ && numIters > 0) {
1035  // If zero iterations, then input X is output X.
1036  X.putScalar (zero);
1037  }
1038  MV R (B.getMap (), B.getNumVectors (), false);
1039  MV P (B.getMap (), B.getNumVectors (), false);
1040  MV Z (B.getMap (), B.getNumVectors (), false);
1041  ST alpha, beta;
1042  for (int i = 0; i < numIters; ++i) {
1043  computeResidual (R, B, A, X); // R = B - A*X
1044  solve (Z, D_inv, R); // z = D_inv * R, that is, D \ R.
1045  if (i == 0) {
1046  P = Z;
1047  alpha = two / d;
1048  } else {
1049  //beta = (c * alpha / two)^2;
1050  //const ST sqrtBeta = c * alpha / two;
1051  //beta = sqrtBeta * sqrtBeta;
1052  beta = alpha * (c/two) * (c/two);
1053  alpha = one / (d - beta);
1054  P.update (one, Z, beta); // P = Z + beta*P
1055  }
1056  X.update (alpha, P, one); // X = X + alpha*P
1057  // If we compute the residual here, we could either do R = B -
1058  // A*X, or R = R - alpha*A*P. Since we choose the former, we
1059  // can move the computeResidual call to the top of the loop.
1060  }
1061 }
1062 
1063 template<class ScalarType, class MV>
1066  Teuchos::Array<MT> norms (X.getNumVectors ());
1067  X.normInf (norms());
1068  return *std::max_element (norms.begin (), norms.end ());
1069 }
1070 
1071 template<class ScalarType, class MV>
1072 void
1074 ifpackApplyImpl (const op_type& A,
1075  const MV& B,
1076  MV& X,
1077  const int numIters,
1078  const ST lambdaMax,
1079  const ST lambdaMin,
1080  const ST eigRatio,
1081  const V& D_inv)
1082 {
1083 #ifdef HAVE_TEUCHOS_DEBUG
1084 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1085  using std::cerr;
1086  using std::endl;
1087  cerr << "\\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1088  cerr << "\\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1089 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1090 #endif // HAVE_TEUCHOS_DEBUG
1091 
1092  if (numIters <= 0) {
1093  return;
1094  }
1095  const ST one = Teuchos::as<ST> (1);
1096  const ST two = Teuchos::as<ST> (2);
1097 
1098  // Quick solve when the matrix A is the identity.
1099  if (lambdaMin == one && lambdaMax == lambdaMin) {
1100  solve (X, D_inv, B);
1101  return;
1102  }
1103 
1104  // Initialize coefficients
1105  const ST alpha = lambdaMax / eigRatio;
1106  const ST beta = Teuchos::as<ST> (1.1) * lambdaMax;
1107  const ST delta = two / (beta - alpha);
1108  const ST theta = (beta + alpha) / two;
1109  const ST s1 = theta * delta;
1110 
1111 #ifdef HAVE_TEUCHOS_DEBUG
1112 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1113  cerr << "alpha = " << alpha << endl
1114  << "beta = " << beta << endl
1115  << "delta = " << delta << endl
1116  << "theta = " << theta << endl
1117  << "s1 = " << s1 << endl;
1118 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1119 #endif // HAVE_TEUCHOS_DEBUG
1120 
1121  // Fetch cached temporary vectors.
1122  Teuchos::RCP<MV> V_ptr, W_ptr;
1123  makeTempMultiVectors (V_ptr, W_ptr, B);
1124 
1125  // mfh 28 Jan 2013: We write V1 instead of V, so as not to confuse
1126  // the multivector V with the typedef V (for Tpetra::Vector).
1127  //MV V1 (B.getMap (), B.getNumVectors (), false);
1128  //MV W (B.getMap (), B.getNumVectors (), false);
1129  MV& V1 = *V_ptr;
1130  MV& W = *W_ptr;
1131 
1132 #ifdef HAVE_TEUCHOS_DEBUG
1133 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1134  cerr << "Iteration " << 1 << ":" << endl
1135  << "- \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1136 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1137 #endif // HAVE_TEUCHOS_DEBUG
1138 
1139  // Special case for the first iteration.
1140  if (! zeroStartingSolution_) {
1141  computeResidual (V1, B, A, X); // V1 = B - A*X
1142 
1143 #ifdef HAVE_TEUCHOS_DEBUG
1144 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1145  cerr << "- \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1146 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1147 #endif // HAVE_TEUCHOS_DEBUG
1148 
1149  solve (W, one/theta, D_inv, V1); // W = (1/theta)*D_inv*(B-A*X)
1150 
1151 #ifdef HAVE_TEUCHOS_DEBUG
1152 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1153  cerr << "- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1154 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1155 #endif // HAVE_TEUCHOS_DEBUG
1156 
1157  X.update (one, W, one); // X = X + W
1158  } else {
1159  solve (W, one/theta, D_inv, B); // W = (1/theta)*D_inv*B
1160 
1161 #ifdef HAVE_TEUCHOS_DEBUG
1162 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1163  cerr << "- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1164 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1165 #endif // HAVE_TEUCHOS_DEBUG
1166 
1167  Tpetra::deep_copy(X, W); // X = 0 + W
1168  }
1169 #ifdef HAVE_TEUCHOS_DEBUG
1170 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1171  cerr << "- \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1172 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1173 #endif // HAVE_TEUCHOS_DEBUG
1174 
1175  // The rest of the iterations.
1176  ST rhok = one / s1;
1177  ST rhokp1, dtemp1, dtemp2;
1178  for (int deg = 1; deg < numIters; ++deg) {
1179 
1180 #ifdef HAVE_TEUCHOS_DEBUG
1181 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1182  cerr << "Iteration " << deg+1 << ":" << endl;
1183  cerr << "- \\|D\\|_{\\infty} = " << D_->normInf () << endl;
1184  cerr << "- \\|B\\|_{\\infty} = " << maxNormInf (B) << endl;
1185  cerr << "- \\|A\\|_{\\text{frob}} = " << A_->getFrobeniusNorm () << endl;
1186  cerr << "- rhok = " << rhok << endl;
1187  V1.putScalar (STS::zero ());
1188 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1189 #endif // HAVE_TEUCHOS_DEBUG
1190 
1191  computeResidual (V1, B, A, X); // V1 = B - A*X
1192 
1193 #ifdef HAVE_TEUCHOS_DEBUG
1194 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1195  cerr << "- \\|B - A*X\\|_{\\infty} = " << maxNormInf (V1) << endl;
1196 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1197 #endif // HAVE_TEUCHOS_DEBUG
1198 
1199  rhokp1 = one / (two * s1 - rhok);
1200  dtemp1 = rhokp1 * rhok;
1201  dtemp2 = two * rhokp1 * delta;
1202  rhok = rhokp1;
1203 
1204 #ifdef HAVE_TEUCHOS_DEBUG
1205 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1206  cerr << "- dtemp1 = " << dtemp1 << endl
1207  << "- dtemp2 = " << dtemp2 << endl;
1208 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1209 #endif // HAVE_TEUCHOS_DEBUG
1210 
1211  W.scale (dtemp1);
1212  W.elementWiseMultiply (dtemp2, D_inv, V1, one);
1213  X.update (one, W, one);
1214 
1215 #ifdef HAVE_TEUCHOS_DEBUG
1216 #ifdef IFPACK_DETAILS_CHEBYSHEV_DEBUG
1217  cerr << "- \\|W\\|_{\\infty} = " << maxNormInf (W) << endl;
1218  cerr << "- \\|X\\|_{\\infty} = " << maxNormInf (X) << endl;
1219 #endif // IFPACK_DETAILS_CHEBYSHEV_DEBUG
1220 #endif // HAVE_TEUCHOS_DEBUG
1221  }
1222 }
1223 
1224 template<class ScalarType, class MV>
1228  const V& D_inv,
1229  const int numIters,
1230  V& x)
1231 {
1232  const ST zero = Teuchos::as<ST> (0);
1233  const ST one = Teuchos::as<ST> (1);
1234  ST lambdaMax = zero;
1235  ST RQ_top, RQ_bottom, norm;
1236 
1237  V y (A.getRangeMap ());
1238  norm = x.norm2 ();
1239  TEUCHOS_TEST_FOR_EXCEPTION(
1240  norm == zero, std::runtime_error,
1241  "Ifpack2::Chebyshev::powerMethodWithInitGuess: "
1242  "The initial guess has zero norm. This could be either because Tpetra::"
1243  "Vector::randomize() filled the vector with zeros (if that was used to "
1244  "compute the initial guess), or because the norm2 method has a bug. The "
1245  "first is not impossible, but unlikely.");
1246 
1247  //std::cerr << "Original norm1(x): " << x.norm1 () << ", norm2(x): " << norm << std::endl;
1248  x.scale (one / norm);
1249  //std::cerr << "norm1(x.scale(one/norm)): " << x.norm1 () << std::endl;
1250  for (int iter = 0; iter < numIters; ++iter) {
1251  A.apply (x, y);
1252  //std::cerr << "norm1(y) before: " << y.norm1 ();
1253  solve (y, D_inv, y);
1254  //std::cerr << ", norm1(y) after: " << y.norm1 () << std::endl;
1255  RQ_top = y.dot (x);
1256  RQ_bottom = x.dot (x);
1257  //std::cerr << "RQ_top: " << RQ_top << ", RQ_bottom: " << RQ_bottom << std::endl;
1258  lambdaMax = RQ_top / RQ_bottom;
1259  norm = y.norm2 ();
1260  if (norm == zero) { // Return something reasonable.
1261  return zero;
1262  }
1263  x.update (one / norm, y, zero);
1264  }
1265  return lambdaMax;
1266 }
1267 
1268 template<class ScalarType, class MV>
1271 powerMethod (const op_type& A, const V& D_inv, const int numIters)
1272 {
1273  const ST zero = Teuchos::as<ST> (0);
1274  V x (A.getDomainMap ());
1275  x.randomize ();
1276 
1277  ST lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1278  // mfh 07 Jan 2015: Taking the real part here is only a concession
1279  // to the compiler, so that this class can build with ScalarType =
1280  // std::complex<T>. Our Chebyshev implementation only works with
1281  // real, symmetric positive definite matrices. The right thing to
1282  // do would be what Belos does, which is provide a partial
1283  // specialization for ScalarType = std::complex<T> with a stub
1284  // implementation (that builds, but whose constructor throws).
1285  if (STS::real (lambdaMax) < STS::real (zero)) {
1286  // Max eigenvalue estimate is negative. Perhaps we got unlucky
1287  // with the random initial guess. Try again with a different (but
1288  // still random) initial guess. Only try again once, so that the
1289  // run time is bounded.
1290  x.randomize ();
1291  lambdaMax = powerMethodWithInitGuess (A, D_inv, numIters, x);
1292  }
1293  return lambdaMax;
1294 }
1295 
1296 template<class ScalarType, class MV>
1297 Teuchos::RCP<const typename Chebyshev<ScalarType, MV>::row_matrix_type>
1299  return A_;
1300 }
1301 
1302 template<class ScalarType, class MV>
1303 bool
1306  // Technically, this is true, because the matrix must be symmetric.
1307  return true;
1308 }
1309 
1310 template<class ScalarType, class MV>
1311 void
1313 makeTempMultiVectors (Teuchos::RCP<MV>& V1,
1314  Teuchos::RCP<MV>& W,
1315  const MV& X)
1316 {
1317  if (V_.is_null ()) {
1318  V_ = Teuchos::rcp (new MV (X.getMap (), X.getNumVectors (), false));
1319  }
1320  //W must be initialized to zero when it is used as a multigrid smoother.
1321  if (W_.is_null ()) {
1322  W_ = Teuchos::rcp (new MV (X.getMap (), X.getNumVectors (), true));
1323  }
1324  V1 = V_;
1325  W = W_;
1326 }
1327 
1328 template<class ScalarType, class MV>
1329 std::string
1331 description () const {
1332  std::ostringstream oss;
1333  // YAML requires quoting the key in this case, to distinguish
1334  // key's colons from the colon that separates key from value.
1335  oss << "\"Ifpack2::Details::Chebyshev\":"
1336  << "{"
1337  << "degree: " << numIters_
1338  << ", lambdaMax: " << lambdaMaxForApply_
1339  << ", alpha: " << eigRatioForApply_
1340  << ", lambdaMin: " << lambdaMinForApply_
1341  << "}";
1342  return oss.str();
1343 }
1344 
1345 template<class ScalarType, class MV>
1346 void
1348 describe (Teuchos::FancyOStream& out,
1349  const Teuchos::EVerbosityLevel verbLevel) const
1350 {
1351  using Teuchos::TypeNameTraits;
1352  using std::endl;
1353 
1354  const Teuchos::EVerbosityLevel vl =
1355  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
1356  if (vl > Teuchos::VERB_NONE) {
1357  if (vl == Teuchos::VERB_LOW) {
1358  out << description () << endl;
1359  } else { // vl > Teuchos::VERB_LOW
1360  // YAML requires quoting the key in this case, to distinguish
1361  // key's colons from the colon that separates key from value.
1362  out << "\"Ifpack2::Details::Chebyshev\":" << endl;
1363  Teuchos::OSTab tab1 (Teuchos::rcpFromRef (out));
1364  out << "Template parameters:" << endl;
1365  {
1366  Teuchos::OSTab tab2 (Teuchos::rcpFromRef (out));
1367  out << "ScalarType: \"" << TypeNameTraits<ScalarType>::name () << "\"" << endl
1368  << "MV: \"" << TypeNameTraits<MV>::name () << "\"" << endl;
1369  }
1370  // "Computed parameters" literally means "parameters whose
1371  // values were computed by compute()."
1372  out << endl << "Computed parameters:" << endl;
1373  {
1374  Teuchos::OSTab tab2 (Teuchos::rcpFromRef (out));
1375  // Users might want to see the values in the computed inverse
1376  // diagonal, so we print them out at the highest verbosity.
1377  out << "D_: ";
1378  if (D_.is_null ()) {
1379  out << "unset" << endl;
1380  } else if (vl <= Teuchos::VERB_HIGH) {
1381  out << "set" << endl;
1382  } else { // D_ not null and vl > Teuchos::VERB_HIGH
1383  out << endl;
1384  {
1385  Teuchos::OSTab tab3 (Teuchos::rcpFromRef (out));
1386  D_->describe (out, vl);
1387  }
1388  out << endl;
1389  }
1390  // V_ and W_ are scratch space; their values are irrelevant.
1391  // All that matters is whether or not they have been set.
1392  out << "V_: " << (V_.is_null () ? "unset" : "set") << endl
1393  << "W_: " << (W_.is_null () ? "unset" : "set") << endl
1394  << "computedLambdaMax_: " << computedLambdaMax_ << endl
1395  << "computedLambdaMin_: " << computedLambdaMin_ << endl
1396  << "lambdaMaxForApply_: " << lambdaMaxForApply_ << endl
1397  << "lambdaMinForApply_: " << lambdaMinForApply_ << endl
1398  << "eigRatioForApply_: " << eigRatioForApply_ << endl;
1399  }
1400  out << "User parameters:" << endl;
1401  {
1402  Teuchos::OSTab tab2 (Teuchos::rcpFromRef (out));
1403  out << "userInvDiag_: ";
1404  if (userInvDiag_.is_null ()) {
1405  out << "unset" << endl;
1406  } else if (vl <= Teuchos::VERB_HIGH) {
1407  out << "set" << endl;
1408  } else { // userInvDiag_ not null and vl > Teuchos::VERB_HIGH
1409  out << endl;
1410  {
1411  Teuchos::OSTab tab3 (Teuchos::rcpFromRef (out));
1412  userInvDiag_->describe (out, vl);
1413  }
1414  out << endl;
1415  }
1416  out << "userLambdaMax_: " << userLambdaMax_ << endl
1417  << "userLambdaMin_: " << userLambdaMin_ << endl
1418  << "userEigRatio_: " << userEigRatio_ << endl
1419  << "numIters_: " << numIters_ << endl
1420  << "eigMaxIters_: " << eigMaxIters_ << endl
1421  << "zeroStartingSolution_: " << zeroStartingSolution_ << endl
1422  << "assumeMatrixUnchanged_: " << assumeMatrixUnchanged_ << endl
1423  << "textbookAlgorithm_: " << textbookAlgorithm_ << endl
1424  << "computeMaxResNorm_: " << computeMaxResNorm_ << endl;
1425  }
1426  out << endl;
1427  }
1428  }
1429 }
1430 
1431 } // namespace Details
1432 } // namespace Ifpack2
1433 
1434 #define IFPACK2_DETAILS_CHEBYSHEV_INSTANT(S,LO,GO,N) \
1435  template class Ifpack2::Details::Chebyshev< S, Tpetra::MultiVector<S, LO, GO, N> >;
1436 
1437 #endif // IFPACK2_DETAILS_CHEBYSHEV_DEF_HPP
void setParameters(Teuchos::ParameterList &plist)
Set (or reset) parameters.
Definition: Ifpack2_Details_Chebyshev_def.hpp:316
MT apply(const MV &B, MV &X)
Solve Ax=b for x with Chebyshev iteration with left diagonal scaling.
Definition: Ifpack2_Details_Chebyshev_def.hpp:825
Tpetra::Operator< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > op_type
Specialization of Tpetra::Operator.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:142
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print a description of the Chebyshev solver to out.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1348
void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Set the matrix.
Definition: Ifpack2_Details_Chebyshev_def.hpp:652
Teuchos::RCP< const row_matrix_type > getMatrix() const
Get the matrix given to the constructor.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1298
Chebyshev(Teuchos::RCP< const row_matrix_type > A)
Definition: Ifpack2_Details_Chebyshev_def.hpp:266
Ifpack2 implementation details.
void print(std::ostream &out)
Print instance data to the given output stream.
Definition: Ifpack2_Details_Chebyshev_def.hpp:872
STS::magnitudeType MT
The type of the absolute value of a ScalarType.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:137
ScalarType ST
The type of entries in the matrix and vectors.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:133
Tpetra::RowMatrix< typename MV::scalar_type, typename MV::local_ordinal_type, typename MV::global_ordinal_type, typename MV::node_type > row_matrix_type
Specialization of Tpetra::RowMatrix.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:147
Left-scaled Chebyshev iteration.
Definition: Ifpack2_Details_Chebyshev_decl.hpp:127
Declaration of Chebyshev implementation.
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:63
std::string description() const
A single-line description of the Chebyshev solver.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1331
bool hasTransposeApply() const
Whether it&#39;s possible to apply the transpose of this operator.
Definition: Ifpack2_Details_Chebyshev_def.hpp:1305
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
void compute()
(Re)compute the left scaling D_inv, and estimate min and max eigenvalues of D_inv * A...
Definition: Ifpack2_Details_Chebyshev_def.hpp:665