Near nullspace lost in fieldsplit
Hi all Because of a combination of settings, our code passes through this line: https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/field... i.e. the matrices associated with each of the sub-KSPs of a fieldsplit are destroyed and then re-created later. The thing is that one of these destroyed matrices had a near nullspace attached, which is lost because the new matrix does not have it anymore. Is this a bug or are we missing something? I'll need some time to come up with a MWE showing how the near nullspace is lost, but if its needed I can take a stab at it. -- jeremy [https://gitlab.com/uploads/-/system/project/avatar/13882401/PETSc_RBG-logo.png]<https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/fieldsplit.c?ref_type=heads#L692> src/ksp/pc/impls/fieldsplit/fieldsplit.c · main · PETSc / petsc · GitLab<https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/fieldsplit.c?ref_type=heads#L692> PETSc, pronounced PET-see (the S is silent), is a suite of data structures and routines for the scalable (parallel) solution of scientific applications modeled by partial differential equations. gitlab.com
On Fri, Feb 9, 2024 at 7:04 AM Jeremy Theler (External) < [email protected]> wrote:
Hi all
Because of a combination of settings, our code passes through this line:
https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/field...
i.e. the matrices associated with each of the sub-KSPs of a fieldsplit are destroyed and then re-created later. The thing is that one of these destroyed matrices had a near nullspace attached, which is lost because the new matrix does not have it anymore.
Is this a bug or are we missing something?
I just want to get a clear picture. You create a PCFIELDSPLIT, set it up, then pull out the matrices and attach a nullspace before the solve. At a later time, you start another solve with this PC, and it has the DIFFERENT_NONZERO_PATTERN flag, so it recreates these matrices and loses your attached nullspace. First, does the matrix really change? Second, I had the same problem. I added callbacks to DM which allow a nullspace to be automatically attached if you extract a certain subfield. Are you using a DM? Thanks, Matt
I'll need some time to come up with a MWE showing how the near nullspace is lost, but if its needed I can take a stab at it.
-- jeremy
<https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/fieldsplit.c?ref_type=heads#L692> src/ksp/pc/impls/fieldsplit/fieldsplit.c · main · PETSc / petsc · GitLab <https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/fieldsplit.c?ref_type=heads#L692> PETSc, pronounced PET-see (the S is silent), is a suite of data structures and routines for the scalable (parallel) solution of scientific applications modeled by partial differential equations. gitlab.com
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
Because of a combination of settings, our code passes through this line:
https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/field...
i.e. the matrices associated with each of the sub-KSPs of a fieldsplit are destroyed and then re-created later. The thing is that one of these destroyed matrices had a near nullspace attached, which is lost because the new matrix does > not have it anymore.
Is this a bug or are we missing something?
I just want to get a clear picture. You create a PCFIELDSPLIT, set it up, then pull out the matrices and attach a nullspace before the solve.
We need to solve an SNES. We use dmplex so we have the jacobian allocated before starting the solve. At setup time we 1. define the PC of the KSP of the SNES to be fieldsplit 2. define the fields with ISes 3. call PCSetup() to create the sub-KSPs 4. retrieve the matrix attached to the sub-KSP that needs the near nullspace and attach it to that matrix
At a later time, you start another solve with this PC, and it has the DIFFERENT_NONZERO_PATTERN flag, so it recreates these matrices and loses your attached nullspace.
At a later time, in the jacobian evaluation we populate the global matrix (i.e. not the matrices attached to each sub-KSPs) and then we set dirichlet bcs with MatZeroRowsColumns() on that same global matrix. For some reason, in serial the near nullspace is not lost but in parallel the call to MatZeroRowsColumns() does change the non-zero structure (even though the manual says it does not) and then the code goes through that line 692 in fieldsplit.c and the near nullspace is lost.
First, does the matrix really change?
Well, the matrix during setup is not filled in, just allocated. The thing is that if we set MAT_KEEP_NONZERO_PATTERN to true with MatSetOption() before setting the dirichlet BCs, then the near nullspace is not lost (because the code does not go through line 692 of fieldsplit.c). So there are (at least) two issues: 1. Code going through line 692 looses 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
Second, I had the same problem. I added callbacks to DM which allow a nullspace to be automatically attached if you extract a certain subfield. Are you using a DM?
Yes. Can you give us an example? Regards -- jeremy
1. Code going through line 692 looses 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 (unless MAT_KEEP_NONZERO_PATTERN is used) MatZeroRowsColumns() manual page states: Unlike `MatZeroRows()` this does not change the nonzero structure of the matrix, it merely zeros those entries in the matrix. MatZeroRowsColumns_MPIAIJ() has the code /* only change matrix nonzero state if pattern was allowed to be changed */ if (!((Mat_SeqAIJ *)(l->A->data))->keepnonzeropattern) { PetscObjectState state = l->A->nonzerostate + l->B->nonzerostate; PetscCall(MPIU_Allreduce(&state, &A->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)A))); } The if() test is simply an optimization to avoid the reduction if keepnonzeropattern is true. In your first run (when MAT_KEEP_NONZERO_PATTERN is not used) the keepnonzeropattern is not set so the A matrix nonzerostate is updated using the nonzerostate of the two submatrices. But that A nonzero state will only change if one of the two submatrices nonzerostate changes. In the second run the A nonzerostate cannot be changed. MatZeroRowsColumns_SeqAIJ() has the code if (diag != 0.0) { PetscCall(MatMissingDiagonal_SeqAIJ(A, &missing, &d)); if (missing) { for (i = 0; i < N; i++) { if (rows[i] >= A->cmap->N) continue; PetscCheck(!a->nonew || rows[i] < d, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry in row %" PetscInt_FMT " (%" PetscInt_FMT ")", d, rows[i]); PetscCall(MatSetValues_SeqAIJ(A, 1, &rows[i], 1, &rows[i], &diag, INSERT_VALUES)); which adds a nonzero to the matrix to fill in a missing diagonal. This will change the nonzerostate of the matrix, but in a different way than keepnonzeropattern flag is for. So, do you have missing diagonals in your matrix? Bug 1 - documentation The documentation statement "Unlike `MatZeroRows()` this does not change the nonzero structure of the matrix, it merely zeros those entries in the matrix." appears to be incorrect for matrices missing diagonal entries since new nonzeros are added (to fill in the diagonal). The documentation should really say "Unlike `MatZeroRows()` this routine cannot remove the zeroed entries from the nonzero structure of the matrix; in other words setting the option `MAT_KEEP_NONZERO_PATTERN to PETSC_FALSE has no effect on this routine. Bug 2 - The short circuit if (!((Mat_SeqAIJ *)(l->A->data))->keepnonzeropattern) is wrong because that flag is meaningless for this operation. The check should be changed (I think) to nonew instead of keepnonzeropattern. I will prepare a bug fix later today. For now you can just use the MAT_KEEP_NONZERO_PATTERN option for you code to work. Barry
On Feb 9, 2024, at 8:00 AM, Jeremy Theler (External) <[email protected]> wrote:
Because of a combination of settings, our code passes through this line:
https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/field...
i.e. the matrices associated with each of the sub-KSPs of a fieldsplit are destroyed and then re-created later. The thing is that one of these destroyed matrices had a near nullspace attached, which is lost because the new matrix does > not have it anymore.
Is this a bug or are we missing something?
I just want to get a clear picture. You create a PCFIELDSPLIT, set it up, then pull out the matrices and attach a nullspace before the solve.
We need to solve an SNES. We use dmplex so we have the jacobian allocated before starting the solve. At setup time we
1. define the PC of the KSP of the SNES to be fieldsplit 2. define the fields with ISes 3. call PCSetup() to create the sub-KSPs 4. retrieve the matrix attached to the sub-KSP that needs the near nullspace and attach it to that matrix
At a later time, you start another solve with this PC, and it has the DIFFERENT_NONZERO_PATTERN flag, so it recreates these matrices and loses your attached nullspace.
At a later time, in the jacobian evaluation we populate the global matrix (i.e. not the matrices attached to each sub-KSPs) and then we set dirichlet bcs with MatZeroRowsColumns() on that same global matrix. For some reason, in serial the near nullspace is not lost but in parallel the call to MatZeroRowsColumns() does change the non-zero structure (even though the manual says it does not) and then the code goes through that line 692 in fieldsplit.c and the near nullspace is lost.
First, does the matrix really change?
Well, the matrix during setup is not filled in, just allocated. The thing is that if we set MAT_KEEP_NONZERO_PATTERN to true with MatSetOption() before setting the dirichlet BCs, then the near nullspace is not lost (because the code does not go through line 692 of fieldsplit.c).
So there are (at least) two issues:
1. Code going through line 692 looses 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
Second, I had the same problem. I added callbacks to DM which allow a nullspace to be automatically attached if you extract a certain subfield. Are you using a DM?
Yes. Can you give us an example?
Regards -- jeremy
The bug fix for 2 is availabel in https://gitlab.com/petsc/petsc/-/merge_requests/7279
On Feb 9, 2024, at 10:50 AM, Barry Smith <[email protected]> wrote:
1. Code going through line 692 looses 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 (unless MAT_KEEP_NONZERO_PATTERN is used)
MatZeroRowsColumns() manual page states:
Unlike `MatZeroRows()` this does not change the nonzero structure of the matrix, it merely zeros those entries in the matrix.
MatZeroRowsColumns_MPIAIJ() has the code
/* only change matrix nonzero state if pattern was allowed to be changed */ if (!((Mat_SeqAIJ *)(l->A->data))->keepnonzeropattern) { PetscObjectState state = l->A->nonzerostate + l->B->nonzerostate; PetscCall(MPIU_Allreduce(&state, &A->nonzerostate, 1, MPIU_INT64, MPI_SUM, PetscObjectComm((PetscObject)A))); }
The if() test is simply an optimization to avoid the reduction if keepnonzeropattern is true. In your first run (when MAT_KEEP_NONZERO_PATTERN is not used) the keepnonzeropattern is not set so the A matrix nonzerostate is updated using the nonzerostate of the two submatrices. But that A nonzero state will only change if one of the two submatrices nonzerostate changes. In the second run the A nonzerostate cannot be changed.
MatZeroRowsColumns_SeqAIJ() has the code
if (diag != 0.0) { PetscCall(MatMissingDiagonal_SeqAIJ(A, &missing, &d)); if (missing) { for (i = 0; i < N; i++) { if (rows[i] >= A->cmap->N) continue; PetscCheck(!a->nonew || rows[i] < d, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONGSTATE, "Matrix is missing diagonal entry in row %" PetscInt_FMT " (%" PetscInt_FMT ")", d, rows[i]); PetscCall(MatSetValues_SeqAIJ(A, 1, &rows[i], 1, &rows[i], &diag, INSERT_VALUES));
which adds a nonzero to the matrix to fill in a missing diagonal. This will change the nonzerostate of the matrix, but in a different way than keepnonzeropattern flag is for.
So, do you have missing diagonals in your matrix?
Bug 1 - documentation The documentation statement "Unlike `MatZeroRows()` this does not change the nonzero structure of the matrix, it merely zeros those entries in the matrix." appears to be incorrect for matrices missing diagonal entries since new nonzeros are added (to fill in the diagonal). The documentation should really say
"Unlike `MatZeroRows()` this routine cannot remove the zeroed entries from the nonzero structure of the matrix; in other words setting the option `MAT_KEEP_NONZERO_PATTERN to PETSC_FALSE has no effect on this routine.
Bug 2 - The short circuit if (!((Mat_SeqAIJ *)(l->A->data))->keepnonzeropattern) is wrong because that flag is meaningless for this operation. The check should be changed (I think) to nonew instead of keepnonzeropattern.
I will prepare a bug fix later today. For now you can just use the MAT_KEEP_NONZERO_PATTERN option for you code to work.
Barry
On Feb 9, 2024, at 8:00 AM, Jeremy Theler (External) <[email protected]> wrote:
Because of a combination of settings, our code passes through this line:
https://gitlab.com/petsc/petsc/-/blob/main/src/ksp/pc/impls/fieldsplit/field...
i.e. the matrices associated with each of the sub-KSPs of a fieldsplit are destroyed and then re-created later. The thing is that one of these destroyed matrices had a near nullspace attached, which is lost because the new matrix does > not have it anymore.
Is this a bug or are we missing something?
I just want to get a clear picture. You create a PCFIELDSPLIT, set it up, then pull out the matrices and attach a nullspace before the solve.
We need to solve an SNES. We use dmplex so we have the jacobian allocated before starting the solve. At setup time we
1. define the PC of the KSP of the SNES to be fieldsplit 2. define the fields with ISes 3. call PCSetup() to create the sub-KSPs 4. retrieve the matrix attached to the sub-KSP that needs the near nullspace and attach it to that matrix
At a later time, you start another solve with this PC, and it has the DIFFERENT_NONZERO_PATTERN flag, so it recreates these matrices and loses your attached nullspace.
At a later time, in the jacobian evaluation we populate the global matrix (i.e. not the matrices attached to each sub-KSPs) and then we set dirichlet bcs with MatZeroRowsColumns() on that same global matrix. For some reason, in serial the near nullspace is not lost but in parallel the call to MatZeroRowsColumns() does change the non-zero structure (even though the manual says it does not) and then the code goes through that line 692 in fieldsplit.c and the near nullspace is lost.
First, does the matrix really change?
Well, the matrix during setup is not filled in, just allocated. The thing is that if we set MAT_KEEP_NONZERO_PATTERN to true with MatSetOption() before setting the dirichlet BCs, then the near nullspace is not lost (because the code does not go through line 692 of fieldsplit.c).
So there are (at least) two issues:
1. Code going through line 692 looses 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
Second, I had the same problem. I added callbacks to DM which allow a nullspace to be automatically attached if you extract a certain subfield. Are you using a DM?
Yes. Can you give us an example?
Regards -- jeremy
Hi Barry
The bug fix for 2 is availabel in https://gitlab.com/petsc/petsc/-/merge_requests/7279
Note that our code goes through MatZeroRowsColumns_MPIAIJ() not through MatZeroRows_MPIAIJ() so this fix does not change anything for the case I mentioned. -- jeremy
On Mon, Feb 12, 2024 at 7:10 AM Jeremy Theler (External) < [email protected]> wrote:
Hi Barry
The bug fix for 2 is availabel in https://gitlab.com/petsc/petsc/-/merge_requests/7279
Note that our code goes through MatZeroRowsColumns_MPIAIJ() not through MatZeroRows_MPIAIJ() so this fix does not change anything for the case I mentioned.
The other fixes are in https://gitlab.com/petsc/petsc/-/merge_requests/7279 Thanks, Matt
-- jeremy
-- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
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))); }
On Feb 12, 2024, at 7:02 AM, Jeremy Theler (External) <[email protected]> wrote:
Hi Barry
The bug fix for 2 is availabel in https://gitlab.com/petsc/petsc/-/merge_requests/7279
Note that our code goes through MatZeroRowsColumns_MPIAIJ() not through MatZeroRows_MPIAIJ() so this fix does not change anything for the case I mentioned.
-- jeremy
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: 1. reads a global matrix A from a file. 2. reads two ISs from another file 3. sets up a field split PC 4. attaches a near nullspace (read from another file) to one of the sub-KSPs 5. 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
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>
The new code is in https://gitlab.com/petsc/petsc/-/merge_requests/7293 and retains the null space on the submatrices for both MatZeroRows() and MatZeroRowsAndColumns() regardless of changes to the nonzero structure of the matrix. Barry
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>
The new code is in https://gitlab.com/petsc/petsc/-/merge_requests/7293 and retains the null space on the submatrices for both MatZeroRows() and MatZeroRowsAndColumns() regardless of changes to the nonzero structure of the matrix.
Thank you very much Barry. It works as expected. -- jeremy
participants (3)
-
Barry Smith -
Jeremy Theler (External) -
Matthew Knepley