Dune Core Modules (2.8.0)

ilusubdomainsolver.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_ILUSUBDOMAINSOLVER_HH
4#define DUNE_ISTL_ILUSUBDOMAINSOLVER_HH
5
6#include <map>
9#include "matrix.hh"
10#include <cmath>
11#include <cstdlib>
12
13namespace Dune {
14
33 template<class M, class X, class Y>
35 public:
37 typedef typename std::remove_const<M>::type matrix_type;
39 typedef X domain_type;
41 typedef Y range_type;
42
49 virtual void apply (X& v, const Y& d) =0;
50
51 virtual ~ILUSubdomainSolver()
52 {}
53
54 protected:
60 template<class S>
61 std::size_t copyToLocalMatrix(const M& A, S& rowset);
62
64 // for ILUN
66 };
67
74 template<class M, class X, class Y>
76 : public ILUSubdomainSolver<M,X,Y>{
77 public:
79 typedef typename std::remove_const<M>::type matrix_type;
80 typedef typename std::remove_const<M>::type rilu_type;
82 typedef X domain_type;
84 typedef Y range_type;
85
86
91 void apply (X& v, const Y& d)
92 {
93 ILU::blockILUBacksolve(this->ILU,v,d);
94 }
102 template<class S>
103 void setSubMatrix(const M& A, S& rowset);
104
105 };
106
107 template<class M, class X, class Y>
108 class ILUNSubdomainSolver
109 : public ILUSubdomainSolver<M,X,Y>{
110 public:
112 typedef typename std::remove_const<M>::type matrix_type;
113 typedef typename std::remove_const<M>::type rilu_type;
115 typedef X domain_type;
117 typedef Y range_type;
118
123 void apply (X& v, const Y& d)
124 {
125 ILU::blockILUBacksolve(RILU,v,d);
126 }
127
135 template<class S>
136 void setSubMatrix(const M& A, S& rowset);
137
138 private:
142 rilu_type RILU;
143 };
144
145
146
147 template<class M, class X, class Y>
148 template<class S>
149 std::size_t ILUSubdomainSolver<M,X,Y>::copyToLocalMatrix(const M& A, S& rowSet)
150 {
151 // Calculate consecutive indices for local problem
152 // while perserving the ordering
153 typedef typename M::size_type size_type;
154 typedef std::map<typename S::value_type,size_type> IndexMap;
155 typedef typename IndexMap::iterator IMIter;
156 IndexMap indexMap;
157 IMIter guess = indexMap.begin();
158 size_type localIndex=0;
159
160 typedef typename S::const_iterator SIter;
161 for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
162 rowIdx!= rowEnd; ++rowIdx, ++localIndex)
163 guess = indexMap.insert(guess,
164 std::make_pair(*rowIdx,localIndex));
165
166
167 // Build Matrix for local subproblem
168 ILU.setSize(rowSet.size(),rowSet.size());
169 ILU.setBuildMode(matrix_type::row_wise);
170
171 // Create sparsity pattern
172 typedef typename matrix_type::CreateIterator CIter;
173 CIter rowCreator = ILU.createbegin();
174 std::size_t offset=0;
175 for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
176 rowIdx!= rowEnd; ++rowIdx, ++rowCreator) {
177 // See which row entries are in our subset and add them to
178 // the sparsity pattern
179 guess = indexMap.begin();
180
181 for(typename matrix_type::ConstColIterator col=A[*rowIdx].begin(),
182 endcol=A[*rowIdx].end(); col != endcol; ++col) {
183 // search for the entry in the row set
184 guess = indexMap.find(col.index());
185 if(guess!=indexMap.end()) {
186 // add local index to row
187 rowCreator.insert(guess->second);
188 offset=std::max(offset,(std::size_t)std::abs((int)(guess->second-rowCreator.index())));
189 }
190 }
191
192 }
193
194 // Insert the matrix values for the local problem
195 typename matrix_type::iterator iluRow=ILU.begin();
196
197 for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
198 rowIdx!= rowEnd; ++rowIdx, ++iluRow) {
199 // See which row entries are in our subset and add them to
200 // the sparsity pattern
201 typename matrix_type::ColIterator localCol=iluRow->begin();
202 for(typename matrix_type::ConstColIterator col=A[*rowIdx].begin(),
203 endcol=A[*rowIdx].end(); col != endcol; ++col) {
204 // search for the entry in the row set
205 guess = indexMap.find(col.index());
206 if(guess!=indexMap.end()) {
207 // set local value
208 (*localCol)=(*col);
209 ++localCol;
210 }
211 }
212 }
213 return offset;
214 }
215
216
217 template<class M, class X, class Y>
218 template<class S>
220 {
221 this->copyToLocalMatrix(A,rowSet);
222 ILU::blockILU0Decomposition(this->ILU);
223 }
224
225 template<class M, class X, class Y>
226 template<class S>
227 void ILUNSubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
228 {
229 std::size_t offset=copyToLocalMatrix(A,rowSet);
230 RILU.setSize(rowSet.size(),rowSet.size(), (1+2*offset)*rowSet.size());
231 RILU.setBuildMode(matrix_type::row_wise);
232 ILU::blockILUDecomposition(this->ILU, (offset+1)/2, RILU);
233 }
234
236} // end name space DUNE
237
238
239#endif
Exact subdomain solver using ILU(p) with appropriate p.
Definition: ilusubdomainsolver.hh:76
X domain_type
The domain type of the preconditioner.
Definition: ilusubdomainsolver.hh:82
Y range_type
The range type of the preconditioner.
Definition: ilusubdomainsolver.hh:84
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: ilusubdomainsolver.hh:79
void apply(X &v, const Y &d)
Apply the subdomain solver.
Definition: ilusubdomainsolver.hh:91
base class encapsulating common algorithms of ILU0SubdomainSolver and ILUNSubdomainSolver.
Definition: ilusubdomainsolver.hh:34
matrix_type ILU
The ILU0 decomposition of the matrix, or the local matrix.
Definition: ilusubdomainsolver.hh:65
X domain_type
The domain type of the preconditioner.
Definition: ilusubdomainsolver.hh:39
Y range_type
The range type of the preconditioner.
Definition: ilusubdomainsolver.hh:41
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: ilusubdomainsolver.hh:37
virtual void apply(X &v, const Y &d)=0
Apply the subdomain solver.
Traits for type conversions and type information.
std::size_t copyToLocalMatrix(const M &A, S &rowset)
Copy the local part of the global matrix to ILU.
Definition: ilusubdomainsolver.hh:149
void setSubMatrix(const M &A, S &rowset)
Set the data of the local problem.
Definition: ilusubdomainsolver.hh:219
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
A dynamic dense block matrix class.
Dune namespace.
Definition: alignedallocator.hh:11
Define general preconditioner interface.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)