Ifpack2 Templated Preconditioning Package Version 1.0
Loading...
Searching...
No Matches
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
54namespace Ifpack2 {
55
56template<class MatrixType>
58Chebyshev (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
75template<class MatrixType>
77}
78
79
80template<class MatrixType>
81void 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
91template<class MatrixType>
92void
93Chebyshev<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
100template<class MatrixType>
101void
103{
104 impl_.setZeroStartingSolution(zeroStartingSolution);
105}
106
107template<class MatrixType>
108Teuchos::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
120template<class MatrixType>
121Teuchos::RCP<const typename Chebyshev<MatrixType>::row_matrix_type>
123getMatrix() const {
124 return impl_.getMatrix ();
125}
126
127
128template<class MatrixType>
129Teuchos::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> >
134getCrsMatrix() 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
141template<class MatrixType>
142Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>
144getDomainMap () 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
155template<class MatrixType>
156Teuchos::RCP<const typename Chebyshev<MatrixType>::map_type>
158getRangeMap () 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
169template<class MatrixType>
171 return impl_.hasTransposeApply ();
172}
173
174
175template<class MatrixType>
177 return NumInitialize_;
178}
179
180
181template<class MatrixType>
183 return NumCompute_;
184}
185
186
187template<class MatrixType>
189 return NumApply_;
190}
191
192
193template<class MatrixType>
195 return InitializeTime_;
196}
197
198
199template<class MatrixType>
201 return ComputeTime_;
202}
203
204
205template<class MatrixType>
207 return ApplyTime_;
208}
209
210
211template<class MatrixType>
213 return ComputeFlops_;
214}
215
216
217template<class MatrixType>
219 return ApplyFlops_;
220}
221
222template<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->getLocalNumRows() + A->getLocalNumEntries();
231}
232
233
234
235template<class MatrixType>
236void
238apply (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
274template<class MatrixType>
275void
277applyMat (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
295template<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
310template<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
337template <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
365template <class MatrixType>
367describe (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
421template<class MatrixType>
422void
424applyImpl (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
500template<class MatrixType>
501typename 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
Diagonally scaled Chebyshev iteration for Tpetra sparse matrices.
Definition: Ifpack2_Chebyshev_decl.hpp:208
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:223
double getInitializeTime() const
The total time spent in all calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:194
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
std::string description() const
A simple one-line description of this object.
Definition: Ifpack2_Chebyshev_def.hpp:338
Chebyshev(const Teuchos::RCP< const row_matrix_type > &A)
Constructor.
Definition: Ifpack2_Chebyshev_def.hpp:58
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
Teuchos::RCP< const map_type > getRangeMap() const
The Tpetra::Map representing the range of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:158
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
void initialize()
Initialize the preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:296
Teuchos::RCP< const row_matrix_type > getMatrix() const
The matrix for which this is a preconditioner.
Definition: Ifpack2_Chebyshev_def.hpp:123
int getNumApply() const
The total number of successful calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:188
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:220
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:229
size_t getNodeSmootherComplexity() const
Get a rough estimate of cost per iteration.
Definition: Ifpack2_Chebyshev_def.hpp:223
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Chebyshev_decl.hpp:217
virtual void setMatrix(const Teuchos::RCP< const row_matrix_type > &A)
Change the matrix to be preconditioned.
Definition: Ifpack2_Chebyshev_def.hpp:81
int getNumCompute() const
The total number of successful calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:182
void setParameters(const Teuchos::ParameterList &params)
Set (or reset) parameters.
Definition: Ifpack2_Chebyshev_def.hpp:93
double getComputeFlops() const
The total number of floating-point operations taken by all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:212
double getComputeTime() const
The total time spent in all calls to compute().
Definition: Ifpack2_Chebyshev_def.hpp:200
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
void setZeroStartingSolution(bool zeroStartingSolution)
Set this preconditioner's parameters.
Definition: Ifpack2_Chebyshev_def.hpp:102
Teuchos::RCP< const map_type > getDomainMap() const
The Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:144
int getNumInitialize() const
The total number of successful calls to initialize().
Definition: Ifpack2_Chebyshev_def.hpp:176
virtual ~Chebyshev()
Destructor.
Definition: Ifpack2_Chebyshev_def.hpp:76
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
The communicator over which the matrix is distributed.
Definition: Ifpack2_Chebyshev_def.hpp:109
double getApplyFlops() const
The total number of floating-point operations taken by all calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:218
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
MatrixType::scalar_type getLambdaMaxForApply() const
The estimate of the maximum eigenvalue used in the apply().
Definition: Ifpack2_Chebyshev_def.hpp:501
bool hasTransposeApply() const
Whether it's possible to apply the transpose of this operator.
Definition: Ifpack2_Chebyshev_def.hpp:170
double getApplyTime() const
The total time spent in all calls to apply().
Definition: Ifpack2_Chebyshev_def.hpp:206
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:74