adding SNESSetLinearSolve()

[Lisandro Dalcin]

On 10/29/07, Barry Smith <bsmith at mcs.anl.gov> wrote:
>   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!!!!!

In practice, the only way to wrap any solver package is
KSPPREONLY+PCSHELL. And you actually implement the linear solve in PC,
not KSP.

This combination work great in many, many (almot all?) situations.
At the cost of lossing some features (perhaps not commonly used). Have
you ever tried to make Eisenstat&Walker working using a
KSPPREONLY+PCSHELL combination, where your Shell PC is actually using
a inner KSP? Of course, you can do it, but the way is not immediate.

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

Understood.

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

You completelly convinded me with this argument ;-)

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

In short, I need to solve two linear systems with the same matrix
using any PETSc KSP+PC combination and different rhs vectors.  This is
for managing the update of a global dof (coupled with all other dofs)
evolving in the nonlinear iteration.

But forget my use case. Special cases aren't special enough to break the rules.

> > - PCSHELL to currently have some deficiencies.
>   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*.

Great!! I asked for changing this many, many time ago, and my proposal
was not accepted!! I had to find my way to pass-over this problem in
petsc4py.

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

Well, PETSc already have some public stuff that could be seen as
hacks, symply because practicality beats purity. But I have to agree
that the way to go is extending/remplementing things in the right
place, and not by adding random stuff.

I'll try to rethink the way a user can implement custom linear solves
on KSP, and let you know my results.

If PETSc internals were implemented in C++, I would simply do:

MyKSP: public KSP {
void solve(Vec b, Vec x) {
   /* do something */
   KSP::solve(b1,x1)
   /* do more*/
   KSP::solve(b2, x2)
  /* fill result */
  x = ...
}
}

But with the current implementation in C, this is possible but not so
easy to do.



-- 
Lisandro Dalcín
---------------
Centro Internacional de Métodos Computacionales en Ingeniería (CIMEC)
Instituto de Desarrollo Tecnológico para la Industria Química (INTEC)
Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET)
PTLC - Güemes 3450, (3000) Santa Fe, Argentina
Tel/Fax: +54-(0)342-451.1594