adding SNESSetLinearSolve()

[Barry Smith]

On Mon, 29 Oct 2007, Lisandro Dalcin wrote:

> Barry, why do you really resist this idea? Do you believe it is an
> example of 'too many ways to do it'? All I am asking for is let users
> replace the call to KSPSolve for a custom call, wich in fact CAN
> internally call KSPSolve after getting the KSP from SNES.
>

  I strongly resist this. Back in the bad old days, 1995?, there was
some pressure to "allow SNES to use any linear solver package, not just
PETSc's". WELL, the whole idea then and now is PETSc KSP is suppose to
be a powerful enough wrapper to wrap ANY linear solver package; BUT
now you want to impose a wrapper class below SNES that can call any
solver package. This is obscene! KSP IS SUPPOSE TO BE EXACTLY THAT
WRAPPER CLASS!!!!! If it is not good enough we should fix it, not hack
a wrapper below SNES. SNES uses KSP, and ONLY KSP to solve linear
systems (of course since KSP/PC is a wrapper class it can do whatever
you want). If we just willy-nilly shove in extra functions randomly
in PETSc, without regard to the overall design of PETSc we end up
with Trilinos.

> Yes, Barry, I know we can do this, but this trick is not fully
                                             ^^^^^^^

  In absolutely NO WAY is this a "trick". This approach is one
of the most basic philosphical foundations of PETSc. Just adding a
SNESSetLinearSolver() is a trick.

> fuctional, or you have to make a lot of more trickery:
> 
> - You have to explicitely use a KSPPREONLY, thus you cannot reuse this
> ksp in your 'LinearSolve' routine. Then, you cannot get EW method work
> in a easy way.

   Can you please describe your linear solver algorithm in more detail
using the following notation: 

   As is the SPARSE matrix of size n
   A is of size n+1 and looks like
             ( As      stuff_1 )
        A =  (stuff_2  something)
  stuff_1 is n by 1 and stuff_2 is 1 by n.

  Newton's method requires solving A x = -F(u) where u is a vector
of size n+1 as is x.

  Is the algorithm: solve A x = -F(u) with GMRES using the preconditioner
XXXXXXXXXX? Or is it something else? Please describe XXXXXX or the something 
else.

> 
> - PCSHELL to currently have some deficiencies. For example, there is
> not a easy way to get matrix operators inside user-defined PCSHELL
> routines, unless you implement PCPreSolve, which gives access to KSP.

  As was discussed last week in email with Victor Eijkhout, the PCSHELL
needs to be fixed so that the first argument to the PCApply_yours() etc
is the PC, NOT the void*.

> And PCSHELL does not have a SetFromOptions() feature.

  This is an oversight and trivially fixed.
> 

  Don't use a poor implementation of one of PETSc's classes
to justify a horrible hack to PETSc.

  Best regards,

   Barry



> 
> On 10/29/07, Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> >   SNESGetKSP(
> >   KSPSetType(ksp,KSPPREONLY)
> >   KSPGetPC(
> >   PCSetType(pc,PCSHELL);
> >   PCShellSetApply(pc,YourSolver);
> >
> >
> > On Mon, 29 Oct 2007, Lisandro Dalcin wrote:
> >
> > > Some time ago I made a request for adding SNESSetLinearSolve(), in
> > > order to let users set a custom linear solve routine to SNES. I was
> > > not clear to me if this idea was finally accepted, so I want to ask
> > > again.
> > >
> > > I really needs this feature in the near future, I I would like this to
> > > go in the next petsc release. My use case is rather simple: I'm
> > > helping a coworker to use SNES for solving a nonlinear optimization
> > > problem related to mesh movement (that is, move nodes retaining
> > > element connectivities). This problem uses a regularized functional,
> > > wich introduces a global degree of freedom coupled with all others
> > > dofs. To avoid assembling a sparse parallel AIJ matrix with a dense
> > > column and row, we endup needing to solve two linear systems with the
> > > same matrix for actually solving the Hessian of the regularized
> > > functional.
> > >
> > >
> > >
> > >
> >
> >
> 
> 
>