Dune Core Modules (2.6.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 <iostream>
9#include <iomanip>
10#include <string>
11#include <set>
12#include <map>
13#include <vector>
14
16#include "istlexception.hh"
17
22namespace Dune {
23
28 class MatrixBlockError : public virtual Dune::FMatrixError {
29 public:
30 int r, c;
31 };
32
33
35 template<class M>
37 {
38 // iterator types
39 typedef typename M::RowIterator rowiterator;
40 typedef typename M::ColIterator coliterator;
41 typedef typename M::block_type block;
42
43 // implement left looking variant with stored inverse
44 rowiterator endi=A.end();
45 for (rowiterator i=A.begin(); i!=endi; ++i)
46 {
47 // coliterator is diagonal after the following loop
48 coliterator endij=(*i).end(); // end of row i
49 coliterator ij;
50
51 // eliminate entries left of diagonal; store L factor
52 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
53 {
54 // find A_jj which eliminates A_ij
55 coliterator jj = A[ij.index()].find(ij.index());
56
57 // compute L_ij = A_jj^-1 * A_ij
58 (*ij).rightmultiply(*jj);
59
60 // modify row
61 coliterator endjk=A[ij.index()].end(); // end of row j
62 coliterator jk=jj; ++jk;
63 coliterator ik=ij; ++ik;
64 while (ik!=endij && jk!=endjk)
65 if (ik.index()==jk.index())
66 {
67 block B(*jk);
68 B.leftmultiply(*ij);
69 *ik -= B;
70 ++ik; ++jk;
71 }
72 else
73 {
74 if (ik.index()<jk.index())
75 ++ik;
76 else
77 ++jk;
78 }
79 }
80
81 // invert pivot and store it in A
82 if (ij.index()!=i.index())
83 DUNE_THROW(ISTLError,"diagonal entry missing");
84 try {
85 (*ij).invert(); // compute inverse of diagonal block
86 }
87 catch (Dune::FMatrixError & e) {
88 DUNE_THROW(MatrixBlockError, "ILU failed to invert matrix block A["
89 << i.index() << "][" << ij.index() << "]" << e.what();
90 th__ex.r=i.index(); th__ex.c=ij.index(););
91 }
92 }
93 }
94
96 template<class M, class X, class Y>
97 void bilu_backsolve (const M& A, X& v, const Y& d)
98 {
99 // iterator types
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;
104
105 // lower triangular solve
106 rowiterator endi=A.end();
107 for (rowiterator i=A.begin(); i!=endi; ++i)
108 {
109 dblock rhs(d[i.index()]);
110 for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
111 (*j).mmv(v[j.index()],rhs);
112 v[i.index()] = rhs; // Lii = I
113 }
114
115 // upper triangular solve
116 rowiterator rendi=A.beforeBegin();
117 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
118 {
119 vblock rhs(v[i.index()]);
120 coliterator j;
121 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
122 (*j).mmv(v[j.index()],rhs);
123 v[i.index()] = 0;
124 (*j).umv(rhs,v[i.index()]); // diagonal stores inverse!
125 }
126 }
127
128
129
130 // recursive function template to access first entry of a matrix
131 template<class M>
132 typename M::field_type& firstmatrixelement (M& A)
133 {
134 return firstmatrixelement(*(A.begin()->begin()));
135 }
136
137 template<class K, int n, int m>
138 K& firstmatrixelement (FieldMatrix<K,n,m>& A)
139 {
140 return A[0][0];
141 }
142
143 template<class K>
144 K& firstmatrixelement (FieldMatrix<K,1,1>& A)
145 {
146 return A[0][0];
147 }
148
149
156 template<class M>
157 void bilu_decomposition (const M& A, int n, M& ILU)
158 {
159 // iterator types
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;
167
168 // symbolic factorization phase, store generation number in first matrix element
169 crowiterator endi=A.end();
170 createiterator ci=ILU.createbegin();
171 for (crowiterator i=A.begin(); i!=endi; ++i)
172 {
173 // std::cout << "in row " << i.index() << std::endl;
174 map rowpattern; // maps column index to generation
175
176 // initialize pattern with row of A
177 for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
178 rowpattern[j.index()] = 0;
179
180 // eliminate entries in row which are to the left of the diagonal
181 for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
182 {
183 if ((*ik).second<n)
184 {
185 // std::cout << " eliminating " << i.index() << "," << (*ik).first
186 // << " level " << (*ik).second << std::endl;
187
188 coliterator endk = ILU[(*ik).first].end(); // end of row k
189 coliterator kj = ILU[(*ik).first].find((*ik).first); // diagonal in k
190 for (++kj; kj!=endk; ++kj) // row k eliminates in row i
191 {
192 // we misuse the storage to store an int. If the field_type is std::complex, we have to access the real part
193 // starting from C++11, we can use std::real to always return a real value, even if it is double/float
194 using std::real;
195 int generation = (int) real( firstmatrixelement(*kj) );
196 if (generation<n)
197 {
198 mapiterator ij = rowpattern.find(kj.index());
199 if (ij==rowpattern.end())
200 {
201 //std::cout << " new entry " << i.index() << "," << kj.index()
202 // << " level " << (*ik).second+1 << std::endl;
203
204 rowpattern[kj.index()] = generation+1;
205 }
206 }
207 }
208 }
209 }
210
211 // create row
212 for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
213 ci.insert((*ik).first);
214 ++ci; // now row i exist
215
216 // write generation index into entries
217 coliterator endILUij = ILU[i.index()].end();;
218 for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
219 firstmatrixelement(*ILUij) = (K) rowpattern[ILUij.index()];
220 }
221
222 // printmatrix(std::cout,ILU,"ilu pattern","row",10,2);
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
255 namespace ILU {
256
257 template <class 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 > 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
371
373 template<class CRS, class mblock, class X, class Y>
374 void bilu_backsolve (const CRS& lower,
375 const CRS& upper,
376 const std::vector< mblock >& inv,
377 X& v, const Y& d)
378 {
379 // iterator types
380 typedef typename Y :: block_type dblock;
381 typedef typename X :: block_type vblock;
382 typedef typename X :: size_type size_type ;
383
384 const size_type iEnd = lower.rows();
385 const size_type lastRow = iEnd - 1;
386 if( iEnd != upper.rows() )
387 {
388 DUNE_THROW(ISTLError,"ILU::bilu_backsolve: lower and upper rows must be the same");
389 }
390
391 // lower triangular solve
392 for( size_type i=0; i<iEnd; ++ i )
393 {
394 dblock rhs( d[ i ] );
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 {
400 lower.values_[ col ].mmv( v[ lower.cols_[ col ] ], rhs );
401 }
402
403 v[ i ] = rhs; // Lii = I
404 }
405
406 // upper triangular solve
407 for( size_type i=0; i<iEnd; ++ i )
408 {
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 ];
413
414 for( size_type col = rowI; col < rowINext; ++ col )
415 {
416 upper.values_[ col ].mmv( v[ upper.cols_[ col ] ], rhs );
417 }
418
419 // apply inverse and store result
420 inv[ i ].mv( rhs, vBlock);
421 }
422 }
423
424 } // end namespace ILU
425
426
429} // end namespace
430
431#endif
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:146
derive error class from the base class in common
Definition: istlexception.hh:16
RowIterator beforeBegin()
Definition: matrix.hh:629
RowIterator beforeEnd()
Definition: matrix.hh:622
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:615
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:609
Implements a matrix constructed from a given type representing a field and compile-time given number ...
const char * what() const noexcept override
output internal message buffer
Definition: exceptions.cc:35
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
void bilu_backsolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition: ilu.hh:97
void bilu_decomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:157
void bilu0_decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:36
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
Dune namespace.
Definition: alignedallocator.hh:10
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)