[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.jamond at cea.fr Phone: +336.78.18.18.25


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



On 27 Jun 2025, at 3:21 PM, Barry Smith <bsmith at 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 <pierre at joliv.et> wrote:



On 27 Jun 2025, at 1:33 PM, Barry Smith <bsmith at 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.JAMOND at 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.jamond at cea.fr>@cea.fr<mailto:olivier.jamond at cea.fr> Phone: +336.78.18.18.25


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