dune-istl  2.2.0
novlpschwarz.hh
Go to the documentation of this file.
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set ts=4 sw=2 et sts=2:
00003 #ifndef DUNE_NOVLPSCHWARZ_HH
00004 #define DUNE_NOVLPSCHWARZ_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 NonoverlappingSchwarzOperator : 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     typedef typename C::PIS PIS;
00074     typedef typename C::RI RI;
00075     typedef typename RI::RemoteIndexList RIL;
00076     typedef typename RI::const_iterator RIIterator;
00077     typedef typename RIL::const_iterator RILIterator;
00078     typedef typename M::ConstColIterator ColIterator;
00079     typedef typename M::ConstRowIterator RowIterator;
00080     typedef std::multimap<int,int> MM;
00081     typedef std::multimap<int,std::pair<int,RILIterator> > RIMap;
00082     typedef typename RIMap::iterator RIMapit;
00083     
00084     enum {
00086       category=SolverCategory::nonoverlapping
00087     };
00088 
00096     NonoverlappingSchwarzOperator (const matrix_type& A, const communication_type& com) 
00097       : _A_(A), communication(com), buildcomm(true)
00098     {}
00099 
00101     virtual void apply (const X& x, Y& y) const
00102     {
00103       y = 0;
00104       novlp_op_apply(x,y,1);
00105       communication.addOwnerCopyToOwnerCopy(y,y); 
00106     }
00107 
00109     virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
00110     {
00111       // only apply communication to alpha*A*x to make it consistent,
00112       // y already has to be consistent. 
00113       Y y1(y);
00114       y = 0;
00115       novlp_op_apply(x,y,alpha);
00116       communication.addOwnerCopyToOwnerCopy(y,y);
00117       y += y1;
00118     }
00119     
00121     virtual const matrix_type& getmat () const
00122     {
00123       return _A_;
00124     }
00125 
00126     void novlp_op_apply (const X& x, Y& y, field_type alpha) const
00127     { 
00128       //get index sets
00129       const PIS& pis=communication.indexSet();
00130       const RI& ri = communication.remoteIndices();
00131       
00132       // at the beginning make a multimap "bordercontribution".
00133       // process has i and j as border dofs but is not the owner
00134       // => only contribute to Ax if i,j is in bordercontribution     
00135       if (buildcomm == true){
00136       
00137         // set up mask vector
00138         if (mask.size()!=static_cast<typename std::vector<double>::size_type>(x.size())){
00139           mask.resize(x.size());
00140           for (typename std::vector<double>::size_type i=0; i<mask.size(); i++) 
00141             mask[i] = 1; 
00142           for (typename PIS::const_iterator i=pis.begin(); i!=pis.end(); ++i)
00143             if (i->local().attribute()==OwnerOverlapCopyAttributeSet::copy)
00144               mask[i->local().local()] = 0;
00145             else if (i->local().attribute()==OwnerOverlapCopyAttributeSet::overlap)
00146               mask[i->local().local()] = 2;
00147         }
00148         
00149         for (MM::iterator iter = bordercontribution.begin();
00150              iter != bordercontribution.end(); ++iter)
00151           bordercontribution.erase(iter);
00152         std::map<int,int> owner; //key: local index i, value: process, that owns i
00153         RIMap rimap;
00154         
00155         // for each local index make multimap rimap: 
00156         // key: local index i, data: pair of process that knows i and pointer to RI entry
00157         for (RowIterator i = _A_.begin(); i != _A_.end(); ++i)
00158           if (mask[i.index()] == 0)
00159             for (RIIterator remote = ri.begin(); remote != ri.end(); ++remote){
00160               RIL& ril = *(remote->second.first);
00161               for (RILIterator rindex = ril.begin(); rindex != ril.end(); ++rindex)
00162                 if (rindex->attribute() != OwnerOverlapCopyAttributeSet::overlap)
00163                   if (rindex->localIndexPair().local().local() == i.index()){
00164                     rimap.insert
00165                       (std::make_pair(i.index(),
00166                                       std::pair<int,RILIterator>(remote->first, rindex)));
00167                     if(rindex->attribute()==OwnerOverlapCopyAttributeSet::owner)
00168                       owner.insert(std::make_pair(i.index(),remote->first));
00169                   }
00170             }
00171         
00172         int iowner = 0; 
00173         for (RowIterator i = _A_.begin(); i != _A_.end(); ++i){
00174           if (mask[i.index()] == 0){
00175             std::map<int,int>::iterator it = owner.find(i.index());
00176             iowner = it->second;
00177             std::pair<RIMapit, RIMapit> foundiit = rimap.equal_range(i.index());
00178             for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j){
00179               if (mask[j.index()] == 0){
00180                 bool flag = true;
00181                 for (RIMapit foundi = foundiit.first; foundi != foundiit.second; ++foundi){ 
00182                   std::pair<RIMapit, RIMapit> foundjit = rimap.equal_range(j.index());
00183                   for (RIMapit foundj = foundjit.first; foundj != foundjit.second; ++foundj)
00184                     if (foundj->second.first == foundi->second.first)
00185                       if (foundj->second.second->attribute() == OwnerOverlapCopyAttributeSet::owner
00186                           || foundj->second.first == iowner
00187                           || foundj->second.first  < communication.communicator().rank()){
00188                         flag = false;
00189                         continue;
00190                       }
00191                   if (flag == false)
00192                     continue;
00193                 }
00194                 // don“t contribute to Ax if 
00195                 // 1. the owner of j has i as interior/border dof
00196                 // 2. iowner has j as interior/border dof 
00197                 // 3. there is another process with smaller rank that has i and j 
00198                 // as interor/border dofs
00199                 // if the owner of j does not have i as interior/border dof, 
00200                 // it will not be taken into account
00201                 if (flag==true)
00202                   bordercontribution.insert(std::pair<int,int>(i.index(),j.index()));
00203               }
00204             }
00205           }
00206         }
00207         buildcomm = false;
00208       }
00209       
00210       //compute alpha*A*x nonoverlapping case
00211       for (RowIterator i = _A_.begin(); i != _A_.end(); ++i){ 
00212         if (mask[i.index()] == 0){ 
00213           //dof doesn't belong to process but is border (not ghost)
00214           for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j){
00215             if (mask[j.index()] == 1) //j is owner => then sum entries
00216               (*j).usmv(alpha,x[j.index()],y[i.index()]);
00217             else if (mask[j.index()] == 0){
00218               std::pair<MM::iterator, MM::iterator> itp = 
00219                 bordercontribution.equal_range(i.index());
00220               for (MM::iterator it = itp.first; it != itp.second; ++it) 
00221                 if ((*it).second == (int)j.index())
00222                   (*j).usmv(alpha,x[j.index()],y[i.index()]); 
00223             }
00224           }
00225         }
00226         else if (mask[i.index()] == 1){
00227           for (ColIterator j = _A_[i.index()].begin(); j != _A_[i.index()].end(); ++j)
00228             if (mask[j.index()] != 2)
00229               (*j).usmv(alpha,x[j.index()],y[i.index()]);
00230         }
00231       }
00232     }
00233     
00234   private:
00235     const matrix_type& _A_;
00236     const communication_type& communication;
00237     mutable bool buildcomm;
00238     mutable std::vector<double> mask;
00239     mutable std::multimap<int,int>  bordercontribution;
00240   };
00241   
00253   template<class X, class C>
00254   class NonoverlappingSchwarzScalarProduct : public ScalarProduct<X>
00255   {
00256   public:
00258     typedef X domain_type;
00260     typedef typename X::field_type field_type;
00262     typedef C communication_type;
00263     
00265     enum {category=SolverCategory::nonoverlapping};
00266 
00271     NonoverlappingSchwarzScalarProduct (const communication_type& com)
00272       : communication(com)
00273     {}
00274 
00279     virtual field_type dot (const X& x, const X& y)
00280     {
00281       field_type result;
00282       communication.dot(x,y,result);
00283       return result;
00284     }
00285 
00289     virtual double norm (const X& x)
00290     {
00291       return communication.norm(x);
00292     }
00293     
00296     void make_consistent (X& x) const
00297     {
00298       communication.copyOwnerToAll(x,x);
00299     }
00300     
00301   private:
00302     const communication_type& communication;
00303   };
00304 
00305   template<class X, class C>
00306   struct ScalarProductChooser<X,C,SolverCategory::nonoverlapping>
00307   {
00309     typedef NonoverlappingSchwarzScalarProduct<X,C> ScalarProduct;
00311     typedef C communication_type;
00312 
00313     enum{
00315       solverCategory=SolverCategory::nonoverlapping
00316     };
00317 
00318     static ScalarProduct* construct(const communication_type& comm)
00319     {
00320       return new ScalarProduct(comm);
00321     }
00322   };
00323 
00324   namespace Amg
00325   {
00326     template<class T> class ConstructionTraits;
00327   }
00328 
00338   template<class C, class P>
00339   class NonoverlappingBlockPreconditioner
00340     : public Dune::Preconditioner<typename P::domain_type,typename P::range_type> {
00341     friend class Amg::ConstructionTraits<NonoverlappingBlockPreconditioner<C,P> >;
00342   public:
00344     typedef typename P::domain_type domain_type;
00346     typedef typename P::range_type range_type;
00348     typedef C communication_type;
00349 
00350     // define the category
00351     enum {
00353       category=SolverCategory::nonoverlapping};
00354     
00362     NonoverlappingBlockPreconditioner (P& prec, const communication_type& c)
00363       : preconditioner(prec), communication(c)
00364     {}
00365     
00371     virtual void pre (domain_type& x, range_type& b) 
00372     {
00373       preconditioner.pre(x,b);
00374     }
00375 
00381     virtual void apply (domain_type& v, const range_type& d)
00382     {
00383       // block preconditioner equivalent to WrappedPreconditioner from 
00384       // pdelab/backend/ovlpistsolverbackend.hh, 
00385       // but not to BlockPreconditioner from schwarz.hh
00386       preconditioner.apply(v,d);
00387       communication.addOwnerCopyToOwnerCopy(v,v);
00388     }
00389 
00395     virtual void post (domain_type& x) 
00396     {
00397       preconditioner.post(x);
00398     }
00399 
00400   private:
00402     P& preconditioner;
00403 
00405     const communication_type& communication;
00406   };
00407 
00410 } // end namespace
00411 
00412 #endif