dune-istl  2.8.0
ilu.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_ILU_HH
4 #define DUNE_ISTL_ILU_HH
5 
6 #include <cmath>
7 #include <complex>
8 #include <map>
9 #include <vector>
10 
11 #include <dune/common/fmatrix.hh>
12 #include <dune/common/deprecated.hh>
13 #include <dune/common/scalarvectorview.hh>
14 #include <dune/common/scalarmatrixview.hh>
15 
16 #include "istlexception.hh"
17 
22 namespace Dune {
23 
28  namespace ILU {
29 
31  template<class M>
33  {
34  // iterator types
35  typedef typename M::RowIterator rowiterator;
36  typedef typename M::ColIterator coliterator;
37  typedef typename M::block_type block;
38 
39  // implement left looking variant with stored inverse
40  rowiterator endi=A.end();
41  for (rowiterator i=A.begin(); i!=endi; ++i)
42  {
43  // coliterator is diagonal after the following loop
44  coliterator endij=(*i).end(); // end of row i
45  coliterator ij;
46 
47  // eliminate entries left of diagonal; store L factor
48  for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
49  {
50  // find A_jj which eliminates A_ij
51  coliterator jj = A[ij.index()].find(ij.index());
52 
53  // compute L_ij = A_jj^-1 * A_ij
54  Impl::asMatrix(*ij).rightmultiply(Impl::asMatrix(*jj));
55 
56  // modify row
57  coliterator endjk=A[ij.index()].end(); // end of row j
58  coliterator jk=jj; ++jk;
59  coliterator ik=ij; ++ik;
60  while (ik!=endij && jk!=endjk)
61  if (ik.index()==jk.index())
62  {
63  block B(*jk);
64  Impl::asMatrix(B).leftmultiply(Impl::asMatrix(*ij));
65  *ik -= B;
66  ++ik; ++jk;
67  }
68  else
69  {
70  if (ik.index()<jk.index())
71  ++ik;
72  else
73  ++jk;
74  }
75  }
76 
77  // invert pivot and store it in A
78  if (ij.index()!=i.index())
79  DUNE_THROW(ISTLError,"diagonal entry missing");
80  try {
81  Impl::asMatrix(*ij).invert(); // compute inverse of diagonal block
82  }
83  catch (Dune::FMatrixError & e) {
84  DUNE_THROW(MatrixBlockError, "ILU failed to invert matrix block A["
85  << i.index() << "][" << ij.index() << "]" << e.what();
86  th__ex.r=i.index(); th__ex.c=ij.index(););
87  }
88  }
89  }
90 
92  template<class M, class X, class Y>
93  void blockILUBacksolve (const M& A, X& v, const Y& d)
94  {
95  // iterator types
96  typedef typename M::ConstRowIterator rowiterator;
97  typedef typename M::ConstColIterator coliterator;
98  typedef typename Y::block_type dblock;
99  typedef typename X::block_type vblock;
100 
101  // lower triangular solve
102  rowiterator endi=A.end();
103  for (rowiterator i=A.begin(); i!=endi; ++i)
104  {
105  // We need to be careful here: Directly using
106  // auto rhs = Impl::asVector(d[ i.index() ]);
107  // is not OK in case this is a proxy. Hence
108  // we first have to copy the value. Notice that
109  // this is still not OK, if the vector type itself returns
110  // proxy references.
111  dblock rhsValue(d[i.index()]);
112  auto&& rhs = Impl::asVector(rhsValue);
113  for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
114  Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
115  Impl::asVector(v[i.index()]) = rhs; // Lii = I
116  }
117 
118  // upper triangular solve
119  rowiterator rendi=A.beforeBegin();
120  for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
121  {
122  // We need to be careful here: Directly using
123  // auto rhs = Impl::asVector(v[ i.index() ]);
124  // is not OK in case this is a proxy. Hence
125  // we first have to copy the value. Notice that
126  // this is still not OK, if the vector type itself returns
127  // proxy references.
128  vblock rhsValue(v[i.index()]);
129  auto&& rhs = Impl::asVector(rhsValue);
130  coliterator j;
131  for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
132  Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
133  auto&& vi = Impl::asVector(v[i.index()]);
134  Impl::asMatrix(*j).mv(rhs,vi); // diagonal stores inverse!
135  }
136  }
137 
138  // recursive function template to access first entry of a matrix
139  template<class M>
140  typename M::field_type& firstMatrixElement (M& A,
141  [[maybe_unused]] typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
142  {
143  return firstMatrixElement(*(A.begin()->begin()));
144  }
145 
146  template<class K>
148  [[maybe_unused]] typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr)
149  {
150  return A;
151  }
152 
153  template<class K, int n, int m>
155  {
156  return A[0][0];
157  }
158 
165  template<class M>
166  void blockILUDecomposition (const M& A, int n, M& ILU)
167  {
168  // iterator types
169  typedef typename M::ColIterator coliterator;
170  typedef typename M::ConstRowIterator crowiterator;
171  typedef typename M::ConstColIterator ccoliterator;
172  typedef typename M::CreateIterator createiterator;
173  typedef typename M::field_type K;
174  typedef std::map<size_t, int> map;
175  typedef typename map::iterator mapiterator;
176 
177  // symbolic factorization phase, store generation number in first matrix element
178  crowiterator endi=A.end();
179  createiterator ci=ILU.createbegin();
180  for (crowiterator i=A.begin(); i!=endi; ++i)
181  {
182  map rowpattern; // maps column index to generation
183 
184  // initialize pattern with row of A
185  for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
186  rowpattern[j.index()] = 0;
187 
188  // eliminate entries in row which are to the left of the diagonal
189  for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
190  {
191  if ((*ik).second<n)
192  {
193  coliterator endk = ILU[(*ik).first].end(); // end of row k
194  coliterator kj = ILU[(*ik).first].find((*ik).first); // diagonal in k
195  for (++kj; kj!=endk; ++kj) // row k eliminates in row i
196  {
197  // we misuse the storage to store an int. If the field_type is std::complex, we have to access the real/abs part
198  // starting from C++11, we can use std::abs to always return a real value, even if it is double/float
199  using std::abs;
200  int generation = (int) Simd::lane(0, abs( firstMatrixElement(*kj) ));
201  if (generation<n)
202  {
203  mapiterator ij = rowpattern.find(kj.index());
204  if (ij==rowpattern.end())
205  {
206  rowpattern[kj.index()] = generation+1;
207  }
208  }
209  }
210  }
211  }
212 
213  // create row
214  for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
215  ci.insert((*ik).first);
216  ++ci; // now row i exist
217 
218  // write generation index into entries
219  coliterator endILUij = ILU[i.index()].end();;
220  for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
221  Simd::lane(0, firstMatrixElement(*ILUij)) = (Simd::Scalar<K>) rowpattern[ILUij.index()];
222  }
223 
224  // copy entries of A
225  for (crowiterator i=A.begin(); i!=endi; ++i)
226  {
227  coliterator ILUij;
228  coliterator endILUij = ILU[i.index()].end();;
229  for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
230  (*ILUij) = 0; // clear row
231  ccoliterator Aij = (*i).begin();
232  ccoliterator endAij = (*i).end();
233  ILUij = ILU[i.index()].begin();
234  while (Aij!=endAij && ILUij!=endILUij)
235  {
236  if (Aij.index()==ILUij.index())
237  {
238  *ILUij = *Aij;
239  ++Aij; ++ILUij;
240  }
241  else
242  {
243  if (Aij.index()<ILUij.index())
244  ++Aij;
245  else
246  ++ILUij;
247  }
248  }
249  }
250 
251  // call decomposition on pattern
253  }
254 
256  template <class B, class Alloc = std::allocator<B>>
257  struct CRS
258  {
259  typedef B block_type;
260  typedef size_t size_type;
261 
262  CRS() : nRows_( 0 ) {}
263 
264  size_type rows() const { return nRows_; }
265 
267  {
268  assert( rows_[ rows() ] != size_type(-1) );
269  return rows_[ rows() ];
270  }
271 
272  void resize( const size_type nRows )
273  {
274  if( nRows_ != nRows )
275  {
276  nRows_ = nRows ;
277  rows_.resize( nRows_+1, size_type(-1) );
278  }
279  }
280 
282  {
283  const size_type needed = values_.size() + nonZeros ;
284  if( values_.capacity() < needed )
285  {
286  const size_type estimate = needed * 1.1;
287  values_.reserve( estimate );
288  cols_.reserve( estimate );
289  }
290  }
291 
292  void push_back( const block_type& value, const size_type index )
293  {
294  values_.push_back( value );
295  cols_.push_back( index );
296  }
297 
298  std::vector< size_type > rows_;
299  std::vector< block_type, Alloc> values_;
300  std::vector< size_type > cols_;
302  };
303 
305  template<class M, class CRS, class InvVector>
306  void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
307  {
308  typedef typename M :: size_type size_type;
309 
310  lower.resize( A.N() );
311  upper.resize( A.N() );
312  inv.resize( A.N() );
313 
314  // lower and upper triangular should store half of non zeros minus diagonal
315  const size_t memEstimate = (A.nonzeroes() - A.N())/2;
316 
317  assert( A.nonzeroes() != 0 );
318  lower.reserveAdditional( memEstimate );
319  upper.reserveAdditional( memEstimate );
320 
321  const auto endi = A.end();
322  size_type row = 0;
323  size_type colcount = 0;
324  lower.rows_[ 0 ] = colcount;
325  for (auto i=A.begin(); i!=endi; ++i, ++row)
326  {
327  const size_type iIndex = i.index();
328 
329  // store entries left of diagonal
330  for (auto j=(*i).begin(); j.index() < iIndex; ++j )
331  {
332  lower.push_back( (*j), j.index() );
333  ++colcount;
334  }
335  lower.rows_[ iIndex+1 ] = colcount;
336  }
337 
338  const auto rendi = A.beforeBegin();
339  row = 0;
340  colcount = 0;
341  upper.rows_[ 0 ] = colcount ;
342 
343  // NOTE: upper and inv store entries in reverse row and col order,
344  // reverse here relative to ILU
345  for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
346  {
347  const auto endij=(*i).beforeBegin(); // end of row i
348 
349  const size_type iIndex = i.index();
350 
351  // store in reverse row order for faster access during backsolve
352  for (auto j=(*i).beforeEnd(); j != endij; --j )
353  {
354  const size_type jIndex = j.index();
355  if( j.index() == iIndex )
356  {
357  inv[ row ] = (*j);
358  break; // assuming consecutive ordering of A
359  }
360  else if ( j.index() >= i.index() )
361  {
362  upper.push_back( (*j), jIndex );
363  ++colcount ;
364  }
365  }
366  upper.rows_[ row+1 ] = colcount;
367  }
368  } // end convertToCRS
369 
370  template<class CRS, class InvVector, class X, class Y>
371  DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use blockILUBacksolve.")
372  void bilu_backsolve (const CRS& lower, const CRS& upper, const InvVector& inv, X& v, const Y& d)
373  { blockILUBacksolve(lower, upper, inv, v, d); }
374 
376  template<class CRS, class InvVector, class X, class Y>
377  void blockILUBacksolve (const CRS& lower,
378  const CRS& upper,
379  const InvVector& inv,
380  X& v, const Y& d)
381  {
382  // iterator types
383  typedef typename Y :: block_type dblock;
384  typedef typename X :: block_type vblock;
385  typedef typename X :: size_type size_type ;
386 
387  const size_type iEnd = lower.rows();
388  const size_type lastRow = iEnd - 1;
389  if( iEnd != upper.rows() )
390  {
391  DUNE_THROW(ISTLError,"ILU::blockILUBacksolve: lower and upper rows must be the same");
392  }
393 
394  // lower triangular solve
395  for( size_type i=0; i<iEnd; ++ i )
396  {
397  dblock rhsValue( d[ i ] );
398  auto&& rhs = Impl::asVector(rhsValue);
399  const size_type rowI = lower.rows_[ i ];
400  const size_type rowINext = lower.rows_[ i+1 ];
401 
402  for( size_type col = rowI; col < rowINext; ++ col )
403  Impl::asMatrix(lower.values_[ col ]).mmv( Impl::asVector(v[ lower.cols_[ col ] ] ), rhs );
404 
405  Impl::asVector(v[ i ]) = rhs; // Lii = I
406  }
407 
408  // upper triangular solve
409  for( size_type i=0; i<iEnd; ++ i )
410  {
411  auto&& vBlock = Impl::asVector(v[ lastRow - i ]);
412  vblock rhsValue ( v[ lastRow - i ] );
413  auto&& rhs = Impl::asVector(rhsValue);
414  const size_type rowI = upper.rows_[ i ];
415  const size_type rowINext = upper.rows_[ i+1 ];
416 
417  for( size_type col = rowI; col < rowINext; ++ col )
418  Impl::asMatrix(upper.values_[ col ]).mmv( Impl::asVector(v[ upper.cols_[ col ] ]), rhs );
419 
420  // apply inverse and store result
421  Impl::asMatrix(inv[ i ]).mv(rhs, vBlock);
422  }
423  }
424 
425  } // end namespace ILU
426 
429  template<class M>
430  DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use ILU::blockILU0Decomposition.")
432 
433  template<class M, class X, class Y>
434  DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use ILU::blockILUBacksolve.")
435  void bilu_backsolve (const M& A, X& v, const Y& d) { ILU::blockILUBacksolve(A, v, d); }
436 
437  template<class M>
438  DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use ILU::firstMatrixElement.")
439  decltype(auto) firstmatrixelement (M& A) { return ILU::firstMatrixElement(A); }
440 
441  template<class M>
442  DUNE_DEPRECATED_MSG("Will be removed after release 2.8. Use ILU::blockILUDecomposition.")
443  void bilu_decomposition (const M& A, int n, M& ilu) { return ILU::blockILUDecomposition(A, n, ilu); }
444 
445 } // end namespace
446 
447 #endif
Col col
Definition: matrixmatrix.hh:349
Definition: allocator.hh:9
decltype(auto) firstmatrixelement(M &A)
Definition: ilu.hh:439
void bilu_backsolve(const M &A, X &v, const Y &d)
Definition: ilu.hh:435
void bilu_decomposition(const M &A, int n, M &ilu)
Definition: ilu.hh:443
void bilu0_decomposition(M &A)
Definition: ilu.hh:431
void convertToCRS(const M &A, CRS &lower, CRS &upper, InvVector &inv)
convert ILU decomposition into CRS format for lower and upper triangular and inverse.
Definition: ilu.hh:306
void blockILUBacksolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition: ilu.hh:93
M::field_type & firstMatrixElement(M &A, [[maybe_unused]] typename std::enable_if_t<!Dune::IsNumber< M >::value > *sfinae=nullptr)
Definition: ilu.hh:140
void blockILU0Decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:32
void bilu_backsolve(const CRS &lower, const CRS &upper, const InvVector &inv, X &v, const Y &d)
Definition: ilu.hh:372
void blockILUDecomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:166
a simple compressed row storage matrix class
Definition: ilu.hh:258
std::vector< size_type > cols_
Definition: ilu.hh:300
size_type nonZeros() const
Definition: ilu.hh:266
void resize(const size_type nRows)
Definition: ilu.hh:272
size_type rows() const
Definition: ilu.hh:264
CRS()
Definition: ilu.hh:262
void reserveAdditional(const size_type nonZeros)
Definition: ilu.hh:281
B block_type
Definition: ilu.hh:259
std::vector< block_type, Alloc > values_
Definition: ilu.hh:299
size_type nRows_
Definition: ilu.hh:301
size_t size_type
Definition: ilu.hh:260
std::vector< size_type > rows_
Definition: ilu.hh:298
void push_back(const block_type &value, const size_type index)
Definition: ilu.hh:292
derive error class from the base class in common
Definition: istlexception.hh:17
Error when performing an operation on a matrix block.
Definition: istlexception.hh:59
int r
Definition: istlexception.hh:61
Definition: matrixutils.hh:25