Andreas, you are right, there seem to be some issues with the pivoting in the LU decomposition. In particular. it would be nice to have full pivoting instead of row pivoting only. Furthermore, the singularity threshold should depend on the data type of the matrix entries (at the moment, it is more or less hard-coded).

But your matrix is very ill-conditioned (condition ~ 10^19) -- it is not a complete surprise it is rejected by the algorithm.

-- Christian and Sven

Sure, I know that it is ill-conditioned and I'm also not surprised. But emerges directly from my physical problem, and I think I should not change the model unnecessarily. What I could do is that I could for example multiply the second row with 10^9, but this would be much more effort than to fix the pivoting (and wouldn't change the conditioning if I'm not mistaken).

To fix the problem properly, I think that the singularity threshold needs to be be dependent on the largest entry that occurs in a column during the decomposition. I just have trouble finding a efficient way to implement this. In the mean time I would like to propose the attached patch to fix the problem, if compiled with -DNDEBUG the singularity detection just breaks down to a comparison with 0, else the singularity thershold is column specific.

oops, the patch doesn't apply properly because I had some local modifications in my git repo. here's a clean patch.

Andreas, you can just change the threshold to 0 in your application.

0) just call Dune::FMatrixPrecision<>::set_singular_limit(0.0);

1) Rescaling is the way to go if the matrix is like in your case, although you should better scale the last line by 1e-9.

2) Full pivoting would be the proper fix, but your suggested fix isn't a fix at all.

If you just force dune to invert your matrix, you A*A^-1 is a lot worse compared to the scaled version (see code attached):

============orig=============
1 0 0 6.53252e-21
9.04729e-24 1 4.69256e-19 -1.55248e-27
0 0 1 0
-3.41061e-13 0 0 1

============scaled=============
1 0 0 0
2.67283e-23 1 1.07997e-19 -1.70084e-18
0 0 1 0
1.12579e-21 7.75205e-18 7.72219e-18 1

Ok, thanks for the hints. I will check whether I implement the rescaling, but this might involve quite a bit of work. (the problem is that this small matrix is on the main diagonal of a BCRS matrix).