[petsc-users] move from KSPSetNullSpace to MatSetNullSpace

[Stephan Kramer]

>> 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