Dune Core Modules (unstable)

ilu.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_ILU_HH
6#define DUNE_ISTL_ILU_HH
7
8#include <sstream>
9
10#include <cmath>
11#include <complex>
12#include <map>
13#include <vector>
14
18
19#include "istlexception.hh"
20
25namespace Dune {
26
31 namespace ILU {
32
34 template<class M>
36 {
37 // iterator types
38 typedef typename M::RowIterator rowiterator;
39 typedef typename M::ColIterator coliterator;
40 typedef typename M::block_type block;
41
42 // implement left looking variant with stored inverse
43 rowiterator endi=A.end();
44 for (rowiterator i=A.begin(); i!=endi; ++i)
45 {
46 // coliterator is diagonal after the following loop
47 coliterator endij=(*i).end(); // end of row i
48 coliterator ij;
49
50 // eliminate entries left of diagonal; store L factor
51 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
52 {
53 // find A_jj which eliminates A_ij
54 coliterator jj = A[ij.index()].find(ij.index());
55
56 // compute L_ij = A_jj^-1 * A_ij
57 Impl::asMatrix(*ij).rightmultiply(Impl::asMatrix(*jj));
58
59 // modify row
60 coliterator endjk=A[ij.index()].end(); // end of row j
61 coliterator jk=jj; ++jk;
62 coliterator ik=ij; ++ik;
63 while (ik!=endij && jk!=endjk)
64 if (ik.index()==jk.index())
65 {
66 block B(*jk);
67 Impl::asMatrix(B).leftmultiply(Impl::asMatrix(*ij));
68 *ik -= B;
69 ++ik; ++jk;
70 }
71 else
72 {
73 if (ik.index()<jk.index())
74 ++ik;
75 else
76 ++jk;
77 }
78 }
79
80 // invert pivot and store it in A
81 if (ij.index()!=i.index())
82 DUNE_THROW(ISTLError,"diagonal entry missing");
83 try {
84 Impl::asMatrix(*ij).invert(); // compute inverse of diagonal block
85 }
86 catch (Dune::FMatrixError & e) {
87 std::ostringstream sstream;
88 sstream << THROWSPEC(MatrixBlockError)
89 << THROWSPEC(MatrixBlockError)
90 << "ILU failed to invert matrix block A["
91 << i.index() << "][" << ij.index() << "]" << e.what();
93 ex.message(sstream.str());
94 ex.r = i.index();
95 ex.c = ij.index();
96 throw ex;
97 }
98 }
99 }
100
102 template<class M, class X, class Y>
103 void blockILUBacksolve (const M& A, X& v, const Y& d)
104 {
105 // iterator types
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;
110
111 // lower triangular solve
112 rowiterator endi=A.end();
113 for (rowiterator i=A.begin(); i!=endi; ++i)
114 {
115 // We need to be careful here: Directly using
116 // auto rhs = Impl::asVector(d[ i.index() ]);
117 // is not OK in case this is a proxy. Hence
118 // we first have to copy the value. Notice that
119 // this is still not OK, if the vector type itself returns
120 // proxy references.
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; // Lii = I
126 }
127
128 // upper triangular solve
129 rowiterator rendi=A.beforeBegin();
130 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
131 {
132 // We need to be careful here: Directly using
133 // auto rhs = Impl::asVector(v[ i.index() ]);
134 // is not OK in case this is a proxy. Hence
135 // we first have to copy the value. Notice that
136 // this is still not OK, if the vector type itself returns
137 // proxy references.
138 vblock rhsValue(v[i.index()]);
139 auto&& rhs = Impl::asVector(rhsValue);
140 coliterator j;
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); // diagonal stores inverse!
145 }
146 }
147
148 // recursive function template to access first entry of a matrix
149 template<class M>
150 typename M::field_type& firstMatrixElement (M& A,
151 [[maybe_unused]] typename std::enable_if_t<!Dune::IsNumber<M>::value>* sfinae = nullptr)
152 {
153 return firstMatrixElement(*(A.begin()->begin()));
154 }
155
156 template<class K>
157 K& firstMatrixElement (K& A,
158 [[maybe_unused]] typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae = nullptr)
159 {
160 return A;
161 }
162
163 template<class K, int n, int m>
164 K& firstMatrixElement (FieldMatrix<K,n,m>& A)
165 {
166 return A[0][0];
167 }
168
175 template<class M>
176 void blockILUDecomposition (const M& A, int n, M& ILU)
177 {
178 // iterator types
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;
186
187 // symbolic factorization phase, store generation number in first matrix element
188 crowiterator endi=A.end();
189 createiterator ci=ILU.createbegin();
190 for (crowiterator i=A.begin(); i!=endi; ++i)
191 {
192 map rowpattern; // maps column index to generation
193
194 // initialize pattern with row of A
195 for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
196 rowpattern[j.index()] = 0;
197
198 // eliminate entries in row which are to the left of the diagonal
199 for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
200 {
201 if ((*ik).second<n)
202 {
203 coliterator endk = ILU[(*ik).first].end(); // end of row k
204 coliterator kj = ILU[(*ik).first].find((*ik).first); // diagonal in k
205 for (++kj; kj!=endk; ++kj) // row k eliminates in row i
206 {
207 // we misuse the storage to store an int. If the field_type is std::complex, we have to access the real/abs part
208 // starting from C++11, we can use std::abs to always return a real value, even if it is double/float
209 using std::abs;
210 int generation = (int) Simd::lane(0, abs( firstMatrixElement(*kj) ));
211 if (generation<n)
212 {
213 mapiterator ij = rowpattern.find(kj.index());
214 if (ij==rowpattern.end())
215 {
216 rowpattern[kj.index()] = generation+1;
217 }
218 }
219 }
220 }
221 }
222
223 // create row
224 for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
225 ci.insert((*ik).first);
226 ++ci; // now row i exist
227
228 // write generation index into entries
229 coliterator endILUij = ILU[i.index()].end();;
230 for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
231 Simd::lane(0, firstMatrixElement(*ILUij)) = (Simd::Scalar<K>) rowpattern[ILUij.index()];
232 }
233
234 // copy entries of A
235 for (crowiterator i=A.begin(); i!=endi; ++i)
236 {
237 coliterator ILUij;
238 coliterator endILUij = ILU[i.index()].end();;
239 for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
240 (*ILUij) = 0; // clear row
241 ccoliterator Aij = (*i).begin();
242 ccoliterator endAij = (*i).end();
243 ILUij = ILU[i.index()].begin();
244 while (Aij!=endAij && ILUij!=endILUij)
245 {
246 if (Aij.index()==ILUij.index())
247 {
248 *ILUij = *Aij;
249 ++Aij; ++ILUij;
250 }
251 else
252 {
253 if (Aij.index()<ILUij.index())
254 ++Aij;
255 else
256 ++ILUij;
257 }
258 }
259 }
260
261 // call decomposition on pattern
263 }
264
266 template <class B, class Alloc = std::allocator<B>>
267 struct CRS
268 {
269 typedef B block_type;
270 typedef size_t size_type;
271
272 CRS() : nRows_( 0 ) {}
273
274 size_type rows() const { return nRows_; }
275
276 size_type nonZeros() const
277 {
278 assert( rows_[ rows() ] != size_type(-1) );
279 return rows_[ rows() ];
280 }
281
282 void resize( const size_type nRows )
283 {
284 if( nRows_ != nRows )
285 {
286 nRows_ = nRows ;
287 rows_.resize( nRows_+1, size_type(-1) );
288 }
289 }
290
291 void reserveAdditional( const size_type nonZeros )
292 {
293 const size_type needed = values_.size() + nonZeros ;
294 if( values_.capacity() < needed )
295 {
296 const size_type estimate = needed * 1.1;
297 values_.reserve( estimate );
298 cols_.reserve( estimate );
299 }
300 }
301
302 void push_back( const block_type& value, const size_type index )
303 {
304 values_.push_back( value );
305 cols_.push_back( index );
306 }
307
308 std::vector< size_type > rows_;
309 std::vector< block_type, Alloc> values_;
310 std::vector< size_type > cols_;
311 size_type nRows_;
312 };
313
315 template<class M, class CRS, class InvVector>
316 void convertToCRS(const M& A, CRS& lower, CRS& upper, InvVector& inv )
317 {
318 typedef typename M :: size_type size_type;
319
320 lower.resize( A.N() );
321 upper.resize( A.N() );
322 inv.resize( A.N() );
323
324 // lower and upper triangular should store half of non zeros minus diagonal
325 const size_t memEstimate = (A.nonzeroes() - A.N())/2;
326
327 assert( A.nonzeroes() != 0 );
328 lower.reserveAdditional( memEstimate );
329 upper.reserveAdditional( memEstimate );
330
331 const auto endi = A.end();
332 size_type row = 0;
333 size_type colcount = 0;
334 lower.rows_[ 0 ] = colcount;
335 for (auto i=A.begin(); i!=endi; ++i, ++row)
336 {
337 const size_type iIndex = i.index();
338
339 // store entries left of diagonal
340 for (auto j=(*i).begin(); j.index() < iIndex; ++j )
341 {
342 lower.push_back( (*j), j.index() );
343 ++colcount;
344 }
345 lower.rows_[ iIndex+1 ] = colcount;
346 }
347
348 const auto rendi = A.beforeBegin();
349 row = 0;
350 colcount = 0;
351 upper.rows_[ 0 ] = colcount ;
352
353 // NOTE: upper and inv store entries in reverse row and col order,
354 // reverse here relative to ILU
355 for (auto i=A.beforeEnd(); i!=rendi; --i, ++ row )
356 {
357 const auto endij=(*i).beforeBegin(); // end of row i
358
359 const size_type iIndex = i.index();
360
361 // store in reverse row order for faster access during backsolve
362 for (auto j=(*i).beforeEnd(); j != endij; --j )
363 {
364 const size_type jIndex = j.index();
365 if( j.index() == iIndex )
366 {
367 inv[ row ] = (*j);
368 break; // assuming consecutive ordering of A
369 }
370 else if ( j.index() >= i.index() )
371 {
372 upper.push_back( (*j), jIndex );
373 ++colcount ;
374 }
375 }
376 upper.rows_[ row+1 ] = colcount;
377 }
378 } // end convertToCRS
379
381 template<class CRS, class InvVector, class X, class Y>
382 void blockILUBacksolve (const CRS& lower,
383 const CRS& upper,
384 const InvVector& inv,
385 X& v, const Y& d)
386 {
387 // iterator types
388 typedef typename Y :: block_type dblock;
389 typedef typename X :: block_type vblock;
390 typedef typename X :: size_type size_type ;
391
392 const size_type iEnd = lower.rows();
393 const size_type lastRow = iEnd - 1;
394 if( iEnd != upper.rows() )
395 {
396 DUNE_THROW(ISTLError,"ILU::blockILUBacksolve: lower and upper rows must be the same");
397 }
398
399 // lower triangular solve
400 for( size_type i=0; i<iEnd; ++ i )
401 {
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 ];
406
407 for( size_type col = rowI; col < rowINext; ++ col )
408 Impl::asMatrix(lower.values_[ col ]).mmv( Impl::asVector(v[ lower.cols_[ col ] ] ), rhs );
409
410 Impl::asVector(v[ i ]) = rhs; // Lii = I
411 }
412
413 // upper triangular solve
414 for( size_type i=0; i<iEnd; ++ i )
415 {
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 ];
421
422 for( size_type col = rowI; col < rowINext; ++ col )
423 Impl::asMatrix(upper.values_[ col ]).mmv( Impl::asVector(v[ upper.cols_[ col ] ]), rhs );
424
425 // apply inverse and store result
426 Impl::asMatrix(inv[ i ]).mv(rhs, vBlock);
427 }
428 }
429
430 } // end namespace ILU
431
434} // end namespace
435
436#endif
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,...)
Definition: exceptions.hh:312
const char * what() const noexcept override
output internal message buffer
Definition: exceptions.cc:37
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)