Ifpack2 Templated Preconditioning Package  Version 1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
Ifpack2_Krylov_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_KRYLOV_DEF_HPP
44 #define IFPACK2_KRYLOV_DEF_HPP
45 
46 #ifdef HAVE_IFPACK2_DEPRECATED_CODE
47 
48 #include "Ifpack2_Chebyshev.hpp"
49 #include "Ifpack2_Heap.hpp"
50 #include "Ifpack2_ILUT.hpp"
51 #include "Ifpack2_Parameters.hpp"
52 #include "Ifpack2_Relaxation.hpp"
53 #include "Ifpack2_RILUK.hpp"
54 
56 #include "BelosBlockCGSolMgr.hpp"
57 
58 #include "Teuchos_Assert.hpp"
59 #include "Teuchos_Time.hpp"
60 
61 #include <iostream>
62 #include <sstream>
63 #include <cmath>
64 
65 
66 namespace Ifpack2 {
67 namespace DeprecatedAndMayDisappearAtAnyTime {
68 
69 template <class MatrixType>
70 Krylov<MatrixType>::
71 Krylov (const Teuchos::RCP<const row_matrix_type>& A) :
72  A_ (A),
73  // Default values
74  iterationType_ ("GMRES"),
75  numIters_ (5),
76  resTol_ (0.001), // FIXME (mfh 17 Jan 2014) Make a function of STS::eps()
77  BlockSize_ (1),
78  ZeroStartingSolution_ (true),
79  PreconditionerType_ (1),
80  // General
81  IsInitialized_ (false),
82  IsComputed_ (false),
83  NumInitialize_ (0),
84  NumCompute_ (0),
85  NumApply_ (0),
86  InitializeTime_ (0.0),
87  ComputeTime_ (0.0),
88  ApplyTime_ (0.0)
89 {}
90 
91 
92 template <class MatrixType>
93 Krylov<MatrixType>::~Krylov() {}
94 
95 
96 template <class MatrixType>
97 void Krylov<MatrixType>::setMatrix (const Teuchos::RCP<const row_matrix_type>& A)
98 {
99  // Check in serial or one-process mode if the matrix is square.
101  ! A.is_null () && A->getComm ()->getSize () == 1 &&
102  A->getNodeNumRows () != A->getNodeNumCols (),
103  std::runtime_error, "Ifpack2::Krylov::setMatrix: If A's communicator only "
104  "contains one process, then A must be square. Instead, you provided a "
105  "matrix A with " << A->getNodeNumRows () << " rows and "
106  << A->getNodeNumCols () << " columns.");
107 
108  // It's legal for A to be null; in that case, you may not call
109  // initialize() until calling setMatrix() with a nonnull input.
110  // Regardless, setting the matrix invalidates any previous
111  // factorization.
112  IsInitialized_ = false;
113  IsComputed_ = false;
114 
115  A_ = A;
116 }
117 
118 
119 template <class MatrixType>
120 void Krylov<MatrixType>::setParameters (const Teuchos::ParameterList& plist)
121 {
122  using Teuchos::as;
126 
127  // FIXME (mfh 12 Sep 2014) Don't rewrite Belos::SolverFactory!!! Use
128  // that instead.
129 
130  ParameterList params = plist;
131 
132  // Get the current parameters' values. We don't assign to the
133  // instance data directly until we've gotten all the parameters.
134  // This ensures "transactional" semantics, so that if attempting to
135  // get some parameter throws an exception, the class' state doesn't
136  // change.
137  magnitude_type resTol = resTol_;
138  int numIters = numIters_;
139  std::string iterType = iterationType_;
140  int blockSize = BlockSize_;
141  bool zeroStartingSolution = ZeroStartingSolution_;
142  int precType = PreconditionerType_;
143 
144  //
145  // Get the "krylov: iteration type" parameter.
146  //
147  // We prefer std::string (name of Belos solver), but allow int
148  // ("enum" value) for backwards compatibility.
149  //
150  bool gotIterType = false;
151  try {
152  iterType = params.get<std::string> ("krylov: iteration type");
153  gotIterType = true;
154  }
155  catch (InvalidParameterName) {
156  gotIterType = true; // the parameter is not there, so don't try to get it anymore
157  }
158  catch (InvalidParameterType) {
159  // Perhaps the user specified it as int.
160  }
161  // If it's not string, it has to be int.
162  // We've already checked whether the name exists.
163  if (! gotIterType) {
164  const int iterTypeInt = params.get<int> ("krylov: iteration type");
165  gotIterType = true;
166 
167  if (iterTypeInt == 1) {
168  iterType = "GMRES";
169  } else if (iterTypeInt == 2) {
170  iterType = "CG";
171  } else {
173  true, std::invalid_argument, "Ifpack2::Krylov::setParameters: Invalid "
174  "\"krylov: iteration type\" value " << iterTypeInt << ". Valid int "
175  "values are 1 (GMRES) and 2 (CG). Please prefer setting this "
176  "parameter as a string (the name of the Krylov solver to use).");
177  }
178  }
179 
180  resTol = params.get ("krylov: residual tolerance", resTol);
181  numIters = params.get ("krylov: number of iterations", numIters);
182  blockSize = params.get ("krylov: block size", blockSize);
183  zeroStartingSolution = params.get ("krylov: zero starting solution",
184  zeroStartingSolution);
185  precType = params.get ("krylov: preconditioner type", precType);
186 
187  // Separate preconditioner parameters into another list
188  //
189  // FIXME (mfh 17 Jan 2014) Inner preconditioner's parameters should
190  // be a sublist, not part of the main list!!!
191  if (PreconditionerType_ == 1) {
192  precParams_.set ("relaxation: sweeps",
193  params.get ("relaxation: sweeps", 1));
194  precParams_.set ("relaxation: damping factor",
195  params.get ("relaxation: damping factor", (scalar_type) 1.0));
196  precParams_.set ("relaxation: min diagonal value",
197  params.get ("relaxation: min diagonal value", STS::one ()));
198  precParams_.set ("relaxation: zero starting solution",
199  params.get ("relaxation: zero starting solution", true));
200  precParams_.set ("relaxation: backward mode",
201  params.get ("relaxation: backward mode", false));
202  }
203  // FIXME (mfh 17 Jan 2014) AdditiveSchwarz's ParameterList no longer
204  // takes parameters for its subdomain solver! You have to pass them
205  // into a sublist.
206  if (PreconditionerType_ == 2 || PreconditionerType_ == 3) {
207  // FIXME (mfh 17 Jan 2014) should be an integer, given how ILUT
208  // works! Furthermore, this parameter does not mean what you
209  // think it means.
210  precParams_.set ("fact: ilut level-of-fill",
211  params.get ("fact: ilut level-of-fill", (double) 1.0));
212  precParams_.set ("fact: iluk level-of-fill",
213  params.get ("fact: iluk level-of-fill", (double) 1.0));
214  // FIXME (mfh 17 Jan 2014) scalar_type or magnitude_type? not
215  // sure, but double is definitely wrong.
216  precParams_.set ("fact: absolute threshold",
217  params.get ("fact: absolute threshold", (double) 0.0));
218  // FIXME (mfh 17 Jan 2014) scalar_type or magnitude_type? not
219  // sure, but double is definitely wrong.
220  precParams_.set ("fact: relative threshold",
221  params.get("fact: relative threshold", (double) 1.0));
222  // FIXME (mfh 17 Jan 2014) scalar_type or magnitude_type? not
223  // sure, but double is definitely wrong.
224  precParams_.set ("fact: relax value",
225  params.get ("fact: relax value", (double) 0.0));
226  }
227 
228  // "Commit" the new values to the instance data.
229  iterationType_ = iterType;
230  numIters_ = numIters;
231  resTol_ = resTol;
232  BlockSize_ = blockSize;
233  ZeroStartingSolution_ = zeroStartingSolution;
234  PreconditionerType_ = precType;
235 }
236 
237 
238 template <class MatrixType>
240 Krylov<MatrixType>::getComm () const {
242  A_.is_null (), std::runtime_error, "Ifpack2::Krylov::getComm: "
243  "The input matrix A is null. Please call setMatrix() with a nonnull "
244  "input matrix before calling this method.");
245  return A_->getComm ();
246 }
247 
248 
249 template <class MatrixType>
251 Krylov<MatrixType>::getMatrix () const {
252  return A_;
253 }
254 
255 
256 template <class MatrixType>
258 Krylov<MatrixType>::getDomainMap () const
259 {
261  A_.is_null (), std::runtime_error, "Ifpack2::Krylov::getDomainMap: "
262  "The input matrix A is null. Please call setMatrix() with a nonnull "
263  "input matrix before calling this method.");
264  return A_->getDomainMap ();
265 }
266 
267 
268 template <class MatrixType>
270 Krylov<MatrixType>::getRangeMap () const
271 {
273  A_.is_null (), std::runtime_error, "Ifpack2::Krylov::getRangeMap: "
274  "The input matrix A is null. Please call setMatrix() with a nonnull "
275  "input matrix before calling this method.");
276  return A_->getRangeMap ();
277 }
278 
279 
280 template <class MatrixType>
281 bool Krylov<MatrixType>::hasTransposeApply () const {
282  // FIXME (mfh 17 Jan 2014) apply() does not currently work with mode
283  // != NO_TRANS, so it's correct to return false here.
284  return false;
285 }
286 
287 
288 template <class MatrixType>
289 int Krylov<MatrixType>::getNumInitialize () const {
290  return NumInitialize_;
291 }
292 
293 
294 template <class MatrixType>
295 int Krylov<MatrixType>::getNumCompute () const {
296  return NumCompute_;
297 }
298 
299 
300 template <class MatrixType>
301 int Krylov<MatrixType>::getNumApply () const {
302  return NumApply_;
303 }
304 
305 
306 template <class MatrixType>
307 double Krylov<MatrixType>::getInitializeTime () const {
308  return InitializeTime_;
309 }
310 
311 
312 template <class MatrixType>
313 double Krylov<MatrixType>::getComputeTime () const {
314  return ComputeTime_;
315 }
316 
317 
318 template <class MatrixType>
319 double Krylov<MatrixType>::getApplyTime () const {
320  return ApplyTime_;
321 }
322 
323 
324 template <class MatrixType>
325 void Krylov<MatrixType>::initialize ()
326 {
328  using Teuchos::RCP;
329  using Teuchos::rcp;
330  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
331  global_ordinal_type, node_type> TMV;
332  typedef Tpetra::Operator<scalar_type, local_ordinal_type,
333  global_ordinal_type, node_type> TOP;
334 
336  A_.is_null (), std::runtime_error, "Ifpack2::Krylov::initialize: "
337  "The input matrix A is null. Please call setMatrix() with a nonnull "
338  "input matrix before calling this method.");
339 
340  // clear any previous allocation
341  IsInitialized_ = false;
342  IsComputed_ = false;
343 
344  Teuchos::Time timer ("initialize");
345  { // The body of code to time
346  Teuchos::TimeMonitor timeMon (timer);
347 
348  // Belos parameter list
349  RCP<ParameterList> belosList = rcp (new ParameterList ("GMRES"));
350  belosList->set ("Maximum Iterations", numIters_);
351  belosList->set ("Convergence Tolerance", resTol_);
352 
353  // FIXME (17 Jan 2014) This whole "preconditioner type" thing is
354  // not how we want Krylov to initialize its inner preconditioner.
355  // Krylov should be initialized like AdditiveSchwarz: the Factory
356  // should create it, in order to avoid circular dependencies.
357 
358  if (PreconditionerType_ == 0) {
359  // no preconditioner
360  }
361  else if (PreconditionerType_==1) {
362  ifpack2_prec_=rcp (new Relaxation<row_matrix_type> (A_));
363  }
364  else if (PreconditionerType_==2) {
365  ifpack2_prec_=rcp (new ILUT<row_matrix_type> (A_));
366  }
367  else if (PreconditionerType_==3) {
368  ifpack2_prec_ = rcp (new RILUK<row_matrix_type> (A_));
369  }
370  else if (PreconditionerType_==4) {
371  ifpack2_prec_ = rcp (new Chebyshev<row_matrix_type> (A_));
372  }
373  if (PreconditionerType_>0) {
374  ifpack2_prec_->initialize();
375  ifpack2_prec_->setParameters(precParams_);
376  }
377  belosProblem_ = rcp (new Belos::LinearProblem<belos_scalar_type,TMV,TOP> ());
378  belosProblem_->setOperator (A_);
379 
380  if (iterationType_ == "GMRES") {
381  belosSolver_ =
382  rcp (new Belos::BlockGmresSolMgr<belos_scalar_type,TMV,TOP> (belosProblem_, belosList));
383  }
384  else {
385  belosSolver_ =
386  rcp (new Belos::BlockCGSolMgr<belos_scalar_type,TMV,TOP> (belosProblem_, belosList));
387  }
388 
389  }
390  IsInitialized_ = true;
391  ++NumInitialize_;
392  InitializeTime_ += timer.totalElapsedTime ();
393 }
394 
395 
396 template <class MatrixType>
397 void Krylov<MatrixType>::compute ()
398 {
400  A_.is_null (), std::runtime_error, "Ifpack2::Krylov::compute: "
401  "The input matrix A is null. Please call setMatrix() with a nonnull "
402  "input matrix before calling this method.");
403 
404  // Don't time the initialize(); that gets timed separately.
405  if (! isInitialized ()) {
406  initialize ();
407  }
408 
409  Teuchos::Time timer ("compute");
410  { // The body of code to time
411  Teuchos::TimeMonitor timeMon (timer);
412  if (PreconditionerType_ > 0) {
413  ifpack2_prec_->compute ();
414  belosProblem_->setLeftPrec (ifpack2_prec_);
415  }
416  }
417  IsComputed_ = true;
418  ++NumCompute_;
419  ComputeTime_ += timer.totalElapsedTime ();
420 }
421 
422 
423 template <class MatrixType>
424 void Krylov<MatrixType>::
425 apply (const Tpetra::MultiVector<typename MatrixType::scalar_type,
426  typename MatrixType::local_ordinal_type,
427  typename MatrixType::global_ordinal_type,
428  typename MatrixType::node_type>& X,
429  Tpetra::MultiVector<typename MatrixType::scalar_type,
430  typename MatrixType::local_ordinal_type,
431  typename MatrixType::global_ordinal_type,
432  typename MatrixType::node_type>& Y,
433  Teuchos::ETransp mode,
434  typename MatrixType::scalar_type alpha,
435  typename MatrixType::scalar_type beta) const
436 {
437  using Teuchos::RCP;
438  using Teuchos::rcp;
439  using Teuchos::rcpFromRef;
440  typedef Tpetra::MultiVector<scalar_type, local_ordinal_type,
441  global_ordinal_type, node_type> MV;
443  ! isComputed (), std::runtime_error,
444  "Ifpack2::Krylov::apply: You must call compute() before you may call apply().");
446  X.getNumVectors () != Y.getNumVectors (), std::invalid_argument,
447  "Ifpack2::Krylov::apply: The MultiVector inputs X and Y do not have the "
448  "same number of columns. X.getNumVectors() = " << X.getNumVectors ()
449  << " != Y.getNumVectors() = " << Y.getNumVectors () << ".");
450 
451  // Catch unimplemented cases: alpha != 1, beta != 0, mode != NO_TRANS.
453  alpha != STS::one (), std::logic_error,
454  "Ifpack2::Krylov::apply: alpha != 1 has not been implemented.");
456  beta != STS::zero (), std::logic_error,
457  "Ifpack2::Krylov::apply: zero != 0 has not been implemented.");
459  mode != Teuchos::NO_TRANS, std::logic_error,
460  "Ifpack2::Krylov::apply: mode != Teuchos::NO_TRANS has not been implemented.");
461 
462  Teuchos::Time timer ("apply");
463  { // The body of code to time
464  Teuchos::TimeMonitor timeMon (timer);
465 
466  // If X and Y are pointing to the same memory location,
467  // we need to create an auxiliary vector, Xcopy
468  RCP<const MV> Xcopy;
469  {
470  auto X_lcl_host = X.template getLocalView<Kokkos::HostSpace> ();
471  auto Y_lcl_host = Y.template getLocalView<Kokkos::HostSpace> ();
472  if (X_lcl_host.data () == Y_lcl_host.data ()) {
473  Xcopy = rcp (new MV (X, Teuchos::Copy));
474  } else {
475  Xcopy = rcpFromRef (X);
476  }
477  }
478 
479  RCP<MV> Ycopy = rcpFromRef (Y);
480  if (ZeroStartingSolution_) {
481  Ycopy->putScalar (STS::zero ());
482  }
483 
484  // Set left and right hand sides for Belos
485  belosProblem_->setProblem (Ycopy, Xcopy);
486  belosSolver_->solve (); // solve the linear system
487  }
488  ++NumApply_;
489  ApplyTime_ += timer.totalElapsedTime ();
490 }
491 
492 
493 template <class MatrixType>
494 std::string Krylov<MatrixType>::description () const
495 {
496  std::ostringstream os;
497 
498  // Output is a valid YAML dictionary in flow style. If you don't
499  // like everything on a single line, you should call describe()
500  // instead.
501  os << "\"Ifpack2::Krylov\": {";
502  if (this->getObjectLabel () != "") {
503  os << "Label: \"" << this->getObjectLabel () << "\", ";
504  }
505  os << "Initialized: " << (isInitialized () ? "true" : "false") << ", "
506  << "Computed: " << (isComputed () ? "true" : "false") << ", ";
507 
508  if (A_.is_null ()) {
509  os << "Matrix: null";
510  }
511  else {
512  os << "Matrix: not null"
513  << ", Global matrix dimensions: ["
514  << A_->getGlobalNumRows () << ", " << A_->getGlobalNumCols () << "]";
515  }
516 
517  os << "}";
518  return os.str ();
519 }
520 
521 
522 template <class MatrixType>
523 void Krylov<MatrixType>::
524 describe (Teuchos::FancyOStream &out,
525  const Teuchos::EVerbosityLevel verbLevel) const
526 {
527  using std::endl;
528  using std::setw;
529  using Teuchos::VERB_DEFAULT;
530  using Teuchos::VERB_NONE;
531  using Teuchos::VERB_LOW;
532  using Teuchos::VERB_MEDIUM;
533  using Teuchos::VERB_HIGH;
534  using Teuchos::VERB_EXTREME;
535 
536  const Teuchos::EVerbosityLevel vl =
537  (verbLevel == VERB_DEFAULT) ? VERB_LOW : verbLevel;
538 
539  if (vl != VERB_NONE) {
540  // describe() always starts with a tab by convention.
541  Teuchos::OSTab tab0 (out);
542  out << "\"Ifpack2::Krylov\":";
543 
544  Teuchos::OSTab tab1 (out);
545  if (this->getObjectLabel () != "") {
546  out << "Label: " << this->getObjectLabel () << endl;
547  }
548  out << "Initialized: " << (isInitialized () ? "true" : "false") << endl
549  << "Computed: " << (isComputed () ? "true" : "false") << endl
550  << "Global number of rows: " << A_->getGlobalNumRows () << endl
551  << "Global number of columns: " << A_->getGlobalNumCols () << endl
552  << "Matrix:";
553  if (A_.is_null ()) {
554  out << " null" << endl;
555  } else {
556  A_->describe (out, vl);
557  }
558  }
559 }
560 
561 } // namespace DeprecatedAndMayDisappearAtAnyTime
562 } // namespace Ifpack2
563 
564 // There's no need to instantiate for CrsMatrix too. All Ifpack2
565 // preconditioners can and should do dynamic casts if they need a type
566 // more specific than RowMatrix.
567 //
568 // In fact, Krylov really doesn't _need_ the RowMatrix methods; it
569 // could very well just rely on the Operator interface. initialize()
570 // need only check whether the domain and range Maps have changed
571 // (which would necessitate reinitializing the Krylov solver).
572 
573 #define IFPACK2_KRYLOV_INSTANT(S,LO,GO,N) \
574  template class Ifpack2::DeprecatedAndMayDisappearAtAnyTime::Krylov< Tpetra::RowMatrix<S, LO, GO, N> >;
575 
576 #endif // HAVE_IFPACK2_DEPRECATED_CODE
577 #endif /* IFPACK2_KRYLOV_DEF_HPP */
T & get(const std::string &name, T def_value)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
TypeTo as(const TypeFrom &t)
bool is_null() const