DUNE PDELab (git)

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