#778 MatrixHelper::detAAT is unstable
Metadata
| Property | Value |
|---|---|
| Reported by | Oliver Sander (oliver.sander@tu-dresden.de) |
| Reported at | Apr 21, 2010 20:53 |
| Type | Bug Report |
| Version | Git (pre2.4) [autotools] |
| Operating System | Unspecified / All |
| Last edited by | Christian Engwer (christi@conan.iwr.uni-heidelberg.de) |
| Last edited at | Nov 1, 2011 17:10 |
Description
For certain almost singular matrices, for example
A[0][0] = 0.099999999999999867;
A[0][1] = -0.010000000000002118;
A[1][0] = 0.099999999999999867;
A[1][1] = -0.0099999999999998979;the method detAAT in genericgeometries/matrix.hh aborts with
testmatrix: ../matrix.hh:243: static void Dune::GenericGeometry::MatrixHelper<Traits>::cholesky_L(const typename Traits::Matrix<n, n>::type&, typename Traits::Matrix<n, n>::type&) [with int n = 2, Traits = Dune::GenericGeometry::DuneCoordTraits<double>]: Assertion `x > FieldType( 0 )' failed.rather than computing something close to zero.
This is not an artificial example. I found it when integrating over remote intersections between two unrelated grids.
A test case has been added in the file genericgeometries/test/testmatrix.cc.