Ifpack2 Templated Preconditioning Package  Version 1.0
Ifpack2_RILUK_decl.hpp
Go to the documentation of this file.
1 /*@HEADER
2 // ***********************************************************************
3 //
4 // Ifpack2: Tempated 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 
45 
46 #ifndef IFPACK2_CRSRILUK_DECL_HPP
47 #define IFPACK2_CRSRILUK_DECL_HPP
48 
51 #include "Tpetra_CrsMatrix_decl.hpp"
52 
53 #include "Ifpack2_ScalingType.hpp"
54 #include "Ifpack2_IlukGraph.hpp"
55 
56 #include <type_traits>
57 
58 namespace Teuchos {
59  class ParameterList; // forward declaration
60 }
61 
62 namespace Ifpack2 {
63 
242 template<class MatrixType>
243 class RILUK:
244  virtual public Ifpack2::Preconditioner<typename MatrixType::scalar_type,
245  typename MatrixType::local_ordinal_type,
246  typename MatrixType::global_ordinal_type,
247  typename MatrixType::node_type>,
248  virtual public Ifpack2::Details::CanChangeMatrix<Tpetra::RowMatrix<typename MatrixType::scalar_type,
249  typename MatrixType::local_ordinal_type,
250  typename MatrixType::global_ordinal_type,
251  typename MatrixType::node_type> >
252 {
253  public:
255  typedef typename MatrixType::scalar_type scalar_type;
256 
258  typedef typename MatrixType::local_ordinal_type local_ordinal_type;
259 
261  typedef typename MatrixType::global_ordinal_type global_ordinal_type;
262 
264  typedef typename MatrixType::node_type node_type;
265 
267  typedef typename Teuchos::ScalarTraits<scalar_type>::magnitudeType magnitude_type;
268 
270  typedef Tpetra::RowMatrix<scalar_type,
271  local_ordinal_type,
272  global_ordinal_type,
273  node_type> row_matrix_type;
274 
275  static_assert(std::is_same<MatrixType, row_matrix_type>::value, "Ifpack2::RILUK: The template parameter MatrixType must be a Tpetra::RowMatrix specialization. Please don't use Tpetra::CrsMatrix (a subclass of Tpetra::RowMatrix) here anymore.");
276 
278  typedef Tpetra::CrsMatrix<scalar_type,
279  local_ordinal_type,
280  global_ordinal_type,
281  node_type> crs_matrix_type;
282 
283  template <class NewMatrixType> friend class RILUK;
284 
288  RILUK (const Teuchos::RCP<const row_matrix_type>& A_in);
289 
297  RILUK (const Teuchos::RCP<const crs_matrix_type>& A_in);
298 
299  private:
302  RILUK (const RILUK<MatrixType> & src);
303 
304  public:
310  template <typename NewMatrixType>
311  Teuchos::RCP< RILUK< NewMatrixType > >
312  clone (const Teuchos::RCP<const NewMatrixType>& A_newnode) const;
313 
315  virtual ~RILUK ();
316 
324  void setParameters (const Teuchos::ParameterList& params);
325 
327  void initialize ();
328 
337  void compute ();
338 
340  bool isInitialized () const {
341  return isInitialized_;
342  }
344  bool isComputed () const {
345  return isComputed_;
346  }
347 
349  int getNumInitialize () const {
350  return numInitialize_;
351  }
353  int getNumCompute () const {
354  return numCompute_;
355  }
357  int getNumApply () const {
358  return numApply_;
359  }
360 
362  double getInitializeTime () const {
363  return initializeTime_;
364  }
366  double getComputeTime () const {
367  return computeTime_;
368  }
370  double getApplyTime () const {
371  return applyTime_;
372  }
373 
375 
376 
399  virtual void
400  setMatrix (const Teuchos::RCP<const row_matrix_type>& A);
401 
403 
405 
407  std::string description () const;
408 
410 
412 
414  Teuchos::RCP<const Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> >
415  getDomainMap () const;
416 
418  Teuchos::RCP<const Tpetra::Map<local_ordinal_type,global_ordinal_type,node_type> >
419  getRangeMap () const;
420 
450  void
451  apply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
452  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
453  Teuchos::ETransp mode = Teuchos::NO_TRANS,
454  scalar_type alpha = Teuchos::ScalarTraits<scalar_type>::one (),
455  scalar_type beta = Teuchos::ScalarTraits<scalar_type>::zero ()) const;
457 
458 private:
480  void
481  multiply (const Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& X,
482  Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>& Y,
483  const Teuchos::ETransp mode = Teuchos::NO_TRANS) const;
484 public:
485 
487  Teuchos::RCP<const row_matrix_type> getMatrix () const;
488 
489  // Attribute access functions
490 
492  magnitude_type getRelaxValue () const { return RelaxValue_; }
493 
495  magnitude_type getAbsoluteThreshold () const { return Athresh_; }
496 
498  magnitude_type getRelativeThreshold () const {return Rthresh_;}
499 
501  int getLevelOfFill () const { return LevelOfFill_; }
502 
504  Tpetra::CombineMode getOverlapMode () {
505  TEUCHOS_TEST_FOR_EXCEPTION(
506  true, std::logic_error, "Ifpack2::RILUK::SetOverlapMode: "
507  "RILUK no longer implements overlap on its own. "
508  "Use RILUK with AdditiveSchwarz if you want overlap.");
509  }
510 
512  Tpetra::global_size_t getGlobalNumEntries () const {
513  return getL ().getGlobalNumEntries () + getU ().getGlobalNumEntries ();
514  }
515 
517  Teuchos::RCP<Ifpack2::IlukGraph<Tpetra::CrsGraph<local_ordinal_type,global_ordinal_type,node_type> > > getGraph () const {
518  return Graph_;
519  }
520 
522  const crs_matrix_type& getL () const;
523 
525  const Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type>&
526  getD () const;
527 
529  const crs_matrix_type& getU () const;
530 
532  Teuchos::RCP<const crs_matrix_type> getCrsMatrix () const;
533 
534 private:
535  typedef Tpetra::MultiVector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> MV;
536  typedef Teuchos::ScalarTraits<scalar_type> STS;
537  typedef Teuchos::ScalarTraits<magnitude_type> STM;
538 
539  void allocate_L_and_U();
540  void initAllValues (const row_matrix_type& A);
541 
547  static Teuchos::RCP<const row_matrix_type>
548  makeLocalFilter (const Teuchos::RCP<const row_matrix_type>& A);
549 
550 protected:
551  typedef Tpetra::Vector<scalar_type,local_ordinal_type,global_ordinal_type,node_type> vec_type;
552 
554  Teuchos::RCP<const row_matrix_type> A_;
555 
557  Teuchos::RCP<Ifpack2::IlukGraph<Tpetra::CrsGraph<local_ordinal_type,
558  global_ordinal_type,
559  node_type> > > Graph_;
568  Teuchos::RCP<const crs_matrix_type> A_local_crs_;
569 
571  Teuchos::RCP<crs_matrix_type> L_;
573  Teuchos::RCP<crs_matrix_type> U_;
575  Teuchos::RCP<vec_type> D_;
576 
577  int LevelOfFill_;
578 
579  bool isAllocated_;
580  bool isInitialized_;
581  bool isComputed_;
582 
583  int numInitialize_;
584  int numCompute_;
585  mutable int numApply_;
586 
587  double initializeTime_;
588  double computeTime_;
589  mutable double applyTime_;
590 
591  magnitude_type RelaxValue_;
592  magnitude_type Athresh_;
593  magnitude_type Rthresh_;
594 };
595 
596 // NOTE (mfh 11 Feb 2015) This used to exist in order to deal with
597 // different behavior of Tpetra::Crs{Graph,Matrix} for
598 // KokkosClassic::ThrustGPUNode. In particular, fillComplete on a
599 // CrsMatrix used to make the graph go away by default, so we had to
600 // pass in a parameter to keep a host copy of the graph. With the new
601 // (Kokkos refactor) version of Tpetra, this problem has gone away.
602 namespace detail {
603  template<class MatrixType, class NodeType>
604  struct setLocalSolveParams{
605  static Teuchos::RCP<Teuchos::ParameterList>
606  setParams (const Teuchos::RCP<Teuchos::ParameterList>& param) {
607  return param;
608  }
609  };
610 } // namespace detail
611 
612 template <class MatrixType>
613 template <typename NewMatrixType>
614 Teuchos::RCP<RILUK<NewMatrixType> >
616 clone (const Teuchos::RCP<const NewMatrixType>& A_newnode) const
617 {
618  using Teuchos::ParameterList;
619  using Teuchos::RCP;
620  using Teuchos::rcp;
621 
622  typedef typename NewMatrixType::scalar_type new_scalar_type;
623  typedef typename NewMatrixType::local_ordinal_type new_local_ordinal_type;
624  typedef typename NewMatrixType::global_ordinal_type new_global_ordinal_type;
625  typedef typename NewMatrixType::node_type new_node_type;
626  typedef Tpetra::RowMatrix<new_scalar_type, new_local_ordinal_type,
627  new_global_ordinal_type, new_node_type> new_row_matrix_type;
628  typedef RILUK<new_row_matrix_type> new_riluk_type;
629 
630  RCP<new_riluk_type> new_riluk = rcp (new new_riluk_type (A_newnode));
631 
632  RCP<ParameterList> plClone = Teuchos::parameterList ();
633  plClone = detail::setLocalSolveParams<NewMatrixType, new_node_type>::setParams (plClone);
634 
635  RCP<new_node_type> new_node = A_newnode->getNode ();
636  new_riluk->L_ = L_->clone (new_node, plClone);
637  new_riluk->U_ = U_->clone (new_node, plClone);
638  new_riluk->D_ = D_->clone (new_node);
639 
640  new_riluk->LevelOfFill_ = LevelOfFill_;
641 
642  new_riluk->isAllocated_ = isAllocated_;
643  new_riluk->isInitialized_ = isInitialized_;
644  new_riluk->isComputed_ = isComputed_;
645 
646  new_riluk->numInitialize_ = numInitialize_;
647  new_riluk->numCompute_ = numCompute_;
648  new_riluk->numApply_ = numApply_;
649 
650  new_riluk->RelaxValue_ = RelaxValue_;
651  new_riluk->Athresh_ = Athresh_;
652  new_riluk->Rthresh_ = Rthresh_;
653 
654  return new_riluk;
655 }
656 
657 } // namespace Ifpack2
658 
659 #endif /* IFPACK2_CRSRILUK_DECL_HPP */
Mix-in interface for preconditioners that can change their matrix after construction.
Definition: Ifpack2_Details_CanChangeMatrix.hpp:93
Declaration and definition of IlukGraph.
magnitude_type getRelaxValue() const
Get RILU(k) relaxation parameter.
Definition: Ifpack2_RILUK_decl.hpp:492
magnitude_type getAbsoluteThreshold() const
Get absolute threshold value.
Definition: Ifpack2_RILUK_decl.hpp:495
Ifpack2::ScalingType enumerable type.
MatrixType::global_ordinal_type global_ordinal_type
The type of global indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:261
Teuchos::RCP< crs_matrix_type > L_
The L (lower triangular) factor of ILU(k).
Definition: Ifpack2_RILUK_decl.hpp:571
Tpetra::RowMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > row_matrix_type
Tpetra::RowMatrix specialization used by this class.
Definition: Ifpack2_RILUK_decl.hpp:273
double getApplyTime() const
Total time in seconds taken by all successful apply() calls for this object.
Definition: Ifpack2_RILUK_decl.hpp:370
ILU(k) factorization of a given Tpetra::RowMatrix.
Definition: Ifpack2_RILUK_decl.hpp:243
Tpetra::CombineMode getOverlapMode()
Get overlap mode type.
Definition: Ifpack2_RILUK_decl.hpp:504
Teuchos::RCP< Ifpack2::IlukGraph< Tpetra::CrsGraph< local_ordinal_type, global_ordinal_type, node_type > > > getGraph() const
Return the Ifpack2::IlukGraph associated with this factored matrix.
Definition: Ifpack2_RILUK_decl.hpp:517
bool isComputed() const
Whether compute() has been called on this object.
Definition: Ifpack2_RILUK_decl.hpp:344
Tpetra::CrsMatrix< scalar_type, local_ordinal_type, global_ordinal_type, node_type > crs_matrix_type
Tpetra::CrsMatrix specialization used by this class for representing L and U.
Definition: Ifpack2_RILUK_decl.hpp:275
MatrixType::node_type node_type
The Node type used by the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:264
MatrixType::local_ordinal_type local_ordinal_type
The type of local indices in the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:258
Teuchos::RCP< const crs_matrix_type > A_local_crs_
The matrix used to to compute ILU(k).
Definition: Ifpack2_RILUK_decl.hpp:568
int getLevelOfFill() const
Get level of fill (the "k" in ILU(k)).
Definition: Ifpack2_RILUK_decl.hpp:501
Teuchos::RCP< Ifpack2::IlukGraph< Tpetra::CrsGraph< local_ordinal_type, global_ordinal_type, node_type > > > Graph_
The ILU(k) graph.
Definition: Ifpack2_RILUK_decl.hpp:559
Interface for all Ifpack2 preconditioners.
Definition: Ifpack2_Preconditioner.hpp:107
magnitude_type getRelativeThreshold() const
Get relative threshold value.
Definition: Ifpack2_RILUK_decl.hpp:498
Teuchos::ScalarTraits< scalar_type >::magnitudeType magnitude_type
The type of the magnitude (absolute value) of a matrix entry.
Definition: Ifpack2_RILUK_decl.hpp:267
Construct a level filled graph for use in computing an ILU(k) incomplete factorization.
Definition: Ifpack2_IlukGraph.hpp:97
Declaration of interface for preconditioners that can change their matrix after construction.
bool isInitialized() const
Whether initialize() has been called on this object.
Definition: Ifpack2_RILUK_decl.hpp:340
Teuchos::RCP< const row_matrix_type > A_
The (original) input matrix for which to compute ILU(k).
Definition: Ifpack2_RILUK_decl.hpp:554
Definition: Ifpack2_Details_Amesos2Wrapper_decl.hpp:63
double getComputeTime() const
Total time in seconds taken by all successful compute() calls for this object.
Definition: Ifpack2_RILUK_decl.hpp:366
Teuchos::RCP< vec_type > D_
The diagonal entries of the ILU(k) factorization.
Definition: Ifpack2_RILUK_decl.hpp:575
int getNumCompute() const
Number of successful compute() calls for this object.
Definition: Ifpack2_RILUK_decl.hpp:353
int getNumApply() const
Number of successful apply() calls for this object.
Definition: Ifpack2_RILUK_decl.hpp:357
Tpetra::global_size_t getGlobalNumEntries() const
Returns the number of nonzero entries in the global graph.
Definition: Ifpack2_RILUK_decl.hpp:512
Preconditioners and smoothers for Tpetra sparse matrices.
Definition: Ifpack2_AdditiveSchwarz_decl.hpp:72
double getInitializeTime() const
Total time in seconds taken by all successful initialize() calls for this object. ...
Definition: Ifpack2_RILUK_decl.hpp:362
Teuchos::RCP< crs_matrix_type > U_
The U (upper triangular) factor of ILU(k).
Definition: Ifpack2_RILUK_decl.hpp:573
int getNumInitialize() const
Number of successful initialize() calls for this object.
Definition: Ifpack2_RILUK_decl.hpp:349
MatrixType::scalar_type scalar_type
The type of the entries of the input MatrixType.
Definition: Ifpack2_RILUK_decl.hpp:255