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