Reference documentation for deal.II version 8.1.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Modules Pages
trilinos_precondition.h
1 // ---------------------------------------------------------------------
2 // @f$Id: trilinos_precondition.h 31932 2013-12-08 02:15:54Z heister @f$
3 //
4 // Copyright (C) 2008 - 2013 by the deal.II authors
5 //
6 // This file is part of the deal.II library.
7 //
8 // The deal.II library is free software; you can use it, redistribute
9 // it, and/or modify it under the terms of the GNU Lesser General
10 // Public License as published by the Free Software Foundation; either
11 // version 2.1 of the License, or (at your option) any later version.
12 // The full text of the license can be found in the file LICENSE at
13 // the top level of the deal.II distribution.
14 //
15 // ---------------------------------------------------------------------
16 
17 #ifndef __deal2__trilinos_precondition_h
18 #define __deal2__trilinos_precondition_h
19 
20 
21 #include <deal.II/base/config.h>
22 
23 #ifdef DEAL_II_WITH_TRILINOS
24 
25 # include <deal.II/base/subscriptor.h>
26 # include <deal.II/base/std_cxx1x/shared_ptr.h>
27 
28 # include <deal.II/lac/trilinos_vector_base.h>
29 # include <deal.II/lac/parallel_vector.h>
30 
31 # ifdef DEAL_II_WITH_MPI
32 # include <Epetra_MpiComm.h>
33 # else
34 # include <Epetra_SerialComm.h>
35 # endif
36 # include <Epetra_Map.h>
37 
38 # include <Teuchos_ParameterList.hpp>
39 # include <Epetra_Operator.h>
40 # include <Epetra_Vector.h>
41 
42 // forward declarations
43 class Ifpack_Preconditioner;
44 class Ifpack_Chebyshev;
45 namespace ML_Epetra
46 {
47  class MultiLevelPreconditioner;
48 }
49 
50 
52 
53 // forward declarations
54 template <typename number> class SparseMatrix;
55 template <typename number> class Vector;
56 class SparsityPattern;
57 
62 namespace TrilinosWrappers
63 {
64  // forward declarations
65  class SparseMatrix;
66  class BlockSparseMatrix;
67  class SolverBase;
68 
78  {
79  public:
83  typedef ::types::global_dof_index size_type;
84 
91  {};
92 
102  PreconditionBase ();
103 
108 
113 
119  void clear ();
120 
124  virtual void vmult (VectorBase &dst,
125  const VectorBase &src) const;
126 
130  virtual void Tvmult (VectorBase &dst,
131  const VectorBase &src) const;
132 
140  virtual void vmult (::Vector<double> &dst,
141  const ::Vector<double> &src) const;
142 
150  virtual void Tvmult (::Vector<double> &dst,
151  const ::Vector<double> &src) const;
152 
159  virtual void vmult (::parallel::distributed::Vector<double> &dst,
160  const ::parallel::distributed::Vector<double> &src) const;
161 
168  virtual void Tvmult (::parallel::distributed::Vector<double> &dst,
169  const ::parallel::distributed::Vector<double> &src) const;
170 
174  DeclException1 (ExcNonMatchingMaps,
175  std::string,
176  << "The sparse matrix the preconditioner is based on "
177  << "uses a map that is not compatible to the one in vector "
178  << arg1
179  << ". Check preconditioner and matrix setup.");
180 
181  friend class SolverBase;
182  friend class PreconditionStokes;
183 
184  protected:
191  std_cxx1x::shared_ptr<Epetra_Operator> preconditioner;
192 
199 #ifdef DEAL_II_WITH_MPI
200  Epetra_MpiComm communicator;
201 #else
202  Epetra_SerialComm communicator;
203 #endif
204 
210  std_cxx1x::shared_ptr<Epetra_Map> vector_distributor;
211  };
212 
213 
232  {
233  public:
234 
261  {
268  AdditionalData (const double omega = 1,
269  const double min_diagonal = 0);
270 
276  double omega;
277 
290  double min_diagonal;
291  };
292 
300  void initialize (const SparseMatrix &matrix,
301  const AdditionalData &additional_data = AdditionalData());
302  };
303 
304 
305 
306 
336  {
337  public:
338 
371  {
383  AdditionalData (const double omega = 1,
384  const double min_diagonal = 0,
385  const unsigned int overlap = 0);
386 
392  double omega;
393 
407  double min_diagonal;
408 
416  unsigned int overlap;
417  };
418 
428  void initialize (const SparseMatrix &matrix,
429  const AdditionalData &additional_data = AdditionalData());
430  };
431 
432 
433 
434 
464  {
465  public:
466 
499  {
511  AdditionalData (const double omega = 1,
512  const double min_diagonal = 0,
513  const unsigned int overlap = 0);
514 
520  double omega;
521 
535  double min_diagonal;
536 
544  unsigned int overlap;
545  };
546 
556  void initialize (const SparseMatrix &matrix,
557  const AdditionalData &additional_data = AdditionalData());
558  };
559 
560 
561 
601  {
602  public:
643  {
663  AdditionalData (const unsigned int ic_fill = 0,
664  const double ic_atol = 0.,
665  const double ic_rtol = 1.,
666  const unsigned int overlap = 0);
667 
687  unsigned int ic_fill;
688 
697  double ic_atol;
698 
706  double ic_rtol;
707 
715  unsigned int overlap;
716  };
717 
725  void initialize (const SparseMatrix &matrix,
726  const AdditionalData &additional_data = AdditionalData());
727  };
728 
729 
730 
770  {
771  public:
812  {
831  AdditionalData (const unsigned int ilu_fill = 0,
832  const double ilu_atol = 0.,
833  const double ilu_rtol = 1.,
834  const unsigned int overlap = 0);
835 
855  unsigned int ilu_fill;
856 
865  double ilu_atol;
866 
874  double ilu_rtol;
875 
883  unsigned int overlap;
884  };
885 
893  void initialize (const SparseMatrix &matrix,
894  const AdditionalData &additional_data = AdditionalData());
895  };
896 
897 
898 
899 
900 
901 
940  {
941  public:
978  {
996  AdditionalData (const double ilut_drop = 0.,
997  const unsigned int ilut_fill = 0,
998  const double ilut_atol = 0.,
999  const double ilut_rtol = 1.,
1000  const unsigned int overlap = 0);
1001 
1008  double ilut_drop;
1009 
1029  unsigned int ilut_fill;
1030 
1039  double ilut_atol;
1040 
1048  double ilut_rtol;
1049 
1057  unsigned int overlap;
1058  };
1059 
1067  void initialize (const SparseMatrix &matrix,
1068  const AdditionalData &additional_data = AdditionalData());
1069  };
1070 
1071 
1072 
1094  {
1095  public:
1102  {
1106  AdditionalData (const unsigned int overlap = 0);
1107 
1108 
1116  unsigned int overlap;
1117  };
1118 
1126  void initialize (const SparseMatrix &matrix,
1127  const AdditionalData &additional_data = AdditionalData());
1128  };
1129 
1130 
1131 
1132 
1133 
1134 
1146  {
1147  public:
1154  {
1158  AdditionalData (const unsigned int degree = 1,
1159  const double max_eigenvalue = 10.,
1160  const double eigenvalue_ratio = 30.,
1161  const double min_eigenvalue = 1.,
1162  const double min_diagonal = 1e-12,
1163  const bool nonzero_starting = false);
1164 
1173  unsigned int degree;
1174 
1183 
1190 
1198 
1206 
1226  };
1227 
1235  void initialize (const SparseMatrix &matrix,
1236  const AdditionalData &additional_data = AdditionalData());
1237  };
1238 
1239 
1240 
1284  {
1285  public:
1286 
1294  {
1301  AdditionalData (const bool elliptic = true,
1302  const bool higher_order_elements = false,
1303  const unsigned int n_cycles = 1,
1304  const bool w_cyle = false,
1305  const double aggregation_threshold = 1e-4,
1306  const std::vector<std::vector<bool> > &constant_modes = std::vector<std::vector<bool> > (1),
1307  const unsigned int smoother_sweeps = 2,
1308  const unsigned int smoother_overlap = 0,
1309  const bool output_details = false);
1310 
1322  bool elliptic;
1323 
1331 
1337  unsigned int n_cycles;
1338 
1344  bool w_cycle;
1345 
1363 
1373  std::vector<std::vector<bool> > constant_modes;
1374 
1393  unsigned int smoother_sweeps;
1394 
1400  unsigned int smoother_overlap;
1401 
1411  };
1412 
1421  void initialize (const SparseMatrix &matrix,
1422  const AdditionalData &additional_data = AdditionalData());
1423 
1444  void initialize (const SparseMatrix &matrix,
1445  const Teuchos::ParameterList &ml_parameters);
1446 
1457  template <typename number>
1458  void initialize (const ::SparseMatrix<number> &deal_ii_sparse_matrix,
1459  const AdditionalData &additional_data = AdditionalData(),
1460  const double drop_tolerance = 1e-13,
1461  const ::SparsityPattern *use_this_sparsity = 0);
1462 
1486  void reinit ();
1487 
1493  void clear ();
1494 
1499  size_type memory_consumption () const;
1500 
1501  private:
1506  std_cxx1x::shared_ptr<SparseMatrix> trilinos_matrix;
1507  };
1508 
1509 
1510 
1519  {
1520  public:
1521 
1525  void vmult (VectorBase &dst,
1526  const VectorBase &src) const;
1527 
1531  void Tvmult (VectorBase &dst,
1532  const VectorBase &src) const;
1533 
1539  void vmult (::Vector<double> &dst,
1540  const ::Vector<double> &src) const;
1541 
1547  void Tvmult (::Vector<double> &dst,
1548  const ::Vector<double> &src) const;
1549 
1556  const ::parallel::distributed::Vector<double> &src) const;
1557 
1564  const ::parallel::distributed::Vector<double> &src) const;
1565  };
1566 
1567 
1568 
1569 // -------------------------- inline and template functions ----------------------
1570 
1571 
1572 #ifndef DOXYGEN
1573 
1574  inline
1575  void
1577  const VectorBase &src) const
1578  {
1579  Assert (dst.vector_partitioner().SameAs(preconditioner->OperatorRangeMap()),
1580  ExcNonMatchingMaps("dst"));
1581  Assert (src.vector_partitioner().SameAs(preconditioner->OperatorDomainMap()),
1582  ExcNonMatchingMaps("src"));
1583 
1584  const int ierr = preconditioner->ApplyInverse (src.trilinos_vector(),
1585  dst.trilinos_vector());
1586  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1587  }
1588 
1589  inline
1590  void
1591  PreconditionBase::Tvmult (VectorBase &dst,
1592  const VectorBase &src) const
1593  {
1594  Assert (dst.vector_partitioner().SameAs(preconditioner->OperatorRangeMap()),
1595  ExcNonMatchingMaps("dst"));
1596  Assert (src.vector_partitioner().SameAs(preconditioner->OperatorDomainMap()),
1597  ExcNonMatchingMaps("src"));
1598 
1599  preconditioner->SetUseTranspose(true);
1600  const int ierr = preconditioner->ApplyInverse (src.trilinos_vector(),
1601  dst.trilinos_vector());
1602  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1603  preconditioner->SetUseTranspose(false);
1604  }
1605 
1606 
1607  // For the implementation of
1608  // the <code>vmult</code>
1609  // function with deal.II data
1610  // structures we note that
1611  // invoking a call of the
1612  // Trilinos preconditioner
1613  // requires us to use Epetra
1614  // vectors as well. We do this
1615  // by providing a view, i.e.,
1616  // feed Trilinos with a
1617  // pointer to the data, so we
1618  // avoid copying the content
1619  // of the vectors during the
1620  // iteration (this function is
1621  // only useful when used in
1622  // serial anyway). In the
1623  // declaration of the right
1624  // hand side, we need to cast
1625  // the source vector (that is
1626  // <code>const</code> in all
1627  // deal.II calls) to
1628  // non-constant value, as this
1629  // is the way Trilinos wants
1630  // to have them.
1631  inline
1633  const ::Vector<double> &src) const
1634  {
1635  AssertDimension (static_cast<TrilinosWrappers::types::int_type>(dst.size()),
1636  preconditioner->OperatorDomainMap().NumMyElements());
1637  AssertDimension (static_cast<TrilinosWrappers::types::int_type>(src.size()),
1638  preconditioner->OperatorRangeMap().NumMyElements());
1639  Epetra_Vector tril_dst (View, preconditioner->OperatorDomainMap(),
1640  dst.begin());
1641  Epetra_Vector tril_src (View, preconditioner->OperatorRangeMap(),
1642  const_cast<double *>(src.begin()));
1643 
1644  const int ierr = preconditioner->ApplyInverse (tril_src, tril_dst);
1645  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1646  }
1647 
1648 
1649  inline
1651  const ::Vector<double> &src) const
1652  {
1653  AssertDimension (static_cast<TrilinosWrappers::types::int_type>(dst.size()),
1654  preconditioner->OperatorDomainMap().NumMyElements());
1655  AssertDimension (static_cast<TrilinosWrappers::types::int_type>(src.size()),
1656  preconditioner->OperatorRangeMap().NumMyElements());
1657  Epetra_Vector tril_dst (View, preconditioner->OperatorDomainMap(),
1658  dst.begin());
1659  Epetra_Vector tril_src (View, preconditioner->OperatorRangeMap(),
1660  const_cast<double *>(src.begin()));
1661 
1662  preconditioner->SetUseTranspose(true);
1663  const int ierr = preconditioner->ApplyInverse (tril_src, tril_dst);
1664  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1665  preconditioner->SetUseTranspose(false);
1666  }
1667 
1668 
1669 
1670  inline
1671  void
1673  const parallel::distributed::Vector<double> &src) const
1674  {
1675  AssertDimension (static_cast<TrilinosWrappers::types::int_type>(dst.local_size()),
1676  preconditioner->OperatorDomainMap().NumMyElements());
1677  AssertDimension (static_cast<TrilinosWrappers::types::int_type>(src.local_size()),
1678  preconditioner->OperatorRangeMap().NumMyElements());
1679  Epetra_Vector tril_dst (View, preconditioner->OperatorDomainMap(),
1680  dst.begin());
1681  Epetra_Vector tril_src (View, preconditioner->OperatorRangeMap(),
1682  const_cast<double *>(src.begin()));
1683 
1684  const int ierr = preconditioner->ApplyInverse (tril_src, tril_dst);
1685  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1686  }
1687 
1688  inline
1689  void
1691  const parallel::distributed::Vector<double> &src) const
1692  {
1693  AssertDimension (static_cast<TrilinosWrappers::types::int_type>(dst.local_size()),
1694  preconditioner->OperatorDomainMap().NumMyElements());
1695  AssertDimension (static_cast<TrilinosWrappers::types::int_type>(src.local_size()),
1696  preconditioner->OperatorRangeMap().NumMyElements());
1697  Epetra_Vector tril_dst (View, preconditioner->OperatorDomainMap(),
1698  dst.begin());
1699  Epetra_Vector tril_src (View, preconditioner->OperatorRangeMap(),
1700  const_cast<double *>(src.begin()));
1701 
1702  preconditioner->SetUseTranspose(true);
1703  const int ierr = preconditioner->ApplyInverse (tril_src, tril_dst);
1704  AssertThrow (ierr == 0, ExcTrilinosError(ierr));
1705  preconditioner->SetUseTranspose(false);
1706  }
1707 
1708 #endif
1709 
1710 }
1711 
1712 
1716 DEAL_II_NAMESPACE_CLOSE
1717 
1718 #endif // DEAL_II_WITH_TRILINOS
1719 
1720 /*---------------------------- trilinos_precondition.h ---------------------------*/
1721 
1722 #endif
1723 /*---------------------------- trilinos_precondition.h ---------------------------*/
AdditionalData(const unsigned int degree=1, const double max_eigenvalue=10., const double eigenvalue_ratio=30., const double min_eigenvalue=1., const double min_diagonal=1e-12, const bool nonzero_starting=false)
#define AssertDimension(dim1, dim2)
Definition: exceptions.h:858
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
AdditionalData(const double omega=1, const double min_diagonal=0, const unsigned int overlap=0)
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
virtual void Tvmult(VectorBase &dst, const VectorBase &src) const
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
#define AssertThrow(cond, exc)
Definition: exceptions.h:362
AdditionalData(const unsigned int ilu_fill=0, const double ilu_atol=0., const double ilu_rtol=1., const unsigned int overlap=0)
AdditionalData(const double omega=1, const double min_diagonal=0)
AdditionalData(const bool elliptic=true, const bool higher_order_elements=false, const unsigned int n_cycles=1, const bool w_cyle=false, const double aggregation_threshold=1e-4, const std::vector< std::vector< bool > > &constant_modes=std::vector< std::vector< bool > >(1), const unsigned int smoother_sweeps=2, const unsigned int smoother_overlap=0, const bool output_details=false)
#define Assert(cond, exc)
Definition: exceptions.h:299
DeclException1(ExcNonMatchingMaps, std::string,<< "The sparse matrix the preconditioner is based on "<< "uses a map that is not compatible to the one in vector "<< arg1<< ". Check preconditioner and matrix setup.")
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void Tvmult(VectorBase &dst, const VectorBase &src) const
std_cxx1x::shared_ptr< Epetra_Map > vector_distributor
AdditionalData(const double ilut_drop=0., const unsigned int ilut_fill=0, const double ilut_atol=0., const double ilut_rtol=1., const unsigned int overlap=0)
const Epetra_Map & vector_partitioner() const
iterator begin()
AdditionalData(const double omega=1, const double min_diagonal=0, const unsigned int overlap=0)
std_cxx1x::shared_ptr< SparseMatrix > trilinos_matrix
void vmult(VectorBase &dst, const VectorBase &src) const
std_cxx1x::shared_ptr< Epetra_Operator > preconditioner
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
size_type local_size() const
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
void initialize(const SparseMatrix &matrix, const AdditionalData &additional_data=AdditionalData())
AdditionalData(const unsigned int ic_fill=0, const double ic_atol=0., const double ic_rtol=1., const unsigned int overlap=0)
virtual void vmult(VectorBase &dst, const VectorBase &src) const
const Epetra_MultiVector & trilinos_vector() const
std::size_t size() const
size_type memory_consumption() const