include/deal.II/lac/trilinos_precondition.h

00001 //---------------------------------------------------------------------------
00002 //    @f$Id: trilinos_precondition.h 25345 2012-03-31 08:37:04Z bangerth @f$
00003 //
00004 //    Copyright (C) 2008, 2010, 2012 by the deal.II authors
00005 //
00006 //    This file is subject to QPL and may not be  distributed
00007 //    without copyright and license information. Please refer
00008 //    to the file deal.II/doc/license.html for the  text  and
00009 //    further information on this license.
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 // forward declarations
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 // forward declarations
00049 template <typename number> class SparseMatrix;
00050 template <typename number> class Vector;
00051 class SparsityPattern;
00052 
00057 namespace TrilinosWrappers
00058 {
00059   // forward declarations
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 // -------------------------- inline and template functions ----------------------
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                                         // For the implementation of
01493                                         // the <code>vmult</code>
01494                                         // function with deal.II data
01495                                         // structures we note that
01496                                         // invoking a call of the
01497                                         // Trilinos preconditioner
01498                                         // requires us to use Epetra
01499                                         // vectors as well. We do this
01500                                         // by providing a view, i.e.,
01501                                         // feed Trilinos with a
01502                                         // pointer to the data, so we
01503                                         // avoid copying the content
01504                                         // of the vectors during the
01505                                         // iteration (this function is
01506                                         // only useful when used in
01507                                         // serial anyway). In the
01508                                         // declaration of the right
01509                                         // hand side, we need to cast
01510                                         // the source vector (that is
01511                                         // <code>const</code> in all
01512                                         // deal.II calls) to
01513                                         // non-constant value, as this
01514                                         // is the way Trilinos wants
01515                                         // to have them.
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 /*----------------------------   trilinos_precondition.h     ---------------------------*/
01566 
01567 #endif
01568 /*----------------------------   trilinos_precondition.h     ---------------------------*/
01569 
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines

deal.II documentation generated on Wed May 23 2012 12:03:30 by doxygen 1.7.3