Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_Chebyshev_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Templated Object-Oriented Algebraic Preconditioner Package
5 // Copyright (2009) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
38 //
39 // ***********************************************************************
40 //@HEADER
41 */
42 
43 #ifndef IFPACK2_CHEBYSHEV_DEF_HPP
44 #define IFPACK2_CHEBYSHEV_DEF_HPP
45 
46 #include "Ifpack2_Parameters.hpp"
47 #include "Teuchos_TimeMonitor.hpp"
48 #include "Tpetra_CrsMatrix.hpp"
49 #include "Teuchos_TypeNameTraits.hpp"
50 #include <iostream>
51 #include <sstream>
52 
53 
54 namespace Ifpack2 {
55 
56 template<class MatrixType>
58 Chebyshev (const Teuchos::RCP<const row_matrix_type>& A)
59  : impl_ (A),
60  IsInitialized_ (false),
61  IsComputed_ (false),
62  NumInitialize_ (0),
63  NumCompute_ (0),
64  NumApply_ (0),
65  InitializeTime_ (0.0),
66  ComputeTime_ (0.0),
67  ApplyTime_ (0.0),
68  ComputeFlops_ (0.0),
69  ApplyFlops_ (0.0)
70 {
71  this->setObjectLabel ("Ifpack2::Chebyshev");
72 }
73 
74 
75 template<class MatrixType>
77 }
78 
79 
80 template<class MatrixType>
81 void Chebyshev<MatrixType>::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
82 {
83  if (A.getRawPtr () != impl_.getMatrix ().getRawPtr ()) {
84  IsInitialized_ = false;
85  IsComputed_ = false;
86  impl_.setMatrix (A);
87  }
88 }
89 
90 
91 template<class MatrixType>
92 void
93 Chebyshev<MatrixType>::setParameters (const Teuchos::ParameterList& List)
94 {
95  // FIXME (mfh 25 Jan 2013) Casting away const is bad here.
96  impl_.setParameters (const_cast<Teuchos::ParameterList&> (List));
97 }
98 
99 
100 template<class MatrixType>
101 void
103 {
104  impl_.setZeroStartingSolution(zeroStartingSolution);
105 }
106 
107 template<class MatrixType>
108 Teuchos::RCP<const Teuchos::Comm<int> >
110 {
111  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
112  TEUCHOS_TEST_FOR_EXCEPTION(
113  A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::getComm: The input "
114  "matrix A is null. Please call setMatrix() with a nonnull input matrix "
115  "before calling this method.");
116  return A->getRowMap ()->getComm ();
117 }
118 
119 
120 template<class MatrixType>
121 Teuchos::RCP<const typename Chebyshev<MatrixType>::row_matrix_type>
123 getMatrix() const {
124  return impl_.getMatrix ();
125 }
126 
127 
128 template<class MatrixType>
129 Teuchos::RCP<const Tpetra::CrsMatrix<typename MatrixType::scalar_type,
130  typename MatrixType::local_ordinal_type,
131  typename MatrixType::global_ordinal_type,
132  typename MatrixType::node_type> >
134 getCrsMatrix() const {
135  typedef Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
136  global_ordinal_type, node_type> crs_matrix_type;
137  return Teuchos::rcp_dynamic_cast<const crs_matrix_type> (impl_.getMatrix ());
138 }
139 
140 
141 template<class MatrixType>
142 Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>
144 getDomainMap () const
145 {
146  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
147  TEUCHOS_TEST_FOR_EXCEPTION(
148  A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::getDomainMap: The "
149  "input matrix A is null. Please call setMatrix() with a nonnull input "
150  "matrix before calling this method.");
151  return A->getDomainMap ();
152 }
153 
154 
155 template<class MatrixType>
156 Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>
158 getRangeMap () const
159 {
160  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
161  TEUCHOS_TEST_FOR_EXCEPTION(
162  A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::getRangeMap: The "
163  "input matrix A is null. Please call setMatrix() with a nonnull input "
164  "matrix before calling this method.");
165  return A->getRangeMap ();
166 }
167 
168 
169 template<class MatrixType>
171  return impl_.hasTransposeApply ();
172 }
173 
174 
175 template<class MatrixType>
177  return NumInitialize_;
178 }
179 
180 
181 template<class MatrixType>
183  return NumCompute_;
184 }
185 
186 
187 template<class MatrixType>
189  return NumApply_;
190 }
191 
192 
193 template<class MatrixType>
195  return InitializeTime_;
196 }
197 
198 
199 template<class MatrixType>
201  return ComputeTime_;
202 }
203 
204 
205 template<class MatrixType>
207  return ApplyTime_;
208 }
209 
210 
211 template<class MatrixType>
213  return ComputeFlops_;
214 }
215 
216 
217 template<class MatrixType>
219  return ApplyFlops_;
220 }
221 
222 template<class MatrixType>
224  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix();
225  TEUCHOS_TEST_FOR_EXCEPTION(
226  A.is_null (), std::runtime_error, "Ifpack2::Chevyshev::getNodeSmootherComplexity: "
227  "The input matrix A is null. Please call setMatrix() with a nonnull "
228  "input matrix, then call compute(), before calling this method.");
229  // Chevyshev costs roughly one apply + one diagonal inverse per iteration
230  return A->getNodeNumRows() + A->getNodeNumEntries();
231 }
232 
233 
234 
235 template<class MatrixType>
236 void
238 apply (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
239  Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
240  Teuchos::ETransp mode,
241  scalar_type alpha,
242  scalar_type beta) const
243 {
244  const std::string timerName ("Ifpack2::Chebyshev::apply");
245  Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
246  if (timer.is_null ()) {
247  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
248  }
249 
250  double startTime = timer->wallTime();
251 
252  // Start timing here.
253  {
254  Teuchos::TimeMonitor timeMon (*timer);
255 
256  // compute() calls initialize() if it hasn't already been called.
257  // Thus, we only need to check isComputed().
258  TEUCHOS_TEST_FOR_EXCEPTION(
259  ! isComputed (), std::runtime_error,
260  "Ifpack2::Chebyshev::apply(): You must call the compute() method before "
261  "you may call apply().");
262  TEUCHOS_TEST_FOR_EXCEPTION(
263  X.getNumVectors () != Y.getNumVectors (), std::runtime_error,
264  "Ifpack2::Chebyshev::apply(): X and Y must have the same number of "
265  "columns. X.getNumVectors() = " << X.getNumVectors() << " != "
266  << "Y.getNumVectors() = " << Y.getNumVectors() << ".");
267  applyImpl (X, Y, mode, alpha, beta);
268  }
269  ++NumApply_;
270  ApplyTime_ += (timer->wallTime() - startTime);
271 }
272 
273 
274 template<class MatrixType>
275 void
277 applyMat (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
278  Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
279  Teuchos::ETransp mode) const
280 {
281  TEUCHOS_TEST_FOR_EXCEPTION(
282  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
283  "Ifpack2::Chebyshev::applyMat: X.getNumVectors() != Y.getNumVectors().");
284 
285  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
286  TEUCHOS_TEST_FOR_EXCEPTION(
287  A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::applyMat: The input "
288  "matrix A is null. Please call setMatrix() with a nonnull input matrix "
289  "before calling this method.");
290 
291  A->apply (X, Y, mode);
292 }
293 
294 
295 template<class MatrixType>
297  // We create the timer, but this method doesn't do anything, so
298  // there is no need to start the timer. The resulting total time
299  // will always be zero.
300  const std::string timerName ("Ifpack2::Chebyshev::initialize");
301  Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
302  if (timer.is_null ()) {
303  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
304  }
305  IsInitialized_ = true;
306  ++NumInitialize_;
307 }
308 
309 
310 template<class MatrixType>
312 {
313  const std::string timerName ("Ifpack2::Chebyshev::compute");
314  Teuchos::RCP<Teuchos::Time> timer = Teuchos::TimeMonitor::lookupCounter (timerName);
315  if (timer.is_null ()) {
316  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
317  }
318 
319  double startTime = timer->wallTime();
320 
321  // Start timing here.
322  {
323  Teuchos::TimeMonitor timeMon (*timer);
324  if (! isInitialized ()) {
325  initialize ();
326  }
327  IsComputed_ = false;
328  impl_.compute ();
329  }
330  IsComputed_ = true;
331  ++NumCompute_;
332 
333  ComputeTime_ += (timer->wallTime() - startTime);
334 }
335 
336 
337 template <class MatrixType>
339  std::ostringstream out;
340 
341  // Output is a valid YAML dictionary in flow style. If you don't
342  // like everything on a single line, you should call describe()
343  // instead.
344  out << "\"Ifpack2::Chebyshev\": {";
345  out << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
346  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
347 
348  out << impl_.description() << ", ";
349 
350  if (impl_.getMatrix ().is_null ()) {
351  out << "Matrix: null";
352  }
353  else {
354  out << "Global matrix dimensions: ["
355  << impl_.getMatrix ()->getGlobalNumRows () << ", "
356  << impl_.getMatrix ()->getGlobalNumCols () << "]"
357  << ", Global nnz: " << impl_.getMatrix ()->getGlobalNumEntries();
358  }
359 
360  out << "}";
361  return out.str ();
362 }
363 
364 
365 template <class MatrixType>
367 describe (Teuchos::FancyOStream& out,
368  const Teuchos::EVerbosityLevel verbLevel) const
369 {
370  using Teuchos::TypeNameTraits;
371  using std::endl;
372 
373  // Default verbosity level is VERB_LOW
374  const Teuchos::EVerbosityLevel vl =
375  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
376 
377  if (vl == Teuchos::VERB_NONE) {
378  return; // print NOTHING, not even the class name
379  }
380 
381  // By convention, describe() starts with a tab.
382  //
383  // This does affect all processes on which it's valid to print to
384  // 'out'. However, it does not actually print spaces to 'out'
385  // unless operator<< gets called, so it's safe to use on all
386  // processes.
387  Teuchos::OSTab tab0 (out);
388  const int myRank = this->getComm ()->getRank ();
389  if (myRank == 0) {
390  // Output is a valid YAML dictionary.
391  // In particular, we quote keys with colons in them.
392  out << "\"Ifpack2::Chebyshev\":" << endl;
393  }
394 
395  Teuchos::OSTab tab1 (out);
396  if (vl >= Teuchos::VERB_LOW && myRank == 0) {
397  out << "Template parameters:" << endl;
398  {
399  Teuchos::OSTab tab2 (out);
400  out << "Scalar: " << TypeNameTraits<scalar_type>::name () << endl
401  << "LocalOrdinal: " << TypeNameTraits<local_ordinal_type>::name () << endl
402  << "GlobalOrdinal: " << TypeNameTraits<global_ordinal_type>::name () << endl
403  << "Device: " << TypeNameTraits<device_type>::name () << endl;
404  }
405  out << "Initialized: " << (isInitialized () ? "true" : "false") << endl
406  << "Computed: " << (isComputed () ? "true" : "false") << endl;
407  impl_.describe (out, vl);
408 
409  if (impl_.getMatrix ().is_null ()) {
410  out << "Matrix: null" << endl;
411  }
412  else {
413  out << "Global matrix dimensions: ["
414  << impl_.getMatrix ()->getGlobalNumRows () << ", "
415  << impl_.getMatrix ()->getGlobalNumCols () << "]" << endl
416  << "Global nnz: " << impl_.getMatrix ()->getGlobalNumEntries() << endl;
417  }
418  }
419 }
420 
421 template<class MatrixType>
422 void
424 applyImpl (const MV& X,
425  MV& Y,
426  Teuchos::ETransp /* mode */,
427  scalar_type alpha,
428  scalar_type beta) const
429 {
430  using Teuchos::ArrayRCP;
431  using Teuchos::as;
432  using Teuchos::RCP;
433  using Teuchos::rcp;
434  using Teuchos::rcp_const_cast;
435  using Teuchos::rcpFromRef;
436 
437  const scalar_type zero = STS::zero();
438  const scalar_type one = STS::one();
439 
440  // Y = beta*Y + alpha*M*X.
441 
442  // If alpha == 0, then we don't need to do Chebyshev at all.
443  if (alpha == zero) {
444  if (beta == zero) { // Obey Sparse BLAS rules; avoid 0*NaN.
445  Y.putScalar (zero);
446  }
447  else {
448  Y.scale (beta);
449  }
450  return;
451  }
452 
453  // If beta != 0, then we need to keep a (deep) copy of the initial
454  // value of Y, so that we can add beta*it to the Chebyshev result at
455  // the end. Usually this method is called with beta == 0, so we
456  // don't have to worry about caching Y_org.
457  RCP<MV> Y_orig;
458  if (beta != zero) {
459  Y_orig = rcp (new MV (Y, Teuchos::Copy));
460  }
461 
462  // If X and Y point to the same memory location, we need to use a
463  // (deep) copy of X (X_copy) as the input MV. Otherwise, just let
464  // X_copy point to X.
465  //
466  // This is hopefully an uncommon use case, so we don't bother to
467  // optimize for it by caching X_copy.
468  RCP<const MV> X_copy;
469  bool copiedInput = false;
470  if (X.aliases(Y)) {
471  X_copy = rcp (new MV (X, Teuchos::Copy));
472  copiedInput = true;
473  } else {
474  X_copy = rcpFromRef (X);
475  }
476 
477  // If alpha != 1, fold alpha into (a deep copy of) X.
478  //
479  // This is an uncommon use case, so we don't bother to optimize for
480  // it by caching X_copy. However, we do check whether we've already
481  // copied X above, to avoid a second copy.
482  if (alpha != one) {
483  RCP<MV> X_copy_nonConst = rcp_const_cast<MV> (X_copy);
484  if (! copiedInput) {
485  X_copy_nonConst = rcp (new MV (X, Teuchos::Copy));
486  copiedInput = true;
487  }
488  X_copy_nonConst->scale (alpha);
489  X_copy = rcp_const_cast<const MV> (X_copy_nonConst);
490  }
491 
492  impl_.apply (*X_copy, Y);
493 
494  if (beta != zero) {
495  Y.update (beta, *Y_orig, one); // Y = beta * Y_orig + 1 * Y
496  }
497 }
498 
499 
500 template<class MatrixType>
501 typename MatrixType::scalar_type Chebyshev<MatrixType>::getLambdaMaxForApply () const {
502  return impl_.getLambdaMaxForApply ();
503 }
504 
505 
506 
507 }//namespace Ifpack2
508 
509 #define IFPACK2_CHEBYSHEV_INSTANT(S,LO,GO,N) \
510  template class Ifpack2::Chebyshev< Tpetra::RowMatrix<S, LO, GO, N> >;
511 
512 #endif // IFPACK2_CHEBYSHEV_DEF_HPP
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to a Teuchos::FancyOStream.
Definition: Ifpack2_Chebyshev_def.hpp:367
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:223
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:217
Teuchos::RCP< const map_type > getRangeMap() const
The Tpetra::Map representing the range of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:158
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Chebyshev_def.hpp:81
void setParameters(const Teuchos::ParameterList &params)
Set (or reset) parameters.
Definition: Ifpack2_Chebyshev_def.hpp:93
double getApplyTime() const
The total time spent in all calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:206
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:229
int getNumCompute() const
The total number of successful calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:182
void apply(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS, scalar_type alpha=Teuchos::ScalarTraits< scalar_type >::one(), scalar_type beta=Teuchos::ScalarTraits< scalar_type >::zero()) const
Apply the preconditioner to X, returning the result in Y.
Definition: Ifpack2_Chebyshev_def.hpp:238
void compute()
(Re)compute the left scaling, and (if applicable) estimate max and min eigenvalues of D_inv * A...
Definition: Ifpack2_Chebyshev_def.hpp:311
Teuchos::RCP< const Tpetra::CrsMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > getCrsMatrix() const
Attempt to return the matrix A as a Tpetra::CrsMatrix.
Definition: Ifpack2_Chebyshev_def.hpp:134
MatrixType::scalar_type getLambdaMaxForApply() const
The estimate of the maximum eigenvalue used in the apply().
Definition: Ifpack2_Chebyshev_def.hpp:501
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix for which this is a preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:123
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:199
Chebyshev(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_Chebyshev_def.hpp:58
void applyMat(const Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &X, Tpetra::MultiVector< scalar_type, local_ordinal_type, global_ordinal_type, node_type > &Y, Teuchos::ETransp mode=Teuchos::NO_TRANS) const
Compute Y = Op(A)*X, where Op(A) is either A, , or .
Definition: Ifpack2_Chebyshev_def.hpp:277
bool hasTransposeApply() const
Whether it&#39;s possible to apply the transpose of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:170
void initialize()
Initialize the preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:296
double getComputeTime() const
The total time spent in all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:200
double getApplyFlops() const
The total number of floating-point operations taken by all calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:218
int getNumInitialize() const
The total number of successful calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:176
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:220
std::string description() const
A simple one-line description of this object.
Definition: Ifpack2_Chebyshev_def.hpp:338
Teuchos::RCP< const map_type > getDomainMap() const
The Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:144
double getInitializeTime() const
The total time spent in all calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:194
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_Chebyshev_def.hpp:223
virtual ~Chebyshev()
Destructor.
Definition: Ifpack2_Chebyshev_def.hpp:76
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition: Ifpack2_Chebyshev_def.hpp:109
int getNumApply() const
The total number of successful calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:188
double getComputeFlops() const
The total number of floating-point operations taken by all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:212
void setZeroStartingSolution(bool zeroStartingSolution)
Set this preconditioner&#39;s parameters.
Definition: Ifpack2_Chebyshev_def.hpp:102