Dune Core Modules (2.4.1)

overlappingschwarz.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
4#define DUNE_ISTL_OVERLAPPINGSCHWARZ_HH
5#include <cassert>
6#include <algorithm>
7#include <functional>
8#include <vector>
9#include <set>
11#include <dune/common/sllist.hh>
12#include <dune/common/unused.hh>
13#include "preconditioners.hh"
14#include "superlu.hh"
15#include "umfpack.hh"
16#include "bvector.hh"
17#include "bcrsmatrix.hh"
18#include "ilusubdomainsolver.hh"
20
21namespace Dune
22{
23
35 template<class M, class X, class TM, class TD, class TA>
36 class SeqOverlappingSchwarz;
37
41 template<class I, class S, class D>
43 {
44 public:
47
48 typedef I InitializerList;
49 typedef typename InitializerList::value_type AtomInitializer;
50 typedef typename AtomInitializer::Matrix Matrix;
51 typedef typename Matrix::const_iterator Iter;
52 typedef typename Matrix::row_type::const_iterator CIter;
53
54 typedef S IndexSet;
55 typedef typename IndexSet::size_type size_type;
56
57 OverlappingSchwarzInitializer(InitializerList& il,
58 const IndexSet& indices,
59 const subdomain_vector& domains);
60
61
62 void addRowNnz(const Iter& row);
63
64 void allocate();
65
66 void countEntries(const Iter& row, const CIter& col) const;
67
68 void calcColstart() const;
69
70 void copyValue(const Iter& row, const CIter& col) const;
71
72 void createMatrix() const;
73 private:
74 class IndexMap
75 {
76 public:
77 typedef typename S::size_type size_type;
78 typedef std::map<size_type,size_type> Map;
79 typedef typename Map::iterator iterator;
80 typedef typename Map::const_iterator const_iterator;
81
82 IndexMap();
83
84 void insert(size_type grow);
85
86 const_iterator find(size_type grow) const;
87
88 iterator find(size_type grow);
89
90 iterator begin();
91
92 const_iterator begin() const;
93
94 iterator end();
95
96 const_iterator end() const;
97
98 private:
99 std::map<size_type,size_type> map_;
100 size_type row;
101 };
102
103
104 typedef typename InitializerList::iterator InitIterator;
105 typedef typename IndexSet::const_iterator IndexIteratur;
106 InitializerList* initializers;
107 const IndexSet *indices;
108 mutable std::vector<IndexMap> indexMaps;
109 const subdomain_vector& domains;
110 };
111
116 {};
117
122 {};
123
129 {};
130
135 template<class M, class X, class Y>
137
138 // Specialization for BCRSMatrix
139 template<class K, int n, class Al, class X, class Y>
140 class DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >
141 {
142 typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> M;
143 public:
145 typedef typename Dune::remove_const<M>::type matrix_type;
146 typedef K field_type;
147 typedef typename Dune::remove_const<M>::type rilu_type;
149 typedef X domain_type;
151 typedef Y range_type;
152
158 {
159 assert(v.size() > 0);
160 assert(v.size() == d.size());
161 assert(A.rows() <= v.size());
162 assert(A.cols() <= v.size());
163 size_t sz = A.rows();
164 v.resize(sz);
165 d.resize(sz);
166 A.solve(v,d);
167 v.resize(v.capacity());
168 d.resize(d.capacity());
169 }
170
178 template<class S>
179 void setSubMatrix(const M& BCRS, S& rowset)
180 {
181 size_t sz = rowset.size();
182 A.resize(sz*n,sz*n);
183 typedef typename S::const_iterator SIter;
184 size_t r = 0, c = 0;
185 for(SIter rowIdx = rowset.begin(), rowEnd=rowset.end();
186 rowIdx!= rowEnd; ++rowIdx, r++)
187 {
188 c = 0;
189 for(SIter colIdx = rowset.begin(), colEnd=rowset.end();
190 colIdx!= colEnd; ++colIdx, c++)
191 {
192 if (BCRS[*rowIdx].find(*colIdx) == BCRS[*rowIdx].end())
193 continue;
194 for (size_t i=0; i<n; i++)
195 {
196 for (size_t j=0; j<n; j++)
197 {
198 A[r*n+i][c*n+j] = BCRS[*rowIdx][*colIdx][i][j];
199 }
200 }
201 }
202 }
203 }
204 private:
205 DynamicMatrix<K> A;
206 };
207
208 template<typename T, bool tag>
209 class OverlappingAssignerHelper
210 {};
211
212 template<typename T>
213 struct OverlappingAssigner : public OverlappingAssignerHelper<T,Dune::StoresColumnCompressed<T>::value >
214 {
215 public:
216 OverlappingAssigner(std::size_t maxlength, const typename T::matrix_type& mat,
217 const typename T::range_type& b, typename T::range_type& x)
218 : OverlappingAssignerHelper<T,Dune::StoresColumnCompressed<T>::value >
219 (maxlength,mat,b,x)
220 {}
221 };
222
223 // specialization for DynamicMatrix
224 template<class K, int n, class Al, class X, class Y>
225 class OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
226 {
227 public:
228 typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> matrix_type;
229 typedef K field_type;
230 typedef Y range_type;
231 typedef typename range_type::block_type block_type;
232 typedef typename matrix_type::size_type size_type;
233
241 OverlappingAssignerHelper(std::size_t maxlength, const BCRSMatrix<FieldMatrix<K,n,n>, Al>& mat_, const X& b_, Y& x_);
242
246 inline
247 void deallocate();
248
252 inline
253 void resetIndexForNextDomain();
254
259 inline
260 DynamicVector<K> & lhs();
261
266 inline
267 DynamicVector<K> & rhs();
268
273 inline
274 void relaxResult(field_type relax);
275
280 void operator()(const size_type& domainIndex);
281
289 inline
290 void assignResult(block_type& res);
291
292 private:
296 const matrix_type* mat;
298 // we need a pointer, because we have to avoid deep copies
299 DynamicVector<field_type> * rhs_;
301 // we need a pointer, because we have to avoid deep copies
302 DynamicVector<field_type> * lhs_;
304 const range_type* b;
306 range_type* x;
308 std::size_t i;
310 std::size_t maxlength_;
311 };
312
313#if HAVE_SUPERLU || HAVE_UMFPACK
314 template<template<class> class S, int n, int m, typename T, typename A>
315 struct OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>, A> >, true>
316 {
317 typedef BCRSMatrix<FieldMatrix<T,n,m>, A> matrix_type;
318 typedef typename S<BCRSMatrix<FieldMatrix<T,n,m>, A> >::range_type range_type;
319 typedef typename range_type::field_type field_type;
320 typedef typename range_type::block_type block_type;
321
322 typedef typename matrix_type::size_type size_type;
323
331 OverlappingAssignerHelper(std::size_t maxlength, const matrix_type& mat,
332 const range_type& b, range_type& x);
338 void deallocate();
339
340 /*
341 * @brief Resets the local index to zero.
342 */
343 void resetIndexForNextDomain();
344
349 field_type* lhs();
350
355 field_type* rhs();
356
361 void relaxResult(field_type relax);
362
367 void operator()(const size_type& domain);
368
376 void assignResult(block_type& res);
377
378 private:
382 const matrix_type* mat;
384 field_type* rhs_;
386 field_type* lhs_;
388 const range_type* b;
390 range_type* x;
392 std::size_t i;
394 std::size_t maxlength_;
395 };
396
397#endif
398
399 template<class M, class X, class Y>
400 class OverlappingAssignerILUBase
401 {
402 public:
403 typedef M matrix_type;
404
405 typedef typename M::field_type field_type;
406
407 typedef typename Y::block_type block_type;
408
409 typedef typename matrix_type::size_type size_type;
417 OverlappingAssignerILUBase(std::size_t maxlength, const M& mat,
418 const Y& b, X& x);
424 void deallocate();
425
429 void resetIndexForNextDomain();
430
435 X& lhs();
436
441 Y& rhs();
442
447 void relaxResult(field_type relax);
448
453 void operator()(const size_type& domain);
454
462 void assignResult(block_type& res);
463
464 private:
468 const M* mat;
470 X* lhs_;
472 Y* rhs_;
474 const Y* b;
476 X* x;
478 size_type i;
479 };
480
481 // specialization for ILU0
482 template<class M, class X, class Y>
483 class OverlappingAssignerHelper<ILU0SubdomainSolver<M,X,Y>, false>
484 : public OverlappingAssignerILUBase<M,X,Y>
485 {
486 public:
494 OverlappingAssignerHelper(std::size_t maxlength, const M& mat,
495 const Y& b, X& x)
496 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
497 {}
498 };
499
500 // specialization for ILUN
501 template<class M, class X, class Y>
502 class OverlappingAssignerHelper<ILUNSubdomainSolver<M,X,Y>,false>
503 : public OverlappingAssignerILUBase<M,X,Y>
504 {
505 public:
513 OverlappingAssignerHelper(std::size_t maxlength, const M& mat,
514 const Y& b, X& x)
515 : OverlappingAssignerILUBase<M,X,Y>(maxlength, mat,b,x)
516 {}
517 };
518
519 template<typename S, typename T>
520 struct AdditiveAdder
521 {};
522
523 template<typename S, typename T, typename A, int n>
524 struct AdditiveAdder<S, BlockVector<FieldVector<T,n>,A> >
525 {
526 typedef typename A::size_type size_type;
527 AdditiveAdder(BlockVector<FieldVector<T,n>,A>& v, BlockVector<FieldVector<T,n>,A>& x,
528 OverlappingAssigner<S>& assigner, const T& relax_);
529 void operator()(const size_type& domain);
530 void axpy();
531
532 private:
533 BlockVector<FieldVector<T,n>,A>* v;
534 BlockVector<FieldVector<T,n>,A>* x;
535 OverlappingAssigner<S>* assigner;
536 T relax;
537 };
538
539 template<typename S,typename T>
540 struct MultiplicativeAdder
541 {};
542
543 template<typename S, typename T, typename A, int n>
544 struct MultiplicativeAdder<S, BlockVector<FieldVector<T,n>,A> >
545 {
546 typedef typename A::size_type size_type;
547 MultiplicativeAdder(BlockVector<FieldVector<T,n>,A>& v, BlockVector<FieldVector<T,n>,A>& x,
548 OverlappingAssigner<S>& assigner_, const T& relax_);
549 void operator()(const size_type& domain);
550 void axpy();
551
552 private:
553 BlockVector<FieldVector<T,n>,A>* x;
554 OverlappingAssigner<S>* assigner;
555 T relax;
556 };
557
567 template<typename T, class X, class S>
569 {};
570
571 template<class X, class S>
573 {
574 typedef AdditiveAdder<S,X> Adder;
575 };
576
577 template<class X, class S>
578 struct AdderSelector<MultiplicativeSchwarzMode,X,S>
579 {
580 typedef MultiplicativeAdder<S,X> Adder;
581 };
582
583 template<class X, class S>
584 struct AdderSelector<SymmetricMultiplicativeSchwarzMode,X,S>
585 {
586 typedef MultiplicativeAdder<S,X> Adder;
587 };
588
600 template<typename T1, typename T2, bool forward>
602 {
603 typedef T1 solver_vector;
604 typedef typename solver_vector::iterator solver_iterator;
605 typedef T2 subdomain_vector;
606 typedef typename subdomain_vector::const_iterator domain_iterator;
607
608 static solver_iterator begin(solver_vector& sv)
609 {
610 return sv.begin();
611 }
612
613 static solver_iterator end(solver_vector& sv)
614 {
615 return sv.end();
616 }
617 static domain_iterator begin(const subdomain_vector& sv)
618 {
619 return sv.begin();
620 }
621
622 static domain_iterator end(const subdomain_vector& sv)
623 {
624 return sv.end();
625 }
626 };
627
628 template<typename T1, typename T2>
629 struct IteratorDirectionSelector<T1,T2,false>
630 {
631 typedef T1 solver_vector;
632 typedef typename solver_vector::reverse_iterator solver_iterator;
633 typedef T2 subdomain_vector;
634 typedef typename subdomain_vector::const_reverse_iterator domain_iterator;
635
636 static solver_iterator begin(solver_vector& sv)
637 {
638 return sv.rbegin();
639 }
640
641 static solver_iterator end(solver_vector& sv)
642 {
643 return sv.rend();
644 }
645 static domain_iterator begin(const subdomain_vector& sv)
646 {
647 return sv.rbegin();
648 }
649
650 static domain_iterator end(const subdomain_vector& sv)
651 {
652 return sv.rend();
653 }
654 };
655
664 template<class T>
666 {
667 typedef T smoother;
668 typedef typename smoother::range_type range_type;
669
670 static void apply(smoother& sm, range_type& v, const range_type& b)
671 {
672 sm.template apply<true>(v, b);
673 }
674 };
675
676 template<class M, class X, class TD, class TA>
678 {
680 typedef typename smoother::range_type range_type;
681
682 static void apply(smoother& sm, range_type& v, const range_type& b)
683 {
684 sm.template apply<true>(v, b);
685 sm.template apply<false>(v, b);
686 }
687 };
688
689 template<class T, bool tag>
690 struct SeqOverlappingSchwarzAssemblerHelper
691 {};
692
693 template<class T>
694 struct SeqOverlappingSchwarzAssembler
695 : public SeqOverlappingSchwarzAssemblerHelper<T,Dune::StoresColumnCompressed<T>::value>
696 {};
697
698
699 template<class K, int n, class Al, class X, class Y>
700 struct SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
701 {
702 typedef BCRSMatrix< FieldMatrix<K,n,n>, Al> matrix_type;
703 template<class RowToDomain, class Solvers, class SubDomains>
704 static std::size_t assembleLocalProblems(const RowToDomain& rowToDomain, const matrix_type& mat,
705 Solvers& solvers, const SubDomains& domains,
706 bool onTheFly);
707 };
708
709 template<template<class> class S, typename T, typename A, int m, int n>
710 struct SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<FieldMatrix<T,m,n>,A> >,true>
711 {
712 typedef BCRSMatrix<FieldMatrix<T,m,n>,A> matrix_type;
713 template<class RowToDomain, class Solvers, class SubDomains>
714 static std::size_t assembleLocalProblems(const RowToDomain& rowToDomain, const matrix_type& mat,
715 Solvers& solvers, const SubDomains& domains,
716 bool onTheFly);
717 };
718
719 template<class M,class X, class Y>
720 struct SeqOverlappingSchwarzAssemblerILUBase
721 {
722 typedef M matrix_type;
723 template<class RowToDomain, class Solvers, class SubDomains>
724 static std::size_t assembleLocalProblems(const RowToDomain& rowToDomain, const matrix_type& mat,
725 Solvers& solvers, const SubDomains& domains,
726 bool onTheFly);
727 };
728
729 template<class M,class X, class Y>
730 struct SeqOverlappingSchwarzAssemblerHelper<ILU0SubdomainSolver<M,X,Y>,false>
731 : public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
732 {};
733
734 template<class M,class X, class Y>
735 struct SeqOverlappingSchwarzAssemblerHelper<ILUNSubdomainSolver<M,X,Y>,false>
736 : public SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>
737 {};
738
749 template<class M, class X, class TM=AdditiveSchwarzMode,
750 class TD=ILU0SubdomainSolver<M,X,X>, class TA=std::allocator<X> >
752 : public Preconditioner<X,X>
753 {
754 public:
758 typedef M matrix_type;
759
763 typedef X domain_type;
764
768 typedef X range_type;
769
776 typedef TM Mode;
777
781 typedef typename X::field_type field_type;
782
784 typedef typename matrix_type::size_type size_type;
785
787 typedef TA allocator;
788
790 typedef std::set<size_type, std::less<size_type>,
791 typename TA::template rebind<size_type>::other>
793
795 typedef std::vector<subdomain_type, typename TA::template rebind<subdomain_type>::other> subdomain_vector;
796
799
801 typedef std::vector<subdomain_list, typename TA::template rebind<subdomain_list>::other > rowtodomain_vector;
802
804 typedef TD slu;
805
807 typedef std::vector<slu, typename TA::template rebind<slu>::other> slu_vector;
808
809 enum {
812 };
813
828 field_type relaxationFactor=1, bool onTheFly_=true);
829
842 field_type relaxationFactor=1, bool onTheFly_=true);
843
849 virtual void pre (X& x, X& b)
850 {
853 }
854
860 virtual void apply (X& v, const X& d);
861
867 virtual void post (X& x)
868 {
870 }
871
872 template<bool forward>
873 void apply(X& v, const X& d);
874
875 private:
876 const M& mat;
877 slu_vector solvers;
878 subdomain_vector subDomains;
879 field_type relax;
880
881 typename M::size_type maxlength;
882
883 bool onTheFly;
884 };
885
886
887
888 template<class I, class S, class D>
890 const IndexSet& idx,
891 const subdomain_vector& domains_)
892 : initializers(&il), indices(&idx), indexMaps(il.size()), domains(domains_)
893 {}
894
895
896 template<class I, class S, class D>
897 void OverlappingSchwarzInitializer<I,S,D>::addRowNnz(const Iter& row)
898 {
899 typedef typename IndexSet::value_type::const_iterator iterator;
900 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
901 (*initializers)[*domain].addRowNnz(row, domains[*domain]);
902 indexMaps[*domain].insert(row.index());
903 }
904 }
905
906 template<class I, class S, class D>
907 void OverlappingSchwarzInitializer<I,S,D>::allocate()
908 {
909 std::for_each(initializers->begin(), initializers->end(),
910 std::mem_fun_ref(&AtomInitializer::allocateMatrixStorage));
911 std::for_each(initializers->begin(), initializers->end(),
912 std::mem_fun_ref(&AtomInitializer::allocateMarker));
913 }
914
915 template<class I, class S, class D>
916 void OverlappingSchwarzInitializer<I,S,D>::countEntries(const Iter& row, const CIter& col) const
917 {
918 typedef typename IndexSet::value_type::const_iterator iterator;
919 for(iterator domain=(*indices)[row.index()].begin(); domain != (*indices)[row.index()].end(); ++domain) {
920 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
921 if(v!= indexMaps[*domain].end()) {
922 (*initializers)[*domain].countEntries(indexMaps[*domain].find(col.index())->second);
923 }
924 }
925 }
926
927 template<class I, class S, class D>
928 void OverlappingSchwarzInitializer<I,S,D>::calcColstart() const
929 {
930 std::for_each(initializers->begin(), initializers->end(),
931 std::mem_fun_ref(&AtomInitializer::calcColstart));
932 }
933
934 template<class I, class S, class D>
935 void OverlappingSchwarzInitializer<I,S,D>::copyValue(const Iter& row, const CIter& col) const
936 {
937 typedef typename IndexSet::value_type::const_iterator iterator;
938 for(iterator domain=(*indices)[row.index()].begin(); domain!= (*indices)[row.index()].end(); ++domain) {
939 typename std::map<size_type,size_type>::const_iterator v = indexMaps[*domain].find(col.index());
940 if(v!= indexMaps[*domain].end()) {
941 assert(indexMaps[*domain].end()!=indexMaps[*domain].find(row.index()));
942 (*initializers)[*domain].copyValue(col, indexMaps[*domain].find(row.index())->second,
943 v->second);
944 }
945 }
946 }
947
948 template<class I, class S, class D>
949 void OverlappingSchwarzInitializer<I,S,D>::createMatrix() const
950 {
951 std::vector<IndexMap>().swap(indexMaps);
952 std::for_each(initializers->begin(), initializers->end(),
953 std::mem_fun_ref(&AtomInitializer::createMatrix));
954 }
955
956 template<class I, class S, class D>
957 OverlappingSchwarzInitializer<I,S,D>::IndexMap::IndexMap()
958 : row(0)
959 {}
960
961 template<class I, class S, class D>
962 void OverlappingSchwarzInitializer<I,S,D>::IndexMap::insert(size_type grow)
963 {
964 assert(map_.find(grow)==map_.end());
965 map_.insert(std::make_pair(grow, row++));
966 }
967
968 template<class I, class S, class D>
969 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
970 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow) const
971 {
972 return map_.find(grow);
973 }
974
975 template<class I, class S, class D>
976 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
977 OverlappingSchwarzInitializer<I,S,D>::IndexMap::find(size_type grow)
978 {
979 return map_.find(grow);
980 }
981
982 template<class I, class S, class D>
983 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
984 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end() const
985 {
986 return map_.end();
987 }
988
989 template<class I, class S, class D>
990 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
991 OverlappingSchwarzInitializer<I,S,D>::IndexMap::end()
992 {
993 return map_.end();
994 }
995
996 template<class I, class S, class D>
997 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::const_iterator
998 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin() const
999 {
1000 return map_.begin();
1001 }
1002
1003 template<class I, class S, class D>
1004 typename OverlappingSchwarzInitializer<I,S,D>::IndexMap::iterator
1005 OverlappingSchwarzInitializer<I,S,D>::IndexMap::begin()
1006 {
1007 return map_.begin();
1008 }
1009
1010 template<class M, class X, class TM, class TD, class TA>
1012 field_type relaxationFactor, bool fly)
1013 : mat(mat_), relax(relaxationFactor), onTheFly(fly)
1014 {
1015 typedef typename rowtodomain_vector::const_iterator RowDomainIterator;
1016 typedef typename subdomain_list::const_iterator DomainIterator;
1017#ifdef DUNE_ISTL_WITH_CHECKING
1018 assert(rowToDomain.size()==mat.N());
1019 assert(rowToDomain.size()==mat.M());
1020
1021 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1022 assert(iter->size()>0);
1023
1024#endif
1025 // calculate the number of domains
1026 size_type domains=0;
1027 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter)
1028 for(DomainIterator d=iter->begin(); d != iter->end(); ++d)
1029 domains=std::max(domains, *d);
1030 ++domains;
1031
1032 solvers.resize(domains);
1033 subDomains.resize(domains);
1034
1035 // initialize subdomains to row mapping from row to subdomain mapping
1036 size_type row=0;
1037 for(RowDomainIterator iter=rowToDomain.begin(); iter != rowToDomain.end(); ++iter, ++row)
1038 for(DomainIterator d=iter->begin(); d != iter->end(); ++d)
1039 subDomains[*d].insert(row);
1040
1041#ifdef DUNE_ISTL_WITH_CHECKING
1042 size_type i=0;
1043 typedef typename subdomain_vector::const_iterator iterator;
1044 for(iterator iter=subDomains.begin(); iter != subDomains.end(); ++iter) {
1045 typedef typename subdomain_type::const_iterator entry_iterator;
1046 Dune::dvverb<<"domain "<<i++<<":";
1047 for(entry_iterator entry = iter->begin(); entry != iter->end(); ++entry) {
1048 Dune::dvverb<<" "<<*entry;
1049 }
1050 Dune::dvverb<<std::endl;
1051 }
1052#endif
1053 maxlength = SeqOverlappingSchwarzAssembler<slu>
1054 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1055 }
1056
1057 template<class M, class X, class TM, class TD, class TA>
1059 const subdomain_vector& sd,
1060 field_type relaxationFactor,
1061 bool fly)
1062 : mat(mat_), solvers(sd.size()), subDomains(sd), relax(relaxationFactor),
1063 onTheFly(fly)
1064 {
1065 typedef typename subdomain_vector::const_iterator DomainIterator;
1066
1067#ifdef DUNE_ISTL_WITH_CHECKING
1068 size_type i=0;
1069
1070 for(DomainIterator d=sd.begin(); d != sd.end(); ++d,++i) {
1071 //std::cout<<i<<": "<<d->size()<<std::endl;
1072 assert(d->size()>0);
1073 typedef typename DomainIterator::value_type::const_iterator entry_iterator;
1074 Dune::dvverb<<"domain "<<i<<":";
1075 for(entry_iterator entry = d->begin(); entry != d->end(); ++entry) {
1076 Dune::dvverb<<" "<<*entry;
1077 }
1078 Dune::dvverb<<std::endl;
1079 }
1080
1081#endif
1082
1083 // Create a row to subdomain mapping
1084 rowtodomain_vector rowToDomain(mat.N());
1085
1086 size_type domainId=0;
1087
1088 for(DomainIterator domain=sd.begin(); domain != sd.end(); ++domain, ++domainId) {
1089 typedef typename subdomain_type::const_iterator iterator;
1090 for(iterator row=domain->begin(); row != domain->end(); ++row)
1091 rowToDomain[*row].push_back(domainId);
1092 }
1093
1094 maxlength = SeqOverlappingSchwarzAssembler<slu>
1095 ::assembleLocalProblems(rowToDomain, mat, solvers, subDomains, onTheFly);
1096 }
1097
1104 template<class M>
1106
1107 template<typename T, typename A, int n, int m>
1109 {
1110 template<class Domain>
1111 static int size(const Domain & d)
1112 {
1113 assert(n==m);
1114 return m*d.size();
1115 }
1116 };
1117
1118 template<class K, int n, class Al, class X, class Y>
1119 template<class RowToDomain, class Solvers, class SubDomains>
1120 std::size_t
1121 SeqOverlappingSchwarzAssemblerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>::
1122 assembleLocalProblems(const RowToDomain& rowToDomain,
1123 const matrix_type& mat,
1124 Solvers& solvers,
1125 const SubDomains& subDomains,
1126 bool onTheFly)
1127 {
1128 DUNE_UNUSED_PARAMETER(onTheFly);
1129 DUNE_UNUSED_PARAMETER(rowToDomain);
1131 DUNE_UNUSED_PARAMETER(solvers);
1132 typedef typename SubDomains::const_iterator DomainIterator;
1133 std::size_t maxlength = 0;
1134
1135 assert(onTheFly);
1136
1137 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1138 maxlength=std::max(maxlength, domain->size());
1139 maxlength*=n;
1140
1141 return maxlength;
1142 }
1143
1144#if HAVE_SUPERLU || HAVE_UMFPACK
1145 template<template<class> class S, typename T, typename A, int m, int n>
1146 template<class RowToDomain, class Solvers, class SubDomains>
1147 std::size_t SeqOverlappingSchwarzAssemblerHelper<S<BCRSMatrix<FieldMatrix<T,m,n>,A> >,true>::assembleLocalProblems(const RowToDomain& rowToDomain,
1148 const matrix_type& mat,
1149 Solvers& solvers,
1150 const SubDomains& subDomains,
1151 bool onTheFly)
1152 {
1153 typedef typename S<BCRSMatrix<FieldMatrix<T,m,n>,A> >::MatrixInitializer MatrixInitializer;
1154 typedef typename std::vector<MatrixInitializer>::iterator InitializerIterator;
1155 typedef typename SubDomains::const_iterator DomainIterator;
1156 typedef typename Solvers::iterator SolverIterator;
1157 std::size_t maxlength = 0;
1158
1159 if(onTheFly) {
1160 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1161 maxlength=std::max(maxlength, domain->size());
1162 maxlength*=mat[0].begin()->N();
1163 }else{
1164 // initialize the initializers
1165 DomainIterator domain=subDomains.begin();
1166
1167 // Create the initializers list.
1168 std::vector<MatrixInitializer> initializers(subDomains.size());
1169
1170 SolverIterator solver=solvers.begin();
1171 for(InitializerIterator initializer=initializers.begin(); initializer!=initializers.end();
1172 ++initializer, ++solver, ++domain) {
1173 solver->getInternalMatrix().N_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1174 solver->getInternalMatrix().M_=SeqOverlappingSchwarzDomainSize<matrix_type>::size(*domain);
1175 //solver->setVerbosity(true);
1176 *initializer=MatrixInitializer(solver->getInternalMatrix());
1177 }
1178
1179 // Set up the supermatrices according to the subdomains
1180 typedef OverlappingSchwarzInitializer<std::vector<MatrixInitializer>,
1181 RowToDomain, SubDomains> Initializer;
1182
1183 Initializer initializer(initializers, rowToDomain, subDomains);
1184 copyToColCompMatrix(initializer, mat);
1185
1186 // Calculate the LU decompositions
1187 std::for_each(solvers.begin(), solvers.end(), std::mem_fun_ref(&S<BCRSMatrix<FieldMatrix<T,m,n>,A> >::decompose));
1188 for (SolverIterator solverIt = solvers.begin(); solverIt != solvers.end(); ++solverIt)
1189 {
1190 assert(solverIt->getInternalMatrix().N() == solverIt->getInternalMatrix().M());
1191 maxlength = std::max(maxlength, solverIt->getInternalMatrix().N());
1192 }
1193 }
1194 return maxlength;
1195 }
1196
1197#endif
1198
1199 template<class M,class X,class Y>
1200 template<class RowToDomain, class Solvers, class SubDomains>
1201 std::size_t SeqOverlappingSchwarzAssemblerILUBase<M,X,Y>::assembleLocalProblems(const RowToDomain& rowToDomain,
1202 const matrix_type& mat,
1203 Solvers& solvers,
1204 const SubDomains& subDomains,
1205 bool onTheFly)
1206 {
1207 DUNE_UNUSED_PARAMETER(rowToDomain);
1208 typedef typename SubDomains::const_iterator DomainIterator;
1209 typedef typename Solvers::iterator SolverIterator;
1210 std::size_t maxlength = 0;
1211
1212 if(onTheFly) {
1213 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end(); ++domain)
1214 maxlength=std::max(maxlength, domain->size());
1215 }else{
1216 // initialize the solvers of the local prolems.
1217 SolverIterator solver=solvers.begin();
1218 for(DomainIterator domain=subDomains.begin(); domain!=subDomains.end();
1219 ++domain, ++solver) {
1220 solver->setSubMatrix(mat, *domain);
1221 maxlength=std::max(maxlength, domain->size());
1222 }
1223 }
1224
1225 return maxlength;
1226
1227 }
1228
1229
1230 template<class M, class X, class TM, class TD, class TA>
1232 {
1234 }
1235
1236 template<class M, class X, class TM, class TD, class TA>
1237 template<bool forward>
1238 void SeqOverlappingSchwarz<M,X,TM,TD,TA>::apply(X& x, const X& b)
1239 {
1240 typedef slu_vector solver_vector;
1241 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::solver_iterator iterator;
1242 typedef typename IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::domain_iterator
1243 domain_iterator;
1244
1245 OverlappingAssigner<TD> assigner(maxlength, mat, b, x);
1246
1249 X v(x); // temporary for the update
1250 v=0;
1251
1252 typedef typename AdderSelector<TM,X,TD >::Adder Adder;
1253 Adder adder(v, x, assigner, relax);
1254
1255 for(; domain != IteratorDirectionSelector<solver_vector,subdomain_vector,forward>::end(subDomains); ++domain) {
1256 //Copy rhs to C-array for SuperLU
1257 std::for_each(domain->begin(), domain->end(), assigner);
1258 assigner.resetIndexForNextDomain();
1259 if(onTheFly) {
1260 // Create the subdomain solver
1261 slu sdsolver;
1262 sdsolver.setSubMatrix(mat, *domain);
1263 // Apply
1264 sdsolver.apply(assigner.lhs(), assigner.rhs());
1265 }else{
1266 solver->apply(assigner.lhs(), assigner.rhs());
1267 ++solver;
1268 }
1269
1270 //Add relaxed correction to from SuperLU to v
1271 std::for_each(domain->begin(), domain->end(), adder);
1272 assigner.resetIndexForNextDomain();
1273
1274 }
1275
1276 adder.axpy();
1277 assigner.deallocate();
1278 }
1279
1280 template<class K, int n, class Al, class X, class Y>
1281 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
1282 ::OverlappingAssignerHelper(std::size_t maxlength, const BCRSMatrix<FieldMatrix<K,n,n>, Al>& mat_,
1283 const X& b_, Y& x_) :
1284 mat(&mat_),
1285 rhs_( new DynamicVector<field_type>(maxlength, 42) ),
1286 lhs_( new DynamicVector<field_type>(maxlength, -42) ),
1287 b(&b_),
1288 x(&x_),
1289 i(0),
1290 maxlength_(maxlength)
1291 {}
1292
1293 template<class K, int n, class Al, class X, class Y>
1294 void
1295 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
1296 ::deallocate()
1297 {
1298 delete rhs_;
1299 delete lhs_;
1300 }
1301
1302 template<class K, int n, class Al, class X, class Y>
1303 void
1304 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
1305 ::resetIndexForNextDomain()
1306 {
1307 i=0;
1308 }
1309
1310 template<class K, int n, class Al, class X, class Y>
1312 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
1313 ::lhs()
1314 {
1315 return *lhs_;
1316 }
1317
1318 template<class K, int n, class Al, class X, class Y>
1320 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
1321 ::rhs()
1322 {
1323 return *rhs_;
1324 }
1325
1326 template<class K, int n, class Al, class X, class Y>
1327 void
1328 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
1329 ::relaxResult(field_type relax)
1330 {
1331 lhs() *= relax;
1332 }
1333
1334 template<class K, int n, class Al, class X, class Y>
1335 void
1336 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
1337 ::operator()(const size_type& domainIndex)
1338 {
1339 lhs() = 0.0;
1340#if 0
1341 //assign right hand side of current domainindex block
1342 for(size_type j=0; j<n; ++j, ++i) {
1343 assert(i<maxlength_);
1344 rhs()[i]=(*b)[domainIndex][j];
1345 }
1346
1347 // loop over all Matrix row entries and calculate defect.
1348 typedef typename matrix_type::ConstColIterator col_iterator;
1349
1350 // calculate defect for current row index block
1351 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1352 block_type tmp(0.0);
1353 (*col).mv((*x)[col.index()], tmp);
1354 i-=n;
1355 for(size_type j=0; j<n; ++j, ++i) {
1356 assert(i<maxlength_);
1357 rhs()[i]-=tmp[j];
1358 }
1359 }
1360#else
1361 //assign right hand side of current domainindex block
1362 for(size_type j=0; j<n; ++j, ++i) {
1363 assert(i<maxlength_);
1364 rhs()[i]=(*b)[domainIndex][j];
1365
1366 // loop over all Matrix row entries and calculate defect.
1367 typedef typename matrix_type::ConstColIterator col_iterator;
1368
1369 // calculate defect for current row index block
1370 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1371 for(size_type k=0; k<n; ++k) {
1372 rhs()[i]-=(*col)[j][k] * (*x)[col.index()][k];
1373 }
1374 }
1375 }
1376#endif
1377 }
1378
1379 template<class K, int n, class Al, class X, class Y>
1380 void
1381 OverlappingAssignerHelper< DynamicMatrixSubdomainSolver< BCRSMatrix< FieldMatrix<K,n,n>, Al>, X, Y >,false>
1382 ::assignResult(block_type& res)
1383 {
1384 // assign the result of the local solve to the global vector
1385 for(size_type j=0; j<n; ++j, ++i) {
1386 assert(i<maxlength_);
1387 res[j]+=lhs()[i];
1388 }
1389 }
1390
1391#if HAVE_SUPERLU || HAVE_UMFPACK
1392
1393 template<template<class> class S, int n, int m, typename T, typename A>
1394 OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>
1395 ::OverlappingAssignerHelper(std::size_t maxlength,
1396 const BCRSMatrix<FieldMatrix<T,n,m>,A>& mat_,
1397 const range_type& b_,
1398 range_type& x_)
1399 : mat(&mat_),
1400 b(&b_),
1401 x(&x_), i(0), maxlength_(maxlength)
1402 {
1403 rhs_ = new field_type[maxlength];
1404 lhs_ = new field_type[maxlength];
1405
1406 }
1407
1408 template<template<class> class S, int n, int m, typename T, typename A>
1409 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::deallocate()
1410 {
1411 delete[] rhs_;
1412 delete[] lhs_;
1413 }
1414
1415 template<template<class> class S, int n, int m, typename T, typename A>
1416 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::operator()(const size_type& domainIndex)
1417 {
1418 //assign right hand side of current domainindex block
1419 // rhs is an array of doubles!
1420 // rhs[starti] = b[domainindex]
1421 for(size_type j=0; j<n; ++j, ++i) {
1422 assert(i<maxlength_);
1423 rhs_[i]=(*b)[domainIndex][j];
1424 }
1425
1426
1427 // loop over all Matrix row entries and calculate defect.
1428 typedef typename matrix_type::ConstColIterator col_iterator;
1429
1430 // calculate defect for current row index block
1431 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1432 block_type tmp;
1433 (*col).mv((*x)[col.index()], tmp);
1434 i-=n;
1435 for(size_type j=0; j<n; ++j, ++i) {
1436 assert(i<maxlength_);
1437 rhs_[i]-=tmp[j];
1438 }
1439
1440 }
1441
1442 }
1443
1444 template<template<class> class S, int n, int m, typename T, typename A>
1445 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::relaxResult(field_type relax)
1446 {
1447 for(size_type j=i+n; i<j; ++i) {
1448 assert(i<maxlength_);
1449 lhs_[i]*=relax;
1450 }
1451 i-=n;
1452 }
1453
1454 template<template<class> class S, int n, int m, typename T, typename A>
1455 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::assignResult(block_type& res)
1456 {
1457 // assign the result of the local solve to the global vector
1458 for(size_type j=0; j<n; ++j, ++i) {
1459 assert(i<maxlength_);
1460 res[j]+=lhs_[i];
1461 }
1462 }
1463
1464 template<template<class> class S, int n, int m, typename T, typename A>
1465 void OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::resetIndexForNextDomain()
1466 {
1467 i=0;
1468 }
1469
1470 template<template<class> class S, int n, int m, typename T, typename A>
1471 typename OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::field_type*
1472 OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::lhs()
1473 {
1474 return lhs_;
1475 }
1476
1477 template<template<class> class S, int n, int m, typename T, typename A>
1478 typename OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::field_type*
1479 OverlappingAssignerHelper<S<BCRSMatrix<FieldMatrix<T,n,m>,A> >,true>::rhs()
1480 {
1481 return rhs_;
1482 }
1483
1484#endif
1485
1486 template<class M, class X, class Y>
1487 OverlappingAssignerILUBase<M,X,Y>::OverlappingAssignerILUBase(std::size_t maxlength,
1488 const M& mat_,
1489 const Y& b_,
1490 X& x_)
1491 : mat(&mat_),
1492 b(&b_),
1493 x(&x_), i(0)
1494 {
1495 rhs_= new Y(maxlength);
1496 lhs_ = new X(maxlength);
1497 }
1498
1499 template<class M, class X, class Y>
1500 void OverlappingAssignerILUBase<M,X,Y>::deallocate()
1501 {
1502 delete rhs_;
1503 delete lhs_;
1504 }
1505
1506 template<class M, class X, class Y>
1507 void OverlappingAssignerILUBase<M,X,Y>::operator()(const size_type& domainIndex)
1508 {
1509 (*rhs_)[i]=(*b)[domainIndex];
1510
1511 // loop over all Matrix row entries and calculate defect.
1512 typedef typename matrix_type::ConstColIterator col_iterator;
1513
1514 // calculate defect for current row index block
1515 for(col_iterator col=(*mat)[domainIndex].begin(); col!=(*mat)[domainIndex].end(); ++col) {
1516 (*col).mmv((*x)[col.index()], (*rhs_)[i]);
1517 }
1518 // Goto next local index
1519 ++i;
1520 }
1521
1522 template<class M, class X, class Y>
1523 void OverlappingAssignerILUBase<M,X,Y>::relaxResult(field_type relax)
1524 {
1525 (*lhs_)[i]*=relax;
1526 }
1527
1528 template<class M, class X, class Y>
1529 void OverlappingAssignerILUBase<M,X,Y>::assignResult(block_type& res)
1530 {
1531 res+=(*lhs_)[i++];
1532 }
1533
1534 template<class M, class X, class Y>
1535 X& OverlappingAssignerILUBase<M,X,Y>::lhs()
1536 {
1537 return *lhs_;
1538 }
1539
1540 template<class M, class X, class Y>
1541 Y& OverlappingAssignerILUBase<M,X,Y>::rhs()
1542 {
1543 return *rhs_;
1544 }
1545
1546 template<class M, class X, class Y>
1547 void OverlappingAssignerILUBase<M,X,Y>::resetIndexForNextDomain()
1548 {
1549 i=0;
1550 }
1551
1552 template<typename S, typename T, typename A, int n>
1553 AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::AdditiveAdder(BlockVector<FieldVector<T,n>,A>& v_,
1555 OverlappingAssigner<S>& assigner_,
1556 const T& relax_)
1557 : v(&v_), x(&x_), assigner(&assigner_), relax(relax_)
1558 {}
1559
1560 template<typename S, typename T, typename A, int n>
1561 void AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::operator()(const size_type& domainIndex)
1562 {
1563 // add the result of the local solve to the current update
1564 assigner->assignResult((*v)[domainIndex]);
1565 }
1566
1567
1568 template<typename S, typename T, typename A, int n>
1569 void AdditiveAdder<S,BlockVector<FieldVector<T,n>,A> >::axpy()
1570 {
1571 // relax the update and add it to the current guess.
1572 x->axpy(relax,*v);
1573 }
1574
1575
1576 template<typename S, typename T, typename A, int n>
1577 MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >
1578 ::MultiplicativeAdder(BlockVector<FieldVector<T,n>,A>& v_,
1579 BlockVector<FieldVector<T,n>,A>& x_,
1580 OverlappingAssigner<S>& assigner_, const T& relax_)
1581 : x(&x_), assigner(&assigner_), relax(relax_)
1582 {
1584 }
1585
1586
1587 template<typename S,typename T, typename A, int n>
1588 void MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >::operator()(const size_type& domainIndex)
1589 {
1590 // add the result of the local solve to the current guess
1591 assigner->relaxResult(relax);
1592 assigner->assignResult((*x)[domainIndex]);
1593 }
1594
1595
1596 template<typename S,typename T, typename A, int n>
1597 void MultiplicativeAdder<S,BlockVector<FieldVector<T,n>,A> >::axpy()
1598 {
1599 // nothing to do, as the corrections already relaxed and added in operator()
1600 }
1601
1602
1604}
1605
1606#endif
Implementation of the BCRSMatrix class.
This file implements a vector space as a tensor product of a given vector space. The number of compon...
A sparse block matrix with compressed row storage.
Definition: bcrsmatrix.hh:413
row_type::ConstIterator ConstColIterator
Const iterator to the entries of a row.
Definition: bcrsmatrix.hh:690
A vector of blocks with memory management.
Definition: bvector.hh:253
Iterator begin()
begin iterator
Definition: densevector.hh:307
size_type size() const
size method
Definition: densevector.hh:296
Iterator end()
end iterator
Definition: densevector.hh:313
Iterator find(size_type i)
return iterator to given element or end()
Definition: densevector.hh:333
Exact subdomain solver using Dune::DynamicMatrix<T>::solve.
Definition: overlappingschwarz.hh:136
size_type capacity() const
Number of elements for which memory has been allocated.
Definition: dynvector.hh:126
A dense n x m matrix.
Definition: fmatrix.hh:67
Index Set Interface base class.
Definition: indexidset.hh:76
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:79
Initializer for SuperLU Matrices representing the subdomains.
Definition: overlappingschwarz.hh:43
D subdomain_vector
The vector type containing the subdomain to row index mapping.
Definition: overlappingschwarz.hh:46
Base class for matrix free definition of preconditioners.
Definition: preconditioner.hh:26
A constant iterator for the SLList.
Definition: sllist.hh:369
A single linked list.
Definition: sllist.hh:42
Sequential overlapping Schwarz preconditioner.
Definition: overlappingschwarz.hh:753
X::field_type field_type
The field type of the preconditioner.
Definition: overlappingschwarz.hh:781
void apply(X &v, const X &d)
Apply one step of the preconditioner to the system A(v)=d.
M matrix_type
The type of the matrix to precondition.
Definition: overlappingschwarz.hh:758
TM Mode
The mode (additive or multiplicative) of the Schwarz method.
Definition: overlappingschwarz.hh:776
X range_type
The range type of the preconditioner.
Definition: overlappingschwarz.hh:768
X domain_type
The domain type of the preconditioner.
Definition: overlappingschwarz.hh:763
TD slu
The type for the subdomain solver in use.
Definition: overlappingschwarz.hh:804
virtual void pre(X &x, X &b)
Prepare the preconditioner.
Definition: overlappingschwarz.hh:849
virtual void post(X &x)
Postprocess the preconditioner.
Definition: overlappingschwarz.hh:867
std::vector< subdomain_type, typename TA::template rebind< subdomain_type >::other > subdomain_vector
The vector type containing the subdomain to row index mapping.
Definition: overlappingschwarz.hh:795
TA allocator
The allocator to use.
Definition: overlappingschwarz.hh:787
SLList< size_type, typename TA::template rebind< size_type >::other > subdomain_list
The type for the row to subdomain mapping.
Definition: overlappingschwarz.hh:798
std::vector< subdomain_list, typename TA::template rebind< subdomain_list >::other > rowtodomain_vector
The vector type containing the row index to subdomain mapping.
Definition: overlappingschwarz.hh:801
std::vector< slu, typename TA::template rebind< slu >::other > slu_vector
The vector type containing subdomain solvers.
Definition: overlappingschwarz.hh:807
matrix_type::size_type size_type
The return type of the size method.
Definition: overlappingschwarz.hh:784
@ category
The category the precondtioner is part of.
Definition: overlappingschwarz.hh:811
std::set< size_type, std::less< size_type >, typename TA::template rebind< size_type >::other > subdomain_type
The type for the subdomain to row index mapping.
Definition: overlappingschwarz.hh:792
RealIterator< const B > const_iterator
iterator class for sequential access
Definition: basearray.hh:195
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: bvector.hh:217
This file implements a dense matrix with dynamic numbers of rows and columns.
virtual void apply(X &v, const X &d)
Apply the precondtioner.
Definition: overlappingschwarz.hh:1231
SeqOverlappingSchwarz(const matrix_type &mat, const subdomain_vector &subDomains, field_type relaxationFactor=1, bool onTheFly_=true)
Construct the overlapping Schwarz method.
Definition: overlappingschwarz.hh:1058
SeqOverlappingSchwarz(const matrix_type &mat, const rowtodomain_vector &rowToDomain, field_type relaxationFactor=1, bool onTheFly_=true)
Definition: overlappingschwarz.hh:1011
DVVerbType dvverb(std::cout)
stream for very verbose output.
Definition: stdstreams.hh:93
Various local subdomain solvers based on ILU for SeqOverlappingSchwarz.
Dune namespace.
Definition: alignment.hh:10
Define general preconditioner interface.
Implements a singly linked list together with the necessary iterators.
Templates characterizing the type of a solver.
template meta program for choosing how to add the correction.
Definition: overlappingschwarz.hh:569
Tag that the tells the schwarz method to be additive.
Definition: overlappingschwarz.hh:116
Helper template meta program for application of overlapping schwarz.
Definition: overlappingschwarz.hh:602
Tag that tells the Schwarz method to be multiplicative.
Definition: overlappingschwarz.hh:122
Helper template meta program for application of overlapping schwarz.
Definition: overlappingschwarz.hh:666
Definition: overlappingschwarz.hh:1105
@ sequential
Category for sequential solvers.
Definition: solvercategory.hh:21
Tag that tells the Schwarz method to be multiplicative and symmetric.
Definition: overlappingschwarz.hh:129
Classes for using SuperLU with ISTL matrices.
Classes for using UMFPack with ISTL matrices.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)