dune-istl
2.2.0
|
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