Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_Hiptmair_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_HIPTMAIR_DEF_HPP
44 #define IFPACK2_HIPTMAIR_DEF_HPP
45 
46 #include "Ifpack2_Details_OneLevelFactory.hpp"
47 #include "Ifpack2_Parameters.hpp"
48 #include "Teuchos_TimeMonitor.hpp"
49 #include "Tpetra_MultiVector.hpp"
50 #include <cmath>
51 #include <iostream>
52 #include <sstream>
53 
54 namespace Ifpack2 {
55 
56 template <class MatrixType>
58 Hiptmair (const Teuchos::RCP<const row_matrix_type>& A,
59  const Teuchos::RCP<const row_matrix_type>& PtAP,
60  const Teuchos::RCP<const row_matrix_type>& P) :
61  A_ (A),
62  PtAP_ (PtAP),
63  P_ (P),
64  // Default values
65  precType1_ ("CHEBYSHEV"),
66  precType2_ ("CHEBYSHEV"),
67  preOrPost_ ("both"),
68  ZeroStartingSolution_ (true),
69  // General
70  IsInitialized_ (false),
71  IsComputed_ (false),
72  NumInitialize_ (0),
73  NumCompute_ (0),
74  NumApply_ (0),
75  InitializeTime_ (0.0),
76  ComputeTime_ (0.0),
77  ApplyTime_ (0.0)
78 {}
79 
80 
81 
82 
83 template <class MatrixType>
85 Hiptmair (const Teuchos::RCP<const row_matrix_type>& A):
86  A_ (A),
87  PtAP_ (),
88  P_ (),
89  // Default values
90  precType1_ ("CHEBYSHEV"),
91  precType2_ ("CHEBYSHEV"),
92  preOrPost_ ("both"),
93  ZeroStartingSolution_ (true),
94  // General
95  IsInitialized_ (false),
96  IsComputed_ (false),
97  NumInitialize_ (0),
98  NumCompute_ (0),
99  NumApply_ (0),
100  InitializeTime_ (0.0),
101  ComputeTime_ (0.0),
102  ApplyTime_ (0.0)
103 {}
104 
105 
106 template <class MatrixType>
108 
109 template <class MatrixType>
110 void Hiptmair<MatrixType>::setParameters (const Teuchos::ParameterList& plist)
111 {
112  using Teuchos::as;
113  using Teuchos::RCP;
114  using Teuchos::ParameterList;
115  using Teuchos::Exceptions::InvalidParameterName;
116  using Teuchos::Exceptions::InvalidParameterType;
117 
118  ParameterList params = plist;
119 
120  // Get the current parameters' values. We don't assign to the
121  // instance data directly until we've gotten all the parameters.
122  // This ensures "transactional" semantics, so that if attempting to
123  // get some parameter throws an exception, the class' state doesn't
124  // change.
125  std::string precType1 = precType1_;
126  std::string precType2 = precType2_;
127  std::string preOrPost = preOrPost_;
128  Teuchos::ParameterList precList1 = precList1_;
129  Teuchos::ParameterList precList2 = precList2_;
130  bool zeroStartingSolution = ZeroStartingSolution_;
131 
132  precType1 = params.get("hiptmair: smoother type 1", precType1);
133  precType2 = params.get("hiptmair: smoother type 2", precType2);
134  precList1 = params.get("hiptmair: smoother list 1", precList1);
135  precList2 = params.get("hiptmair: smoother list 2", precList2);
136  preOrPost = params.get("hiptmair: pre or post", preOrPost);
137  zeroStartingSolution = params.get("hiptmair: zero starting solution",
138  zeroStartingSolution);
139 
140  // Grab the matrices off of the parameter list if we need them
141  // This will intentionally throw if they're not there and we need them
142  if(PtAP_.is_null())
143  PtAP_ = params.get<RCP<row_matrix_type> >("PtAP");
144  if(P_.is_null())
145  P_ = params.get<RCP<row_matrix_type> >("P");
146 
147 
148  // "Commit" the new values to the instance data.
149  precType1_ = precType1;
150  precType2_ = precType2;
151  precList1_ = precList1;
152  precList2_ = precList2;
153  preOrPost_ = preOrPost;
154  ZeroStartingSolution_ = zeroStartingSolution;
155 }
156 
157 
158 template <class MatrixType>
159 Teuchos::RCP<const Teuchos::Comm<int> >
161  TEUCHOS_TEST_FOR_EXCEPTION(
162  A_.is_null (), std::runtime_error, "Ifpack2::Hiptmair::getComm: "
163  "The input matrix A is null. Please call setMatrix() with a nonnull "
164  "input matrix before calling this method.");
165  return A_->getComm ();
166 }
167 
168 
169 template <class MatrixType>
170 Teuchos::RCP<const Tpetra::RowMatrix<typename MatrixType::scalar_type,typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
172  return A_;
173 }
174 
175 
176 template <class MatrixType>
177 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
179 {
180  TEUCHOS_TEST_FOR_EXCEPTION(
181  A_.is_null (), std::runtime_error, "Ifpack2::Hiptmair::getDomainMap: "
182  "The input matrix A is null. Please call setMatrix() with a nonnull "
183  "input matrix before calling this method.");
184  return A_->getDomainMap ();
185 }
186 
187 
188 template <class MatrixType>
189 Teuchos::RCP<const Tpetra::Map<typename MatrixType::local_ordinal_type,typename MatrixType::global_ordinal_type,typename MatrixType::node_type> >
191 {
192  TEUCHOS_TEST_FOR_EXCEPTION(
193  A_.is_null (), std::runtime_error, "Ifpack2::Hiptmair::getRangeMap: "
194  "The input matrix A is null. Please call setMatrix() with a nonnull "
195  "input matrix before calling this method.");
196  return A_->getRangeMap ();
197 }
198 
199 
200 template <class MatrixType>
202  // FIXME (mfh 17 Jan 2014) apply() does not currently work with mode
203  // != NO_TRANS, so it's correct to return false here.
204  return false;
205 }
206 
207 
208 template <class MatrixType>
210  return NumInitialize_;
211 }
212 
213 
214 template <class MatrixType>
216  return NumCompute_;
217 }
218 
219 
220 template <class MatrixType>
222  return NumApply_;
223 }
224 
225 
226 template <class MatrixType>
228  return InitializeTime_;
229 }
230 
231 
232 template <class MatrixType>
234  return ComputeTime_;
235 }
236 
237 
238 template <class MatrixType>
240  return ApplyTime_;
241 }
242 
243 
244 template <class MatrixType>
246 {
247  using Teuchos::ParameterList;
248  using Teuchos::RCP;
249  using Teuchos::rcp;
250 
251  TEUCHOS_TEST_FOR_EXCEPTION(
252  A_.is_null (), std::runtime_error, "Ifpack2::Hiptmair::initialize: "
253  "The input matrix A is null. Please call setMatrix() with a nonnull "
254  "input matrix before calling this method.");
255 
256  // clear any previous allocation
257  IsInitialized_ = false;
258  IsComputed_ = false;
259 
260  Teuchos::Time timer ("initialize");
261  double startTime = timer.wallTime();
262  { // The body of code to time
263  Teuchos::TimeMonitor timeMon (timer);
264 
266 
267  ifpack2_prec1_=factory.create(precType1_,A_);
268  ifpack2_prec1_->initialize();
269  ifpack2_prec1_->setParameters(precList1_);
270 
271  ifpack2_prec2_=factory.create(precType2_,PtAP_);
272  ifpack2_prec2_->initialize();
273  ifpack2_prec2_->setParameters(precList2_);
274 
275  }
276  IsInitialized_ = true;
277  ++NumInitialize_;
278  InitializeTime_ += (timer.wallTime() - startTime);
279 }
280 
281 
282 template <class MatrixType>
284 {
285  TEUCHOS_TEST_FOR_EXCEPTION(
286  A_.is_null (), std::runtime_error, "Ifpack2::Hiptmair::compute: "
287  "The input matrix A is null. Please call setMatrix() with a nonnull "
288  "input matrix before calling this method.");
289 
290  // Don't time the initialize(); that gets timed separately.
291  if (! isInitialized ()) {
292  initialize ();
293  }
294 
295  Teuchos::Time timer ("compute");
296  double startTime = timer.wallTime();
297  { // The body of code to time
298  Teuchos::TimeMonitor timeMon (timer);
299  ifpack2_prec1_->compute();
300  ifpack2_prec2_->compute();
301  }
302  IsComputed_ = true;
303  ++NumCompute_;
304  ComputeTime_ += (timer.wallTime() - startTime);
305 }
306 
307 
308 template <class MatrixType>
310 apply (const Tpetra::MultiVector<typename MatrixType::scalar_type,
311  typename MatrixType::local_ordinal_type,
312  typename MatrixType::global_ordinal_type,
313  typename MatrixType::node_type>& X,
314  Tpetra::MultiVector<typename MatrixType::scalar_type,
315  typename MatrixType::local_ordinal_type,
316  typename MatrixType::global_ordinal_type,
317  typename MatrixType::node_type>& Y,
318  Teuchos::ETransp mode,
319  typename MatrixType::scalar_type alpha,
320  typename MatrixType::scalar_type beta) const
321 {
322  using Teuchos::RCP;
323  using Teuchos::rcp;
324  using Teuchos::rcpFromRef;
325  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
327  TEUCHOS_TEST_FOR_EXCEPTION(
328  ! isComputed (), std::runtime_error,
329  "Ifpack2::Hiptmair::apply: You must call compute() before you may call apply().");
330  TEUCHOS_TEST_FOR_EXCEPTION(
331  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
332  "Ifpack2::Hiptmair::apply: The MultiVector inputs X and Y do not have the "
333  "same number of columns. X.getNumVectors() = " << X.getNumVectors ()
334  << " != Y.getNumVectors() = " << Y.getNumVectors () << ".");
335 
336  // Catch unimplemented cases: alpha != 1, beta != 0, mode != NO_TRANS.
337  TEUCHOS_TEST_FOR_EXCEPTION(
338  alpha != STS::one (), std::logic_error,
339  "Ifpack2::Hiptmair::apply: alpha != 1 has not been implemented.");
340  TEUCHOS_TEST_FOR_EXCEPTION(
341  beta != STS::zero (), std::logic_error,
342  "Ifpack2::Hiptmair::apply: zero != 0 has not been implemented.");
343  TEUCHOS_TEST_FOR_EXCEPTION(
344  mode != Teuchos::NO_TRANS, std::logic_error,
345  "Ifpack2::Hiptmair::apply: mode != Teuchos::NO_TRANS has not been implemented.");
346 
347  Teuchos::Time timer ("apply");
348  double startTime = timer.wallTime();
349  { // The body of code to time
350  Teuchos::TimeMonitor timeMon (timer);
351 
352  // If X and Y are pointing to the same memory location,
353  // we need to create an auxiliary vector, Xcopy
354  RCP<const MV> Xcopy;
355  {
356  if (X.aliases(Y)) {
357  Xcopy = rcp (new MV (X, Teuchos::Copy));
358  } else {
359  Xcopy = rcpFromRef (X);
360  }
361  }
362 
363  RCP<MV> Ycopy = rcpFromRef (Y);
364  if (ZeroStartingSolution_) {
365  Ycopy->putScalar (STS::zero ());
366  }
367 
368  // apply Hiptmair Smoothing
369  applyHiptmairSmoother (*Xcopy, *Ycopy);
370 
371  }
372  ++NumApply_;
373  ApplyTime_ += (timer.wallTime() - startTime);
374 }
375 
376 template <class MatrixType>
378 applyHiptmairSmoother(const Tpetra::MultiVector<typename MatrixType::scalar_type,
379  typename MatrixType::local_ordinal_type,
380  typename MatrixType::global_ordinal_type,
381  typename MatrixType::node_type>& X,
382  Tpetra::MultiVector<typename MatrixType::scalar_type,
383  typename MatrixType::local_ordinal_type,
384  typename MatrixType::global_ordinal_type,
385  typename MatrixType::node_type>& Y) const
386 {
387  using Teuchos::RCP;
388  using Teuchos::rcp;
389  using Teuchos::rcpFromRef;
390  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
391  global_ordinal_type, node_type> MV;
392  const scalar_type ZERO = STS::zero ();
393  const scalar_type ONE = STS::one ();
394 
395  RCP<MV> res1 = rcp (new MV (A_->getRowMap (), X.getNumVectors ()));
396  RCP<MV> vec1 = rcp (new MV (A_->getRowMap (), X.getNumVectors ()));
397  RCP<MV> res2 = rcp (new MV (PtAP_->getRowMap (), X.getNumVectors ()));
398  RCP<MV> vec2 = rcp (new MV (PtAP_->getRowMap (), X.getNumVectors ()));
399 
400  if (preOrPost_ == "pre" || preOrPost_ == "both") {
401  // apply initial relaxation to primary space
402  A_->apply (Y, *res1);
403  res1->update (ONE, X, -ONE);
404  vec1->putScalar (ZERO);
405  ifpack2_prec1_->apply (*res1, *vec1);
406  Y.update (ONE, *vec1, ONE);
407  }
408 
409  // project to auxiliary space and smooth
410  A_->apply (Y, *res1);
411  res1->update (ONE, X, -ONE);
412  P_->apply (*res1, *res2, Teuchos::TRANS);
413  vec2->putScalar (ZERO);
414  ifpack2_prec2_->apply (*res2, *vec2);
415  P_->apply (*vec2, *vec1, Teuchos::NO_TRANS);
416  Y.update (ONE,*vec1,ONE);
417 
418  if (preOrPost_ == "post" || preOrPost_ == "both") {
419  // smooth again on primary space
420  A_->apply (Y, *res1);
421  res1->update (ONE, X, -ONE);
422  vec1->putScalar (ZERO);
423  ifpack2_prec1_->apply (*res1, *vec1);
424  Y.update (ONE, *vec1, ONE);
425  }
426 }
427 
428 template <class MatrixType>
430 {
431  std::ostringstream os;
432 
433  // Output is a valid YAML dictionary in flow style. If you don't
434  // like everything on a single line, you should call describe()
435  // instead.
436  os << "\"Ifpack2::Hiptmair\": {";
437  if (this->getObjectLabel () != "") {
438  os << "Label: \"" << this->getObjectLabel () << "\", ";
439  }
440  os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
441  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
442 
443  if (A_.is_null ()) {
444  os << "Matrix: null";
445  }
446  else {
447  os << "Matrix: not null"
448  << ", Global matrix dimensions: ["
449  << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "]";
450  }
451 
452  os << "}";
453  return os.str ();
454 }
455 
456 
457 template <class MatrixType>
459 describe (Teuchos::FancyOStream &out,
460  const Teuchos::EVerbosityLevel verbLevel) const
461 {
462  using std::endl;
463  using std::setw;
464  using Teuchos::VERB_DEFAULT;
465  using Teuchos::VERB_NONE;
466  using Teuchos::VERB_LOW;
467  using Teuchos::VERB_MEDIUM;
468  using Teuchos::VERB_HIGH;
469  using Teuchos::VERB_EXTREME;
470 
471  const Teuchos::EVerbosityLevel vl =
472  (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
473 
474  if (vl != VERB_NONE) {
475  // describe() always starts with a tab by convention.
476  Teuchos::OSTab tab0 (out);
477  out << "\"Ifpack2::Hiptmair\":";
478 
479  Teuchos::OSTab tab1 (out);
480  if (this->getObjectLabel () != "") {
481  out << "Label: " << this->getObjectLabel () << endl;
482  }
483  out << "Initialized: " << (isInitialized () ? "true" : "false") << endl
484  << "Computed: " << (isComputed () ? "true" : "false") << endl
485  << "Global number of rows: " << A_->getGlobalNumRows () << endl
486  << "Global number of columns: " << A_->getGlobalNumCols () << endl
487  << "Matrix:";
488  if (A_.is_null ()) {
489  out << " null" << endl;
490  } else {
491  A_->describe (out, vl);
492  }
493  }
494 }
495 
496 } // namespace Ifpack2
497 
498 #define IFPACK2_HIPTMAIR_INSTANT(S,LO,GO,N) \
499  template class Ifpack2::Hiptmair< Tpetra::RowMatrix<S, LO, GO, N> >;
500 
501 #endif /* IFPACK2_HIPTMAIR_DEF_HPP */
void initialize()
Do any initialization that depends on the input matrix&#39;s structure.
Definition: Ifpack2_Hiptmair_def.hpp:245
int getNumInitialize() const
Returns the number of calls to Initialize().
Definition: Ifpack2_Hiptmair_def.hpp:209
double getComputeTime() const
Returns the time spent in Compute().
Definition: Ifpack2_Hiptmair_def.hpp:233
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_Hiptmair_decl.hpp:85
Hiptmair(const Teuchos::RCP< const row_matrix_type > &A)
Constructor that takes 1 Tpetra matrix (assumes we&#39;ll get the rest off the parameter list) ...
Definition: Ifpack2_Hiptmair_def.hpp:85
double getApplyTime() const
Returns the time spent in apply().
Definition: Ifpack2_Hiptmair_def.hpp:239
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_Hiptmair_decl.hpp:82
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getRangeMap() const
Tpetra::Map representing the range of this operator.
Definition: Ifpack2_Hiptmair_def.hpp:190
int getNumCompute() const
Returns the number of calls to Compute().
Definition: Ifpack2_Hiptmair_def.hpp:215
void setParameters(const Teuchos::ParameterList &params)
Set the preconditioner&#39;s parameters.
Definition: Ifpack2_Hiptmair_def.hpp:110
Teuchos::RCP< const Tpetra::Map< local_ordinal_type, global_ordinal_type, node_type > > getDomainMap() const
Tpetra::Map representing the domain of this operator.
Definition: Ifpack2_Hiptmair_def.hpp:178
Teuchos::RCP< const Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > > getMatrix() const
Returns a reference to the matrix to be preconditioned.
Definition: Ifpack2_Hiptmair_def.hpp:171
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_Hiptmair_decl.hpp:88
Wrapper for Hiptmair smoothers.
Definition: Ifpack2_Hiptmair_decl.hpp:71
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_Hiptmair_decl.hpp:91
void describe(Teuchos::FancyOStream &out, const Teuchos::EVerbosityLevel verbLevel=Teuchos::Describable::verbLevel_default) const
Print the object with some verbosity level to an FancyOStream object.
Definition: Ifpack2_Hiptmair_def.hpp:459
void compute()
Do any initialization that depends on the input matrix&#39;s values.
Definition: Ifpack2_Hiptmair_def.hpp:283
virtual ~Hiptmair()
Destructor.
Definition: Ifpack2_Hiptmair_def.hpp:107
double getInitializeTime() const
Returns the time spent in Initialize().
Definition: Ifpack2_Hiptmair_def.hpp:227
std::string description() const
Return a simple one-line description of this object.
Definition: Ifpack2_Hiptmair_def.hpp:429
bool hasTransposeApply() const
Whether this object&#39;s apply() method can apply the transpose (or conjugate transpose, if applicable).
Definition: Ifpack2_Hiptmair_def.hpp:201
int getNumApply() const
Returns the number of calls to apply().
Definition: Ifpack2_Hiptmair_def.hpp:221
"Factory" for creating single-level preconditioners.
Definition: Ifpack2_Details_OneLevelFactory_decl.hpp:124
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:73
Teuchos::RCP< prec_type > create(const std::string &precType, const Teuchos::RCP< const row_matrix_type > &matrix) const
Create an instance of Preconditioner given the string name of the preconditioner type.
Definition: Ifpack2_Details_OneLevelFactory_def.hpp:84
Teuchos::RCP< const Teuchos::Comm< int > > getComm() const
Returns the operator&#39;s communicator.
Definition: Ifpack2_Hiptmair_def.hpp:160
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, putting the result in Y.
Definition: Ifpack2_Hiptmair_def.hpp:310