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

  Fande,

What SNES method are you using?  If you use SNESKSPONLY I think it is ok, it
will solve for the norm minimizing least square solution during the one
KSPSolve() and then return.

Yes, if you use SNESNEWTONLS or others though the SNES solver will, as you
say, think that progress has not been made.

I do not like what you propose to do, changing the right hand side of the
system the user provides is a nasty and surprising side effect.

What is your goal? To make it look like the SNES system has had a residual
norm reduction?

We could generalize you question and ask what about solving for nonlinear
problems: find the minimal norm solution of min_x || F(x) - b||. This may or
may not belong in Tao, currently SNES doesn't do any kind of nonlinear least
squares.

Barry

> On Oct 13, 2016, at 5:20 PM, Kong, Fande <fande.k...@inl.gov> wrote:
>
> 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 <knep...@gmail.com> wrote:
> On Thu, Oct 13, 2016 at 9:06 AM, Kong, Fande <fande.k...@inl.gov> wrote:
>
>
> On Wed, Oct 12, 2016 at 10:21 PM, Jed Brown <j...@jedbrown.org> wrote:
> Barry Smith <bsm...@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