3#ifndef DUNE_ISTL_ILU_HH
4#define DUNE_ISTL_ILU_HH
15#include "istlexception.hh"
38 typedef typename M::RowIterator rowiterator;
39 typedef typename M::ColIterator coliterator;
40 typedef typename M::block_type block;
43 rowiterator endi=A.
end();
44 for (rowiterator i=A.
begin(); i!=endi; ++i)
47 coliterator endij=(*i).end();
51 for (ij=(*i).begin(); ij.index()<i.index(); ++ij)
54 coliterator jj = A[ij.index()].find(ij.index());
57 (*ij).rightmultiply(*jj);
60 coliterator endjk=A[ij.index()].
end();
61 coliterator jk=jj; ++jk;
62 coliterator ik=ij; ++ik;
63 while (ik!=endij && jk!=endjk)
64 if (ik.index()==jk.index())
73 if (ik.index()<jk.index())
81 if (ij.index()!=i.index())
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(););
95 template<
class M,
class X,
class Y>
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;
105 rowiterator endi=A.
end();
106 for (rowiterator i=A.
begin(); i!=endi; ++i)
108 dblock rhs(d[i.index()]);
109 for (coliterator j=(*i).begin(); j.index()<i.index(); ++j)
110 (*j).mmv(v[j.index()],rhs);
116 for (rowiterator i=A.
beforeEnd(); i!=rendi; --i)
118 vblock rhs(v[i.index()]);
120 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
121 (*j).mmv(v[j.index()],rhs);
123 (*j).umv(rhs,v[i.index()]);
131 typename M::field_type& firstmatrixelement (M& A)
133 return firstmatrixelement(*(A.
begin()->
begin()));
136 template<
class K,
int n,
int m>
137 K& firstmatrixelement (FieldMatrix<K,n,m>& A)
143 K& firstmatrixelement (FieldMatrix<K,1,1>& A)
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;
168 crowiterator endi=A.
end();
169 createiterator ci=ILU.createbegin();
170 for (crowiterator i=A.
begin(); i!=endi; ++i)
176 for (ccoliterator j=(*i).begin(); j!=(*i).end(); ++j)
177 rowpattern[j.index()] = 0;
180 for (mapiterator ik=rowpattern.begin(); (*ik).first<i.index(); ++ik)
187 coliterator endk = ILU[(*ik).first].end();
188 coliterator kj = ILU[(*ik).first].find((*ik).first);
189 for (++kj; kj!=endk; ++kj)
193 int generation = (int) std::real( firstmatrixelement(*kj) );
196 mapiterator ij = rowpattern.find(kj.index());
197 if (ij==rowpattern.end())
202 rowpattern[kj.index()] = generation+1;
210 for (mapiterator ik=rowpattern.begin(); ik!=rowpattern.end(); ++ik)
211 ci.insert((*ik).first);
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()];
223 for (crowiterator i=A.
begin(); i!=endi; ++i)
226 coliterator endILUij = ILU[i.index()].end();;
227 for (ILUij=ILU[i.index()].begin(); ILUij!=endILUij; ++ILUij)
229 ccoliterator Aij = (*i).begin();
230 ccoliterator endAij = (*i).end();
231 ILUij = ILU[i.index()].begin();
232 while (Aij!=endAij && ILUij!=endILUij)
234 if (Aij.index()==ILUij.index())
241 if (Aij.index()<ILUij.index())
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