3 #ifndef DUNE_ISTL_ILU_HH
4 #define DUNE_ISTL_ILU_HH
15 #include <dune/common/fmatrix.hh>
39 typedef typename M::RowIterator rowiterator;
40 typedef typename M::ColIterator coliterator;
41 typedef typename M::block_type block;
44 rowiterator endi=A.
end();
45 for (rowiterator i=A.
begin(); i!=endi; ++i)
48 coliterator endij=(*i).end();
52 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
55 coliterator jj = A[ij.index()].find(ij.index());
58 (*ij).rightmultiply(*jj);
61 coliterator endjk=A[ij.index()].
end();
62 coliterator jk=jj; ++jk;
63 coliterator ik=ij; ++ik;
64 while (ik!=endij && jk!=endjk)
65 if (ik.index()==jk.index())
74 if (ik.index()<jk.index())
82 if (ij.index()!=i.index())
83 DUNE_THROW(
ISTLError,
"diagonal entry missing");
87 catch (Dune::FMatrixError & e) {
89 << i.index() <<
"][" << ij.index() <<
"]" << e.what();
90 th__ex.
r=i.index(); th__ex.c=ij.index(););
96 template<
class M,
class X,
class Y>
100 typedef typename M::ConstRowIterator rowiterator;
101 typedef typename M::ConstColIterator coliterator;
102 typedef typename Y::block_type dblock;
103 typedef typename X::block_type vblock;
106 rowiterator endi=A.
end();
107 for (rowiterator i=A.
begin(); i!=endi; ++i)
109 dblock rhs(d[i.index()]);
110 for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
111 (*j).mmv(v[j.index()],rhs);
117 for (rowiterator i=A.
beforeEnd(); i!=rendi; --i)
119 vblock rhs(v[i.index()]);
121 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
122 (*j).mmv(v[j.index()],rhs);
124 (*j).umv(rhs,v[i.index()]);
137 template<
class K,
int n,
int m>
160 typedef typename M::ColIterator coliterator;
161 typedef typename M::ConstRowIterator crowiterator;
162 typedef typename M::ConstColIterator ccoliterator;
163 typedef typename M::CreateIterator createiterator;
164 typedef typename M::field_type K;
165 typedef std::map<size_t, int> map;
166 typedef typename map::iterator mapiterator;
169 crowiterator endi=A.
end();
170 createiterator ci=ILU.createbegin();
171 for (crowiterator i=A.
begin(); i!=endi; ++i)
177 for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
178 rowpattern[j.index()] = 0;
181 for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
188 coliterator endk = ILU[(*ik).first].end();
189 coliterator kj = ILU[(*ik).first].find((*ik).first);
190 for (++kj; kj!=endk; ++kj)
198 mapiterator ij = rowpattern.find(kj.index());
199 if (ij==rowpattern.end())
204 rowpattern[kj.index()] = generation+1;
212 for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
213 ci.insert((*ik).first);
217 coliterator endILUij = ILU[i.index()].end();;
218 for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
225 for (crowiterator i=A.
begin(); i!=endi; ++i)
228 coliterator endILUij = ILU[i.index()].end();;
229 for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
231 ccoliterator Aij = (*i).begin();
232 ccoliterator endAij = (*i).end();
233 ILUij = ILU[i.index()].begin();
234 while (Aij!=endAij && ILUij!=endILUij)
236 if (Aij.index()==ILUij.index())
243 if (Aij.index()<ILUij.index())
285 if(
values_.capacity() < needed )
289 cols_.reserve( estimate );
296 cols_.push_back( index );
306 template<
class M,
class CRS,
class InvVector>
309 typedef typename M :: size_type size_type;
316 const size_t memEstimate = (A.nonzeroes() - A.N())/2;
318 assert( A.nonzeroes() != 0 );
322 const auto endi = A.end();
324 size_type colcount = 0;
325 lower.
rows_[ 0 ] = colcount;
326 for (
auto i=A.begin(); i!=endi; ++i, ++row)
328 const size_type iIndex = i.index();
331 for (
auto j=(*i).begin(); j.index() < iIndex; ++j )
336 lower.
rows_[ iIndex+1 ] = colcount;
339 const auto rendi = A.beforeBegin();
342 upper.
rows_[ 0 ] = colcount ;
346 for (
auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
348 const auto endij=(*i).beforeBegin();
350 const size_type iIndex = i.index();
353 for (
auto j=(*i).beforeEnd(); j != endij; --j )
355 const size_type jIndex = j.index();
356 if( j.index() == iIndex )
361 else if ( j.index() >= i.index() )
367 upper.
rows_[ row+1 ] = colcount;
373 template<
class CRS,
class mblock,
class X,
class Y>
376 const std::vector< mblock >& inv,
380 typedef typename Y :: block_type dblock;
381 typedef typename X :: block_type vblock;
382 typedef typename X :: size_type size_type ;
384 const size_type iEnd = lower.
rows();
385 const size_type lastRow = iEnd - 1;
386 if( iEnd != upper.
rows() )
388 DUNE_THROW(
ISTLError,
"ILU::bilu_backsolve: lower and upper rows must be the same");
392 for( size_type i=0; i<iEnd; ++ i )
394 dblock rhs( d[ i ] );
395 const size_type rowI = lower.
rows_[ i ];
396 const size_type rowINext = lower.
rows_[ i+1 ];
398 for( size_type
col = rowI;
col < rowINext; ++
col )
407 for( size_type i=0; i<iEnd; ++ i )
409 vblock& vBlock = v[ lastRow - i ];
410 vblock rhs ( vBlock );
411 const size_type rowI = upper.
rows_[ i ];
412 const size_type rowINext = upper.
rows_[ i+1 ];
414 for( size_type
col = rowI;
col < rowINext; ++
col )
420 inv[ i ].mv( rhs, vBlock);