[petsc-users] Efficient handling of missing diagonal entities

[JAMOND Olivier]

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 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/20250627/9083eec8/attachment.html>