[Dune] [#778] MatrixHelper::detAAT is unstable

[Dune]

THIS IS AN AUTOMATED MESSAGE, DO NOT REPLY.

The following task has a new comment added:

FS#778 - MatrixHelper::detAAT is unstable
User who did this - Martin Nolte (nolte)

----------
I think I don't get your exact point. The determinant in your example is very nearly zero:

gnuplot> print 0.099999999999999867 * (-0.0099999999999998979) - (-0.010000000000002118) * 0.099999999999999867
2.22044604925031e-16

I have to admit, though, that the above formula should be numerically stable in this context.

Now, detAAT first (A * A^T) and then its determinant, i.e.,  the determinant of (computed by gnuplot)
AAT[0][0] = 0.0101
AAT[0][1] = 0.0101
AAT[1][0] = 0.0101
AAT[1][1] = 0.0101
The determinant of this matrix is obviously zero.

Moreover, the matrix has one eigenvalue pretty close to 0 and one eigenvalue around 0.01, i.e., it is quite ill-conditioned. I wonder if the inverse can be computed numerically (with double precision)?

Since this MatrixHelper shall only be used for the Jacobians of geometries, the following questions arise:
(a) is there an algorithm that is numerically more stable (and works in the nonsquare case)?
(b) do we really need to consider geometries where the jacobian is this ill-conditioned?

If someone really must use such ill-conditioned geometries, I think he is well-advised to use some higher precision floating point type for the geometry's ctype.
----------

More information can be found at the following URL:
http://www.dune-project.org/flyspray/index.php?do=details&task_id=778#comment1946

You are receiving this message because you have requested it from the Flyspray bugtracking system.  If you did not expect this message or don't want to receive mails in future, you can change your notification settings at the URL shown above.