solution of linear system using LAPACKFullMatrix - problems with singular matrix

[Simon]

Dear all,

I have to deal with a large number of very small linear systems A*a=b; I need to solve this system for each vertex of my triangulation.

The matrix A is purely determinand by dyadic products of certain quadrature points.

Example: Point<dim> qp =[x,y]; Vector<double> qp_ = [1, x, y, xy, x^2, y^2]^T

-> A = A + qp_*qp_^T

The number of qps is restricted to the qps of those cells which are in the patch. For instance there are 16 qps when using Q2 elements in 2d with reduced integration (QGauss<dim>(2)).

First I decided to actually invert the matrix A,

i.e. A.invert(A)....,

but the computation of the inverse introduced too much numerical errors. Hence I decided to switch to LAPACKFullMatrix, computed LU factorization and finally solved the system via the .solve() member functions.

Switching to LAPACK made my results much better but I still reach a point (fine meshes) where the LU decomposition fails. For nearly each vertex the reciprocal condition number reaches the machine accuracy.

I guess the problem is that fine meshes lead to large differences in the x- and y-coordinates of  the qps (Point<dim> dummy = [x,y] = [10.0 ,0.01]) so that some entries in the A matrix are close to zero and others in turn become very large. IMHO the LU decomposition becomes more demanding if the entries in A differ by several orders. Here is a typical example for A (lowest order e-01 , highest order e+02):

1.600e+01  6.099e+01  3.798e+00  2.327e+02  1.203e+00  1.449e+01
 6.099e+01  2.327e+02  1.449e+01  8.884e+02  4.594e+00  5.532e+01
 3.798e+00  1.449e+01  1.203e+00  5.532e+01  4.292e-01  4.594e+00
 2.327e+02  8.884e+02  5.532e+01  3.395e+03  1.755e+01  2.114e+02
 1.203e+00  4.594e+00  4.292e-01  1.755e+01  1.632e-01  1.640e+00
 1.449e+01  5.532e+01  4.594e+00  2.114e+02  1.640e+00  1.755e+01

My question is if I can give some "massage" to A before calling the LU decomposition? What I mean is something like a preconditioner which are often used when talking about iterative solvers. But when checking the documentation I couldn't find certain functionalities for the Matrix classes.

For me it would be great to know if the A matrix could be "rescued" somehow by certain solution techniques, preconditioning or whatever. Otherwise I have to think about the way of constructing A itself.

Thanks for giving me feedback!

Best

Simon

[Bruno Turcksin]

Simon,

Unfortunately, there is no such thing as a preconditioner for direct solver. The best way to fix your problem is to try to scale x and y so that they have the same order of magnitude. Another way is to try SVD. You may need to check that you get an accurate solution for the finest mesh.

Best,

Bruno

[Simon Wiesheier]

Thank you for your reply!

I scaled the x and y coordinates so that they are both between -1 and 1. Actually this brought the condition numbers down to 1e-04.

I am just a little bit surprised that my results are nearly the same as those obtained with reciprocal conditon numbers smaller than my machine accuracy, i.e. about 1e-18. 

Would you say that one can continue work with reciprocal condition numbers of 1e-04 or is this still too big? 

But anyway, glad to see that my program is no longer aborting because LU fails. :-)

Best

Simon 

--
The deal.II project is located at http://www.dealii.org/
For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
---
You received this message because you are subscribed to a topic in the Google Groups "deal.II User Group" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/dealii/IEhvdDouaBc/unsubscribe.
To unsubscribe from this group and all its topics, send an email to dealii+un...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a67c534a-c7da-472a-99fa-72d02990c10en%40googlegroups.com.

[Wolfgang Bangerth]

On 7/2/21 10:44 AM, Simon Wiesheier wrote:
>
> I scaled the x and y coordinates so that they are both between -1 and 1.
> Actually this brought the condition numbers down to 1e-04.
> I am just a little bit surprised that my results are nearly the same as those
> obtained with reciprocal conditon numbers smaller than my machine accuracy,
> i.e. about 1e-18.
>
> Would you say that one can continue work with reciprocal condition numbers of
> 1e-04 or is this still too big?

A condition number of 1e4 is quite moderate and should not cause any undue
problems.

The monomial basis you chose is known to be quite poorly conditioned. As Bruno
already remarked, a better conditioned basis is to use
1, (x-x0)/dx, (y-y0)/dy, (x-x0)^2/dx^2, ...
where x0,y0 are the center of the cell and dx, dy are the extents of the cell
in x and y direction (or you can choose dx=dy=h if your cells are all
reasonably well behaved). This basis is appropriate if you only need linear
and quadratic terms, but not if you go to higher polynomial degrees. In that
case, it is useful to go with something like Chebyshev or Legendre
polynomials. That's what all higher order finite element implementations are
based on.

Best
W.

--
------------------------------------------------------------------------
Wolfgang Bangerth email: bang...@colostate.edu
www: http://www.math.colostate.edu/~bangerth/