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.
The code that begins
l = ipvt[k - 1];
if (l != k) {
ax = &a[k3 + 1];
ay = &a[6 * l + 1];
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.
Barry
> On Jun 5, 2025, at 7:49 AM, Matthew Knepley <knep...@gmail.com> wrote:
>
> On Thu, Jun 5, 2025 at 2:27 AM Frank Bramkamp <bramk...@nsc.liu.se
> <mailto:bramk...@nsc.liu.se>> wrote:
>> Thanks for the response,
>>
>>
>> I refer to the routine src/mat/impls/baij/seq/dgefa6.c
>>
>>
>> /*
>> Inverts 6 by 6 matrix using gaussian elimination with partial pivoting.
>>
>> Used by the sparse factorization routines in
>> src/mat/impls/baij/seq
>>
>> This is a combination of the Linpack routines
>> dgefa() and dgedi() specialized for a size of 6.
>>
>> */
>>
>>
>> which implements the the computation of an inverse for a dense 6x6 block.
>> I think this routine is used in the ILU numeric factorization to compute the
>> inverse diagonal block
>> when we later apply ILU preconditioning where we need the inverse.
>>
>> And there it looks we use a direct LU solve and not triangular solves ?!
>>
>> In this context I am not sure if something is missing at the end or not
>> regarding the inverse pivoting ?!
>
> We do not think so, but there is a simple test for this. You can construct a
> small block diagonal matrix with
> blocksize 6, and apply it to the vector [1, 2, 3, 4, 5, 6 ...]. Use that
> vector as the RHS for a solve with
> preonly/ILU. That should get the exact inverse and output the original
> vector. If the vector is permuted, then
> we indeed have a bug.
>
> Thanks!
>
> Matt
>
>> Greetings, Frank
>>
>>
>>
>>
>
>
>
> --
> What most experimenters take for granted before they begin their experiments
> is infinitely more interesting than any results to which their experiments
> lead.
> -- Norbert Wiener
>
> https://urldefense.us/v3/__https://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!bMXreiYRM3eT86Jioq3IATcR2EVaQfLh4ml3pmh88QtCtsm-80eT25djN21OkkgGiBm1v8KIhRDmklo_cvaliCI$
>
> <https://urldefense.us/v3/__http://www.cse.buffalo.edu/*knepley/__;fg!!G_uCfscf7eWS!d8vX5tFZOpXshlbVrdeXXKCaT9OQePsnBaR3YBpe5DEoej9o3r1ZB9Rs-RgRm4xVgXdw-FNm8ShhNTEksffZ$>
