dune-istl  2.2.0
overlappingschwarz.hh
Go to the documentation of this file.
00001 // -*- tab-width: 8; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=8 sw=2 sts=2:
00003 #ifndef DUNE_OVERLAPPINGSCHWARZ_HH
00004 #define DUNE_OVERLAPPINGSCHWARZ_HH
00005 #include<cassert>
00006 #include<algorithm>
00007 #include<functional>
00008 #include<vector>
00009 #include<set>
00010 #include<dune/common/dynmatrix.hh>
00011 #include<dune/common/sllist.hh>
00012 #include"preconditioners.hh"
00013 #include"superlu.hh"
00014 #include"bvector.hh"
00015 #include"bcrsmatrix.hh"
00016 #include"ilusubdomainsolver.hh"
00017 
00018 namespace Dune
00019 {
00020 
00032   template<class M, class X, class TM, class TD, class TA>
00033   class SeqOverlappingSchwarz;
00034   
00038   template<class I, class S, class D>
00039   class OverlappingSchwarzInitializer
00040   {
00041   public:
00043     typedef D subdomain_vector;
00044 
00045     typedef I InitializerList;
00046     typedef typename InitializerList::value_type AtomInitializer;
00047     typedef typename AtomInitializer::Matrix Matrix;
00048     typedef typename Matrix::const_iterator Iter;
00049     typedef typename Matrix::row_type::const_iterator CIter;
00050     
00051     typedef S IndexSet;
00052     typedef typename IndexSet::size_type size_type;
00053     
00054     OverlappingSchwarzInitializer(InitializerList& il,
00055                                   const IndexSet& indices,
00056                                   const subdomain_vector& domains);
00057     
00058     
00059     void addRowNnz(const Iter& row);
00060     
00061     void allocate();
00062         
00063     void countEntries(const Iter& row, const CIter& col) const;
00064 
00065     void calcColstart() const;
00066     
00067     void copyValue(const Iter& row, const CIter& col) const;
00068     
00069     void createMatrix() const;
00070   private:
00071     class IndexMap
00072     {
00073     public:
00074       typedef typename S::size_type size_type;
00075       typedef std::map<size_type,size_type> Map;      
00076       typedef typename Map::iterator iterator;
00077       typedef typename Map::const_iterator const_iterator;
00078       
00079       IndexMap();
00080       
00081       void insert(size_type grow);
00082       
00083       const_iterator find(size_type grow) const;
00084       
00085       iterator find(size_type grow);
00086 
00087       iterator begin();
00088       
00089       const_iterator begin()const;
00090 
00091       iterator end();
00092       
00093       const_iterator end() const;
00094             
00095     private:
00096       std::map<size_type,size_type> map_;
00097       size_type row;
00098     };
00099     
00100     
00101     typedef typename InitializerList::iterator InitIterator;
00102     typedef typename IndexSet::const_iterator IndexIteratur;
00103     InitializerList* initializers;
00104     const IndexSet *indices;
00105     mutable std::vector<IndexMap> indexMaps;
00106     const subdomain_vector& domains;
00107   };
00108 
00112   struct AdditiveSchwarzMode
00113   {};
00114 
00118   struct MultiplicativeSchwarzMode
00119   {};
00120   
00125   struct SymmetricMultiplicativeSchwarzMode
00126   {};  
00127 
00132   template<class M, class X, class Y>
00133   class DynamicMatrixSubdomainSolver;
00134 
00135   // Specialization for BCRSMatrix
00136   template<class K, int n, class Al, class X, class Y>
00137   class DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >
00138   {
00139     typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> M;
00140   public:
00142     typedef typename Dune::remove_const<M>::type matrix_type;
00143     typedef K field_type;
00144     typedef typename Dune::remove_const<M>::type rilu_type;
00146     typedef X domain_type;
00148     typedef Y range_type;
00149     
00154     void apply (DynamicVector<field_type>& v, DynamicVector<field_type>& d)
00155     {
00156       assert(v.size() > 0);
00157       assert(v.size() == d.size());
00158       assert(A.rows() <= v.size());
00159       assert(A.cols() <= v.size());
00160       size_t sz = A.rows();
00161       v.resize(sz);
00162       d.resize(sz);
00163       A.solve(v,d);
00164       v.resize(v.capacity());
00165       d.resize(d.capacity());
00166     }
00167 
00175     template<class S>
00176     void setSubMatrix(const M& BCRS, S& rowset)
00177     {
00178       size_t sz = rowset.size();
00179       A.resize(sz*n,sz*n);
00180       typename DynamicMatrix<K>::RowIterator rIt = A.begin();
00181       typedef typename S::const_iterator SIter;
00182       size_t r = 0, c = 0;
00183       for(SIter rowIdx = rowset.begin(), rowEnd=rowset.end();
00184           rowIdx!= rowEnd; ++rowIdx, r++)
00185       {
00186         c = 0;
00187         for(SIter colIdx = rowset.begin(), colEnd=rowset.end(); 
00188             colIdx!= colEnd; ++colIdx, c++)
00189         {
00190           if (BCRS[*rowIdx].find(*colIdx) == BCRS[*rowIdx].end())
00191             continue;
00192           for (size_t i=0; i<n; i++)
00193           {
00194             for (size_t j=0; j<n; j++)
00195             {
00196               A[r*n+i][c*n+j] = BCRS[*rowIdx][*colIdx][i][j];
00197             }
00198           }
00199         }
00200       }
00201     }
00202   private:
00203     DynamicMatrix<K> A;
00204   };
00205 
00206   template<typename T>
00207   struct OverlappingAssigner
00208   {
00209   };
00210     
00211   // specialization for DynamicMatrix
00212   template<class K, int n, class Al, class X, class Y>
00213   class OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
00214   {
00215   public:
00216     typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> matrix_type;
00217     typedef K field_type;
00218     typedef Y range_type;
00219     typedef typename range_type::block_type block_type;
00220     typedef typename matrix_type::size_type size_type;
00221     
00229     OverlappingAssigner(std::size_t maxlength, const BCRSMatrix<FieldMatrix<K,n,n>, Al>& mat_, const X& b_, Y& x_);
00230 
00234     inline
00235     void deallocate();
00236     
00240     inline
00241     void resetIndexForNextDomain();
00242     
00247     inline
00248     DynamicVector<K> & lhs();
00249     
00254     inline
00255     DynamicVector<K> & rhs();
00256     
00261     inline
00262     void relaxResult(field_type relax);
00263     
00268     void operator()(const size_type& domainIndex);
00269     
00277     inline
00278     void assignResult(block_type& res);
00279     
00280   private:
00284     const matrix_type* mat;
00286     // we need a pointer, because we have to avoid deep copies
00287     DynamicVector<field_type> * rhs_;
00289     // we need a pointer, because we have to avoid deep copies
00290     DynamicVector<field_type> * lhs_;
00292     const range_type* b;
00294     range_type* x;
00296     std::size_t i;
00298     std::size_t maxlength_;
00299   };
00300 
00301 #if HAVE_SUPERLU
00302   template<typename T, typename A, int n, int m>
00303   struct OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
00304     {
00305       typedef BCRSMatrix<FieldMatrix<T,n,m>,A> matrix_type;
00306       typedef typename SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> >::range_type range_type;
00307       typedef typename range_type::field_type field_type;
00308       typedef typename range_type::block_type block_type;
00309       
00310       typedef typename matrix_type::size_type size_type;
00311       
00319       OverlappingAssigner(std::size_t maxlength, const BCRSMatrix<FieldMatrix<T,n,m>,A>& mat, 
00320                           const range_type& b, range_type& x);
00326       void deallocate();
00327       
00328       /*
00329        * @brief Resets the local index to zero.
00330        */
00331       void resetIndexForNextDomain();
00332       
00337       field_type* lhs();
00338             
00343       field_type* rhs();
00344       
00349       void relaxResult(field_type relax);
00350 
00355       void operator()(const size_type& domain);
00356 
00364       void assignResult(block_type& res);
00365       
00366     private:
00370       const matrix_type* mat;
00372       field_type* rhs_;
00374       field_type* lhs_;
00376       const range_type* b;
00378       range_type* x;
00380       std::size_t i;
00382       std::size_t maxlength_;
00383     };
00384 
00385 #endif
00386 
00387   template<class M, class X, class Y>
00388   class OverlappingAssignerILUBase
00389   {
00390   public:
00391     typedef M matrix_type;
00392     
00393     typedef typename M::field_type field_type;
00394   
00395     typedef typename Y::block_type block_type;
00396       
00397     typedef typename matrix_type::size_type size_type;
00405     OverlappingAssignerILUBase(std::size_t maxlength, const M& mat, 
00406                         const Y& b, X& x);
00412     void deallocate();
00413     
00417     void resetIndexForNextDomain();
00418     
00423     X& lhs();
00424     
00429     Y& rhs();
00430     
00435     void relaxResult(field_type relax);
00436     
00441     void operator()(const size_type& domain);
00442     
00450     void assignResult(block_type& res);
00451 
00452   private:
00456     const M* mat;
00458     X* lhs_;
00460     Y* rhs_;
00462     const Y* b;
00464     X* x;
00466     size_type i;
00467   };
00468   
00469   // specialization for ILU0
00470   template<class M, class X, class Y>
00471   class OverlappingAssigner<ILU0SubdomainSolver<M,X,Y> >
00472     : public OverlappingAssignerILUBase<M,X,Y>
00473   {
00474   public:
00482     OverlappingAssigner(std::size_t maxlength, const M& mat, 
00483                         const Y& b, X& x)
00484       : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
00485     {}
00486   };
00487 
00488   // specialization for ILUN
00489   template<class M, class X, class Y>
00490   class OverlappingAssigner<ILUNSubdomainSolver<M,X,Y> >
00491     : public OverlappingAssignerILUBase<M,X,Y>
00492   {
00493   public:
00501     OverlappingAssigner(std::size_t maxlength, const M& mat, 
00502                         const Y& b, X& x)
00503       : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
00504     {}
00505   };
00506 
00507     template<typename S, typename T>
00508     struct AdditiveAdder
00509     {
00510     };
00511     
00512     template<typename S, typename T, typename A, int n>
00513     struct AdditiveAdder<S, BlockVector<FieldVector<T,n>,A> >
00514     {
00515       typedef typename A::size_type size_type;
00516       AdditiveAdder(BlockVector<FieldVector<T,n>,A>& v, BlockVector<FieldVector<T,n>,A>& x,
00517                     OverlappingAssigner<S>& assigner, const T& relax_);
00518       void operator()(const size_type& domain);
00519       void axpy();
00520       
00521     private:
00522       BlockVector<FieldVector<T,n>,A>* v;
00523       BlockVector<FieldVector<T,n>,A>* x;
00524       OverlappingAssigner<S>* assigner;
00525       T relax;
00526     };
00527 
00528     template<typename S,typename T>
00529     struct MultiplicativeAdder
00530     {
00531     };
00532     
00533     template<typename S, typename T, typename A, int n>
00534     struct MultiplicativeAdder<S, BlockVector<FieldVector<T,n>,A> >
00535     {
00536       typedef typename A::size_type size_type;
00537       MultiplicativeAdder(BlockVector<FieldVector<T,n>,A>& v, BlockVector<FieldVector<T,n>,A>& x, 
00538                           OverlappingAssigner<S>& assigner_, const T& relax_);
00539       void operator()(const size_type& domain);
00540       void axpy();
00541       
00542     private:
00543       BlockVector<FieldVector<T,n>,A>* x;
00544       OverlappingAssigner<S>* assigner;
00545       T relax;
00546     };
00547 
00557     template<typename T, class X, class S>
00558     struct AdderSelector
00559     {};
00560     
00561   template<class X, class S>
00562     struct AdderSelector<AdditiveSchwarzMode,X,S>
00563     {
00564       typedef AdditiveAdder<S,X> Adder;
00565     };
00566     
00567   template<class X, class S>
00568     struct AdderSelector<MultiplicativeSchwarzMode,X,S>
00569     {
00570       typedef MultiplicativeAdder<S,X> Adder;
00571     };    
00572 
00573   template<class X, class S>
00574     struct AdderSelector<SymmetricMultiplicativeSchwarzMode,X,S>
00575     {
00576       typedef MultiplicativeAdder<S,X> Adder;
00577     };    
00578 
00590     template<typename T1, typename T2, bool forward>
00591     struct IteratorDirectionSelector
00592     {
00593       typedef T1 solver_vector;
00594       typedef typename solver_vector::iterator solver_iterator;
00595       typedef T2 subdomain_vector;
00596       typedef typename subdomain_vector::const_iterator domain_iterator;
00597     
00598       static solver_iterator begin(solver_vector& sv)
00599       {
00600         return sv.begin();
00601       }
00602       
00603       static solver_iterator end(solver_vector& sv)
00604       {
00605         return sv.end();
00606       }
00607       static domain_iterator begin(const subdomain_vector& sv)
00608       {
00609         return sv.begin();
00610       }
00611       
00612       static domain_iterator end(const subdomain_vector& sv)
00613       {
00614         return sv.end();
00615       }
00616     };
00617 
00618     template<typename T1, typename T2>
00619     struct IteratorDirectionSelector<T1,T2,false>
00620     {
00621       typedef T1 solver_vector;
00622       typedef typename solver_vector::reverse_iterator solver_iterator;
00623       typedef T2 subdomain_vector;
00624       typedef typename subdomain_vector::const_reverse_iterator domain_iterator;
00625     
00626       static solver_iterator begin(solver_vector& sv)
00627       {
00628         return sv.rbegin();
00629       }
00630       
00631       static solver_iterator end(solver_vector& sv)
00632       {
00633         return sv.rend();
00634       }
00635       static domain_iterator begin(const subdomain_vector& sv)
00636       {
00637         return sv.rbegin();
00638       }
00639       
00640       static domain_iterator end(const subdomain_vector& sv)
00641       {
00642         return sv.rend();
00643       }
00644     };
00645     
00654     template<class T>
00655     struct SeqOverlappingSchwarzApplier
00656     {
00657       typedef T smoother;
00658       typedef typename smoother::range_type range_type;
00659       
00660       static void apply(smoother& sm, range_type& v, const range_type& b)
00661       {
00662         sm.template apply<true>(v, b);
00663       }
00664     };
00665     
00666     template<class M, class X, class TD, class TA>
00667     struct SeqOverlappingSchwarzApplier<SeqOverlappingSchwarz<M,X,SymmetricMultiplicativeSchwarzMode,TD,TA> >
00668     {
00669       typedef SeqOverlappingSchwarz<M,X,SymmetricMultiplicativeSchwarzMode,TD,TA> smoother;
00670       typedef typename smoother::range_type range_type;
00671       
00672       static void apply(smoother& sm, range_type& v, const range_type& b)
00673       {
00674         sm.template apply<true>(v, b);
00675         sm.template apply<false>(v, b);
00676       }
00677     };
00678 
00679   template<class T>
00680   struct SeqOverlappingSchwarzAssembler
00681   {};
00682   
00683 
00684   template<class K, int n, class Al, class X, class Y>
00685   struct SeqOverlappingSchwarzAssembler< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
00686   {
00687     typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> matrix_type;
00688     template<class RowToDomain, class Solvers, class SubDomains>
00689     static std::size_t assembleLocalProblems(const RowToDomain& rowToDomain, const matrix_type& mat,
00690                                              Solvers& solvers, const SubDomains& domains,
00691                                              bool onTheFly);
00692   };
00693 
00694 #if HAVE_SUPERLU
00695   template<class T>
00696   struct SeqOverlappingSchwarzAssembler<SuperLU<T> >
00697   {
00698     typedef T matrix_type;
00699     template<class RowToDomain, class Solvers, class SubDomains>
00700     static std::size_t assembleLocalProblems(const RowToDomain& rowToDomain, const matrix_type& mat,
00701                                              Solvers& solvers, const SubDomains& domains,
00702                                              bool onTheFly);
00703   };
00704 #endif
00705 
00706   template<class M,class X, class Y>
00707   struct SeqOverlappingSchwarzAssemblerILUBase
00708   {
00709     typedef M matrix_type;
00710     template<class RowToDomain, class Solvers, class SubDomains>
00711     static std::size_t assembleLocalProblems(const RowToDomain& rowToDomain, const matrix_type& mat,
00712                                              Solvers& solvers, const SubDomains& domains,
00713                                              bool onTheFly);
00714   };
00715 
00716   template<class M,class X, class Y>
00717   struct SeqOverlappingSchwarzAssembler<ILU0SubdomainSolver<M,X,Y> >
00718     : public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
00719   {};
00720 
00721   template<class M,class X, class Y>
00722   struct SeqOverlappingSchwarzAssembler<ILUNSubdomainSolver<M,X,Y> >
00723     : public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
00724   {};
00725   
00736   template<class M, class X, class TM=AdditiveSchwarzMode, 
00737            class TD=ILU0SubdomainSolver<M,X,X>, class TA=std::allocator<X> >
00738   class SeqOverlappingSchwarz
00739     : public Preconditioner<X,X>
00740   {
00741   public:
00745     typedef M matrix_type;
00746 
00750     typedef X domain_type;
00751 
00755     typedef X range_type;
00756 
00763     typedef TM Mode;
00764 
00768     typedef typename X::field_type field_type;
00769 
00771     typedef typename matrix_type::size_type size_type;
00772 
00774     typedef TA allocator;
00775         
00777     typedef std::set<size_type, std::less<size_type>, 
00778                      typename TA::template rebind<size_type>::other> 
00779       subdomain_type;
00780 
00782     typedef std::vector<subdomain_type, typename TA::template rebind<subdomain_type>::other> subdomain_vector;
00783 
00785     typedef SLList<size_type, typename TA::template rebind<size_type>::other> subdomain_list;
00786       
00788     typedef std::vector<subdomain_list, typename TA::template rebind<subdomain_list>::other > rowtodomain_vector;
00789 
00791     typedef TD slu;
00792     
00794     typedef std::vector<slu, typename TA::template rebind<slu>::other> slu_vector;
00795 
00796     enum{
00798       category = SolverCategory::sequential
00799         };
00800     
00814     SeqOverlappingSchwarz(const matrix_type& mat, const subdomain_vector& subDomains,
00815                           field_type relaxationFactor=1, bool onTheFly_=true);
00816 
00828     SeqOverlappingSchwarz(const matrix_type& mat, const rowtodomain_vector& rowToDomain,
00829                           field_type relaxationFactor=1, bool onTheFly_=true);
00830                           
00836     virtual void pre (X& x, X& b) {}
00837 
00843     virtual void apply (X& v, const X& d);
00844 
00845     template<bool forward>
00846     void apply(X& v, const X& d);
00847     
00853     virtual void post (X& x) {
00854       Dune::dverb<<" avg nnz over subdomain is "<<nnz<<std::endl;
00855     }
00856   
00857   private:
00858     const M& mat;
00859     slu_vector solvers;
00860     subdomain_vector subDomains;
00861     field_type relax;
00862     
00863     typename M::size_type maxlength;
00864     std::size_t nnz;
00865 
00866     bool onTheFly;
00867   };
00868   
00869 
00870     
00871   template<class I, class S, class D>
00872   OverlappingSchwarzInitializer<I,S,D>::OverlappingSchwarzInitializer(InitializerList& il,
00873                                                                     const IndexSet& idx,
00874                                                                     const subdomain_vector& domains_)
00875     : initializers(&il), indices(&idx), indexMaps(il.size()), domains(domains_)
00876   {
00877   }
00878   
00879   
00880   template<class I, class S, class D>
00881   void OverlappingSchwarzInitializer<I,S,D>::addRowNnz(const Iter& row)
00882   {
00883     typedef typename IndexSet::value_type::const_iterator iterator;
00884     for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain){
00885       (*initializers)[*domain].addRowNnz(row, domains[*domain]);
00886       indexMaps[*domain].insert(row.index());
00887     }
00888   }
00889   
00890   template<class I, class S, class D>
00891   void OverlappingSchwarzInitializer<I,S,D>::allocate()
00892   {
00893     std::for_each(initializers->begin(), initializers->end(),
00894                   std::mem_fun_ref(&AtomInitializer::allocateMatrixStorage));
00895     std::for_each(initializers->begin(), initializers->end(), 
00896                   std::mem_fun_ref(&AtomInitializer::allocateMarker));
00897   }
00898   
00899   template<class I, class S, class D>
00900   void OverlappingSchwarzInitializer<I,S,D>::countEntries(const Iter& row, const CIter& col) const
00901   {
00902     typedef typename IndexSet::value_type::const_iterator iterator;
00903     for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain){
00904       typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
00905       if(v!= indexMaps[*domain].end()){
00906         (*initializers)[*domain].countEntries(indexMaps[*domain].find(col.index())->second);
00907       }
00908     }
00909   }
00910 
00911   template<class I, class S, class D>
00912   void OverlappingSchwarzInitializer<I,S,D>::calcColstart() const
00913   {
00914     std::for_each(initializers->begin(), initializers->end(),
00915                   std::mem_fun_ref(&AtomInitializer::calcColstart));
00916   }
00917 
00918   template<class I, class S, class D>
00919   void OverlappingSchwarzInitializer<I,S,D>::copyValue(const Iter& row, const CIter& col) const
00920   {
00921     typedef typename IndexSet::value_type::const_iterator iterator;
00922     for(iterator domain=(*indices)[row.index()].begin(); domain!= (*indices)[row.index()].end(); ++domain){
00923       typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
00924       if(v!= indexMaps[*domain].end()){
00925         assert(indexMaps[*domain].end()!=indexMaps[*domain].find(row.index()));
00926         (*initializers)[*domain].copyValue(col, indexMaps[*domain].find(row.index())->second, 
00927                                            v->second);
00928       }
00929     }
00930   }
00931     
00932   template<class I, class S, class D>
00933   void OverlappingSchwarzInitializer<I,S,D>::createMatrix() const
00934   {
00935     indexMaps.clear();
00936     indexMaps.swap(std::vector<IndexMap>(indexMaps));
00937     std::for_each(initializers->begin(), initializers->end(),
00938                   std::mem_fun_ref(&AtomInitializer::createMatrix));
00939   }
00940 
00941   template<class I, class S, class D>
00942   OverlappingSchwarzInitializer<I,S,D>::IndexMap::IndexMap()
00943     : row(0)
00944   {}
00945 
00946   template<class I, class S, class D>
00947   void OverlappingSchwarzInitializer<I,S,D>::IndexMap::insert(size_type grow)
00948   {
00949     assert(map_.find(grow)==map_.end());
00950     map_.insert(std::make_pair(grow, row++));
00951   }
00952 
00953   template<class I, class S, class D>
00954   typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator 
00955   OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow) const
00956   {
00957     return map_.find(grow);
00958   }
00959   
00960   template<class I, class S, class D>
00961   typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator 
00962   OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
00963   {
00964     return map_.find(grow);
00965   }
00966 
00967   template<class I, class S, class D>
00968   typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator 
00969   OverlappingSchwarzInitializer<I,S,D>::IndexMap::end() const
00970   {
00971     return map_.end();
00972   }
00973   
00974   template<class I, class S, class D>
00975   typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator 
00976   OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
00977   {
00978     return map_.end();
00979   }
00980 
00981   template<class I, class S, class D>
00982   typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator 
00983   OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin() const
00984   {
00985     return map_.begin();
00986   }
00987   
00988   template<class I, class S, class D>
00989   typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator 
00990   OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
00991   {
00992     return map_.begin();
00993   }
00994 
00995   template<class M, class X, class TM, class TD, class TA>
00996   SeqOverlappingSchwarz<M,X,TM,TD,TA>::SeqOverlappingSchwarz(const matrix_type& mat_, const rowtodomain_vector& rowToDomain,
00997                                                           field_type relaxationFactor, bool fly)
00998     : mat(mat_), relax(relaxationFactor), onTheFly(fly)
00999   {
01000     typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
01001     typedef typename subdomain_list::const_iterator DomainIterator;
01002 #ifdef DUNE_ISTL_WITH_CHECKING
01003     assert(rowToDomain.size()==mat.N());
01004     assert(rowToDomain.size()==mat.M());
01005     
01006     for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
01007       assert(iter->size()>0);
01008     
01009 #endif
01010     // calculate the number of domains
01011     size_type domains=0;
01012     for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
01013       for(DomainIterator d=iter->begin(); d != iter->end(); ++d)
01014         domains=std::max(domains, *d);
01015     ++domains;
01016     
01017     solvers.resize(domains);
01018     subDomains.resize(domains);
01019     
01020     // initialize subdomains to row mapping from row to subdomain mapping
01021     size_type row=0;
01022     for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter, ++row)
01023       for(DomainIterator d=iter->begin(); d != iter->end(); ++d)
01024         subDomains[*d].insert(row);
01025 
01026 #ifdef DUNE_ISTL_WITH_CHECKING
01027     size_type i=0;
01028     typedef typename subdomain_vector::const_iterator iterator;
01029     for(iterator iter=subDomains.begin(); iter != subDomains.end(); ++iter){
01030       typedef typename subdomain_type::const_iterator entry_iterator;
01031       Dune::dvverb<<"domain "<<i++<<":";
01032       for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry){
01033         Dune::dvverb<<" "<<*entry;
01034       }
01035       Dune::dvverb<<std::endl;
01036     }
01037 #endif
01038     maxlength = SeqOverlappingSchwarzAssembler<slu>
01039       ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
01040   }
01041   
01042   template<class M, class X, class TM, class TD, class TA>
01043   SeqOverlappingSchwarz<M,X,TM,TD,TA>::SeqOverlappingSchwarz(const matrix_type& mat_,
01044                                                           const subdomain_vector& sd,
01045                                                           field_type relaxationFactor,
01046                                                           bool fly)
01047     :  mat(mat_), solvers(sd.size()), subDomains(sd), relax(relaxationFactor),
01048        onTheFly(fly)
01049   {
01050     typedef typename subdomain_vector::const_iterator DomainIterator;
01051 
01052 #ifdef DUNE_ISTL_WITH_CHECKING
01053     size_type i=0;
01054     
01055     for(DomainIterator d=sd.begin(); d != sd.end(); ++d,++i){
01056       //std::cout<<i<<": "<<d->size()<<std::endl;
01057       assert(d->size()>0);
01058       typedef typename DomainIterator::value_type::const_iterator entry_iterator;
01059       Dune::dvverb<<"domain "<<i<<":";
01060       for(entry_iterator entry = d->begin(); entry != d->end(); ++entry){
01061         Dune::dvverb<<" "<<*entry;
01062       }
01063       Dune::dvverb<<std::endl;
01064     }
01065     
01066 #endif
01067     
01068     // Create a row to subdomain mapping
01069     rowtodomain_vector rowToDomain(mat.N());
01070 
01071     size_type domainId=0;
01072     
01073     for(DomainIterator domain=sd.begin(); domain != sd.end(); ++domain, ++domainId){
01074       typedef typename subdomain_type::const_iterator iterator;
01075       for(iterator row=domain->begin(); row != domain->end(); ++row)
01076         rowToDomain[*row].push_back(domainId);
01077     }
01078 
01079     maxlength = SeqOverlappingSchwarzAssembler<slu>
01080       ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
01081   }
01082 
01089   template<class M>
01090   struct SeqOverlappingSchwarzDomainSize {};
01091   
01092   template<typename T, typename A, int n, int m>
01093   struct SeqOverlappingSchwarzDomainSize<BCRSMatrix<FieldMatrix<T,n,m>,A > >
01094   {
01095     template<class Domain>
01096     static int size(const Domain & d)
01097     {
01098       assert(n==m);
01099       return m*d.size();
01100     }
01101   };
01102   
01103   template<class K, int n, class Al, class X, class Y>
01104   template<class RowToDomain, class Solvers, class SubDomains>
01105   std::size_t
01106   SeqOverlappingSchwarzAssembler< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >::
01107   assembleLocalProblems(const RowToDomain& rowToDomain, 
01108     const matrix_type& mat,
01109     Solvers& solvers,
01110     const SubDomains& subDomains,
01111     bool onTheFly)
01112   {
01113     typedef typename SubDomains::const_iterator DomainIterator;
01114     std::size_t maxlength = 0;
01115     
01116     assert(onTheFly);
01117     
01118     for(DomainIterator domain=subDomains.begin();domain!=subDomains.end();++domain)
01119       maxlength=std::max(maxlength, domain->size());
01120     maxlength*=n;
01121     
01122     return maxlength;
01123   }
01124 
01125 #if HAVE_SUPERLU
01126   template<class T>
01127   template<class RowToDomain, class Solvers, class SubDomains>
01128   std::size_t SeqOverlappingSchwarzAssembler<SuperLU<T> >::assembleLocalProblems(const RowToDomain& rowToDomain, 
01129                                                                                 const matrix_type& mat,
01130                                                                                 Solvers& solvers,
01131                                                                                 const SubDomains& subDomains,
01132                                                                                 bool onTheFly)
01133   {
01134     typedef typename std::vector<SuperMatrixInitializer<matrix_type> >::iterator InitializerIterator;
01135     typedef typename SubDomains::const_iterator DomainIterator;
01136     typedef typename Solvers::iterator SolverIterator;
01137     std::size_t maxlength = 0;
01138     
01139     if(onTheFly){
01140       for(DomainIterator domain=subDomains.begin();domain!=subDomains.end();++domain)
01141         maxlength=std::max(maxlength, domain->size());
01142       maxlength*=mat[0].begin()->N();
01143     }else{
01144       // initialize the initializers
01145       DomainIterator domain=subDomains.begin();
01146     
01147       // Create the initializers list.
01148       std::vector<SuperMatrixInitializer<matrix_type> > initializers(subDomains.size());    
01149       
01150       SolverIterator solver=solvers.begin();
01151       for(InitializerIterator initializer=initializers.begin(); initializer!=initializers.end();
01152           ++initializer, ++solver, ++domain){
01153         solver->mat.N_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
01154         solver->mat.M_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
01155         //solver->setVerbosity(true);
01156         *initializer=SuperMatrixInitializer<matrix_type>(solver->mat);
01157       }
01158 
01159       // Set up the supermatrices according to the subdomains
01160       typedef OverlappingSchwarzInitializer<std::vector<SuperMatrixInitializer<matrix_type> >,
01161         RowToDomain, SubDomains> Initializer;
01162 
01163       Initializer initializer(initializers, rowToDomain, subDomains);
01164       copyToSuperMatrix(initializer, mat);
01165       if(solvers.size()==1)
01166         assert(solvers[0].mat==mat);
01167     
01168       /*    for(SolverIterator solver=solvers.begin(); solver!=solvers.end(); ++solver)
01169             dPrint_CompCol_Matrix("superlu", &static_cast<SuperMatrix&>(solver->mat)); */
01170       
01171       // Calculate the LU decompositions
01172       std::for_each(solvers.begin(), solvers.end(), std::mem_fun_ref(&SuperLU<matrix_type>::decompose));
01173       for(SolverIterator solver=solvers.begin(); solver!=solvers.end(); ++solver){
01174         assert(solver->mat.N()==solver->mat.M());
01175         maxlength=std::max(maxlength, solver->mat.N());
01176         //writeCompColMatrixToMatlab(static_cast<SuperLUMatrix<M>&>(solver->mat), std::cout);
01177       }
01178     }
01179     return maxlength;
01180   }
01181 
01182 #endif
01183 
01184   template<class M,class X,class Y>
01185   template<class RowToDomain, class Solvers, class SubDomains>
01186   std::size_t SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>::assembleLocalProblems(const RowToDomain& rowToDomain, 
01187                                                                                                  const matrix_type& mat,
01188                                                                                                  Solvers& solvers,
01189                                                                                                  const SubDomains& subDomains,
01190                                                                                                  bool onTheFly)
01191   {
01192     typedef typename SubDomains::const_iterator DomainIterator;
01193     typedef typename Solvers::iterator SolverIterator;
01194     std::size_t maxlength = 0;
01195     
01196     if(onTheFly){
01197       for(DomainIterator domain=subDomains.begin();domain!=subDomains.end();++domain)
01198         maxlength=std::max(maxlength, domain->size());
01199     }else{
01200       // initialize the solvers of the local prolems.
01201       SolverIterator solver=solvers.begin();
01202       for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end();
01203           ++domain, ++solver){
01204         solver->setSubMatrix(mat, *domain);
01205         maxlength=std::max(maxlength, domain->size());
01206       }
01207     }
01208     
01209     return maxlength;
01210       
01211   }
01212 
01213   
01214   template<class M, class X, class TM, class TD, class TA>
01215   void SeqOverlappingSchwarz<M,X,TM,TD,TA>::apply(X& x, const X& b)
01216   {
01217     SeqOverlappingSchwarzApplier<SeqOverlappingSchwarz>::apply(*this, x, b);
01218   }
01219   
01220   template<class M, class X, class TM, class TD, class TA>
01221   template<bool forward>
01222   void SeqOverlappingSchwarz<M,X,TM,TD,TA>::apply(X& x, const X& b)
01223   {
01224     typedef typename X::block_type block;
01225     typedef slu_vector solver_vector;
01226     typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::solver_iterator iterator;
01227     typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::domain_iterator
01228       domain_iterator;
01229     
01230     OverlappingAssigner<TD> assigner(maxlength, mat, b, x);
01231         
01232     domain_iterator domain=IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::begin(subDomains);
01233     iterator solver = IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::begin(solvers);
01234     X v(x); // temporary for the update
01235     v=0;
01236     
01237     typedef typename AdderSelector<TM,X,TD >::Adder Adder;    
01238     Adder adder(v, x, assigner, relax);
01239     
01240     nnz=0;
01241     std::size_t no=0;
01242     for(;domain != IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::end(subDomains); ++domain){
01243       //Copy rhs to C-array for SuperLU
01244       std::for_each(domain->begin(), domain->end(), assigner);
01245       assigner.resetIndexForNextDomain();
01246       if(onTheFly){
01247         // Create the subdomain solver
01248         slu sdsolver;
01249         sdsolver.setSubMatrix(mat, *domain);
01250         // Apply
01251         sdsolver.apply(assigner.lhs(), assigner.rhs());
01252         //nnz+=sdsolver.nnz();
01253       }else{
01254         solver->apply(assigner.lhs(), assigner.rhs());
01255         //nnz+=solver->nnz();
01256         ++solver;
01257       }
01258       ++no;
01259       //Add relaxed correction to from SuperLU to v
01260       std::for_each(domain->begin(), domain->end(), adder);
01261       assigner.resetIndexForNextDomain();
01262       
01263     }
01264     nnz/=no;
01265     
01266     adder.axpy();
01267     assigner.deallocate();
01268   }
01269   
01270   template<class K, int n, class Al, class X, class Y>
01271   OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
01272   ::OverlappingAssigner(std::size_t maxlength, const BCRSMatrix<FieldMatrix<K,n,n>, Al>& mat_, 
01273     const X& b_, Y& x_) :
01274     mat(&mat_), 
01275     rhs_( new DynamicVector<field_type>(maxlength, 42) ),
01276     lhs_( new DynamicVector<field_type>(maxlength, -42) ),
01277     b(&b_),
01278     x(&x_),
01279     i(0),
01280     maxlength_(maxlength)
01281   {}
01282 
01283   template<class K, int n, class Al, class X, class Y>
01284   void
01285   OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
01286   ::deallocate()
01287   {
01288     delete rhs_;
01289     delete lhs_;
01290   }
01291     
01292   template<class K, int n, class Al, class X, class Y>
01293   void
01294   OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
01295   ::resetIndexForNextDomain()
01296   {
01297     i=0;
01298   }
01299   
01300   template<class K, int n, class Al, class X, class Y>
01301   DynamicVector<K> &
01302   OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
01303   ::lhs()
01304   {
01305     return *lhs_;
01306   }
01307     
01308   template<class K, int n, class Al, class X, class Y>
01309   DynamicVector<K> &   
01310   OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
01311   ::rhs()
01312   {
01313     return *rhs_;
01314   }
01315     
01316   template<class K, int n, class Al, class X, class Y>
01317   void
01318   OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
01319   ::relaxResult(field_type relax)
01320   {
01321     lhs() *= relax;
01322   }
01323   
01324   template<class K, int n, class Al, class X, class Y>
01325   void
01326   OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
01327   ::operator()(const size_type& domainIndex)
01328   {
01329     lhs() = 0.0;
01330 #if 0
01331     //assign right hand side of current domainindex block
01332     for(size_type j=0; j<n; ++j, ++i){
01333       assert(i<maxlength_);
01334       rhs()[i]=(*b)[domainIndex][j];
01335     }
01336     
01337     // loop over all Matrix row entries and calculate defect.
01338     typedef typename matrix_type::ConstColIterator col_iterator;
01339     
01340     // calculate defect for current row index block
01341     for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col){
01342       block_type tmp(0.0);
01343       (*col).mv((*x)[col.index()], tmp);
01344       i-=n;
01345       for(size_type j=0; j<n; ++j, ++i){
01346         assert(i<maxlength_);
01347         rhs()[i]-=tmp[j];
01348       }
01349     }
01350 #else
01351     //assign right hand side of current domainindex block
01352     for(size_type j=0; j<n; ++j, ++i){
01353       assert(i<maxlength_);
01354       rhs()[i]=(*b)[domainIndex][j];
01355       
01356       // loop over all Matrix row entries and calculate defect.
01357       typedef typename matrix_type::ConstColIterator col_iterator;
01358       
01359       // calculate defect for current row index block
01360       for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col){
01361         for(size_type k=0; k<n; ++k){
01362           rhs()[i]-=(*col)[j][k] * (*x)[col.index()][k];
01363         }
01364       }
01365     }
01366 #endif 
01367   }
01368 
01369   template<class K, int n, class Al, class X, class Y>
01370   void
01371   OverlappingAssigner< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y > >
01372   ::assignResult(block_type& res)
01373   {
01374     // assign the result of the local solve to the global vector
01375     for(size_type j=0; j<n; ++j, ++i){
01376       assert(i<maxlength_);
01377       res[j]+=lhs()[i];
01378     }
01379   }
01380 
01381 #if HAVE_SUPERLU
01382 
01383   template<typename T, typename A, int n, int m>
01384   OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >
01385   ::OverlappingAssigner(std::size_t maxlength, 
01386                         const BCRSMatrix<FieldMatrix<T,n,m>,A>& mat_, 
01387                         const range_type& b_,
01388                         range_type& x_)
01389     : mat(&mat_), 
01390       b(&b_), 
01391       x(&x_), i(0), maxlength_(maxlength)
01392   {
01393     rhs_ = new field_type[maxlength];
01394     lhs_ = new field_type[maxlength];
01395     
01396   }
01397   
01398   template<typename T, typename A, int n, int m>
01399   void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::deallocate()
01400   {
01401     delete[] rhs_;
01402     delete[] lhs_;
01403   }
01404   
01405   template<typename T, typename A, int n, int m>
01406   void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::operator()(const size_type& domainIndex)
01407   {
01408     //assign right hand side of current domainindex block
01409     // rhs is an array of doubles!
01410     // rhs[starti] = b[domainindex]
01411     for(size_type j=0; j<n; ++j, ++i){
01412       assert(i<maxlength_);
01413       rhs_[i]=(*b)[domainIndex][j];
01414     }
01415     
01416     
01417     // loop over all Matrix row entries and calculate defect.
01418     typedef typename matrix_type::ConstColIterator col_iterator;
01419 
01420     // calculate defect for current row index block
01421     for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col){
01422         block_type tmp;
01423         (*col).mv((*x)[col.index()], tmp);
01424         i-=n;
01425         for(size_type j=0; j<n; ++j, ++i){
01426           assert(i<maxlength_);
01427           rhs_[i]-=tmp[j];
01428         }
01429         
01430     }
01431     
01432   }
01433   template<typename T, typename A, int n, int m>
01434   void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::relaxResult(field_type relax)
01435   {
01436     for(size_type j=i+n; i<j; ++i){
01437       assert(i<maxlength_);
01438       lhs_[i]*=relax;
01439     }
01440     i-=n;
01441   }
01442   
01443   template<typename T, typename A, int n, int m>
01444   void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::assignResult(block_type& res)
01445   {
01446     // assign the result of the local solve to the global vector
01447     for(size_type j=0; j<n; ++j, ++i){
01448       assert(i<maxlength_);
01449       res[j]+=lhs_[i];
01450     }
01451   }
01452   
01453   template<typename T, typename A, int n, int m>
01454   void OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::resetIndexForNextDomain()
01455   {
01456     i=0;
01457   }
01458   
01459   template<typename T, typename A, int n, int m>
01460   typename OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::field_type* 
01461   OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::lhs()
01462   {
01463     return lhs_;
01464   }
01465   
01466   template<typename T, typename A, int n, int m>
01467   typename OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::field_type* 
01468   OverlappingAssigner<SuperLU<BCRSMatrix<FieldMatrix<T,n,m>,A> > >::rhs()
01469   {
01470     return rhs_;
01471   }
01472 
01473 #endif
01474 
01475   template<class M, class X, class Y>
01476   OverlappingAssignerILUBase<M,X,Y>::OverlappingAssignerILUBase(std::size_t maxlength, 
01477                                                                         const M& mat_,
01478                                                                         const Y& b_,
01479                                                                         X& x_)
01480     : mat(&mat_), 
01481       b(&b_), 
01482       x(&x_), i(0)
01483   {
01484     rhs_= new Y(maxlength);
01485     lhs_ = new X(maxlength);
01486   }
01487   
01488   template<class M, class X, class Y>
01489   void OverlappingAssignerILUBase<M,X,Y>::deallocate()
01490   {
01491     delete rhs_;
01492     delete lhs_;
01493   }
01494   
01495   template<class M, class X, class Y>
01496   void OverlappingAssignerILUBase<M,X,Y>::operator()(const size_type& domainIndex)
01497   {
01498     (*rhs_)[i]=(*b)[domainIndex];
01499     
01500      // loop over all Matrix row entries and calculate defect.
01501     typedef typename matrix_type::ConstColIterator col_iterator;
01502 
01503     // calculate defect for current row index block
01504     for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col){
01505       (*col).mmv((*x)[col.index()], (*rhs_)[i]);
01506     }
01507     // Goto next local index
01508     ++i;
01509   }
01510   
01511   template<class M, class X, class Y>
01512   void OverlappingAssignerILUBase<M,X,Y>::relaxResult(field_type relax)
01513   {
01514     (*lhs_)[i]*=relax;
01515   }
01516   
01517   template<class M, class X, class Y>
01518   void OverlappingAssignerILUBase<M,X,Y>::assignResult(block_type& res)
01519   {
01520     res+=(*lhs_)[i++];
01521   }
01522 
01523   template<class M, class X, class Y>
01524   X& OverlappingAssignerILUBase<M,X,Y>::lhs()
01525   {
01526     return *lhs_;
01527   }
01528   
01529   template<class M, class X, class Y>
01530   Y& OverlappingAssignerILUBase<M,X,Y>::rhs()
01531   {
01532     return *rhs_;
01533   }
01534 
01535   template<class M, class X, class Y>
01536   void OverlappingAssignerILUBase<M,X,Y>::resetIndexForNextDomain()
01537   {
01538     i=0;
01539   }
01540   
01541   template<typename S, typename T, typename A, int n>
01542   AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::AdditiveAdder(BlockVector<FieldVector<T,n>,A>& v_, 
01543                                                                    BlockVector<FieldVector<T,n>,A>& x_, 
01544                                                                    OverlappingAssigner<S>& assigner_, 
01545                                                                    const T& relax_)
01546     : v(&v_), x(&x_), assigner(&assigner_), relax(relax_)
01547   {}
01548 
01549   template<typename S, typename T, typename A, int n>
01550   void AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::operator()(const size_type& domainIndex)
01551   {
01552     // add the result of the local solve to the current update
01553     assigner->assignResult((*v)[domainIndex]);
01554   }
01555 
01556  
01557   template<typename S, typename T, typename A, int n>
01558   void AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::axpy()
01559   {
01560     // relax the update and add it to the current guess.
01561     x->axpy(relax,*v);
01562   }
01563 
01564  
01565   template<typename S, typename T, typename A, int n>
01566   MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >
01567   ::MultiplicativeAdder(BlockVector<FieldVector<T,n>,A>& v_, 
01568                         BlockVector<FieldVector<T,n>,A>& x_, 
01569                         OverlappingAssigner<S>& assigner_, const T& relax_)
01570     : x(&x_), assigner(&assigner_), relax(relax_)
01571   {}
01572 
01573   
01574   template<typename S,typename T, typename A, int n>
01575   void MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >::operator()(const size_type& domainIndex)
01576   {
01577     // add the result of the local solve to the current guess
01578     assigner->relaxResult(relax);
01579     assigner->assignResult((*x)[domainIndex]);
01580   }
01581 
01582   
01583   template<typename S,typename T, typename A, int n>
01584   void MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >::axpy()
01585   {
01586     // nothing to do, as the corrections already relaxed and added in operator()
01587   }
01588   
01589   
01591 }
01592 
01593 #endif