[petsc-users] Is precondition works for ill-conditioned jacobian matrix

[Gong Ding]

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 <gongding at cn.cogenda.com> 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<petsc-users-bounces at mcs.anl.gov>on behalf of Mark 
>>> Adams<mfadams at lbl.gov>
>>> *Sent:*Thursday, September 14, 2023 5:35 AM
>>> *To:*Gong Ding<gongding at cn.cogenda.com>
>>> *Cc:*petsc-users at mcs.anl.gov<petsc-users at mcs.anl.gov>
>>> *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 <gongding at cn.cogenda.com> 
>>> 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
>>>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20230915/492e9fd0/attachment-0001.html>