Thank you for the code. A) By default MatZeroRows() does change the nonzero structure of the matrix B) PCFIELDSPLIT loses the null spaces attached to the submatrices if the nonzero structure of the matrix changes. For the example code if one sets MatSetOption(A,MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); then using MatZeroRows() no longer changes the nonzero structure and thus the code behaves similarly for both MatZeroRows() and MatZeroRowsColumns(). Should B lose the null spaces when the nonzero structure changes? The documentation for MatSetNullSpace() does not explicitly state that the attached null space will remain with the matrix for its entire life unless the user calls MatSetNullspace again, but it implicitly implies that. Thus, I think the B should not lose the attached null spaces. I will change the code to preserve the null spaces in PCFIELDSPLIT when the nonzero structure changes. Going back to your two points 1. Code going through line 692 loses the near nullspace of the matrices attached to the sub-KSPs 2. The call to MatZeroRowsColumns() changes then non-zero structure for MPIAIJ but not for SEQAIJ My new branch will prevent the loss of the nullspace in 1. My previous fix (that is now in main) fixes the bug where MatZeroRowsColumns() in parallel thought it changed the nonzero structure (while it actually did not change the nonzero structure). Note: Independent of 1 and 2 most likely when you use MatZeroRows() in your code as you describe it you will want to use MatSetOption(A,MAT_KEEP_NONZERO_PATTERN, PETSC_TRUE); since the code will likely be a bit more efficient and will have less memory churn.
On Feb 13, 2024, at 7:12 AM, Jeremy Theler (External) <[email protected]> wrote:
Hi Barry
7279 does change the code for MatZeroRowsColumns_MPIAIJ(). But perhaps that does not resolve the problem you are seeing? If that is the case we will need a reproducible example so we can determine exactly what else is happening in your code to cause the difficulties.
Here is the diff for MatZeroRowsColumns_MPIAIJ()
@@ -1026,7 +1023,7 @@ static PetscErrorCode MatZeroRowsColumns_MPIAIJ(Mat A, PetscInt N, const PetscIn
PetscCall(PetscFree(lrows));
/* only change matrix nonzero state if pattern was allowed to be changed */ - if (!((Mat_SeqAIJ *)(l->A->data))->keepnonzeropattern) { + if (!((Mat_SeqAIJ *)(l->A->data))->nonew) { PetscObjectState state = l->A->nonzerostate + l->B->nonzerostate; PetscCall(MPIU_Allreduce(&state, &A->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)A))); }
Fair enough, I might have overlooked this single line. But this change does not fix our issue. The if condition is still true. I'll need some time to come up with a reproducible example because it involves setting up the PC of the KSP of an SNES of a TS before doing TSSolve() with the DM-allocated jacobian.
However, regarding bug #1, I do have a MWE. It is not CI-friendly. I tried to modify ex47.c to illustrate the issue but could not make it work. Anyway, the attached tarball: reads a global matrix A from a file. reads two ISs from another file sets up a field split PC attaches a near nullspace (read from another file) to one of the sub-KSPs sets "dirichlet BCs" with MatZeroRowsColumns() or MatZeroRows() depending if -columns was given in the command line
The issue comes when calling MatZeroRowsColumns() or MatZeroRows(). In the first case, the near nullspace is not lost but it is in the second:
$ make lost $ ./lost -columns -ksp_view | grep "near null" has attached near null space has attached near null space $ mpiexec -n 2 ./lost -columns -ksp_view | grep "near null" has attached near null space has attached near null space $ ./lost -ksp_view | grep "near null" $ mpiexec -n 2 ./lost -ksp_view | grep "near null" $
When using MatZeroRows(), the code passes through fieldsplit.c:692 which looses the near nullspace.
Note that the original issue we see in our code is that there is a difference between the serial and parallel implementation of MatZeroRowsColumns(). We loose the near nullspace only in parallel but not in serial because of a combination of bugs #1 and #2.
-- jeremy
<lost.tar.gz>