[petsc-users] BLOCK INVERSE FOR ILU

[Barry Smith]

   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 <knepley at gmail.com> wrote:
> 
> On Thu, Jun 5, 2025 at 2:27 AM Frank Bramkamp <bramkamp at nsc.liu.se <mailto:bramkamp at 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$>

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20250605/5fe55829/attachment-0001.html>