[petsc-users] Algorithms to remove null spaces in a singular system

[Kong, Fande]

One more question.

Suppose that we are solving the singular linear system Ax = b. N(A) is the
null space of A, and N(A^T) is the null space of the transpose of A.

The linear system is solved using SNES, that is, F(x) = Ax-b = Ax -b_r -
b_n. Here b_n  in N(A^T),  and b_r in R(A).  During each nonlinear
iteration, a linear system A \delta x = F(x) is solved. N(A) is applied to
Krylov space  during the linear iterating. Before the actual solve
"(*ksp->ops->solve)(ksp)" for \delta x,  a temporary copy of F(x) is made,
F_tmp. N(A^T) is applied to F_tmp. We will get a \delta x.  F(x+\delta x )
= A(x+\delta x)-b_r - b_n.

F(x+\delta x ) always contain the vector b_n, and then the algorithm never
converges because the normal of F is at least 1.

Should we apply N(A^T) to F instead of F_tmp so that b_n can be removed
from F?

MatGetTransposeNullSpace(pmat,&nullsp);
if (nullsp) {
   VecDuplicate(ksp->vec_rhs,&btmp);
   VecCopy(ksp->vec_rhs,btmp);
   MatNullSpaceRemove(nullsp,btmp);
   vec_rhs      = ksp->vec_rhs;
   ksp->vec_rhs = btmp;
}

should  be changed to

MatGetTransposeNullSpace(pmat,&nullsp);
if (nullsp) {
   MatNullSpaceRemove(nullsp,ksp->vec_rhs);
}
???

Or other solutions to this issue?


Fande Kong,





On Thu, Oct 13, 2016 at 8:23 AM, Matthew Knepley <knepley at gmail.com> wrote:

> On Thu, Oct 13, 2016 at 9:06 AM, Kong, Fande <fande.kong at inl.gov> wrote:
>
>>
>>
>> On Wed, Oct 12, 2016 at 10:21 PM, Jed Brown <jed at jedbrown.org> wrote:
>>
>>> Barry Smith <bsmith at mcs.anl.gov> writes:
>>> >   I would make that a separate routine that the users would call first.
>>>
>>> We have VecMDot and VecMAXPY.  I would propose adding
>>>
>>>   VecQR(PetscInt nvecs,Vec *vecs,PetscScalar *R);
>>>
>>> (where R can be NULL).
>>>
>>
>> What does R mean here?
>>
>
> It means the coefficients of the old basis vectors in the new basis.
>
>   Matt
>
>
>> If nobody working on this, I will be going to take a try.
>>
>> Fande,
>>
>>
>>>
>>> Does anyone use the "Vecs" type?
>>>
>>
>>
>
>
> --
> 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
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20161013/0fd010fc/attachment.html>