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: [email protected] Phone: +336.78.18.18.25
