1 #ifndef DUNE_ISTL_ILDL_HH
2 #define DUNE_ISTL_ILDL_HH
4 #include <dune/common/scalarvectorview.hh>
5 #include <dune/common/scalarmatrixview.hh>
21 template<
class K,
int m,
int n >
24 for(
int i = 0; i < m; ++i )
26 for(
int j = 0; j < n; ++j )
28 for(
int k = 0; k < n; ++k )
29 A[ i ][ j ] -= B[ i ][ k ] * CT[ j ][ k ];
36 typename std::enable_if_t<Dune::IsNumber<K>::value>* sfinae =
nullptr )
41 template<
class Matrix >
43 typename std::enable_if_t<!Dune::IsNumber<Matrix>::value>* sfinae =
nullptr )
45 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
48 auto &&B_i = B[ i.index() ];
49 const auto ikend = B_i.
end();
50 for(
auto j = A_i.begin(), jend = A_i.end(); j != jend; ++j )
53 auto &&CT_j = CT[ j.index() ];
54 const auto jkend = CT_j.
end();
55 for(
auto ik = B_i.begin(), jk = CT_j.begin(); (ik != ikend) && (jk != jkend); )
57 if( ik.index() == jk.index() )
62 else if( ik.index() < jk.index() )
85 template<
class Matrix >
88 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
92 auto ij = A_i.begin();
93 for( ; ij.index() < i.index(); ++ij )
96 auto &&A_j = A[ ij.index() ];
100 auto ik = A_i.
begin();
101 auto jk = A_j.begin();
102 while( (ik != ij) && (jk.index() < ij.index()) )
104 if( ik.index() == jk.index() )
109 else if( ik.index() < jk.index() )
116 if( ij.index() != i.index() )
117 DUNE_THROW(
ISTLError,
"diagonal entry missing" );
121 for(
auto ik = A_i.begin(); ik != ij; ++ik )
124 const auto &A_k = A[ ik.index() ];
127 Impl::asMatrix(A_ik).rightmultiply( Impl::asMatrix(*A_k.find( ik.index() )) );
132 Impl::asMatrix(A_ii).invert();
134 catch(
const Dune::FMatrixError &e )
136 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() );
146 template<
class Matrix,
class X,
class Y >
150 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
152 const auto &A_i = *i;
153 v[ i.index() ] = d[ i.index() ];
154 for(
auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
156 auto&& vi = Impl::asVector( v[ i.index() ] );
157 Impl::asMatrix(*ij).mmv(Impl::asVector( v[ ij.index() ] ), vi);
162 if( isLowerTriangular )
166 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
168 const auto &A_i = *i;
169 const auto ii = A_i.beforeEnd();
170 assert( ii.index() == i.index() );
177 auto rhsValue = v[ i.index() ];
178 auto&& rhs = Impl::asVector(rhsValue);
179 auto&& vi = Impl::asVector( v[ i.index() ] );
180 Impl::asMatrix(*ii).mv(rhs, vi);
187 for(
auto i = A.
begin(), iend = A.
end(); i != iend; ++i )
189 const auto &A_i = *i;
190 const auto ii = A_i.find( i.index() );
191 assert( ii.index() == i.index() );
198 auto rhsValue = v[ i.index() ];
199 auto&& rhs = Impl::asVector(rhsValue);
200 auto&& vi = Impl::asVector( v[ i.index() ] );
201 Impl::asMatrix(*ii).mv(rhs, vi);
209 const auto &A_i = *i;
210 for(
auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
212 auto&& vij = Impl::asVector( v[ ij.index() ] );
213 Impl::asMatrix(*ij).mmtv(Impl::asVector( v[ i.index() ] ), vij);
220 #endif // #ifndef DUNE_ISTL_ILDL_HH