# SNES for bordered block matrix

Even with one process, I don’t get a deterministic behavior.

`````` finding initial point
initial point found. moving to continuation part
starting loop: 1/5. Current nu = 0.1. Predicted nu = 0.09973
0 SNES Function norm 2.805428653365e-07
nu = 0.09973
1 SNES Function norm 1.695141030744e-07
nu = 0.0997301
2 SNES Function norm 6.278092520575e-13
nu = 0.0997308
3 SNES Function norm 1.354552924947e-14
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 3
starting loop: 2/5. Current nu = 0.0997308. Predicted nu = 0.0994608
0 SNES Function norm 1.891297254561e-07
nu = 0.0994608
1 SNES Function norm 7.815465447968e-13
nu = 0.0994616
2 SNES Function norm 1.375679948496e-14
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 2
starting loop: 3/5. Current nu = 0.0994616. Predicted nu = 0.0991916
0 SNES Function norm 2.822903166950e-07
nu = 0.0991916
1 SNES Function norm 1.709241490889e-07
nu = 0.0991917
2 SNES Function norm 6.395819877626e-13
nu = 0.0991924
3 SNES Function norm 1.349931255097e-14
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 3
``````

vs.

`````` finding initial point
initial point found. moving to continuation part
starting loop: 1/5. Current nu = 0.1. Predicted nu = 0.09973
0 SNES Function norm 2.805428653365e-07
nu = 0.09973
1 SNES Function norm 1.695141030744e-07
nu = 0.0997301
2 SNES Function norm 6.278092520575e-13
nu = 0.0997308
3 SNES Function norm 1.354552924947e-14
Nonlinear solve converged due to CONVERGED_SNORM_RELATIVE iterations 3
starting loop: 2/5. Current nu = 0.0997308. Predicted nu = 0.0994608
0 SNES Function norm 2.081688271474e-03
nu = 0.0994608
1 SNES Function norm 5.759118650479e-09
nu = 0.0994141
2 SNES Function norm 1.409543210835e-14
Nonlinear solve converged due to CONVERGED_FNORM_RELATIVE iterations 2
starting loop: 3/5. Current nu = 0.0994141. Predicted nu = 0.0991441
0 SNES Function norm 1.900550339798e-07
nu = 0.0991441
Nonlinear solve did not converge due to DIVERGED_LINE_SEARCH iterations 0
starting loop: 4/5. Current nu = 0.0991441. Predicted nu = 0.0988741
0 SNES Function norm 3.809037608138e-07
``````

If you find where it starts to behave differently, I could try to make sense of this, but I cannot debug the full code right now, sorry.

Interesting. Thank you. Appreciate your help. Will follow up if I find any clues.

Ah, OK, I see at least one simple error. You are calling `J(yq,yq)`, but `yq` is the result of a `KSPSolve()`, i.e., it’s in PETSc numbering. For `J(.,.)` you need to feed vectors in FreeFEM numbering. So, either:

• replace `KSPSolve()` by `^-1` or add a call to `ChangeNumbering(..., inverse = true, exchange = true)`, or,
• compute the scalar product yourself with a vector in PETSc numbering (`tmp = yq' * yq; mpiAllReduce(...);`).

Thank you for pointing out that issue. This does not solve the non-determinism I’m seeing, but I think it must be related to some similar bug I have involving a mix between PETSc and FreeFEM numbering. I noticed a couple other mix-ups as well, and I think I’m getting close to a working code. I should be able to figure it out from here, but thank you so much for all your help! Hi there Pierre, I finally found the time to get back to this, and I got both the “simple” 2x2 block system (N+1)x(N+1) and the more complex 3x3 block system (N+N+1)x(N+N+1) to work. For once, I don’t have any problems or questions for you, and instead just wanted to report back to say thank you, and all is well! Here’s the route I took:

In the 2x2 case, it was relatively straighforward using `-pc_fieldsplit_type schur` as you showed me. Here we have a block structure like

``````Mat JlPM(JlP.n, mpirank == 0 ? 1 : 0), yqPM(yqP.n, mpirank == 0 ? 1 : 0);
Ja = [[J,JlPM],[yqPM',-1.0]];
``````

which can be solved using

``````set(Ja, sparams = "-ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self ", setup = 1);
set(J, sparams = " -fieldsplit_0_ksp_type preonly -fieldsplit_0_pc_type lu ", prefix = "fieldsplit_0_", parent = Ja);
``````

However, the 3x3 case was a bit trickier. Eventually, I gave up on dealing with it directly as a 3x3 and did it like this:

``````Mat JJ = [[J, 0],[H, J]];
Mat JlPM(2*JlP.n, mpirank == 0 ? 1 : 0), qmPM(2*qmP.n, mpirank == 0 ? 1 : 0);
Mat JJa = [[JJ,JlPM],[qqmPM',0]];
``````

with

``````  func real[int] blockLU(real[int]& qqm) {
PetscScalar[int] out1, out2, in1 = qqm(0:J.n-1), in2 = qqm(J.n:2*J.n-1);
KSPSolve(J, in1, out1);
MatMult(H, out1, out2);
in2 -= out2;
KSPSolve(J, in2, out2);
PetscScalar[int] outPETSc = [out1, out2];
return outPETSc;
}
set(JJa, sparams = "-ksp_type preonly -pc_type fieldsplit -pc_fieldsplit_type schur -pc_fieldsplit_schur_precondition self ", setup = 1);
set(JJ, sparams = "-fieldsplit_0_ksp_type preonly -fieldsplit_0_pc_type shell ", precon = blockLU, prefix = "fieldsplit_0_", parent = JJa);
set(J, sparams = "-fieldsplit_ksp_type preonly -fieldsplit_pc_type lu", parent = JJ);
``````

Thanks again!

Eheh, a nest of nest is something completely valid which can be useful to breakdown complexity, as seen right here. Nicely done!

1 Like