Skip to content
Snippets Groups Projects
Commit 05c0a63a authored by Markus Blatt's avatar Markus Blatt
Browse files

[release] Do not assume a symmetric sparsity pattern in SymmetricDependency.

Previously, we assumed that if a_{ij} is stored in the sparse matrix
then a_{ji} must be stored also and used mat[i][j] to access. If the entry
was not stored then an exception somewhere in basearray.hh was thrown.
Unfortunately knowing the cause for this exception is insider knowledge,
seldomly leaked ;). In addition there is now reason to throw an expection here.

With this commit we the find method on the row to search for the entry a_{ji}.
If it is not present we simply treat it as being zero.
parent 2e9dda17
No related branches found
No related tags found
No related merge requests found
......@@ -1409,7 +1409,14 @@ namespace Dune
inline void SymmetricDependency<M,N>::examine(const ColIter& col)
{
real_type eij = norm_(*col);
real_type eji = norm_(matrix_->operator[](col.index())[row_]);
typename Matrix::ConstColIterator opposite_entry =
matrix_->operator[](col.index()).find(row_);
if ( opposite_entry == matrix_->operator[](col.index()).end() )
{
// Consider this a weak connection we disregard.
return;
}
real_type eji = norm_(*opposite_entry);
// skip positive offdiagonals if norm preserves sign of them.
if(!N::is_sign_preserving || eij<0 || eji<0)
......@@ -1423,15 +1430,21 @@ namespace Dune
inline void SymmetricDependency<M,N>::examine(G& graph, const typename G::EdgeIterator& edge, const ColIter& col)
{
real_type eij = norm_(*col);
real_type eji = norm_(matrix_->operator[](col.index())[row_]);
typename Matrix::ConstColIterator opposite_entry =
matrix_->operator[](col.index()).find(row_);
if ( opposite_entry == matrix_->operator[](col.index()).end() )
{
// Consider this as a weak connection we disregard.
return;
}
real_type eji = norm_(*opposite_entry);
// skip positve offdiagonals if norm preserves sign of them.
if(!N::is_sign_preserving || (eij<0 || eji<0))
if(eji / norm_(matrix_->operator[](edge.target())[edge.target()]) *
eij/ diagonal_ > alpha() * maxValue_) {
edge.properties().setDepends();
edge.properties().setInfluences();
typename G::EdgeProperties& other = graph.getEdgeProperties(edge.target(), edge.source());
other.setInfluences();
other.setDepends();
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment