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.

## Advertising

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 > lead. > -- Norbert Wiener >