dune-istl issueshttps://gitlab.dune-project.org/core/dune-istl/issues2017-09-22T04:19:36Zhttps://gitlab.dune-project.org/core/dune-istl/issues/21Unexpected behaviour of BCRSMatrix implict buildmode2017-09-22T04:19:36ZCarsten GrĂ¤serUnexpected behaviour of BCRSMatrix implict buildmodeThe behaviour of the implicit build mode is in several ways unexpected and does not fit its documentation:
* (a) In contrast to the documentation, the allocated overflow area is not used while filling the matrix in favour of a separate `std::map`.
* (b) The allocated overflow area does always contain `4*avg` more entries than expected.
* (c) The given `avg` value is not really related to average row sizes. If you give the exact average, `compres()` may fail if the overflow size is zero.
* (d) The allocated overflow area is in fact used to balance the mismatch between the given `avg` value and rows having more entries than this.
* (e) `compress()` may fail depending on the row order: The following will fail for `row=0` but succeed for `row=5`
```cpp
using M = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
M m(6, 6, 1, 0.0, M::implicit);
for(std::size_t i=0; i<6; ++i)
m.entry(row,i) = i;
m.compress();
```
Explanation for this behaviour: When allocating storage, this is partitioned into the overflow area followed by `avg` reserved entries per row.
When more than `avg` entries are used in a row, they are stored in an additional `std::map`. On `compress()` all entries are shifted forward towards the first
free position, inserting intermediate entries from the overflow map when necessary. This explains, why early rows cannot have more than `avg`+`overflow`
entries, while later ones may consume remaining space left by preceding less populated rows.
Due to (a) there is a possible fix for several of those problems. Before compressing forward with insertion of overflow entries, one can compress backwards
without this insertion. This would solve (c)-(e) at the price of an additional shifting operation. Even more, this would allow to compute the number of non-zeros in the
first sweep and than allocate new storage if the overflow was to small instead of throwing an exception and letting the user assemble all entries again with a hopefully better guess.The behaviour of the implicit build mode is in several ways unexpected and does not fit its documentation:
* (a) In contrast to the documentation, the allocated overflow area is not used while filling the matrix in favour of a separate `std::map`.
* (b) The allocated overflow area does always contain `4*avg` more entries than expected.
* (c) The given `avg` value is not really related to average row sizes. If you give the exact average, `compres()` may fail if the overflow size is zero.
* (d) The allocated overflow area is in fact used to balance the mismatch between the given `avg` value and rows having more entries than this.
* (e) `compress()` may fail depending on the row order: The following will fail for `row=0` but succeed for `row=5`
```cpp
using M = Dune::BCRSMatrix<Dune::FieldMatrix<double,1,1>>;
M m(6, 6, 1, 0.0, M::implicit);
for(std::size_t i=0; i<6; ++i)
m.entry(row,i) = i;
m.compress();
```
Explanation for this behaviour: When allocating storage, this is partitioned into the overflow area followed by `avg` reserved entries per row.
When more than `avg` entries are used in a row, they are stored in an additional `std::map`. On `compress()` all entries are shifted forward towards the first
free position, inserting intermediate entries from the overflow map when necessary. This explains, why early rows cannot have more than `avg`+`overflow`
entries, while later ones may consume remaining space left by preceding less populated rows.
Due to (a) there is a possible fix for several of those problems. Before compressing forward with insertion of overflow entries, one can compress backwards
without this insertion. This would solve (c)-(e) at the price of an additional shifting operation. Even more, this would allow to compute the number of non-zeros in the
first sweep and than allocate new storage if the overflow was to small instead of throwing an exception and letting the user assemble all entries again with a hopefully better guess.https://gitlab.dune-project.org/core/dune-istl/issues/7Rethinking matrix classes and interfaces2017-09-28T17:14:29ZElias PippingRethinking matrix classes and interfacesThis task is supposed to collect and organise criticism and suggestions that have accumulated in flyspray/FS#1669.
So far, it has been suggested that (I'm making a few additions myself and do not both to properly trace back suggestions to their respective authors)
* `IdentityMatrix` (dune-common) and `ScaledIdentityMatrix` (dune-istl) are too close in functionality. Suggestions include:
1. Removing `IdentityMatrix` and replacing any instances with `ScaledIdentityMatrix(1)`.
2. Replacing `ScaledIdentityMatrix` with a more general template class `ScaledIdentityMatrix<MatrixClass>`, so that current uses of `ScaledIdentityMatrix` would be replaced with `ScaledIdentityMatrix<IdentityMatrix>`.
* The CRTP approach with the `DenseMatrix` class appears rather complicated, leading to issues such as `flyspray/FS#1712`.
* Sparse `blocklevel=0` Matrix classes (`DiagonalMatrix`, `ScaledIdentityMatrix`, `IdentityMatrix`, potentially more?) could be plugged into general-purpose implementations of functions like `FieldMatrix::rightmultiply` with optimal arithmetic efficiency if they provided iterators.
* Symmetric matrices could be implemented such that they take less storage (in 3D, 6 scalars instead of 9) and/or perform fewer arithmetic operations.
Update: added symmetric matrices.This task is supposed to collect and organise criticism and suggestions that have accumulated in flyspray/FS#1669.
So far, it has been suggested that (I'm making a few additions myself and do not both to properly trace back suggestions to their respective authors)
* `IdentityMatrix` (dune-common) and `ScaledIdentityMatrix` (dune-istl) are too close in functionality. Suggestions include:
1. Removing `IdentityMatrix` and replacing any instances with `ScaledIdentityMatrix(1)`.
2. Replacing `ScaledIdentityMatrix` with a more general template class `ScaledIdentityMatrix<MatrixClass>`, so that current uses of `ScaledIdentityMatrix` would be replaced with `ScaledIdentityMatrix<IdentityMatrix>`.
* The CRTP approach with the `DenseMatrix` class appears rather complicated, leading to issues such as `flyspray/FS#1712`.
* Sparse `blocklevel=0` Matrix classes (`DiagonalMatrix`, `ScaledIdentityMatrix`, `IdentityMatrix`, potentially more?) could be plugged into general-purpose implementations of functions like `FieldMatrix::rightmultiply` with optimal arithmetic efficiency if they provided iterators.
* Symmetric matrices could be implemented such that they take less storage (in 3D, 6 scalars instead of 9) and/or perform fewer arithmetic operations.
Update: added symmetric matrices.https://gitlab.dune-project.org/core/dune-istl/issues/2add indexed iterator ranges2017-09-22T04:19:36ZChristian Engweradd indexed iterator rangescurrently we can not use ranged-based for in dune-istl, because many algorithms have to know the index of an entry in the sparse matrix. This index is not the same as the `std::distance(begin,it)` and can not be computed using standard iterators.
I suggest to add a wrapper type, with the following interface:
```
template<typename T>
struct IndexedValue {
std::size_t index() const;
T & value();
const T & value() const;
};
```
This will allow us to implement a special `index_range(Matrix)` or `index_range(Vector)` which yields `IndexedValue` instead of the value directly.
These new ranges and new types can be implemented completely outside the current istl containers.
The `BCRSMatrix::umv` might now look like:
```
void umv (const X& x, Y& y) const
{
for (const auto & row : indexed_range(*this))
for (const auto & entry : indexed_range(row.value()))
entry.value().umv(x[entry.index()],y[row.index()]);
}
```
instead of
```
void umv (const X& x, Y& y) const
{
ConstRowIterator endi=end();
for (ConstRowIterator i=begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
(*j).umv(x[j.index()],y[i.index()]);
}
}
```currently we can not use ranged-based for in dune-istl, because many algorithms have to know the index of an entry in the sparse matrix. This index is not the same as the `std::distance(begin,it)` and can not be computed using standard iterators.
I suggest to add a wrapper type, with the following interface:
```
template<typename T>
struct IndexedValue {
std::size_t index() const;
T & value();
const T & value() const;
};
```
This will allow us to implement a special `index_range(Matrix)` or `index_range(Vector)` which yields `IndexedValue` instead of the value directly.
These new ranges and new types can be implemented completely outside the current istl containers.
The `BCRSMatrix::umv` might now look like:
```
void umv (const X& x, Y& y) const
{
for (const auto & row : indexed_range(*this))
for (const auto & entry : indexed_range(row.value()))
entry.value().umv(x[entry.index()],y[row.index()]);
}
```
instead of
```
void umv (const X& x, Y& y) const
{
ConstRowIterator endi=end();
for (ConstRowIterator i=begin(); i!=endi; ++i)
{
ConstColIterator endj = (*i).end();
for (ConstColIterator j=(*i).begin(); j!=endj; ++j)
(*j).umv(x[j.index()],y[i.index()]);
}
}
```