Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_Chebyshev_def.hpp
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated 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"
50 #include <iostream>
51 #include <sstream>
52 
53 
54 namespace Ifpack2 {
55 
56 template<class MatrixType>
57 Chebyshev<MatrixType>::
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>
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
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>
103 {
104  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
106  A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::getComm: The input "
107  "matrix A is null. Please call setMatrix() with a nonnull input matrix "
108  "before calling this method.");
109  return A->getRowMap ()->getComm ();
110 }
111 
112 
113 template<class MatrixType>
116 getMatrix() const {
117  return impl_.getMatrix ();
118 }
119 
120 
121 template<class MatrixType>
122 Teuchos::RCP<const Tpetra::CrsMatrix<typename MatrixType::scalar_type,
123  typename MatrixType::local_ordinal_type,
124  typename MatrixType::global_ordinal_type,
125  typename MatrixType::node_type> >
127 getCrsMatrix() const {
128  typedef Tpetra::CrsMatrix<scalar_type, local_ordinal_type,
129  global_ordinal_type, node_type> crs_matrix_type;
130  return Teuchos::rcp_dynamic_cast<const crs_matrix_type> (impl_.getMatrix ());
131 }
132 
133 
134 template<class MatrixType>
137 getDomainMap () const
138 {
139  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
141  A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::getDomainMap: The "
142  "input matrix A is null. Please call setMatrix() with a nonnull input "
143  "matrix before calling this method.");
144  return A->getDomainMap ();
145 }
146 
147 
148 template<class MatrixType>
151 getRangeMap () const
152 {
153  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
155  A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::getRangeMap: The "
156  "input matrix A is null. Please call setMatrix() with a nonnull input "
157  "matrix before calling this method.");
158  return A->getRangeMap ();
159 }
160 
161 
162 template<class MatrixType>
164  return impl_.hasTransposeApply ();
165 }
166 
167 
168 template<class MatrixType>
170  return NumInitialize_;
171 }
172 
173 
174 template<class MatrixType>
176  return NumCompute_;
177 }
178 
179 
180 template<class MatrixType>
182  return NumApply_;
183 }
184 
185 
186 template<class MatrixType>
188  return InitializeTime_;
189 }
190 
191 
192 template<class MatrixType>
194  return ComputeTime_;
195 }
196 
197 
198 template<class MatrixType>
200  return ApplyTime_;
201 }
202 
203 
204 template<class MatrixType>
206  return ComputeFlops_;
207 }
208 
209 
210 template<class MatrixType>
212  return ApplyFlops_;
213 }
214 
215 
216 template<class MatrixType>
217 void
219 apply (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
220  Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
221  Teuchos::ETransp mode,
222  scalar_type alpha,
223  scalar_type beta) const
224 {
225  const std::string timerName ("Ifpack2::Chebyshev::apply");
227  if (timer.is_null ()) {
228  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
229  }
230 
231  // Start timing here.
232  {
233  Teuchos::TimeMonitor timeMon (*timer);
234 
235  // compute() calls initialize() if it hasn't already been called.
236  // Thus, we only need to check isComputed().
238  ! isComputed (), std::runtime_error,
239  "Ifpack2::Chebyshev::apply(): You must call the compute() method before "
240  "you may call apply().");
242  X.getNumVectors () != Y.getNumVectors (), std::runtime_error,
243  "Ifpack2::Chebyshev::apply(): X and Y must have the same number of "
244  "columns. X.getNumVectors() = " << X.getNumVectors() << " != "
245  << "Y.getNumVectors() = " << Y.getNumVectors() << ".");
246  applyImpl (X, Y, mode, alpha, beta);
247  }
248  ++NumApply_;
249 
250  // timer->totalElapsedTime() returns the total time over all timer
251  // calls. Thus, we use = instead of +=.
252  ApplyTime_ = timer->totalElapsedTime ();
253 }
254 
255 
256 template<class MatrixType>
257 void
259 applyMat (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
260  Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
261  Teuchos::ETransp mode) const
262 {
264  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
265  "Ifpack2::Chebyshev::applyMat: X.getNumVectors() != Y.getNumVectors().");
266 
267  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
269  A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::applyMat: The input "
270  "matrix A is null. Please call setMatrix() with a nonnull input matrix "
271  "before calling this method.");
272 
273  A->apply (X, Y, mode);
274 }
275 
276 
277 template<class MatrixType>
279  // We create the timer, but this method doesn't do anything, so
280  // there is no need to start the timer. The resulting total time
281  // will always be zero.
282  const std::string timerName ("Ifpack2::Chebyshev::initialize");
284  if (timer.is_null ()) {
285  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
286  }
287  IsInitialized_ = true;
288  ++NumInitialize_;
289 }
290 
291 
292 template<class MatrixType>
294 {
295  const std::string timerName ("Ifpack2::Chebyshev::compute");
297  if (timer.is_null ()) {
298  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
299  }
300 
301  // Start timing here.
302  {
303  Teuchos::TimeMonitor timeMon (*timer);
304  if (! isInitialized ()) {
305  initialize ();
306  }
307  IsComputed_ = false;
308  impl_.compute ();
309  }
310  IsComputed_ = true;
311  ++NumCompute_;
312 
313  // timer->totalElapsedTime() returns the total time over all timer
314  // calls. Thus, we use = instead of +=.
315  ComputeTime_ = timer->totalElapsedTime ();
316 }
317 
318 
319 template <class MatrixType>
321  std::ostringstream out;
322 
323  // Output is a valid YAML dictionary in flow style. If you don't
324  // like everything on a single line, you should call describe()
325  // instead.
326  out << "\"Ifpack2::Chebyshev\": {";
327  out << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
328  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
329 
330  out << impl_.description() << ", ";
331 
332  if (impl_.getMatrix ().is_null ()) {
333  out << "Matrix: null";
334  }
335  else {
336  out << "Global matrix dimensions: ["
337  << impl_.getMatrix ()->getGlobalNumRows () << ", "
338  << impl_.getMatrix ()->getGlobalNumCols () << "]"
339  << ", Global nnz: " << impl_.getMatrix ()->getGlobalNumEntries();
340  }
341 
342  out << "}";
343  return out.str ();
344 }
345 
346 
347 template <class MatrixType>
350  const Teuchos::EVerbosityLevel verbLevel) const
351 {
353  using std::endl;
354 
355  // Default verbosity level is VERB_LOW
356  const Teuchos::EVerbosityLevel vl =
357  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
358 
359  if (vl == Teuchos::VERB_NONE) {
360  return; // print NOTHING, not even the class name
361  }
362 
363  // By convention, describe() starts with a tab.
364  //
365  // This does affect all processes on which it's valid to print to
366  // 'out'. However, it does not actually print spaces to 'out'
367  // unless operator<< gets called, so it's safe to use on all
368  // processes.
369  Teuchos::OSTab tab0 (out);
370  const int myRank = this->getComm ()->getRank ();
371  if (myRank == 0) {
372  // Output is a valid YAML dictionary.
373  // In particular, we quote keys with colons in them.
374  out << "\"Ifpack2::Chebyshev\":" << endl;
375  }
376 
377  Teuchos::OSTab tab1 (out);
378  if (vl >= Teuchos::VERB_LOW && myRank == 0) {
379  out << "Template parameters:" << endl;
380  {
381  Teuchos::OSTab tab2 (out);
382  out << "Scalar: " << TypeNameTraits<scalar_type>::name () << endl
383  << "LocalOrdinal: " << TypeNameTraits<local_ordinal_type>::name () << endl
384  << "GlobalOrdinal: " << TypeNameTraits<global_ordinal_type>::name () << endl
385  << "Device: " << TypeNameTraits<device_type>::name () << endl;
386  }
387  out << "Initialized: " << (isInitialized () ? "true" : "false") << endl
388  << "Computed: " << (isComputed () ? "true" : "false") << endl;
389  impl_.describe (out, vl);
390 
391  if (impl_.getMatrix ().is_null ()) {
392  out << "Matrix: null" << endl;
393  }
394  else {
395  out << "Global matrix dimensions: ["
396  << impl_.getMatrix ()->getGlobalNumRows () << ", "
397  << impl_.getMatrix ()->getGlobalNumCols () << "]" << endl
398  << "Global nnz: " << impl_.getMatrix ()->getGlobalNumEntries() << endl;
399  }
400  }
401 }
402 
403 template<class MatrixType>
404 void
406 applyImpl (const MV& X,
407  MV& Y,
408  Teuchos::ETransp mode,
409  scalar_type alpha,
410  scalar_type beta) const
411 {
412  using Teuchos::ArrayRCP;
413  using Teuchos::as;
414  using Teuchos::RCP;
415  using Teuchos::rcp;
416  using Teuchos::rcp_const_cast;
417  using Teuchos::rcpFromRef;
418 
419  const scalar_type zero = STS::zero();
420  const scalar_type one = STS::one();
421 
422  // Y = beta*Y + alpha*M*X.
423 
424  // If alpha == 0, then we don't need to do Chebyshev at all.
425  if (alpha == zero) {
426  if (beta == zero) { // Obey Sparse BLAS rules; avoid 0*NaN.
427  Y.putScalar (zero);
428  }
429  else {
430  Y.scale (beta);
431  }
432  return;
433  }
434 
435  // If beta != 0, then we need to keep a (deep) copy of the initial
436  // value of Y, so that we can add beta*it to the Chebyshev result at
437  // the end. Usually this method is called with beta == 0, so we
438  // don't have to worry about caching Y_org.
439  RCP<MV> Y_orig;
440  if (beta != zero) {
441  Y_orig = rcp (new MV (Y, Teuchos::Copy));
442  }
443 
444  // If X and Y point to the same memory location, we need to use a
445  // (deep) copy of X (X_copy) as the input MV. Otherwise, just let
446  // X_copy point to X.
447  //
448  // This is hopefully an uncommon use case, so we don't bother to
449  // optimize for it by caching X_copy.
450  RCP<const MV> X_copy;
451  bool copiedInput = false;
452  {
453  auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
454  auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
455  if (X_lcl_host.ptr_on_device () == Y_lcl_host.ptr_on_device ()) {
456  X_copy = rcp (new MV (X, Teuchos::Copy));
457  copiedInput = true;
458  } else {
459  X_copy = rcpFromRef (X);
460  }
461  }
462 
463  // If alpha != 1, fold alpha into (a deep copy of) X.
464  //
465  // This is an uncommon use case, so we don't bother to optimize for
466  // it by caching X_copy. However, we do check whether we've already
467  // copied X above, to avoid a second copy.
468  if (alpha != one) {
469  RCP<MV> X_copy_nonConst = rcp_const_cast<MV> (X_copy);
470  if (! copiedInput) {
471  X_copy_nonConst = rcp (new MV (X, Teuchos::Copy));
472  copiedInput = true;
473  }
474  X_copy_nonConst->scale (alpha);
475  X_copy = rcp_const_cast<const MV> (X_copy_nonConst);
476  }
477 
478  impl_.apply (*X_copy, Y);
479 
480  if (beta != zero) {
481  Y.update (beta, *Y_orig, one); // Y = beta * Y_orig + 1 * Y
482  }
483 }
484 
485 
486 template<class MatrixType>
487 typename MatrixType::scalar_type Chebyshev<MatrixType>::getLambdaMaxForApply () const {
488  return impl_.getLambdaMaxForApply ();
489 }
490 
491 
492 
493 }//namespace Ifpack2
494 
495 #define IFPACK2_CHEBYSHEV_INSTANT(S,LO,GO,N) \
496  template class Ifpack2::Chebyshev< Tpetra::RowMatrix<S, LO, GO, N> >;
497 
498 #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:349
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
VERB_DEFAULT
Teuchos::RCP< const map_type > getRangeMap() const
The Tpetra::Map representing the range of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:151
VERB_LOW
VERB_NONE
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Chebyshev_def.hpp:81
static RCP< Time > getNewCounter(const std::string &name)
static RCP< Time > lookupCounter(const std::string &name)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
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:199
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:175
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:219
void compute()
(Re)compute the left scaling, and (if applicable) estimate max and min eigenvalues of D_inv * A...
Definition: Ifpack2_Chebyshev_def.hpp:293
bool is_null() const
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:127
MatrixType::scalar_type getLambdaMaxForApply() const
The estimate of the maximum eigenvalue used in the apply().
Definition: Ifpack2_Chebyshev_def.hpp:487
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix for which this is a preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:116
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:199
LinearOp zero(const VectorSpace &vs)
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:259
bool hasTransposeApply() const
Whether it&#39;s possible to apply the transpose of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:163
void initialize()
Initialize the preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:278
double getComputeTime() const
The total time spent in all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:193
double totalElapsedTime(bool readCurrentTime=false) const
double getApplyFlops() const
The total number of floating-point operations taken by all calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:211
int getNumInitialize() const
The total number of successful calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:169
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:320
Teuchos::RCP< const map_type > getDomainMap() const
The Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:137
TypeTo as(const TypeFrom &t)
T * getRawPtr() const
double getInitializeTime() const
The total time spent in all calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:187
virtual ~Chebyshev()
Destructor.
Definition: Ifpack2_Chebyshev_def.hpp:76
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition: Ifpack2_Chebyshev_def.hpp:102
int getNumApply() const
The total number of successful calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:181
double getComputeFlops() const
The total number of floating-point operations taken by all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:205