00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012 #ifndef __deal2__trilinos_precondition_h
00013 #define __deal2__trilinos_precondition_h
00014
00015
00016 #include <deal.II/base/config.h>
00017
00018 #ifdef DEAL_II_USE_TRILINOS
00019
00020 # include <deal.II/base/subscriptor.h>
00021 # include <deal.II/base/std_cxx1x/shared_ptr.h>
00022
00023 # include <deal.II/lac/trilinos_vector_base.h>
00024 # include <deal.II/lac/parallel_vector.h>
00025
00026 # ifdef DEAL_II_COMPILER_SUPPORTS_MPI
00027 # include <Epetra_MpiComm.h>
00028 # else
00029 # include <Epetra_SerialComm.h>
00030 # endif
00031 # include <Epetra_Map.h>
00032
00033 # include <Teuchos_ParameterList.hpp>
00034 # include <Epetra_Operator.h>
00035 # include <Epetra_Vector.h>
00036
00037
00038 class Ifpack_Preconditioner;
00039 class Ifpack_Chebyshev;
00040 namespace ML_Epetra
00041 {
00042 class MultiLevelPreconditioner;
00043 }
00044
00045
00046 DEAL_II_NAMESPACE_OPEN
00047
00048
00049 template <typename number> class SparseMatrix;
00050 template <typename number> class Vector;
00051 class SparsityPattern;
00052
00057 namespace TrilinosWrappers
00058 {
00059
00060 class SparseMatrix;
00061 class BlockSparseMatrix;
00062 class SolverBase;
00063
00072 class PreconditionBase : public Subscriptor
00073 {
00074 public:
00075
00081 struct AdditionalData
00082 {};
00083
00093 PreconditionBase ();
00094
00098 PreconditionBase (const PreconditionBase &);
00099
00103 ~PreconditionBase ();
00104
00110 void clear ();
00111
00115 void vmult (VectorBase &dst,
00116 const VectorBase &src) const;
00117
00125 void vmult (::Vector<double> &dst,
00126 const ::Vector<double> &src) const;
00127
00134 void vmult (::parallel::distributed::Vector<double> &dst,
00135 const ::parallel::distributed::Vector<double> &src) const;
00136
00140 DeclException1 (ExcNonMatchingMaps,
00141 std::string,
00142 << "The sparse matrix the preconditioner is based on "
00143 << "uses a map that is not compatible to the one in vector "
00144 << arg1
00145 << ". Check preconditioner and matrix setup.");
00146
00147 friend class SolverBase;
00148 friend class PreconditionStokes;
00149
00150 protected:
00157 std_cxx1x::shared_ptr<Epetra_Operator> preconditioner;
00158
00165 #ifdef DEAL_II_COMPILER_SUPPORTS_MPI
00166 Epetra_MpiComm communicator;
00167 #else
00168 Epetra_SerialComm communicator;
00169 #endif
00170
00176 std_cxx1x::shared_ptr<Epetra_Map> vector_distributor;
00177 };
00178
00179
00197 class PreconditionJacobi : public PreconditionBase
00198 {
00199 public:
00200
00226 struct AdditionalData
00227 {
00234 AdditionalData (const double omega = 1,
00235 const double min_diagonal = 0);
00236
00242 double omega;
00243
00256 double min_diagonal;
00257 };
00258
00266 void initialize (const SparseMatrix &matrix,
00267 const AdditionalData &additional_data = AdditionalData());
00268 };
00269
00270
00271
00272
00301 class PreconditionSSOR : public PreconditionBase
00302 {
00303 public:
00304
00336 struct AdditionalData
00337 {
00349 AdditionalData (const double omega = 1,
00350 const double min_diagonal = 0,
00351 const unsigned int overlap = 0);
00352
00358 double omega;
00359
00373 double min_diagonal;
00374
00382 unsigned int overlap;
00383 };
00384
00394 void initialize (const SparseMatrix &matrix,
00395 const AdditionalData &additional_data = AdditionalData());
00396 };
00397
00398
00399
00400
00429 class PreconditionSOR : public PreconditionBase
00430 {
00431 public:
00432
00464 struct AdditionalData
00465 {
00477 AdditionalData (const double omega = 1,
00478 const double min_diagonal = 0,
00479 const unsigned int overlap = 0);
00480
00486 double omega;
00487
00501 double min_diagonal;
00502
00510 unsigned int overlap;
00511 };
00512
00522 void initialize (const SparseMatrix &matrix,
00523 const AdditionalData &additional_data = AdditionalData());
00524 };
00525
00526
00527
00566 class PreconditionIC : public PreconditionBase
00567 {
00568 public:
00608 struct AdditionalData
00609 {
00629 AdditionalData (const unsigned int ic_fill = 0,
00630 const double ic_atol = 0.,
00631 const double ic_rtol = 1.,
00632 const unsigned int overlap = 0);
00633
00653 unsigned int ic_fill;
00654
00663 double ic_atol;
00664
00672 double ic_rtol;
00673
00681 unsigned int overlap;
00682 };
00683
00691 void initialize (const SparseMatrix &matrix,
00692 const AdditionalData &additional_data = AdditionalData());
00693 };
00694
00695
00696
00735 class PreconditionILU : public PreconditionBase
00736 {
00737 public:
00777 struct AdditionalData
00778 {
00797 AdditionalData (const unsigned int ilu_fill = 0,
00798 const double ilu_atol = 0.,
00799 const double ilu_rtol = 1.,
00800 const unsigned int overlap = 0);
00801
00821 unsigned int ilu_fill;
00822
00831 double ilu_atol;
00832
00840 double ilu_rtol;
00841
00849 unsigned int overlap;
00850 };
00851
00859 void initialize (const SparseMatrix &matrix,
00860 const AdditionalData &additional_data = AdditionalData());
00861 };
00862
00863
00864
00865
00866
00867
00905 class PreconditionILUT : public PreconditionBase
00906 {
00907 public:
00943 struct AdditionalData
00944 {
00962 AdditionalData (const double ilut_drop = 0.,
00963 const unsigned int ilut_fill = 0,
00964 const double ilut_atol = 0.,
00965 const double ilut_rtol = 1.,
00966 const unsigned int overlap = 0);
00967
00974 double ilut_drop;
00975
00995 unsigned int ilut_fill;
00996
01005 double ilut_atol;
01006
01014 double ilut_rtol;
01015
01023 unsigned int overlap;
01024 };
01025
01033 void initialize (const SparseMatrix &matrix,
01034 const AdditionalData &additional_data = AdditionalData());
01035 };
01036
01037
01038
01059 class PreconditionBlockwiseDirect : public PreconditionBase
01060 {
01061 public:
01067 struct AdditionalData
01068 {
01072 AdditionalData (const unsigned int overlap = 0);
01073
01074
01082 unsigned int overlap;
01083 };
01084
01092 void initialize (const SparseMatrix &matrix,
01093 const AdditionalData &additional_data = AdditionalData());
01094 };
01095
01096
01097
01098
01099
01100
01111 class PreconditionChebyshev : public PreconditionBase
01112 {
01113 public:
01119 struct AdditionalData
01120 {
01124 AdditionalData (const unsigned int degree = 1,
01125 const double max_eigenvalue = 10.,
01126 const double eigenvalue_ratio = 30.,
01127 const double min_eigenvalue = 1.,
01128 const double min_diagonal = 1e-12,
01129 const bool nonzero_starting = false);
01130
01139 unsigned int degree;
01140
01148 double max_eigenvalue;
01149
01155 double eigenvalue_ratio;
01156
01163 double min_eigenvalue;
01164
01171 double min_diagonal;
01172
01191 bool nonzero_starting;
01192 };
01193
01201 void initialize (const SparseMatrix &matrix,
01202 const AdditionalData &additional_data = AdditionalData());
01203 };
01204
01205
01206
01249 class PreconditionAMG : public PreconditionBase
01250 {
01251 public:
01252
01253 struct AdditionalData
01254 {
01261 AdditionalData (const bool elliptic = true,
01262 const bool higher_order_elements = false,
01263 const unsigned int n_cycles = 1,
01264 const bool w_cyle = false,
01265 const double aggregation_threshold = 1e-4,
01266 const std::vector<std::vector<bool> > &constant_modes = std::vector<std::vector<bool> > (1),
01267 const unsigned int smoother_sweeps = 2,
01268 const unsigned int smoother_overlap = 0,
01269 const bool output_details = false);
01270
01282 bool elliptic;
01283
01290 bool higher_order_elements;
01291
01297 unsigned int n_cycles;
01298
01304 bool w_cycle;
01305
01322 double aggregation_threshold;
01323
01333 std::vector<std::vector<bool> > constant_modes;
01334
01353 unsigned int smoother_sweeps;
01354
01360 unsigned int smoother_overlap;
01361
01370 bool output_details;
01371 };
01372
01381 void initialize (const SparseMatrix &matrix,
01382 const AdditionalData &additional_data = AdditionalData());
01383
01404 void initialize (const SparseMatrix &matrix,
01405 const Teuchos::ParameterList &ml_parameters);
01406
01417 template <typename number>
01418 void initialize (const ::SparseMatrix<number> &deal_ii_sparse_matrix,
01419 const AdditionalData &additional_data = AdditionalData(),
01420 const double drop_tolerance = 1e-13,
01421 const ::SparsityPattern *use_this_sparsity = 0);
01422
01446 void reinit ();
01447
01453 void clear ();
01454
01459 std::size_t memory_consumption () const;
01460
01461 private:
01466 std_cxx1x::shared_ptr<SparseMatrix> trilinos_matrix;
01467 };
01468
01469
01470
01471
01472
01473
01474 #ifndef DOXYGEN
01475
01476 inline
01477 void
01478 PreconditionBase::vmult (VectorBase &dst,
01479 const VectorBase &src) const
01480 {
01481 Assert (dst.vector_partitioner().SameAs(preconditioner->OperatorRangeMap()),
01482 ExcNonMatchingMaps("dst"));
01483 Assert (src.vector_partitioner().SameAs(preconditioner->OperatorDomainMap()),
01484 ExcNonMatchingMaps("src"));
01485
01486 const int ierr = preconditioner->ApplyInverse (src.trilinos_vector(),
01487 dst.trilinos_vector());
01488 AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01489 }
01490
01491
01492
01493
01494
01495
01496
01497
01498
01499
01500
01501
01502
01503
01504
01505
01506
01507
01508
01509
01510
01511
01512
01513
01514
01515
01516 inline
01517 void PreconditionBase::vmult (::Vector<double> &dst,
01518 const ::Vector<double> &src) const
01519 {
01520 AssertDimension (static_cast<int>(dst.size()),
01521 preconditioner->OperatorDomainMap().NumMyElements());
01522 AssertDimension (static_cast<int>(src.size()),
01523 preconditioner->OperatorRangeMap().NumMyElements());
01524 Epetra_Vector tril_dst (View, preconditioner->OperatorDomainMap(),
01525 dst.begin());
01526 Epetra_Vector tril_src (View, preconditioner->OperatorRangeMap(),
01527 const_cast<double*>(src.begin()));
01528
01529 const int ierr = preconditioner->ApplyInverse (tril_src, tril_dst);
01530 AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01531 }
01532
01533
01534
01535 inline
01536 void
01537 PreconditionBase::vmult (parallel::distributed::Vector<double> &dst,
01538 const parallel::distributed::Vector<double> &src) const
01539 {
01540 AssertDimension (static_cast<int>(dst.local_size()),
01541 preconditioner->OperatorDomainMap().NumMyElements());
01542 AssertDimension (static_cast<int>(src.local_size()),
01543 preconditioner->OperatorRangeMap().NumMyElements());
01544 Epetra_Vector tril_dst (View, preconditioner->OperatorDomainMap(),
01545 dst.begin());
01546 Epetra_Vector tril_src (View, preconditioner->OperatorRangeMap(),
01547 const_cast<double*>(src.begin()));
01548
01549 const int ierr = preconditioner->ApplyInverse (tril_src, tril_dst);
01550 AssertThrow (ierr == 0, ExcTrilinosError(ierr));
01551 }
01552
01553 #endif
01554
01555 }
01556
01557
01561 DEAL_II_NAMESPACE_CLOSE
01562
01563 #endif // DEAL_II_USE_TRILINOS
01564
01565
01566
01567 #endif
01568
01569