The matrix has a bad condition number, but not singular. It comes from
real physical problem and the floating zone is weekly controlled by
remote boundary condition.
Yes, I am also afraid that with 64 bit floating number, the matrix is
numerically singular since the construction of jacaobian has already
lost precision.
Anyway, I can build the jacobian at 128 bit precision and then truncate
to 64 bit. Our solution x and function f can also be evaluated at 128
bit precision.
The main purpose is, always do LU factorization at 64 bit for
performance issue.
Method 1 tries to "precondition" a direct solver. I don't know if possible.
Method 2 wants to use post refinement to improve the accuracy of a
direct solver. Theoretically, I think it should work.
On 2023/9/15 01:37, Barry Smith wrote:
Method 1 and 2 are unlikely to work.
It sounds like your matrix is (in exact precision) singular, but
using 128 bit floats allows a "stable" factorization to go through
giving you a descent direction for Newton.
I think you really need to fix the singularity at the modeling
level, it is not robust to fix it at the numerical algorithm level. If
you know the exact form of the null spaces you can use
MatSetNullSpace() but you cannot use a direct solver for the system
since the factorization will still see the singular matrix.
Barry
On Sep 14, 2023, at 12:30 PM, Gong Ding <[email protected]> wrote:
The physical problem itself is ill-conditioned since there are
floating regions in the simulation domain.
I use MUMPS as 64 bit LU solver, and a special improved SuperLU as
128 bit LU solver (https://github.com/cogenda/superlu, added float128
support).
Although 128 bit solver works, it is 10x slower.
I'd like to try, if jacobian can be processed under 64 bit precision
while keeps the Newton iteration convergence.
Method 1:
Use a block inversion of the main diagonal of jacobian as
preconditioner (or ILU? ). Then factorize M*J.
Both the precondition matrix and jacobian matrix are 64 bit.
Method 2:
Do a 64 bit LU factorization of jacobian matrix, and use the
factorization result as a preconditioner for higher precision krylov
solver (such as iterative refinement)
On 2023/9/14 23:05, Zhang, Hong wrote:
Gong Ding,
When you use a LU solver, the preconditioner M = inv(LU) = inv (J)
on theory. I suspect your jacobian evaluation by 64bit might be
inaccurate. What LU solver did you use? Run your code with option
'-snes_view -snes_monitor -ksp_monitor' and compare the displays.
Hong
------------------------------------------------------------------------
*From:*petsc-users<[email protected]>on behalf of Mark
Adams<[email protected]>
*Sent:*Thursday, September 14, 2023 5:35 AM
*To:*Gong Ding<[email protected]>
*Cc:*[email protected]<[email protected]>
*Subject:*Re: [petsc-users] Is precondition works for
ill-conditioned jacobian matrix
I would first verify that you are happy with the solution that works.
Next, I would worry about losing accuracy in computing M*J, but you
could try it and search for any related work. There may be some tricks.
And MUMPS is good at high accuracy, you might try that and if it
fails look at the MUMPS docs for any flags for high-accuracy.
Good luck,
Mark
On Thu, Sep 14, 2023 at 5:35 AM Gong Ding <[email protected]>
wrote:
Hi all
I find such a nonlinear problem, the jacobian matrix is ill
conditioned.
Solve the jacobian matrix by 64bit LU solver, the Newton method
failed
to convergence.
However, when solve the jacobian matrix by 128bit LU solver , Newton
iteration will convergence.
I think this phenomena indicate that , the jacobian matrix is ill
conditioned.
The question is, if I do a precondition as M*J*dx = -M*f(x),
here M is
the precondition matrix, . then I solve the matrix A=M*J by a LU
solver.
Can I expect that solve A=M*J has a better precision result that
help
the convergence of Newton iteration?
Gong Ding
