Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
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"
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 template<class MatrixType>
217  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix();
219  A.is_null (), std::runtime_error, "Ifpack2::Chevyshev::getNodeSmootherComplexity: "
220  "The input matrix A is null. Please call setMatrix() with a nonnull "
221  "input matrix, then call compute(), before calling this method.");
222  // Chevyshev costs roughly one apply + one diagonal inverse per iteration
223  return A->getNodeNumRows() + A->getNodeNumEntries();
224 }
225 
226 
227 
228 template<class MatrixType>
229 void
231 apply (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
232  Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
233  Teuchos::ETransp mode,
234  scalar_type alpha,
235  scalar_type beta) const
236 {
237  const std::string timerName ("Ifpack2::Chebyshev::apply");
239  if (timer.is_null ()) {
240  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
241  }
242 
243  // Start timing here.
244  {
245  Teuchos::TimeMonitor timeMon (*timer);
246 
247  // compute() calls initialize() if it hasn't already been called.
248  // Thus, we only need to check isComputed().
250  ! isComputed (), std::runtime_error,
251  "Ifpack2::Chebyshev::apply(): You must call the compute() method before "
252  "you may call apply().");
254  X.getNumVectors () != Y.getNumVectors (), std::runtime_error,
255  "Ifpack2::Chebyshev::apply(): X and Y must have the same number of "
256  "columns. X.getNumVectors() = " << X.getNumVectors() << " != "
257  << "Y.getNumVectors() = " << Y.getNumVectors() << ".");
258  applyImpl (X, Y, mode, alpha, beta);
259  }
260  ++NumApply_;
261 
262  // timer->totalElapsedTime() returns the total time over all timer
263  // calls. Thus, we use = instead of +=.
264  ApplyTime_ = timer->totalElapsedTime ();
265 }
266 
267 
268 template<class MatrixType>
269 void
271 applyMat (const Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& X,
272  Tpetra::MultiVector<scalar_type, local_ordinal_type, global_ordinal_type, node_type>& Y,
273  Teuchos::ETransp mode) const
274 {
276  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
277  "Ifpack2::Chebyshev::applyMat: X.getNumVectors() != Y.getNumVectors().");
278 
279  Teuchos::RCP<const row_matrix_type> A = impl_.getMatrix ();
281  A.is_null (), std::runtime_error, "Ifpack2::Chebyshev::applyMat: The input "
282  "matrix A is null. Please call setMatrix() with a nonnull input matrix "
283  "before calling this method.");
284 
285  A->apply (X, Y, mode);
286 }
287 
288 
289 template<class MatrixType>
291  // We create the timer, but this method doesn't do anything, so
292  // there is no need to start the timer. The resulting total time
293  // will always be zero.
294  const std::string timerName ("Ifpack2::Chebyshev::initialize");
296  if (timer.is_null ()) {
297  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
298  }
299  IsInitialized_ = true;
300  ++NumInitialize_;
301 }
302 
303 
304 template<class MatrixType>
306 {
307  const std::string timerName ("Ifpack2::Chebyshev::compute");
309  if (timer.is_null ()) {
310  timer = Teuchos::TimeMonitor::getNewCounter (timerName);
311  }
312 
313  // Start timing here.
314  {
315  Teuchos::TimeMonitor timeMon (*timer);
316  if (! isInitialized ()) {
317  initialize ();
318  }
319  IsComputed_ = false;
320  impl_.compute ();
321  }
322  IsComputed_ = true;
323  ++NumCompute_;
324 
325  // timer->totalElapsedTime() returns the total time over all timer
326  // calls. Thus, we use = instead of +=.
327  ComputeTime_ = timer->totalElapsedTime ();
328 }
329 
330 
331 template <class MatrixType>
333  std::ostringstream out;
334 
335  // Output is a valid YAML dictionary in flow style. If you don't
336  // like everything on a single line, you should call describe()
337  // instead.
338  out << "\"Ifpack2::Chebyshev\": {";
339  out << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
340  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
341 
342  out << impl_.description() << ", ";
343 
344  if (impl_.getMatrix ().is_null ()) {
345  out << "Matrix: null";
346  }
347  else {
348  out << "Global matrix dimensions: ["
349  << impl_.getMatrix ()->getGlobalNumRows () << ", "
350  << impl_.getMatrix ()->getGlobalNumCols () << "]"
351  << ", Global nnz: " << impl_.getMatrix ()->getGlobalNumEntries();
352  }
353 
354  out << "}";
355  return out.str ();
356 }
357 
358 
359 template <class MatrixType>
362  const Teuchos::EVerbosityLevel verbLevel) const
363 {
365  using std::endl;
366 
367  // Default verbosity level is VERB_LOW
368  const Teuchos::EVerbosityLevel vl =
369  (verbLevel == Teuchos::VERB_DEFAULT) ? Teuchos::VERB_LOW : verbLevel;
370 
371  if (vl == Teuchos::VERB_NONE) {
372  return; // print NOTHING, not even the class name
373  }
374 
375  // By convention, describe() starts with a tab.
376  //
377  // This does affect all processes on which it's valid to print to
378  // 'out'. However, it does not actually print spaces to 'out'
379  // unless operator<< gets called, so it's safe to use on all
380  // processes.
381  Teuchos::OSTab tab0 (out);
382  const int myRank = this->getComm ()->getRank ();
383  if (myRank == 0) {
384  // Output is a valid YAML dictionary.
385  // In particular, we quote keys with colons in them.
386  out << "\"Ifpack2::Chebyshev\":" << endl;
387  }
388 
389  Teuchos::OSTab tab1 (out);
390  if (vl >= Teuchos::VERB_LOW && myRank == 0) {
391  out << "Template parameters:" << endl;
392  {
393  Teuchos::OSTab tab2 (out);
394  out << "Scalar: " << TypeNameTraits<scalar_type>::name () << endl
395  << "LocalOrdinal: " << TypeNameTraits<local_ordinal_type>::name () << endl
396  << "GlobalOrdinal: " << TypeNameTraits<global_ordinal_type>::name () << endl
397  << "Device: " << TypeNameTraits<device_type>::name () << endl;
398  }
399  out << "Initialized: " << (isInitialized () ? "true" : "false") << endl
400  << "Computed: " << (isComputed () ? "true" : "false") << endl;
401  impl_.describe (out, vl);
402 
403  if (impl_.getMatrix ().is_null ()) {
404  out << "Matrix: null" << endl;
405  }
406  else {
407  out << "Global matrix dimensions: ["
408  << impl_.getMatrix ()->getGlobalNumRows () << ", "
409  << impl_.getMatrix ()->getGlobalNumCols () << "]" << endl
410  << "Global nnz: " << impl_.getMatrix ()->getGlobalNumEntries() << endl;
411  }
412  }
413 }
414 
415 template<class MatrixType>
416 void
418 applyImpl (const MV& X,
419  MV& Y,
420  Teuchos::ETransp /* mode */,
421  scalar_type alpha,
422  scalar_type beta) const
423 {
424  using Teuchos::ArrayRCP;
425  using Teuchos::as;
426  using Teuchos::RCP;
427  using Teuchos::rcp;
428  using Teuchos::rcp_const_cast;
429  using Teuchos::rcpFromRef;
430 
431  const scalar_type zero = STS::zero();
432  const scalar_type one = STS::one();
433 
434  // Y = beta*Y + alpha*M*X.
435 
436  // If alpha == 0, then we don't need to do Chebyshev at all.
437  if (alpha == zero) {
438  if (beta == zero) { // Obey Sparse BLAS rules; avoid 0*NaN.
439  Y.putScalar (zero);
440  }
441  else {
442  Y.scale (beta);
443  }
444  return;
445  }
446 
447  // If beta != 0, then we need to keep a (deep) copy of the initial
448  // value of Y, so that we can add beta*it to the Chebyshev result at
449  // the end. Usually this method is called with beta == 0, so we
450  // don't have to worry about caching Y_org.
451  RCP<MV> Y_orig;
452  if (beta != zero) {
453  Y_orig = rcp (new MV (Y, Teuchos::Copy));
454  }
455 
456  // If X and Y point to the same memory location, we need to use a
457  // (deep) copy of X (X_copy) as the input MV. Otherwise, just let
458  // X_copy point to X.
459  //
460  // This is hopefully an uncommon use case, so we don't bother to
461  // optimize for it by caching X_copy.
462  RCP<const MV> X_copy;
463  bool copiedInput = false;
464  {
465 #ifdef KOKKOS_ENABLE_DEPRECATED_CODE
466  auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
467  auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
468 #else
469  auto X_lcl_host = X.getLocalViewHost ();
470  auto Y_lcl_host = Y.getLocalViewHost ();
471 #endif
472  if (X_lcl_host.data () == Y_lcl_host.data ()) {
473  X_copy = rcp (new MV (X, Teuchos::Copy));
474  copiedInput = true;
475  } else {
476  X_copy = rcpFromRef (X);
477  }
478  }
479 
480  // If alpha != 1, fold alpha into (a deep copy of) X.
481  //
482  // This is an uncommon use case, so we don't bother to optimize for
483  // it by caching X_copy. However, we do check whether we've already
484  // copied X above, to avoid a second copy.
485  if (alpha != one) {
486  RCP<MV> X_copy_nonConst = rcp_const_cast<MV> (X_copy);
487  if (! copiedInput) {
488  X_copy_nonConst = rcp (new MV (X, Teuchos::Copy));
489  copiedInput = true;
490  }
491  X_copy_nonConst->scale (alpha);
492  X_copy = rcp_const_cast<const MV> (X_copy_nonConst);
493  }
494 
495  impl_.apply (*X_copy, Y);
496 
497  if (beta != zero) {
498  Y.update (beta, *Y_orig, one); // Y = beta * Y_orig + 1 * Y
499  }
500 }
501 
502 
503 template<class MatrixType>
504 typename MatrixType::scalar_type Chebyshev<MatrixType>::getLambdaMaxForApply () const {
505  return impl_.getLambdaMaxForApply ();
506 }
507 
508 
509 
510 }//namespace Ifpack2
511 
512 #define IFPACK2_CHEBYSHEV_INSTANT(S,LO,GO,N) \
513  template class Ifpack2::Chebyshev< Tpetra::RowMatrix<S, LO, GO, N> >;
514 
515 #endif // IFPACK2_CHEBYSHEV_DEF_HPP
double getApplyFlops() const
The total number of floating-point operations taken by all calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:211
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::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
double getApplyTime() const
The total time spent in all calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:199
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 getComputeTime() const
The total time spent in all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:193
int getNumApply() const
The total number of successful calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:181
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:229
void compute()
(Re)compute the left scaling, and (if applicable) estimate max and min eigenvalues of D_inv * A...
Definition: Ifpack2_Chebyshev_def.hpp:305
bool hasTransposeApply() const
Whether it&#39;s possible to apply the transpose of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:163
T * getRawPtr() const
double getComputeFlops() const
The total number of floating-point operations taken by all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:205
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:199
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:231
std::string description() const
A simple one-line description of this object.
Definition: Ifpack2_Chebyshev_def.hpp:332
void initialize()
Initialize the preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:290
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition: Ifpack2_Chebyshev_def.hpp:102
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:220
int getNumInitialize() const
The total number of successful calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:169
double getInitializeTime() const
The total time spent in all calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:187
MatrixType::scalar_type getLambdaMaxForApply() const
The estimate of the maximum eigenvalue used in the apply().
Definition: Ifpack2_Chebyshev_def.hpp:504
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_Chebyshev_def.hpp:216
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:361
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix for which this is a preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:116
TypeTo as(const TypeFrom &t)
double totalElapsedTime(bool readCurrentTime=false) const
Teuchos::RCP< const map_type > getDomainMap() const
The Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:137
virtual ~Chebyshev()
Destructor.
Definition: Ifpack2_Chebyshev_def.hpp:76
Teuchos::RCP< const map_type > getRangeMap() const
The Tpetra::Map representing the range of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:151
int getNumCompute() const
The total number of successful calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:175
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:271
bool is_null() const