3 #ifndef DUNE_ISTL_ILU_HH
4 #define DUNE_ISTL_ILU_HH
11 #include <dune/common/fmatrix.hh>
12 #include <dune/common/scalarvectorview.hh>
13 #include <dune/common/scalarmatrixview.hh>
37 typedef typename M::RowIterator rowiterator;
38 typedef typename M::ColIterator coliterator;
39 typedef typename M::block_type block;
42 rowiterator endi=A.
end();
43 for (rowiterator i=A.
begin(); i!=endi; ++i)
46 coliterator endij=(*i).end();
50 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
53 coliterator jj = A[ij.index()].find(ij.index());
56 Impl::asMatrix(*ij).rightmultiply(Impl::asMatrix(*jj));
59 coliterator endjk=A[ij.index()].
end();
60 coliterator jk=jj; ++jk;
61 coliterator ik=ij; ++ik;
62 while (ik!=endij && jk!=endjk)
63 if (ik.index()==jk.index())
66 Impl::asMatrix(B).leftmultiply(Impl::asMatrix(*ij));
72 if (ik.index()<jk.index())
80 if (ij.index()!=i.index())
81 DUNE_THROW(
ISTLError,
"diagonal entry missing");
83 Impl::asMatrix(*ij).invert();
85 catch (Dune::FMatrixError & e) {
87 << i.index() <<
"][" << ij.index() <<
"]" << e.what();
88 th__ex.
r=i.index(); th__ex.c=ij.index(););
94 template<
class M,
class X,
class Y>
98 typedef typename M::ConstRowIterator rowiterator;
99 typedef typename M::ConstColIterator coliterator;
100 typedef typename Y::block_type dblock;
101 typedef typename X::block_type vblock;
104 rowiterator endi=A.
end();
105 for (rowiterator i=A.
begin(); i!=endi; ++i)
113 dblock rhsValue(d[i.index()]);
114 auto&& rhs = Impl::asVector(rhsValue);
115 for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
116 Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
117 Impl::asVector(v[i.index()]) = rhs;
122 for (rowiterator i=A.
beforeEnd(); i!=rendi; --i)
130 vblock rhsValue(v[i.index()]);
131 auto&& rhs = Impl::asVector(rhsValue);
133 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
134 Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
135 auto&& vi = Impl::asVector(v[i.index()]);
136 Impl::asMatrix(*j).mv(rhs,vi);
144 typename M::field_type&
firstmatrixelement (M& A,
typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae =
nullptr)
155 template<
class K,
int n,
int m>
172 typedef typename M::ColIterator coliterator;
173 typedef typename M::ConstRowIterator crowiterator;
174 typedef typename M::ConstColIterator ccoliterator;
175 typedef typename M::CreateIterator createiterator;
176 typedef typename M::field_type K;
177 typedef std::map<size_t, int> map;
178 typedef typename map::iterator mapiterator;
181 crowiterator endi=A.
end();
182 createiterator ci=ILU.createbegin();
183 for (crowiterator i=A.
begin(); i!=endi; ++i)
188 for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
189 rowpattern[j.index()] = 0;
192 for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
196 coliterator endk = ILU[(*ik).first].end();
197 coliterator kj = ILU[(*ik).first].find((*ik).first);
198 for (++kj; kj!=endk; ++kj)
206 mapiterator ij = rowpattern.find(kj.index());
207 if (ij==rowpattern.end())
209 rowpattern[kj.index()] = generation+1;
217 for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
218 ci.insert((*ik).first);
222 coliterator endILUij = ILU[i.index()].end();;
223 for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
224 Simd::lane(0,
firstmatrixelement(*ILUij)) = (Simd::Scalar<K>) rowpattern[ILUij.index()];
228 for (crowiterator i=A.
begin(); i!=endi; ++i)
231 coliterator endILUij = ILU[i.index()].end();;
232 for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
234 ccoliterator Aij = (*i).begin();
235 ccoliterator endAij = (*i).end();
236 ILUij = ILU[i.index()].begin();
237 while (Aij!=endAij && ILUij!=endILUij)
239 if (Aij.index()==ILUij.index())
246 if (Aij.index()<ILUij.index())
260 template <
class B,
class Alloc = std::allocator<B>>
288 if(
values_.capacity() < needed )
292 cols_.reserve( estimate );
299 cols_.push_back( index );
309 template<
class M,
class CRS,
class InvVector>
312 typedef typename M :: size_type size_type;
319 const size_t memEstimate = (A.nonzeroes() - A.N())/2;
321 assert( A.nonzeroes() != 0 );
325 const auto endi = A.end();
327 size_type colcount = 0;
328 lower.
rows_[ 0 ] = colcount;
329 for (
auto i=A.begin(); i!=endi; ++i, ++row)
331 const size_type iIndex = i.index();
334 for (
auto j=(*i).begin(); j.index() < iIndex; ++j )
339 lower.
rows_[ iIndex+1 ] = colcount;
342 const auto rendi = A.beforeBegin();
345 upper.
rows_[ 0 ] = colcount ;
349 for (
auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
351 const auto endij=(*i).beforeBegin();
353 const size_type iIndex = i.index();
356 for (
auto j=(*i).beforeEnd(); j != endij; --j )
358 const size_type jIndex = j.index();
359 if( j.index() == iIndex )
364 else if ( j.index() >= i.index() )
370 upper.
rows_[ row+1 ] = colcount;
376 template<
class CRS,
class InvVector,
class X,
class Y>
379 const InvVector& inv,
383 typedef typename Y :: block_type dblock;
384 typedef typename X :: block_type vblock;
385 typedef typename X :: size_type size_type ;
387 const size_type iEnd = lower.
rows();
388 const size_type lastRow = iEnd - 1;
389 if( iEnd != upper.
rows() )
391 DUNE_THROW(
ISTLError,
"ILU::bilu_backsolve: lower and upper rows must be the same");
395 for( size_type i=0; i<iEnd; ++ i )
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 ];
402 for( size_type
col = rowI;
col < rowINext; ++
col )
403 Impl::asMatrix(lower.
values_[
col ]).mmv( Impl::asVector(v[ lower.
cols_[
col ] ] ), rhs );
405 Impl::asVector(v[ i ]) = rhs;
409 for( size_type i=0; i<iEnd; ++ i )
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 ];
417 for( size_type
col = rowI;
col < rowINext; ++
col )
418 Impl::asMatrix(upper.
values_[
col ]).mmv( Impl::asVector(v[ upper.
cols_[
col ] ]), rhs );
421 Impl::asMatrix(inv[ i ]).mv(rhs, vBlock);
Col col
Definition: matrixmatrix.hh:349
M::field_type & firstmatrixelement(M &A, typename std::enable_if_t<!Dune::IsNumber< M >::value > *sfinae=nullptr)
Definition: ilu.hh:144
void bilu_backsolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition: ilu.hh:95
void bilu_decomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:169
void bilu0_decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:34
Definition: allocator.hh:7
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:310
void bilu_backsolve(const CRS &lower, const CRS &upper, const InvVector &inv, X &v, const Y &d)
LU backsolve with stored inverse in CRS format for lower and upper triangular.
Definition: ilu.hh:377
int c
Definition: ilu.hh:28
int r
Definition: ilu.hh:28
std::vector< size_type > cols_
Definition: ilu.hh:304
size_type nonZeros() const
Definition: ilu.hh:270
void resize(const size_type nRows)
Definition: ilu.hh:276
size_type rows() const
Definition: ilu.hh:268
CRS()
Definition: ilu.hh:266
void reserveAdditional(const size_type nonZeros)
Definition: ilu.hh:285
B block_type
Definition: ilu.hh:263
std::vector< block_type, Alloc > values_
Definition: ilu.hh:303
size_type nRows_
Definition: ilu.hh:305
size_t size_type
Definition: ilu.hh:264
std::vector< size_type > rows_
Definition: ilu.hh:302
void push_back(const block_type &value, const size_type index)
Definition: ilu.hh:296
derive error class from the base class in common
Definition: istlexception.hh:16
RowIterator beforeBegin()
Definition: matrix.hh:630
RowIterator beforeEnd()
Definition: matrix.hh:623
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:616
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:610
Definition: matrixutils.hh:26