<html>
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<style type="text/css" style="display:none;"> P {margin-top:0;margin-bottom:0;} </style>
</head>
<body dir="ltr">
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);" class="elementToProof">
Gong,</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);" class="elementToProof ContentPasted0 ContentPasted1">
Iterative refinement might help, i.e., more than one MatSolve ( LU x = b) is applied. If you use superlu_dist, this operation can be activated at runtime by the option '-pc_type lu -pc_factor_mat_solver_type superlu_dist -mat_superlu_dist_iterrefine' (run your
code withe '-help' for available options).</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);" class="elementToProof ContentPasted0 ContentPasted1">
<br>
</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);" class="elementToProof ContentPasted0 ContentPasted1">
Using mumps, the iterative refinement can be activated by </div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);" class="elementToProof ContentPasted0 ContentPasted1 ContentPasted2">
-mat_mumps_icntl_10: <now 0 : formerly 0>: ICNTL(10): max num of refinements (None)<br>
</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);" class="elementToProof ContentPasted0 ContentPasted1 ContentPasted2">
i.e., '<span style="display: inline !important; background-color: rgb(255, 255, 255);" class="ContentPasted3">-mat_mumps_icntl_10 1' applies one refinement. </span></div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);" class="elementToProof ContentPasted0 ContentPasted1">
<br>
</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);" class="elementToProof ContentPasted0 ContentPasted1">
I suspect the LU solver would be called again if one MatSolve() does not reach the required tolerance. Check your runs with '-snes_monitor -ksp_monitor' to get more understanding.</div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);" class="elementToProof ContentPasted0 ContentPasted1">
Hong</div>
<div id="appendonsend"></div>
<div style="font-family: Calibri, Arial, Helvetica, sans-serif; font-size: 12pt; color: rgb(0, 0, 0);">
<br>
</div>
<hr tabindex="-1" style="display:inline-block; width:98%">
<div id="divRplyFwdMsg" dir="ltr"><font face="Calibri, sans-serif" style="font-size: 11pt; color: rgb(0, 0, 0);"><b>From:</b> petsc-users <petsc-users-bounces@mcs.anl.gov> on behalf of Gong Ding <gongding@cn.cogenda.com><br>
<b>Sent:</b> Thursday, September 14, 2023 10:29 PM<br>
<b>To:</b> Barry Smith <bsmith@petsc.dev><br>
<b>Cc:</b> petsc-users@mcs.anl.gov <petsc-users@mcs.anl.gov><br>
<b>Subject:</b> Re: [petsc-users] Is precondition works for ill-conditioned jacobian matrix</font>
<div> </div>
</div>
<div>
<p>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.</p>
<p>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.</p>
<p>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.
<br>
</p>
<p>The main purpose is, always do LU factorization at 64 bit for performance issue.
<br>
</p>
<p>Method 1 tries to "precondition" a direct solver. I don't know if possible.<br>
</p>
<p>Method 2 wants to use post refinement to improve the accuracy of a direct solver. Theoretically, I think it should work.
<br>
</p>
<p><br>
</p>
<div class="x_moz-cite-prefix">On 2023/9/15 01:37, Barry Smith wrote:<br>
</div>
<blockquote type="cite">
<div><br>
</div>
Method 1 and 2 are unlikely to work.
<div><br>
</div>
<div> 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.</div>
<div><br>
</div>
<div> 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.</div>
<div><br>
</div>
<div> Barry</div>
<div><br>
<div><br>
<blockquote type="cite">
<div>On Sep 14, 2023, at 12:30 PM, Gong Ding <a href="mailto:gongding@cn.cogenda.com" class="x_moz-txt-link-rfc2396E OWAAutoLink" id="OWAc10da529-1841-d5a5-9f25-955328adeeeb">
<gongding@cn.cogenda.com></a> wrote:</div>
<br class="x_Apple-interchange-newline">
<div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
The physical problem itself is ill-conditioned since there are floating regions in the simulation domain.<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
I use MUMPS as 64 bit LU solver, and a special improved SuperLU as 128 bit LU solver (<a href="https://github.com/cogenda/superlu" data-auth="NotApplicable" class="x_moz-txt-link-freetext OWAAutoLink" id="OWA23e575c7-5c53-c66d-e0bb-7839464d376b">https://github.com/cogenda/superlu</a>,
added float128 support).<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
Although 128 bit solver works, it is 10x slower.<span class="x_Apple-converted-space"> </span><br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
I'd like to try, if jacobian can be processed under 64 bit precision while keeps the Newton iteration convergence.<span class="x_Apple-converted-space"> </span><br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
Method 1:<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
Use a block inversion of the main diagonal of jacobian as preconditioner (or ILU? ). Then factorize M*J.<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
Both the precondition matrix and jacobian matrix are 64 bit.<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
Method 2:<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
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)<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
<br>
</div>
<div style="margin-top:0px; margin-bottom:0px; font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
<br>
</div>
<div class="x_moz-cite-prefix" style="font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; text-align:start; text-indent:0px; text-transform:none; white-space:normal; word-spacing:0px; text-decoration:none">
On 2023/9/14 23:05, Zhang, Hong wrote:<br>
</div>
<blockquote type="cite" style="font-family:Helvetica; font-size:18px; font-style:normal; font-variant-caps:normal; font-weight:400; letter-spacing:normal; orphans:auto; text-align:start; text-indent:0px; text-transform:none; white-space:normal; widows:auto; word-spacing:0px; text-decoration:none">
<div class="x_elementToProof" style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">
<span class="x_ContentPasted0" style="color: rgb(36, 36, 36); background-color: rgb(255, 255, 255);">Gong Ding,</span><br>
</div>
<div class="x_elementToProof" style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">
<span class="x_ContentPasted0" style="color: rgb(36, 36, 36); background-color: rgb(255, 255, 255);">When you use a LU solver, the preconditioner M = inv(LU) = inv (J) on theory. I suspect your jacobian evaluation by<span class="x_ContentPasted1" style="display: inline !important; background-color: rgb(255, 255, 255);"> 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.</span></span></div>
<div class="x_elementToProof" style="font-family:Calibri,Arial,Helvetica,sans-serif; font-size:12pt">
<span class="x_ContentPasted0" style="color: rgb(36, 36, 36); background-color: rgb(255, 255, 255);"><span class="x_ContentPasted1" style="display: inline !important; background-color: rgb(255, 255, 255);">Hong</span></span></div>
<hr tabindex="-1" style="display:inline-block; width:925.109375px">
<div id="x_divRplyFwdMsg" dir="ltr"><font face="Calibri, sans-serif" style="font-size:11pt"><b>From:</b><span class="x_Apple-converted-space"> </span>petsc-users<span class="x_Apple-converted-space"> </span><a href="mailto:petsc-users-bounces@mcs.anl.gov" class="x_moz-txt-link-rfc2396E OWAAutoLink" id="OWAa0cbaf7a-6d47-eb66-9865-f1bd1bb4e1ff"><petsc-users-bounces@mcs.anl.gov></a><span class="x_Apple-converted-space"> </span>on
behalf of Mark Adams<span class="x_Apple-converted-space"> </span><a href="mailto:mfadams@lbl.gov" class="x_moz-txt-link-rfc2396E OWAAutoLink" id="OWAe879ac8f-0ec9-39eb-aad6-41e34ced48d1"><mfadams@lbl.gov></a><br>
<b>Sent:</b><span class="x_Apple-converted-space"> </span>Thursday, September 14, 2023 5:35 AM<br>
<b>To:</b><span class="x_Apple-converted-space"> </span>Gong Ding<span class="x_Apple-converted-space"> </span><a href="mailto:gongding@cn.cogenda.com" class="x_moz-txt-link-rfc2396E OWAAutoLink" id="OWA7ae7145e-44f8-f51e-2c85-ae92b1191564"><gongding@cn.cogenda.com></a><br>
<b>Cc:</b><span class="x_Apple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov" class="x_moz-txt-link-abbreviated x_moz-txt-link-freetext OWAAutoLink" id="OWA99056779-0670-673e-19f7-424b2e2fe869">petsc-users@mcs.anl.gov</a><span class="x_Apple-converted-space"> </span><a href="mailto:petsc-users@mcs.anl.gov" class="x_moz-txt-link-rfc2396E OWAAutoLink" id="OWAa2d25aa1-40ed-eb15-adba-bd9f323d7c30"><petsc-users@mcs.anl.gov></a><br>
<b>Subject:</b><span class="x_Apple-converted-space"> </span>Re: [petsc-users] Is precondition works for ill-conditioned jacobian matrix</font>
<div> </div>
</div>
<div>
<div dir="ltr">I would first verify that you are happy with the solution that works.
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div>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.</div>
<div><br>
</div>
<div>Good luck,</div>
<div>Mark</div>
</div>
<br>
<div class="x_x_gmail_quote">
<div dir="ltr" class="x_x_gmail_attr">On Thu, Sep 14, 2023 at 5:35 AM Gong Ding <<a href="mailto:gongding@cn.cogenda.com" class="x_moz-txt-link-freetext OWAAutoLink" id="OWA18b13257-ee69-8cf8-43c8-7b01fe82ed62">gongding@cn.cogenda.com</a>> wrote:<br>
</div>
<blockquote class="x_x_gmail_quote" style="margin:0px
0px 0px 0.8ex; border-left-width:1px; border-left-style:solid; border-left-color:rgb(204,204,204); padding-left:1ex">
Hi all<br>
<br>
I find such a nonlinear problem, the jacobian matrix is ill conditioned.<br>
<br>
Solve the jacobian matrix by 64bit LU solver, the Newton method failed<span class="x_Apple-converted-space"> </span><br>
to convergence.<br>
<br>
However, when solve the jacobian matrix by 128bit LU solver , Newton<span class="x_Apple-converted-space"> </span><br>
iteration will convergence.<br>
<br>
I think this phenomena indicate that , the jacobian matrix is ill<span class="x_Apple-converted-space"> </span><br>
conditioned.<br>
<br>
<br>
The question is, if I do a precondition as M*J*dx = -M*f(x), here M is<span class="x_Apple-converted-space"> </span><br>
the precondition matrix, . then I solve the matrix A=M*J by a LU solver.<br>
<br>
Can I expect that solve A=M*J has a better precision result that help<span class="x_Apple-converted-space"> </span><br>
the convergence of Newton iteration?<br>
<br>
Gong Ding</blockquote>
</div>
</div>
</blockquote>
</div>
</blockquote>
</div>
<br>
</div>
</blockquote>
</div>
</body>
</html>