3#ifndef DUNE_ISTL_ILDL_HH
4#define DUNE_ISTL_ILDL_HH
25 template<
class K,
int m,
int n >
28 for(
int i = 0; i < m; ++i )
30 for(
int j = 0; j < n; ++j )
32 for(
int k = 0; k < n; ++k )
33 A[ i ][ j ] -= B[ i ][ k ] * CT[ j ][ k ];
39 inline static void bildl_subtractBCT (
const K &B,
const K &CT, K &A,
45 template<
class Matrix >
46 inline static void bildl_subtractBCT (
const Matrix &B,
const Matrix &CT, Matrix &A,
49 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
52 auto &&B_i = B[ i.index() ];
53 const auto ikend = B_i.end();
54 for(
auto j = A_i.begin(), jend = A_i.end(); j != jend; ++j )
57 auto &&CT_j = CT[ j.index() ];
58 const auto jkend = CT_j.end();
59 for(
auto ik = B_i.begin(), jk = CT_j.begin(); (ik != ikend) && (jk != jkend); )
61 if( ik.index() == jk.index() )
63 bildl_subtractBCT( *ik, *jk, A_ij );
66 else if( ik.index() < jk.index() )
89 template<
class Matrix >
92 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
96 auto ij = A_i.begin();
97 for( ; ij.index() < i.index(); ++ij )
100 auto &&A_j = A[ ij.index() ];
104 auto ik = A_i.
begin();
105 auto jk = A_j.begin();
106 while( (ik != ij) && (jk.index() < ij.index()) )
108 if( ik.index() == jk.index() )
110 bildl_subtractBCT(*ik, *jk, A_ij);
113 else if( ik.index() < jk.index() )
120 if( ij.index() != i.index() )
125 for(
auto ik = A_i.begin(); ik != ij; ++ik )
128 const auto &A_k = A[ ik.index() ];
131 Impl::asMatrix(A_ik).rightmultiply( Impl::asMatrix(*A_k.find( ik.index() )) );
132 bildl_subtractBCT( B, A_ik, A_ii );
136 Impl::asMatrix(A_ii).invert();
140 std::ostringstream sstream;
142 <<
"ILDL failed to invert matrix block A[" << i.index() <<
"][" << ij.index() <<
"]" << e.
what();
157 template<
class Matrix,
class X,
class Y >
158 inline void bildl_backsolve (
const Matrix &A, X &v,
const Y &d,
bool isLowerTriangular =
false )
161 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
163 const auto &A_i = *i;
164 v[ i.index() ] = d[ i.index() ];
165 for(
auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
167 auto&& vi = Impl::asVector( v[ i.index() ] );
168 Impl::asMatrix(*ij).mmv(Impl::asVector( v[ ij.index() ] ), vi);
173 if( isLowerTriangular )
177 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
179 const auto &A_i = *i;
180 const auto ii = A_i.beforeEnd();
181 assert( ii.index() == i.index() );
188 auto rhsValue = v[ i.index() ];
189 auto&& rhs = Impl::asVector(rhsValue);
190 auto&& vi = Impl::asVector( v[ i.index() ] );
191 Impl::asMatrix(*ii).mv(rhs, vi);
198 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
200 const auto &A_i = *i;
201 const auto ii = A_i.find( i.index() );
202 assert( ii.index() == i.index() );
209 auto rhsValue = v[ i.index() ];
210 auto&& rhs = Impl::asVector(rhsValue);
211 auto&& vi = Impl::asVector( v[ i.index() ] );
212 Impl::asMatrix(*ii).mv(rhs, vi);
220 const auto &A_i = *i;
221 for(
auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
223 auto&& vij = Impl::asVector( v[ ij.index() ] );
224 Impl::asMatrix(*ij).mmtv(Impl::asVector( v[ i.index() ] ), vij);
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
A generic dynamic dense matrix.
Definition: matrix.hh:561
RowIterator beforeBegin()
Definition: matrix.hh:630
RowIterator beforeEnd()
Definition: matrix.hh:623
RowIterator end()
Get iterator to one beyond last row.
Definition: matrix.hh:616
RowIterator begin()
Get iterator to first row.
Definition: matrix.hh:610
void message(const std::string &msg)
store string in internal message buffer
Definition: exceptions.cc:32
const char * what() const noexcept override
output internal message buffer
Definition: exceptions.cc:37
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
The incomplete LU factorization kernels.
Dune namespace.
Definition: alignedallocator.hh:13
void bildl_decompose(Matrix &A)
compute ILDL decomposition of a symmetric matrix A
Definition: ildl.hh:90
Implements a scalar matrix view wrapper around an existing scalar.
Implements a scalar vector view wrapper around an existing scalar.
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition: typetraits.hh:194