Dune Core Modules (2.9.0)

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