[petsc-users] A matrix-free Jacobian for SNES using MatShell
Barry Smith
bsmith at mcs.anl.gov
Tue Aug 31 13:27:49 CDT 2010
On Aug 31, 2010, at 2:17 PM, Jeremy Roberts wrote:
> To implement the safe method---i.e. copying the unknown vector u within FOO---what is proper generic form of FOO?
>
> I guess I'm not clear exactly how FOO and jac_shell (the J*x routine) are linked. It appears they aren't, given I can use an empty subroutine, but how then is a local copy of u within FOO used in for J(u) in J(U)*x of jac_shell?
Don't ask me. You wrote the code, I have no idea how you are accessing u. But you need to compute J(u)*x so you need to use u.
subroutine yourfoo(snes,u,A,B,flg,ierr)
common /mycommonblock/ v
call veccopy(u,v,ierr)
return
and create the v upfront in your main program. Of course, you don't need to use common blocks but that is one way of maintaining this global variables.
Barry
>
>
> On Tue, Aug 31, 2010 at 1:14 PM, Barry Smith <bsmith at mcs.anl.gov> wrote:
>
> On Aug 31, 2010, at 1:03 PM, Jeremy Roberts wrote:
>
> >
> > Hi Barry,
> >
> > Thanks for the quick response.
> >
> > I have the unknown vector u globally accessible and so forgo use of the user-defined context (I tried to get it to work---can it be used in Fortran?).
> >
> > I used
> >
> > call SNESSetJacobian( snes, Jshell, Jshell, FOO, PETSC_NULL_OBJECT, ierr )
> >
> > where FOO is an empty function having the correct argument types. It seems to work now, and I get the exact same answer using an explicit Jacobian matrix and the shell version. Is there a way to get around needing FOO at all?
> >
>
> If it works great. But I how you know what the "unknown vector u " that is "globally accessible" in SNES is? How do you know it is the same each time a new Jacobian is needed? The only save way is to copy the vector passed into the FOO funcion.
>
> Barry
>
>
>
> > Jeremy
> >
> >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20100831/4ed932eb/attachment.htm>
More information about the petsc-users
mailing list