<html><head><meta http-equiv="content-type" content="text/html; charset=utf-8"></head><body style="overflow-wrap: break-word; -webkit-nbsp-mode: space; line-break: after-white-space;"><div><br></div> I do urge you to run a test as Matt suggests to detect possible bugs in our code and satisfy yourself the code is correct or incorrect.<div><br></div><div> The code that begins</div><div><br></div><div><div> l = ipvt[k - 1];</div><div> if (l != k) {</div><div> ax = &a[k3 + 1];</div><div> ay = &a[6 * l + 1];</div><div><br></div><div>is attempting to take into account the pivoting. This chunk of code may be incorrect, of course. Note that there is similar code for block size 5 etc that may also be buggy I guess.</div><div><br></div><div> Barry</div><div><br></div><div><br><blockquote type="cite"><div>On Jun 5, 2025, at 7:49 AM, Matthew Knepley <knepley@gmail.com> wrote:</div><br class="Apple-interchange-newline"><div><div dir="ltr"><div dir="ltr">On Thu, Jun 5, 2025 at 2:27 AM Frank Bramkamp <<a href="mailto:bramkamp@nsc.liu.se">bramkamp@nsc.liu.se</a>> wrote:</div><div class="gmail_quote gmail_quote_container"><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">Thanks for the response,<br>
<br>
<br>
I refer to the routine src/mat/impls/baij/seq/dgefa6.c<br>
<br>
<br>
/*<br>
Inverts 6 by 6 matrix using gaussian elimination with partial pivoting.<br>
<br>
Used by the sparse factorization routines in<br>
src/mat/impls/baij/seq<br>
<br>
This is a combination of the Linpack routines<br>
dgefa() and dgedi() specialized for a size of 6.<br>
<br>
*/<br>
<br>
<br>
which implements the the computation of an inverse for a dense 6x6 block.<br>
I think this routine is used in the ILU numeric factorization to compute the inverse diagonal block<br>
when we later apply ILU preconditioning where we need the inverse.<br>
<br>
And there it looks we use a direct LU solve and not triangular solves ?!<br>
<br>
In this context I am not sure if something is missing at the end or not regarding the inverse pivoting ?!<br></blockquote><div><br></div><div>We do not think so, but there is a simple test for this. You can construct a small block diagonal matrix with</div><div>blocksize 6, and apply it to the vector [1, 2, 3, 4, 5, 6 ...]. Use that vector as the RHS for a solve with</div><div>preonly/ILU. That should get the exact inverse and output the original vector. If the vector is permuted, then</div><div>we indeed have a bug.</div><div><br></div><div> Thanks!</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
Greetings, Frank<br>
<br>
<br>
<br>
<br>
</blockquote></div><div><br clear="all"></div><div><br></div><span class="gmail_signature_prefix">-- </span><br><div dir="ltr" class="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d8vX5tFZOpXshlbVrdeXXKCaT9OQePsnBaR3YBpe5DEoej9o3r1ZB9Rs-RgRm4xVgXdw-FNm8ShhNTEksffZ$" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div></div></div>
</div></blockquote></div><br></div></body></html>