1 #ifndef DUNE_ISTL_ILDL_HH
2 #define DUNE_ISTL_ILDL_HH
19 template<
class K,
int m,
int n >
22 for(
int i = 0; i < m; ++i )
24 for(
int j = 0; j < n; ++j )
26 for(
int k = 0; k < n; ++k )
27 A[ i ][ j ] -= B[ i ][ k ] * CT[ j ][ k ];
32 template<
class Matrix >
35 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
38 auto &&B_i = B[ i.index() ];
39 const auto ikend = B_i.
end();
40 for(
auto j = A_i.begin(), jend = A_i.end(); j != jend; ++j )
43 auto &&CT_j = CT[ j.index() ];
44 const auto jkend = CT_j.
end();
45 for(
auto ik = B_i.begin(), jk = CT_j.begin(); (ik != ikend) && (jk != jkend); )
47 if( ik.index() == jk.index() )
52 else if( ik.index() < jk.index() )
75 template<
class Matrix >
78 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
82 auto ij = A_i.begin();
83 for( ; ij.index() < i.index(); ++ij )
86 auto &&A_j = A[ ij.index() ];
90 auto ik = A_i.
begin();
91 auto jk = A_j.begin();
92 while( (ik != ij) && (jk.index() < ij.index()) )
94 if( ik.index() == jk.index() )
99 else if( ik.index() < jk.index() )
106 if( ij.index() != i.index() )
107 DUNE_THROW(
ISTLError,
"diagonal entry missing" );
111 for(
auto ik = A_i.begin(); ik != ij; ++ik )
114 const auto &A_k = A[ ik.index() ];
117 A_ik.rightmultiply( *A_k.find( ik.index() ) );
124 catch(
const Dune::FMatrixError &e )
126 DUNE_THROW(
MatrixBlockError,
"ILDL failed to invert matrix block A[" << i.index() <<
"][" << ij.index() <<
"]" << e.what(); th__ex.
r = i.index(); th__ex.c = ij.index() );
136 template<
class Matrix,
class X,
class Y >
140 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
142 const auto &A_i = *i;
143 v[ i.index() ] = d[ i.index() ];
144 for(
auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
145 (*ij).mmv( v[ ij.index() ], v[ i.index() ] );
149 if( isLowerTriangular )
153 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
155 const auto &A_i = *i;
156 const auto ii = A_i.beforeEnd();
157 assert( ii.index() == i.index() );
158 auto rhs = v[ i.index() ];
159 ii->mv( rhs, v[ i.index() ] );
166 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
168 const auto &A_i = *i;
169 const auto ii = A_i.find( i.index() );
170 assert( ii.index() == i.index() );
171 auto rhs = v[ i.index() ];
172 ii->mv( rhs, v[ i.index() ] );
180 const auto &A_i = *i;
181 for(
auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
182 (*ij).mmtv( v[ i.index() ], v[ ij.index() ] );
188 #endif // #ifndef DUNE_ISTL_ILDL_HH