dune-istl  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>
7 #include <dune/common/typetraits.hh>
9 #include "matrix.hh"
10 #include <cmath>
11 #include <cstdlib>
12 
13 namespace 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 
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>
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>
219  void ILU0SubdomainSolver<M,X,Y>::setSubMatrix(const M& A, S& rowSet)
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
A dynamic dense block matrix class.
Define general preconditioner interface.
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:227
void setSubMatrix(const M &A, S &rowset)
Set the data of the local problem.
Definition: ilusubdomainsolver.hh:219
Col col
Definition: matrixmatrix.hh:349
Definition: allocator.hh:9
void blockILUBacksolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition: ilu.hh:93
void blockILU0Decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:32
void blockILUDecomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:166
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
virtual ~ILUSubdomainSolver()
Definition: ilusubdomainsolver.hh:51
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.
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 rilu_type
Definition: ilusubdomainsolver.hh:80
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
Definition: ilusubdomainsolver.hh:109
X domain_type
The domain type of the preconditioner.
Definition: ilusubdomainsolver.hh:115
std::remove_const< M >::type matrix_type
The matrix type the preconditioner is for.
Definition: ilusubdomainsolver.hh:112
std::remove_const< M >::type rilu_type
Definition: ilusubdomainsolver.hh:113
void apply(X &v, const Y &d)
Apply the subdomain solver.
Definition: ilusubdomainsolver.hh:123
Y range_type
The range type of the preconditioner.
Definition: ilusubdomainsolver.hh:117