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); 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
