5#ifndef DUNE_ISTL_ILU_HH
6#define DUNE_ISTL_ILU_HH
19#include "istlexception.hh"
38 typedef typename M::RowIterator rowiterator;
39 typedef typename M::ColIterator coliterator;
40 typedef typename M::block_type block;
43 rowiterator endi=A.end();
44 for (rowiterator i=A.begin(); i!=endi; ++i)
47 coliterator endij=(*i).end();
51 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
54 coliterator jj = A[ij.index()].find(ij.index());
57 Impl::asMatrix(*ij).rightmultiply(Impl::asMatrix(*jj));
60 coliterator endjk=A[ij.index()].end();
61 coliterator jk=jj; ++jk;
62 coliterator ik=ij; ++ik;
63 while (ik!=endij && jk!=endjk)
64 if (ik.index()==jk.index())
67 Impl::asMatrix(B).leftmultiply(Impl::asMatrix(*ij));
73 if (ik.index()<jk.index())
81 if (ij.index()!=i.index())
84 Impl::asMatrix(*ij).invert();
87 std::ostringstream sstream;
90 <<
"ILU failed to invert matrix block A["
91 << i.index() <<
"][" << ij.index() <<
"]" << e.what();
102 template<
class M,
class X,
class Y>
106 typedef typename M::ConstRowIterator rowiterator;
107 typedef typename M::ConstColIterator coliterator;
108 typedef typename Y::block_type dblock;
109 typedef typename X::block_type vblock;
112 rowiterator endi=A.end();
113 for (rowiterator i=A.begin(); i!=endi; ++i)
121 dblock rhsValue(d[i.index()]);
122 auto&& rhs = Impl::asVector(rhsValue);
123 for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
124 Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
125 Impl::asVector(v[i.index()]) = rhs;
129 rowiterator rendi=A.beforeBegin();
130 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
138 vblock rhsValue(v[i.index()]);
139 auto&& rhs = Impl::asVector(rhsValue);
141 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
142 Impl::asMatrix(*j).mmv(Impl::asVector(v[j.index()]),rhs);
143 auto&& vi = Impl::asVector(v[i.index()]);
144 Impl::asMatrix(*j).mv(rhs,vi);
150 typename M::field_type& firstMatrixElement (M& A,
153 return firstMatrixElement(*(A.begin()->begin()));
157 K& firstMatrixElement (K& A,
163 template<
class K,
int n,
int m>
164 K& firstMatrixElement (FieldMatrix<K,n,m>& A)
179 typedef typename M::ColIterator coliterator;
180 typedef typename M::ConstRowIterator crowiterator;
181 typedef typename M::ConstColIterator ccoliterator;
182 typedef typename M::CreateIterator createiterator;
183 typedef typename M::field_type K;
184 typedef std::map<size_t, int> map;
185 typedef typename map::iterator mapiterator;
188 crowiterator endi=A.end();
189 createiterator ci=ILU.createbegin();
190 for (crowiterator i=A.begin(); i!=endi; ++i)
195 for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
196 rowpattern[j.index()] = 0;
199 for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
203 coliterator endk = ILU[(*ik).first].end();
204 coliterator kj = ILU[(*ik).first].find((*ik).first);
205 for (++kj; kj!=endk; ++kj)
210 int generation = (int)
Simd::lane(0, abs( firstMatrixElement(*kj) ));
213 mapiterator ij = rowpattern.find(kj.index());
214 if (ij==rowpattern.end())
216 rowpattern[kj.index()] = generation+1;
224 for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
225 ci.insert((*ik).first);
229 coliterator endILUij = ILU[i.index()].end();;
230 for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
235 for (crowiterator i=A.begin(); i!=endi; ++i)
238 coliterator endILUij = ILU[i.index()].end();;
239 for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
241 ccoliterator Aij = (*i).begin();
242 ccoliterator endAij = (*i).end();
243 ILUij = ILU[i.index()].begin();
244 while (Aij!=endAij && ILUij!=endILUij)
246 if (Aij.index()==ILUij.index())
253 if (Aij.index()<ILUij.index())
266 template <
class B,
class Alloc = std::allocator<B>>
269 typedef B block_type;
270 typedef size_t size_type;
272 CRS() : nRows_( 0 ) {}
274 size_type rows()
const {
return nRows_; }
276 size_type nonZeros()
const
278 assert( rows_[ rows() ] != size_type(-1) );
279 return rows_[ rows() ];
282 void resize(
const size_type nRows )
284 if( nRows_ != nRows )
287 rows_.resize( nRows_+1, size_type(-1) );
291 void reserveAdditional(
const size_type nonZeros )
293 const size_type needed = values_.size() + nonZeros ;
294 if( values_.capacity() < needed )
296 const size_type estimate = needed * 1.1;
297 values_.reserve( estimate );
298 cols_.reserve( estimate );
302 void push_back(
const block_type& value,
const size_type index )
304 values_.push_back( value );
305 cols_.push_back( index );
308 std::vector< size_type > rows_;
309 std::vector< block_type, Alloc> values_;
310 std::vector< size_type > cols_;
315 template<
class M,
class CRS,
class InvVector>
318 typedef typename M :: size_type size_type;
320 lower.resize( A.N() );
321 upper.resize( A.N() );
325 const size_t memEstimate = (A.nonzeroes() - A.N())/2;
327 assert( A.nonzeroes() != 0 );
328 lower.reserveAdditional( memEstimate );
329 upper.reserveAdditional( memEstimate );
331 const auto endi = A.end();
333 size_type colcount = 0;
334 lower.rows_[ 0 ] = colcount;
335 for (
auto i=A.begin(); i!=endi; ++i, ++row)
337 const size_type iIndex = i.index();
340 for (
auto j=(*i).begin(); j.index() < iIndex; ++j )
342 lower.push_back( (*j), j.index() );
345 lower.rows_[ iIndex+1 ] = colcount;
348 const auto rendi = A.beforeBegin();
351 upper.rows_[ 0 ] = colcount ;
355 for (
auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
357 const auto endij=(*i).beforeBegin();
359 const size_type iIndex = i.index();
362 for (
auto j=(*i).beforeEnd(); j != endij; --j )
364 const size_type jIndex = j.index();
365 if( j.index() == iIndex )
370 else if ( j.index() >= i.index() )
372 upper.push_back( (*j), jIndex );
376 upper.rows_[ row+1 ] = colcount;
381 template<
class CRS,
class InvVector,
class X,
class Y>
384 const InvVector& inv,
388 typedef typename Y :: block_type dblock;
389 typedef typename X :: block_type vblock;
390 typedef typename X :: size_type size_type ;
392 const size_type iEnd = lower.rows();
393 const size_type lastRow = iEnd - 1;
394 if( iEnd != upper.rows() )
400 for( size_type i=0; i<iEnd; ++ i )
402 dblock rhsValue( d[ i ] );
403 auto&& rhs = Impl::asVector(rhsValue);
404 const size_type rowI = lower.rows_[ i ];
405 const size_type rowINext = lower.rows_[ i+1 ];
407 for( size_type col = rowI; col < rowINext; ++ col )
408 Impl::asMatrix(lower.values_[ col ]).mmv( Impl::asVector(v[ lower.cols_[ col ] ] ), rhs );
410 Impl::asVector(v[ i ]) = rhs;
414 for( size_type i=0; i<iEnd; ++ i )
416 auto&& vBlock = Impl::asVector(v[ lastRow - i ]);
417 vblock rhsValue ( v[ lastRow - i ] );
418 auto&& rhs = Impl::asVector(rhsValue);
419 const size_type rowI = upper.rows_[ i ];
420 const size_type rowINext = upper.rows_[ i+1 ];
422 for( size_type col = rowI; col < rowINext; ++ col )
423 Impl::asMatrix(upper.values_[ col ]).mmv( Impl::asVector(v[ upper.cols_[ col ] ]), rhs );
426 Impl::asMatrix(inv[ i ]).mv(rhs, vBlock);
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:126
derive error class from the base class in common
Definition: istlexception.hh:19
Error when performing an operation on a matrix block.
Definition: istlexception.hh:52
Implements a matrix constructed from a given type representing a field and compile-time given number ...
void message(const std::string &msg)
store string in internal message buffer
Definition: exceptions.cc:32
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
decltype(auto) lane(std::size_t l, V &&v)
Extract an element of a SIMD type.
Definition: interface.hh:324
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
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:316
void blockILUBacksolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition: ilu.hh:103
void blockILU0Decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:35
void blockILUDecomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:176
Dune namespace.
Definition: alignedallocator.hh:13
Implements a scalar matrix view wrapper around an existing scalar.
Implements a scalar vector view wrapper around an existing scalar.
a simple compressed row storage matrix class
Definition: ilu.hh:268
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194