Re: [petsc-users] Efficient handling of missing diagonal entities

[JAMOND Olivier]

Hi Pierre and Barry,


Thanks for your answers and sorry for my reaction time (a bit overwhelmed...).


Yes I confirm that I already tried `MAT_FORCE_DIAGONAL_ENTRIES` following a 
Pierre suggestion, but it did not help...


I understand from your messages that this option indeed needs a fix. Do you 
think that such a fix for this option could be envisaged in a next future in 
petsc's roadmap?


Many thanks,

_________________________________________
Olivier Jamond
Research Engineer
French Atomic Energy and Alternative Energies Commission
DES/ISAS/DM2S/SEMT/DYN
91191 Gif sur Yvette, Cedex, France

Email: olivier.jam...@cea.fr Phone: +336.78.18.18.25


________________________________
From: petsc-users <petsc-users-boun...@mcs.anl.gov> on behalf of Pierre Jolivet 
<pie...@joliv.et>
Sent: Friday, June 27, 2025 3:50:29 PM
To: Barry Smith
Cc: petsc-users@mcs.anl.gov
Subject: Re: [petsc-users] Efficient handling of missing diagonal entities



On 27 Jun 2025, at 3:21 PM, Barry Smith <bsm...@petsc.dev> wrote:



  Because I completely forgot that this option existed, and the LLM didn't save 
me from embarrassing myself.

  I see that this option sets mat->force_diagonals, but this variable is never 
used in the mat assembly routines, meaning it will not help in this situation.

  Presumably, MatAssemblyXXX_YYY() could/should be fixed to respect this flag?

Yes, Olivier asked me the same question previously, I told him that this option 
should probably be revamped because it’s there but I don’t think it’s doing its 
job.

Thanks,
Pierre

Then it would help the Olivier.

  Barry




On Jun 27, 2025, at 7:49 AM, Pierre Jolivet <pie...@joliv.et> wrote:



On 27 Jun 2025, at 1:33 PM, Barry Smith <bsm...@petsc.dev> wrote:



   Handling empty diagonal entries on matrices is often problematic, just as 
you describe.

   I suggest placing explicit zeros on the diagonal first before providing the 
other entries, which might be the cleanest and most efficient approach. So have 
each MPI rank loop over its local rows and call MatSetValue() for each diagonal 
entry and then continue with your other MatSetValues(). Do not call 
MatAssemblyBegin/End() after you have provided the zeros on the diagonal just 
chug straight into setting the other values.

   Barry

  As you observed, trying to add the zero entries in the matrix after it is 
assembled is terribly inefficient and not the way to go.

   I've considered adding a matrix option to force zero entries on the 
diagonal, but I never completed my consideration. For example,   
MatSetOption(A, MAT_NONEMPTY_DIAGONAL,PETSC_TRUE);

Why would you need another option when there is already 
MAT_FORCE_DIAGONAL_ENTRIES?

Thanks,
Pierre

and when this option is set, MatAssemblyBegin fills up any empty diagonal 
entries automatically.



On Jun 27, 2025, at 6:26 AM, JAMOND Olivier <olivier.jam...@cea.fr> wrote:


Hello,

I am working on a PDE solver which uses petsc to solve its sparse distributed 
linear systems. I am mainly dealing with MPIAIJ matrices.

In some situations, it may happen that the matrices considered does not have 
non-zero term on the diagonal. For instance I work on a case which have a 
stokes like saddle-point structure (in a MPIAIJ, not a MATNEST):

[A Bt][U]=[F]
[B 0 ][L] [0]

I do not insert null terms in the zero block.

In some cases, I use the function `MatZeroRowsColumns` to handle "Dirichlet" 
boundary conditions. In this particular case, I apply Dirichlet BCs only on 
dofs of "U". But I get an error `Matrix is missing diagonal entry in row X` 
from the function `MatZeroRowsColumns`, where X is a row related to "L".

My first question is: is it normal that I get an error for a missing diagonal 
in the function `MatZeroRowsColumns`entry for a dof that is not involved in the 
list of dofs that I pass to `MatZeroRowsColumns`?

I then tried to make my code to detect that there are some missing diagonal 
entries, and add an explicit zero to them. My code which adds the missing 
diagonal entries looks like what follows. This is certainly not the best way to 
do that, as in my test case about ~80% of the total computation time is spent 
in this piece of code (more precisely in `MatSetValue(D, k, k, 0., 
ADD_VALUES)`).
So my second question is: what would be the most efficient way to detect the 
missing diagonal entries, and ad explicit zeros on the diagonal at these places?

Many thanks,
Olivier

    ...

    MatAssemblyBegin(A, MAT_FINAL_ASSEMBLY);
    MatAssemblyEnd(A, MAT_FINAL_ASSEMBLY);


    Mat D;
    MatGetDiagonalBlock(A, &D);

    PetscBool missing;
    MatMissingDiagonal(D, &missing, NULL);

    if (missing) {

      IS missingDiagEntryRows;
      MatFindZeroDiagonals(D, &missingDiagEntryRows)

      PetscInt size;
      ISGetLocalSize(missingDiagEntryRows, &size);
      const PetscInt *ptr;
      ISGetIndices(missingDiagEntryRows, &ptr);

      for (Index i = 0; i < size; ++i) {
        PetscInt k = ptr[i];
        MatSetValue(D, k, k, 0., ADD_VALUES);
      }
      MatAssemblyBegin(D, MAT_FINAL_ASSEMBLY);
      MatAssemblyEnd(D, MAT_FINAL_ASSEMBLY);

      ISRestoreIndices(missingDiagEntryRows, &ptr);
    }



_________________________________________
Olivier Jamond
Research Engineer
French Atomic Energy and Alternative Energies Commission
DES/ISAS/DM2S/SEMT/DYN
91191 Gif sur Yvette, Cedex, France

Email: 
olivier.jamond<mailto:olivier.jam...@cea.fr>@cea.fr<mailto:olivier.jam...@cea.fr>
 Phone: +336.78.18.18.25