Dune Core Modules (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
15
16#include "istlexception.hh"
17
22namespace 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>
147 K& firstMatrixElement (K& A,
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>
154 K& firstMatrixElement (FieldMatrix<K,n,m>& A)
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
266 size_type nonZeros() const
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
281 void reserveAdditional( const size_type nonZeros )
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_;
301 size_type nRows_;
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.")
431 void bilu0_decomposition (M& A) { ILU::blockILU0Decomposition(A); }
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
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:151
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
Definition of the DUNE_DEPRECATED macro for the case that config.h is not available.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
#define DUNE_DEPRECATED_MSG(text)
Mark some entity as deprecated.
Definition: deprecated.hh:32
const char * what() const noexcept override
output internal message buffer
Definition: exceptions.cc:35
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
decltype(auto) lane(std::size_t l, V &&v)
Extract an element of a SIMD type.
Definition: interface.hh:322
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:233
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
void blockILU0Decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:32
void blockILUDecomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:166
Dune namespace.
Definition: alignedallocator.hh:11
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
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:258
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)