dune-istl issueshttps://gitlab.dune-project.org/core/dune-istl/-/issues2020-01-27T11:17:00Zhttps://gitlab.dune-project.org/core/dune-istl/-/issues/88Use custom exceptions for solverfactory.2020-01-27T11:17:00ZMarkus BlattUse custom exceptions for solverfactory.We are currently throwing rather generic exceptions if something is wrong with the configuration.
Unfortunately, in user simulation code this might result in the algorithm misjudging this as convergence problems (e.g. in a Newton metho...We are currently throwing rather generic exceptions if something is wrong with the configuration.
Unfortunately, in user simulation code this might result in the algorithm misjudging this as convergence problems (e.g. in a Newton method) and will run the whole simulation.https://gitlab.dune-project.org/core/dune-istl/-/issues/84Inconsistencies in naming conventions in dune-istl2020-01-21T10:30:22ZTimo KochInconsistencies in naming conventions in dune-istlClasses in dune-istl seem to follow the CamelCase naming scheme used throughout Dune. However free functions and member functions, as well as data members and aliases, seem to be an arbitrary mix between something like STL-style naming c...Classes in dune-istl seem to follow the CamelCase naming scheme used throughout Dune. However free functions and member functions, as well as data members and aliases, seem to be an arbitrary mix between something like STL-style naming conventions (all lowercase or lowercase+underscore) and "Dune naming conventions" (camelCase). Differences even occur in the same scope and in the same file (without visible reason).
For example:
* `dune/istl/io.hh` contains the free functions
```c++
void recursive_printvector(..);
void printvector (...);
inline void fill_row (...);
void print_row (...);
void printmatrix (...);
void printSparseMatrix(...);
void writeMatrixToMatlabHelper(...);
void writeMatrixToMatlab(...);
void writeVectorToMatlabHelper(...);
void writeVectorToMatlab(...);
```
* `dune/istl/common/registry.hh` (addressed in !350) contains a macro (usually all uppercase in Dune)
```c++
#define registry_put(Tag, id, ...) ...
```
and the free functions
```c++
auto registry_get(...);
int addRegistryToFactory(...);
```
* In `dune/istl/bcrs.hh`, `Dune::BCRSMatrix` has the typedefs
```c++
typedef Imp::CompressedBlockVectorWindow<B,A> row_type;
typedef ::Dune::CompressionStatistics<size_type> CompressionStatistics;
```
and the public data member
```c++
static constexpr unsigned int blocklevel = ...;
```
the public member functions
```c++
beforeBegin ();
beforeEnd ();
void setBuildMode(BuildMode bm);
void setSize();
void setImplicitBuildModeParameters();
CreateIterator createbegin ();
CreateIterator createend ();
void setrowsize ();
void getrowsize ();
void incrementrowsize ();
void addindex ();
void setIndices(...);
real_type frobenius_norm2 () const;
real_type infinity_norm() const;
N ();
M ();
nonzeroes();
BuildStage buildStage() const;
BuildMode buildMode() const;
```
and some private member functions like
```c++
void allocateData();
void implicit_allocate(...);
```
This mix looks unprofessional/ugly to me and makes it harder for users to guess the correct name of a function. It also makes the code harder to read because you wonder what the difference between two functions is that have different naming scheme ("something should be different right? because it looks different..."). I think I listed some rather surprising cases above.
Does it make sense to unify the naming scheme? Deprecate old names and promote new names? Or is that not bothering anyone else (enough)? Which naming scheme would be appropriate (I would guess camelCase for functions to be consistent with the rest of Dune. However there are also interfaces in dune-common like `DenseVector::two_norm` which are not camelCase although they are Dune-specific types. I'm willing to do some cleanup work if someone cares and the objective/rules are clear.
There is has been a related discussion about some of the dune-common types in https://gitlab.dune-project.org/core/dune-common/issues/83. So there might be already some rules that were decided on that I don't know about...https://gitlab.dune-project.org/core/dune-istl/-/issues/80matrixredisttest leaks memory2019-12-17T09:35:11ZOliver Sanderoliver.sander@tu-dresden.dematrixredisttest leaks memoryWhen running `matrixredisttest` with the valgrind leak checker, it reports an error somewhere in `CollectiveCommunication::allreduce`:
==16702== 1,344 bytes in 1 blocks are definitely lost in loss record 65 of 66
==16702== at...When running `matrixredisttest` with the valgrind leak checker, it reports an error somewhere in `CollectiveCommunication::allreduce`:
==16702== 1,344 bytes in 1 blocks are definitely lost in loss record 65 of 66
==16702== at 0x483577F: malloc (vg_replace_malloc.c:299)
==16702== by 0x5051C3F: ompi_op_create_user (in /usr/lib/x86_64-linux-gnu/openmpi/lib/libmpi.so.40.10.3)
==16702== by 0x5089CE8: PMPI_Op_create (in /usr/lib/x86_64-linux-gnu/openmpi/lib/libmpi.so.40.10.3)
==16702== by 0x1CF0B8: Dune::Generic_MPI_Op<Dune::bigunsignedint<56>, Dune::Max<Dune::bigunsignedint<56> > >::get() (mpicollectivecommunication.hh:45)
==16702== by 0x1C2875: int Dune::CollectiveCommunication<ompi_communicator_t*>::allreduce<Dune::Max<Dune::bigunsignedint<56> >, Dune::bigunsignedint<56> >(Dune::bigunsignedint<56> const*, Dune::bigunsignedint<56>*, int) const (mpicollectivecommunication.hh:329)
==16702== by 0x1B383C: Dune::bigunsignedint<56> Dune::CollectiveCommunication<ompi_communicator_t*>::max<Dune::bigunsignedint<56> >(Dune::bigunsignedint<56> const&) const (mpicollectivecommunication.hh:228)
==16702== by 0x1B0F50: void Dune::fillIndexSetHoles<Dune::Amg::MatrixGraph<Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>, std::allocator<Dune::FieldMatrix<double, 1, 1> > > >, Dune::bigunsignedint<56>, int>(Dune::Amg::MatrixGraph<Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>, std::allocator<Dune::FieldMatrix<double, 1, 1> > > > const&, Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<56>, int>&) (repartition.hh:115)
==16702== by 0x1A8589: bool Dune::graphRepartition<Dune::Amg::MatrixGraph<Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>, std::allocator<Dune::FieldMatrix<double, 1, 1> > > >, Dune::bigunsignedint<56>, int>(Dune::Amg::MatrixGraph<Dune::BCRSMatrix<Dune::FieldMatrix<double, 1, 1>, std::allocator<Dune::FieldMatrix<double, 1, 1> > > > const&, Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<56>, int>&, int, std::shared_ptr<Dune::OwnerOverlapCopyCommunication<Dune::bigunsignedint<56>, int> >&, Dune::RedistributeInterface&, bool) (repartition.hh:1275)
==16702== by 0x1A5995: int testRepart<Dune::FieldMatrix<double, 1, 1> >(int, int) (matrixredisttest.cc:66)
==16702== by 0x19DA57: main (matrixredisttest.cc:138)https://gitlab.dune-project.org/core/dune-istl/-/issues/74GlobalIndex of ParallelIndexSet used out of specification.2019-07-10T12:33:11ZSimon PraetoriusGlobalIndex of ParallelIndexSet used out of specification.The `GlobalIndex` in `ParallelIndexSet` is by documentation a type that "[...] has to provide at least a operator< for sorting". In the documentation of the class it is named a *global id* and it is stated:
```
[...] the global id might...The `GlobalIndex` in `ParallelIndexSet` is by documentation a type that "[...] has to provide at least a operator< for sorting". In the documentation of the class it is named a *global id* and it is stated:
```
[...] the global id might not be consecutive but definitely is persistent.
```
So, it is a persistent (over grid changes) partially ordered id type. (Similar to the global entity id type of the grid)
While this is already not enough for the `ParallelIndexSet` and `GlobalLookupIndexSet` implementation (I think, all comparison operators are used somewhere in the code), it gets worse in dune-istl code that uses the `ParallelIndexSet`.
Example: `fillIndexSetHoles()` from `repartition.hh`:
- a maximal GlobalIndex is determined by first initializing a global index with 0 (not possible in general) and then pairwise `std::max`. Why not using the already ordered list `ParallelIndexSet::localIndices_`?
- to determine the next free global index, the maximum is then incremented (using `operator++`) and shifted (using `operator+`).
Both points are not guaranteed to work by the `ParallelIndexSet`.
So, here some cleanup or at least better specification of the requirements is necessary!https://gitlab.dune-project.org/core/dune-istl/-/issues/42Better documentation for copyCopyToAll::copyCopyToAll2018-07-12T07:47:33ZOliver Sanderoliver.sander@tu-dresden.deBetter documentation for copyCopyToAll::copyCopyToAllThe meaning of copyCopyToAll::copyCopyToAll could benefit from some improvement. Witness Christian's mail:
Dear all,
can someone enlighten me?
What is the meaning of copyCopyToAll::copyCopyToAll(...)?
It seems strange to me. We star...The meaning of copyCopyToAll::copyCopyToAll could benefit from some improvement. Witness Christian's mail:
Dear all,
can someone enlighten me?
What is the meaning of copyCopyToAll::copyCopyToAll(...)?
It seems strange to me. We start with a partition of indices according
to owner and then we might have copies of indices owned by an other
process.
Now, copyCopyToAll is supposed to communicate those indices tagged as
copy to all other processes having this index. But what is this
supposed to do?
a) Assume that the copies don't contain consistent values. Now every
process is sending around his own value. Depending on which message
arrives last the receiving process gets an arbitrary value. This
seems to be wrong and we consequently we should not be allowed to
call copyCopyToAll, if the copies are not consistent.
b) assume we have consistent copies. Why should be communicate at all?
The copies are sent around, but we retrieve the value that we
previously sent to an other process.
It seems I misunderstood something dramatically.
Best
Christian
PS: the two functions
OwnerOverlapCopyCommunication::addOwnerOverlapToAll and
OwnerOverlapCopyCommunication::addOwnerCopyToAll are not used in ISTL
at all. Is anybody using them outside of the core modules?
And Markus' answer:
Hi,
On Wed, Dec 20, 2017 at 09:35:51PM +0100, Christian Engwer wrote:
> can someone enlighten me?
>
> What is the meaning of copyCopyToAll::copyCopyToAll(...)?
>
exactly what the name says it sends the values marked as copy to all
others.
> It seems strange to me. We start with a partition of indices according
> to owner and then we might have copies of indices owned by an other
> process.
>
> Now, copyCopyToAll is supposed to communicate those indices tagged as
> copy to all other processes having this index. But what is this
> supposed to do?
> a) Assume that the copies don't contain consistent values. Now every
> process is sending around his own value. Depending on which message
> arrives last the receiving process gets an arbitrary value. This
> seems to be wrong and we consequently we should not be allowed to
> call copyCopyToAll, if the copies are not consistent.
Well, you are missing some patterns. Just assume the partitioning is
not static, but you repartition from N to n processors (n>N). For an
entity that was previously marked as copy on two adjacent processors
but will only be on one process after the repartitioning you might
need this. Maybe the entity had less visible neighbours on each
process before the repartitioning than before. If you e.g. repartition
a matrix, then the some row entries might even need to up added up.
Concerning "not allowed". This is C++ and you are allowed to shoot
yourself in the food.
> b) assume we have consistent copies. Why should be communicate at all?
> The copies are sent around, but we retrieve the value that we
> previously sent to an other process.
>
Again: Due to repartitioning to fewer processors.
Markus
> It seems I misunderstood something dramatically.
>
> Best
> Christian
>
> PS: the two functions
> OwnerOverlapCopyCommunication::addOwnerOverlapToAll and
> OwnerOverlapCopyCommunication::addOwnerCopyToAll are not used in ISTL
> at all. Is anybody using them outside of the core modules?
>
The were used in downstream modules in the past and if somebody
decides to do additive domain decomposition solvers the might be
again.https://gitlab.dune-project.org/core/dune-istl/-/issues/22Allocator template parameter of classes without memory management2018-05-16T08:58:47ZCarsten Gräsergraeser@math.fau.deAllocator template parameter of classes without memory managementSeveral classes in bvector.hh and basearray.hh have an allocator template parameter but don't do any memory management. This is true for
`base_array_unmanaged`, `base_array_window`, `compressed_base_array_unmanaged`, `block_vector_unma...Several classes in bvector.hh and basearray.hh have an allocator template parameter but don't do any memory management. This is true for
`base_array_unmanaged`, `base_array_window`, `compressed_base_array_unmanaged`, `block_vector_unmanaged`, `BlockVectorWindow`, `compressed_block_vector_unmanaged`, `CompressedBlockVectorWindow`.https://gitlab.dune-project.org/core/dune-istl/-/issues/47Tests should use CMAKE_GUARD to show up as skipped.2018-04-19T12:25:07ZJö Fahlkejorrit@jorrit.deTests should use CMAKE_GUARD to show up as skipped.Many tests in ISTL are inside some condition, e.g.
```cmake
if(cond)
dune_add_test(...)
endif()
```
Instead, they should make use of `CMAKE_GUARD`:
```cmake
dune_add_test(... CMAKE_GUARD cond)
```
This would have the benefit that tho...Many tests in ISTL are inside some condition, e.g.
```cmake
if(cond)
dune_add_test(...)
endif()
```
Instead, they should make use of `CMAKE_GUARD`:
```cmake
dune_add_test(... CMAKE_GUARD cond)
```
This would have the benefit that those tests will show up as skipped so you know that you're not testing everything.https://gitlab.dune-project.org/core/dune-istl/-/issues/30Bug / RFC: Expected behaviour of BCRSMatrix::axpy()2017-11-10T16:01:24ZElias PippingBug / RFC: Expected behaviour of BCRSMatrix::axpy()You can currently put all kinds of matrices into a `BCRSMatrix`. When you're assembling a lumped mass matrix, maybe you'll want to use a `ScaledIdentityMatrix` instead of a `FieldMatrix` -- saves space and matrix-vector multiplication at...You can currently put all kinds of matrices into a `BCRSMatrix`. When you're assembling a lumped mass matrix, maybe you'll want to use a `ScaledIdentityMatrix` instead of a `FieldMatrix` -- saves space and matrix-vector multiplication at the block level will be cheaper.
The `BCRSMatrix` class also has a function `axpy` which does `m1 += c * m2`, i.e. it adds the `c`-fold multiple of `m2` to `m1`. Patterns need to be compatible, naturally, and the way it is implemented, it hands over the work to the `axpy` method of the blocks inside the loop.
Now here's the problem: Some matrix classes do not have an `axpy` function. These include `DiagonalMatrix` and `ScaledIdentityMatrix` (`IdentityMatrix` is already deprecated so I won't bother with pointing out issues that class had). So when you try `m1.axpy(c, m2)` with blocks that are instances of the `ScaledIdentityMatrix` class, you'll get an error at compile-time because there is no `ScaledIdentityMatrix::axpy` while `m1.axpy(c, m2)` will compile and work as expected when the blocks are of type `FieldMatrix` instead.
My question is: Is all of that intentional? Should it be that way? I would expect `BCRSMatrix<ScaledIdentityMatrix>::axpy(scalar, matrix)` to work when matrix is itself of type `BCRSMatrix<ScaledIdentityMatrix>`. Should `BCRSMatrix::axpy` fail differently? Should there be a `ScaledIdentityMatrix::axpy()`?
Here's the code that I used for testing
```c++
#include "config.h"
#include <dune/common/diagonalmatrix.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/dynmatrix.hh>
#include <dune/common/identitymatrix.hh>
#include <dune/istl/bcrsmatrix.hh>
#include <dune/istl/scaledidmatrix.hh>
using namespace Dune;
template <class M>
struct MatrixProvider {
M static get();
};
template <>
struct MatrixProvider<Dune::DynamicMatrix<double>> {
Dune::DynamicMatrix<double> static get() { return {{1, 2}, {3, 4}}; }
};
template <>
struct MatrixProvider<Dune::FieldMatrix<double, 2, 2>> {
Dune::FieldMatrix<double, 2, 2> static get() { return {{1, 5}, {9, 16}}; }
};
template <>
struct MatrixProvider<Dune::IdentityMatrix<double, 2>> {
Dune::IdentityMatrix<double, 2> static get() { return {}; }
};
template <>
struct MatrixProvider<Dune::DiagonalMatrix<double, 2>> {
Dune::DiagonalMatrix<double, 2> static get() { return {25}; }
};
template <>
struct MatrixProvider<Dune::ScaledIdentityMatrix<double, 2>> {
Dune::ScaledIdentityMatrix<double, 2> static get() { return {11}; }
};
template <class T>
void test() {
using M = BCRSMatrix<T>;
T const a = MatrixProvider<T>::get();
M m(2, 2, BCRSMatrix<T>::random);
m.setrowsize(0, 2);
m.setrowsize(1, 2);
m.endrowsizes();
m.addindex(0, 0);
m.addindex(0, 1);
m.addindex(1, 0);
m.addindex(1, 1);
m.endindices();
for (auto &row : m)
for (auto &col : row)
col = a;
M m2 = m;
m2.axpy(2.5, m);
}
int main() {
test<FieldMatrix<double, 2, 2>>();
test<DynamicMatrix<double>>();
//test<DiagonalMatrix<double, 2>>(); // has no axpy
//test<IdentityMatrix<double, 2>>(); // has no blocklevel
//test<ScaledIdentityMatrix<double, 2>>(); // has no axpy
}
```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 suggestion...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/21Unexpected behaviour of BCRSMatrix implict buildmode2017-09-22T04:19:36ZCarsten Gräsergraeser@math.fau.deUnexpected 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 `...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/23Performance of `compressed_base_array_unmanaged::operator[]`2017-09-22T04:19:36ZTobias MalkmusPerformance of `compressed_base_array_unmanaged::operator[]`In one of my applications the operator[] access to the array entries takes up to 5 % of total run-time.
If i comment the index check out i gain this 5% run-time improvement.
I would suggest the two options:
1) use preprocessor dir...In one of my applications the operator[] access to the array entries takes up to 5 % of total run-time.
If i comment the index check out i gain this 5% run-time improvement.
I would suggest the two options:
1) use preprocessor directives to (de-)activate the checks:
```
if (lb == j+n || *lb != i)
DUNE_THROW(ISTLError,"index "<<i<<" not in compressed array");
```
2) implement `at(size_type i)` which includes the index checks and operator[] returns without checking.
The first option seems to be a little more 'dune-istl' like and less invasive, while the second one is similar to the behavior of a `std::vector`.