dune-istl  2.2.0
schwarz.hh
Go to the documentation of this file.
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=4 sw=2 sts=2:
00003 #ifndef DUNE_SCHWARZ_HH
00004 #define DUNE_SCHWARZ_HH
00005 
00006 #include <iostream>              // for input/output to shell
00007 #include <fstream>               // for input/output to files
00008 #include <vector>                // STL vector class
00009 #include <sstream>
00010 
00011 #include <cmath>                // Yes, we do some math here
00012 
00013 #include <dune/common/timer.hh>
00014 
00015 #include"io.hh"
00016 #include"bvector.hh"
00017 #include"vbvector.hh"
00018 #include"bcrsmatrix.hh"
00019 #include"io.hh"
00020 #include"gsetc.hh"
00021 #include"ilu.hh"
00022 #include"operators.hh"
00023 #include"solvers.hh"
00024 #include"preconditioners.hh"
00025 #include"scalarproducts.hh"
00026 #include"owneroverlapcopy.hh"
00027 
00028 namespace Dune {
00029 
00058   template<class M, class X, class Y, class C>
00059   class OverlappingSchwarzOperator : public AssembledLinearOperator<M,X,Y>
00060   {
00061   public:
00063         typedef M matrix_type;
00065         typedef X domain_type;
00067         typedef Y range_type;
00069         typedef typename X::field_type field_type;
00071     typedef C communication_type;
00072     
00073         enum {
00075           category=SolverCategory::overlapping
00076         };
00077 
00085         OverlappingSchwarzOperator (const matrix_type& A, const communication_type& com) 
00086           : _A_(A), communication(com) 
00087         {}
00088 
00090         virtual void apply (const X& x, Y& y) const
00091         {
00092           y = 0;
00093           _A_.umv(x,y); // result is consistent on interior+border
00094           communication.project(y); // we want this here to avoid it before the preconditioner
00095                                     // since there d is const!
00096         }
00097 
00099         virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
00100         {
00101           _A_.usmv(alpha,x,y); // result is consistent on interior+border
00102           communication.project(y); // we want this here to avoid it before the preconditioner
00103                                     // since there d is const!
00104         }
00105 
00107         virtual const matrix_type& getmat () const
00108         {
00109           return _A_;
00110         }
00111 
00112   private:
00113         const matrix_type& _A_;
00114         const communication_type& communication;
00115   };
00116 
00128   template<class X, class C>
00129   class OverlappingSchwarzScalarProduct : public ScalarProduct<X>
00130   {
00131   public:
00133         typedef X domain_type;
00135         typedef typename X::field_type field_type;
00137     typedef C communication_type;
00138     
00140         enum {category=SolverCategory::overlapping};
00141 
00146         OverlappingSchwarzScalarProduct (const communication_type& com)
00147           : communication(com)
00148         {}
00149 
00154         virtual field_type dot (const X& x, const X& y)
00155         {
00156           field_type result;
00157           communication.dot(x,y,result);
00158           return result;
00159         }
00160 
00164         virtual double norm (const X& x)
00165         {
00166           return communication.norm(x);
00167         }
00168 
00169   private:
00170         const communication_type& communication;
00171   };
00172 
00173   template<class X, class C>
00174   struct ScalarProductChooser<X,C,SolverCategory::overlapping>
00175   {
00177     typedef OverlappingSchwarzScalarProduct<X,C> ScalarProduct;
00179     typedef C communication_type;
00180 
00181     enum{
00183       solverCategory=SolverCategory::overlapping
00184         };
00185 
00186     static ScalarProduct* construct(const communication_type& comm)
00187     {
00188       return new ScalarProduct(comm);
00189     }
00190   };
00191 
00198 
00199   template<class M, class X, class Y, class C>
00200   class ParSSOR : public Preconditioner<X,Y> {
00201   public:
00203     typedef M matrix_type;
00205     typedef X domain_type;
00207     typedef Y range_type;
00209     typedef typename X::field_type field_type;
00211     typedef C communication_type;
00212     
00213     // define the category
00214     enum {
00216       category=SolverCategory::overlapping};
00217     
00227     ParSSOR (const matrix_type& A, int n, field_type w, const communication_type& c)
00228       : _A_(A), _n(n), _w(w), communication(c)
00229     {   }
00230 
00236     virtual void pre (X& x, Y& b) 
00237         {
00238           communication.copyOwnerToAll(x,x); // make dirichlet values consistent
00239         }
00240 
00246     virtual void apply (X& v, const Y& d)
00247     {
00248       for (int i=0; i<_n; i++){
00249                 bsorf(_A_,v,d,_w);
00250                 bsorb(_A_,v,d,_w);
00251       }
00252           communication.copyOwnerToAll(v,v);
00253     }
00254 
00260     virtual void post (X& x) {}
00261 
00262   private:
00264     const matrix_type& _A_;
00266     int _n;
00268     field_type _w;
00270         const communication_type& communication;
00271   };
00272 
00273   namespace Amg
00274   {
00275     template<class T> class ConstructionTraits;
00276   }
00277 
00286   template<class X, class Y, class C, class T=Preconditioner<X,Y> >
00287   class BlockPreconditioner : public Preconditioner<X,Y> {
00288     friend class Amg::ConstructionTraits<BlockPreconditioner<X,Y,C,T> >;
00289   public:
00291     typedef X domain_type;
00293     typedef Y range_type;
00295     typedef typename X::field_type field_type;
00297     typedef C communication_type;
00298 
00299     // define the category
00300     enum {
00302       category=SolverCategory::overlapping};
00303     
00311     BlockPreconditioner (T& p, const communication_type& c)
00312       : preconditioner(p), communication(c)
00313     {   }
00314 
00320     virtual void pre (X& x, Y& b) 
00321         {
00322           communication.copyOwnerToAll(x,x); // make dirichlet values consistent
00323           preconditioner.pre(x,b);
00324         }
00325 
00331     virtual void apply (X& v, const Y& d)
00332     {
00333           preconditioner.apply(v,d);
00334           communication.copyOwnerToAll(v,v);
00335     }
00336 
00337     template<bool forward>
00338     void apply (X& v, const Y& d)
00339     {
00340       preconditioner.template apply<forward>(v,d);
00341       communication.copyOwnerToAll(v,v);
00342     }
00343     
00349     virtual void post (X& x) 
00350         {
00351           preconditioner.post(x);
00352         }
00353 
00354   private:
00356         T& preconditioner;
00357 
00359         const communication_type& communication;
00360   };
00361 
00364 } // end namespace
00365 
00366 #endif