Sorry if we're talking completely cross-purposes here: but the pcksp solve
(that doesn't actually have a nullspace) is inside a solve that does have a
nullspace. If what you are saying is that applying the nullspace inside the
pcksp solve does not affect the outer solve, I can only see that if the
difference in outcome of the pcksp solve (with or without the nullspace)
would be inside the nullspace, because such a difference would be projected
out anyway straight after the pcksp solve. I don't see that this is true
however. My reasoning is the following:
Say the original PCKSP solve (that doesn't apply the nullspace) gets
passed some residual r from the outer solve, so it solves:
Pz=r
where P is our mass matrix approximation to the system in the outer solve,
and P is full rank. Now if we do remove the nullspace during the pcksp
solve, effectively we change the system to the following:
NM^{-1}P z = NM^{-1}r
where M^{-1} is the preconditioner inside the pcksp solve (assuming
left-preconditiong here), and N is the projection operator that projects
out the nullspace. This system is now rank-deficient, where I can add:
z -> z + P^{-1}M n
for arbitrary n in the nullspace.
So not only is the possible difference between solving with or without a
nullspace not found in that nullspace, but worse, I've ended up with a
preconditioner that's rank deficient in the outer solve.
This analysis does not make sense to me. The whole point of having all this
nullspace stuff is so that we do not get
a component of the null space in the solution. The way we avoid these in
Krylov methods is to project them out,
which is what happens when you attach it. Thus, your PCKSP solution 'z'
will not have any 'n' component. Your
outer solve will also not have any 'n' component, so I do not see where the
inconsistency comes in.
The nullspace component does indeed get projected out in the pcksp
solve. However that does not mean it gives the right answer - the right
answer would be the same answer I would get before (when we weren't
subtracting the nullspace in the pcksp) modulo the nullspace. For instance
if my initial residual r is exactly Mn I would get the answer 0 instead
of P^{-1}Mn - or equally correct would have been P^{-1}Mn with
the nullspace component projected out afterwards, but not 0. So the
inconsistency is
in having a preconditioner that's not full-rank and whose nullspace
is different than the nullspace of the real matrix in the outer solve.
Some more observations:
* if I switch the preconditioner inside the pcksp to be a right-preconditioner
(-ksp_ksp_pc_side right) I get back the correct answer - which is kind of what
I expect as in that case it should just give the "correct" pcksp answer
with the nullspace projected out
* if instead of gmres+sor for the pcksp solve, I select cg+sor (which is
what I actually want), the pcksp solve fails with an indefinite pc. Both cases
produce the same wrong answer in the outer solve (both using default
left-preconditioning)
Hope this makes sense
Cheers
Stephan
