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 <https://dealii.org/current/doxygen/deal.II/classLAPACKFullMatrix.html#aa855f1fb4ccc0d3d263382d8545bbc6a>. You may need to check that you get an accurate solution for the finest mesh.
Best, Bruno On Tuesday, June 29, 2021 at 6:23:32 AM UTC-4 Simon wrote: > 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 > -- 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 the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to [email protected]. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/a67c534a-c7da-472a-99fa-72d02990c10en%40googlegroups.com.
