"Fossies" - the Fresh Open Source Software Archive  

Source code changes of the file "dune/istl/ildl.hh" between
dune-istl-2.6.0.tar.gz and dune-istl-2.7.0.tar.gz

About: dune-istl - DUNE (Distributed and Unified Numerics Environment) is a modular toolbox for solving partial differential equations (PDEs) with grid-based methods: DUNE Iterative Solver Template Library.

ildl.hh  (dune-istl-2.6.0):ildl.hh  (dune-istl-2.7.0)
#ifndef DUNE_ISTL_ILDL_HH #ifndef DUNE_ISTL_ILDL_HH
#define DUNE_ISTL_ILDL_HH #define DUNE_ISTL_ILDL_HH
#include <dune/common/scalarvectorview.hh>
#include <dune/common/scalarmatrixview.hh>
#include "ilu.hh" #include "ilu.hh"
/** /**
* \file * \file
* *
* \brief Incomplete LDL decomposition * \brief Incomplete LDL decomposition
* \author Martin Nolte * \author Martin Nolte
**/ **/
namespace Dune namespace Dune
skipping to change at line 32 skipping to change at line 34
for( int i = 0; i < m; ++i ) for( int i = 0; i < m; ++i )
{ {
for( int j = 0; j < n; ++j ) for( int j = 0; j < n; ++j )
{ {
for( int k = 0; k < n; ++k ) for( int k = 0; k < n; ++k )
A[ i ][ j ] -= B[ i ][ k ] * CT[ j ][ k ]; A[ i ][ j ] -= B[ i ][ k ] * CT[ j ][ k ];
} }
} }
} }
template< class K >
inline static void bildl_subtractBCT ( const K &B, const K &CT, K &A,
typename std::enable_if_t<Dune::IsNumbe
r<K>::value>* sfinae = nullptr )
{
A -= B * CT;
}
template< class Matrix > template< class Matrix >
inline static void bildl_subtractBCT ( const Matrix &B, const Matrix &CT, Matr inline static void bildl_subtractBCT ( const Matrix &B, const Matrix &CT, Matr
ix &A ) ix &A,
typename std::enable_if_t<!Dune::IsNumb
er<Matrix>::value>* sfinae = nullptr )
{ {
for( auto i = A.begin(), iend = A.end(); i != iend; ++i ) for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
{ {
auto &&A_i = *i; auto &&A_i = *i;
auto &&B_i = B[ i.index() ]; auto &&B_i = B[ i.index() ];
const auto ikend = B_i.end(); const auto ikend = B_i.end();
for( auto j = A_i.begin(), jend = A_i.end(); j != jend; ++j ) for( auto j = A_i.begin(), jend = A_i.end(); j != jend; ++j )
{ {
auto &&A_ij = *j; auto &&A_ij = *j;
auto &&CT_j = CT[ j.index() ]; auto &&CT_j = CT[ j.index() ];
skipping to change at line 115 skipping to change at line 125
DUNE_THROW( ISTLError, "diagonal entry missing" ); DUNE_THROW( ISTLError, "diagonal entry missing" );
// update diagonal and multiply A_ij by D_j^{-1} // update diagonal and multiply A_ij by D_j^{-1}
auto &&A_ii = *ij; auto &&A_ii = *ij;
for( auto ik = A_i.begin(); ik != ij; ++ik ) for( auto ik = A_i.begin(); ik != ij; ++ik )
{ {
auto &&A_ik = *ik; auto &&A_ik = *ik;
const auto &A_k = A[ ik.index() ]; const auto &A_k = A[ ik.index() ];
auto B = A_ik; auto B = A_ik;
A_ik.rightmultiply( *A_k.find( ik.index() ) ); Impl::asMatrix(A_ik).rightmultiply( Impl::asMatrix(*A_k.find( ik.index() )) );
bildl_subtractBCT( B, A_ik, A_ii ); bildl_subtractBCT( B, A_ik, A_ii );
} }
try try
{ {
A_ii.invert(); Impl::asMatrix(A_ii).invert();
} }
catch( const Dune::FMatrixError &e ) catch( const Dune::FMatrixError &e )
{ {
DUNE_THROW( MatrixBlockError, "ILDL failed to invert matrix block A[" << i.index() << "][" << ij.index() << "]" << e.what(); th__ex.r = i.index(); th__e x.c = ij.index() ); DUNE_THROW( MatrixBlockError, "ILDL failed to invert matrix block A[" << i.index() << "][" << ij.index() << "]" << e.what(); th__ex.r = i.index(); th__e x.c = ij.index() );
} }
} }
} }
// bildl_backsolve // bildl_backsolve
// --------------- // ---------------
template< class Matrix, class X, class Y > template< class Matrix, class X, class Y >
inline void bildl_backsolve ( const Matrix &A, X &v, const Y &d, bool isLowerT riangular = false ) inline void bildl_backsolve ( const Matrix &A, X &v, const Y &d, bool isLowerT riangular = false )
{ {
// solve L v = d, note: Lii = I // solve L v = d, note: Lii = I
for( auto i = A.begin(), iend = A.end(); i != iend; ++i ) for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
{ {
const auto &A_i = *i; const auto &A_i = *i;
v[ i.index() ] = d[ i.index() ]; v[ i.index() ] = d[ i.index() ];
for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij ) for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
(*ij).mmv( v[ ij.index() ], v[ i.index() ] ); {
auto&& vi = Impl::asVector( v[ i.index() ] );
Impl::asMatrix(*ij).mmv(Impl::asVector( v[ ij.index() ] ), vi);
}
} }
// solve D w = v, note: diagonal stores Dii^{-1} // solve D w = v, note: diagonal stores Dii^{-1}
if( isLowerTriangular ) if( isLowerTriangular )
{ {
// The matrix is lower triangular, so the diagonal entry is the // The matrix is lower triangular, so the diagonal entry is the
// last one in each row. // last one in each row.
for( auto i = A.begin(), iend = A.end(); i != iend; ++i ) for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
{ {
const auto &A_i = *i; const auto &A_i = *i;
const auto ii = A_i.beforeEnd(); const auto ii = A_i.beforeEnd();
assert( ii.index() == i.index() ); assert( ii.index() == i.index() );
auto rhs = v[ i.index() ]; // We need to be careful here: Directly using
ii->mv( rhs, v[ i.index() ] ); // auto rhs = Impl::asVector(v[ i.index() ]);
// is not OK in case this is a proxy. Hence
// we first have to copy the value. Notice that
// this is still not OK, if the vector type itself returns
// proxy references.
auto rhsValue = v[ i.index() ];
auto&& rhs = Impl::asVector(rhsValue);
auto&& vi = Impl::asVector( v[ i.index() ] );
Impl::asMatrix(*ii).mv(rhs, vi);
} }
} }
else else
{ {
// Without assumptions on the sparsity pattern we have to search // Without assumptions on the sparsity pattern we have to search
// for the diagonal entry in each row. // for the diagonal entry in each row.
for( auto i = A.begin(), iend = A.end(); i != iend; ++i ) for( auto i = A.begin(), iend = A.end(); i != iend; ++i )
{ {
const auto &A_i = *i; const auto &A_i = *i;
const auto ii = A_i.find( i.index() ); const auto ii = A_i.find( i.index() );
assert( ii.index() == i.index() ); assert( ii.index() == i.index() );
auto rhs = v[ i.index() ]; // We need to be careful here: Directly using
ii->mv( rhs, v[ i.index() ] ); // auto rhs = Impl::asVector(v[ i.index() ]);
// is not OK in case this is a proxy. Hence
// we first have to copy the value. Notice that
// this is still not OK, if the vector type itself returns
// proxy references.
auto rhsValue = v[ i.index() ];
auto&& rhs = Impl::asVector(rhsValue);
auto&& vi = Impl::asVector( v[ i.index() ] );
Impl::asMatrix(*ii).mv(rhs, vi);
} }
} }
// solve L^T v = w, note: only L is stored // solve L^T v = w, note: only L is stored
// note: we perform the operation column-wise from right to left // note: we perform the operation column-wise from right to left
for( auto i = A.beforeEnd(), iend = A.beforeBegin(); i != iend; --i ) for( auto i = A.beforeEnd(), iend = A.beforeBegin(); i != iend; --i )
{ {
const auto &A_i = *i; const auto &A_i = *i;
for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij ) for( auto ij = A_i.begin(); ij.index() < i.index(); ++ij )
(*ij).mmtv( v[ i.index() ], v[ ij.index() ] ); {
auto&& vij = Impl::asVector( v[ ij.index() ] );
Impl::asMatrix(*ij).mmtv(Impl::asVector( v[ i.index() ] ), vij);
}
} }
} }
} // namespace Dune } // namespace Dune
#endif // #ifndef DUNE_ISTL_ILDL_HH #endif // #ifndef DUNE_ISTL_ILDL_HH
 End of changes. 9 change blocks. 
10 lines changed or deleted 44 lines changed or added

Home  |  About  |  Features  |  All  |  Newest  |  Dox  |  Diffs  |  RSS Feeds  |  Screenshots  |  Comments  |  Imprint  |  Privacy  |  HTTP(S)