dune-istl  2.2.1
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_ILUSUBDOMAIN_HH
4 #define DUNE_ISTL_ILUSUBDOMAIN_HH
5 
6 #include<map>
7 #include<dune/common/typetraits.hh>
8 #include"matrix.hh"
9 #include<cmath>
10 #include<cstdlib>
11 
12 namespace Dune{
13 
32  template<class M, class X, class Y>
34  public:
38  typedef X domain_type;
40  typedef Y range_type;
41 
48  virtual void apply (X& v, const Y& d) =0;
49 
51  {}
52 
53  protected:
59  template<class S>
60  std::size_t copyToLocalMatrix(const M& A, S& rowset);
61 
63  // for ILUN
65  };
66 
73  template<class M, class X, class Y>
75  : public ILUSubdomainSolver<M,X,Y>{
76  public:
81  typedef X domain_type;
83  typedef Y range_type;
84 
85 
90  void apply (X& v, const Y& d)
91  {
92  bilu_backsolve(this->ILU,v,d);
93  }
101  template<class S>
102  void setSubMatrix(const M& A, S& rowset);
103 
104  };
105 
106  template<class M, class X, class Y>
108  : public ILUSubdomainSolver<M,X,Y>{
109  public:
114  typedef X domain_type;
116  typedef Y range_type;
117 
122  void apply (X& v, const Y& d)
123  {
124  bilu_backsolve(RILU,v,d);
125  }
126 
134  template<class S>
135  void setSubMatrix(const M& A, S& rowset);
136 
137 private:
141  rilu_type RILU;
142  };
143 
144 
145 
146  template<class M, class X, class Y>
147  template<class S>
148  std::size_t ILUSubdomainSolver<M,X,Y>::copyToLocalMatrix(const M& A, S& rowSet)
149  {
150  // Calculate consecutive indices for local problem
151  // while perserving the ordering
152  typedef typename M::size_type size_type;
153  typedef std::map<typename S::value_type,size_type> IndexMap;
154  typedef typename IndexMap::iterator IMIter;
155  IndexMap indexMap;
156  IMIter guess = indexMap.begin();
157  size_type localIndex=0;
158 
159  typedef typename S::const_iterator SIter;
160  for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
161  rowIdx!= rowEnd; ++rowIdx, ++localIndex)
162  guess = indexMap.insert(guess,
163  std::make_pair(*rowIdx,localIndex));
164 
165 
166  // Build Matrix for local subproblem
167  ILU.setSize(rowSet.size(),rowSet.size());
168  ILU.setBuildMode(matrix_type::row_wise);
169 
170  // Create sparsity pattern
171  typedef typename matrix_type::CreateIterator CIter;
172  CIter rowCreator = ILU.createbegin();
173  typedef typename S::const_iterator RIter;
174  std::size_t offset=0;
175  for(SIter rowIdx = rowSet.begin(), rowEnd=rowSet.end();
176  rowIdx!= rowEnd; ++rowIdx, ++rowCreator){
177  // See wich 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 wich 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  bilu0_decomposition(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  bilu_decomposition(this->ILU, (offset+1)/2, RILU);
233  }
234 
236 } // end name space DUNE
237 
238 
239 #endif