Dune Core Modules (2.4.1)

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
15#include "istlexception.hh"
16
21namespace Dune {
22
27 class MatrixBlockError : public virtual Dune::FMatrixError {
28 public:
29 int r, c;
30 };
31
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 (*ij).rightmultiply(*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 B.leftmultiply(*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 (*ij).invert(); // compute inverse of diagonal block
85 }
86 catch (Dune::FMatrixError & e) {
87 DUNE_THROW(MatrixBlockError, "ILU failed to invert matrix block A["
88 << i.index() << "][" << ij.index() << "]" << e.what();
89 th__ex.r=i.index(); th__ex.c=ij.index(););
90 }
91 }
92 }
93
95 template<class M, class X, class Y>
96 void bilu_backsolve (const M& A, X& v, const Y& d)
97 {
98 // iterator types
99 typedef typename M::ConstRowIterator rowiterator;
100 typedef typename M::ConstColIterator coliterator;
101 typedef typename Y::block_type dblock;
102 typedef typename X::block_type vblock;
103
104 // lower triangular solve
105 rowiterator endi=A.end();
106 for (rowiterator i=A.begin(); i!=endi; ++i)
107 {
108 dblock rhs(d[i.index()]);
109 for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
110 (*j).mmv(v[j.index()],rhs);
111 v[i.index()] = rhs; // Lii = I
112 }
113
114 // upper triangular solve
115 rowiterator rendi=A.beforeBegin();
116 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
117 {
118 vblock rhs(v[i.index()]);
119 coliterator j;
120 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
121 (*j).mmv(v[j.index()],rhs);
122 v[i.index()] = 0;
123 (*j).umv(rhs,v[i.index()]); // diagonal stores inverse!
124 }
125 }
126
127
128
129 // recursive function template to access first entry of a matrix
130 template<class M>
131 typename M::field_type& firstmatrixelement (M& A)
132 {
133 return firstmatrixelement(*(A.begin()->begin()));
134 }
135
136 template<class K, int n, int m>
137 K& firstmatrixelement (FieldMatrix<K,n,m>& A)
138 {
139 return A[0][0];
140 }
141
142 template<class K>
143 K& firstmatrixelement (FieldMatrix<K,1,1>& A)
144 {
145 return A[0][0];
146 }
147
148
155 template<class M>
156 void bilu_decomposition (const M& A, int n, M& ILU)
157 {
158 // iterator types
159 typedef typename M::ColIterator coliterator;
160 typedef typename M::ConstRowIterator crowiterator;
161 typedef typename M::ConstColIterator ccoliterator;
162 typedef typename M::CreateIterator createiterator;
163 typedef typename M::field_type K;
164 typedef std::map<size_t, int> map;
165 typedef typename map::iterator mapiterator;
166
167 // symbolic factorization phase, store generation number in first matrix element
168 crowiterator endi=A.end();
169 createiterator ci=ILU.createbegin();
170 for (crowiterator i=A.begin(); i!=endi; ++i)
171 {
172 // std::cout << "in row " << i.index() << std::endl;
173 map rowpattern; // maps column index to generation
174
175 // initialize pattern with row of A
176 for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
177 rowpattern[j.index()] = 0;
178
179 // eliminate entries in row which are to the left of the diagonal
180 for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
181 {
182 if ((*ik).second<n)
183 {
184 // std::cout << " eliminating " << i.index() << "," << (*ik).first
185 // << " level " << (*ik).second << std::endl;
186
187 coliterator endk = ILU[(*ik).first].end(); // end of row k
188 coliterator kj = ILU[(*ik).first].find((*ik).first); // diagonal in k
189 for (++kj; kj!=endk; ++kj) // row k eliminates in row i
190 {
191 // we misuse the storage to store an int. If the field_type is std::complex, we have to access the real part
192 // starting from C++11, we can use std::real to always return a real value, even if it is double/float
193 int generation = (int) std::real( firstmatrixelement(*kj) );
194 if (generation<n)
195 {
196 mapiterator ij = rowpattern.find(kj.index());
197 if (ij==rowpattern.end())
198 {
199 //std::cout << " new entry " << i.index() << "," << kj.index()
200 // << " level " << (*ik).second+1 << std::endl;
201
202 rowpattern[kj.index()] = generation+1;
203 }
204 }
205 }
206 }
207 }
208
209 // create row
210 for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
211 ci.insert((*ik).first);
212 ++ci; // now row i exist
213
214 // write generation index into entries
215 coliterator endILUij = ILU[i.index()].end();;
216 for (coliterator ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
217 firstmatrixelement(*ILUij) = (K) rowpattern[ILUij.index()];
218 }
219
220 // printmatrix(std::cout,ILU,"ilu pattern","row",10,2);
221
222 // copy entries of A
223 for (crowiterator i=A.begin(); i!=endi; ++i)
224 {
225 coliterator ILUij;
226 coliterator endILUij = ILU[i.index()].end();;
227 for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
228 (*ILUij) = 0; // clear row
229 ccoliterator Aij = (*i).begin();
230 ccoliterator endAij = (*i).end();
231 ILUij = ILU[i.index()].begin();
232 while (Aij!=endAij && ILUij!=endILUij)
233 {
234 if (Aij.index()==ILUij.index())
235 {
236 *ILUij = *Aij;
237 ++Aij; ++ILUij;
238 }
239 else
240 {
241 if (Aij.index()<ILUij.index())
242 ++Aij;
243 else
244 ++ILUij;
245 }
246 }
247 }
248
249 // call decomposition on pattern
251 }
252
253
256} // end namespace
257
258#endif
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:171
derive error class from the base class in common
Definition: istlexception.hh:16
RowIterator beforeBegin()
Definition: matrix.hh:99
RowIterator beforeEnd()
Definition: matrix.hh:92
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:85
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:79
iterator begin()
begin iterator
Definition: basearray.hh:160
Implements a matrix constructed from a given type representing a field and compile-time given number ...
const std::string & what() const
output internal message buffer
Definition: exceptions.hh:199
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
void bilu_backsolve(const M &A, X &v, const Y &d)
LU backsolve with stored inverse.
Definition: ilu.hh:96
void bilu_decomposition(const M &A, int n, M &ILU)
Definition: ilu.hh:156
void bilu0_decomposition(M &A)
compute ILU decomposition of A. A is overwritten by its decomposition
Definition: ilu.hh:35
Dune namespace.
Definition: alignment.hh:10
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 13, 23:29, 2024)