[petsc-users] Matrix-free method in PETSc

[Matthew Knepley]

On Mon, Feb 17, 2020 at 8:56 AM Yuyun Yang <yyang85 at stanford.edu> wrote:

> Hello,
>
> I actually have a question about the usage of DMDA since I'm quite new to
> this. I wonder if the DMDA suite of functions can be directly called on
> vectors created from VecCreate? Or the vectors have to be formed by
> DMDACreateGlobalVector?
>

Most things work the same. About the only thing that is different is that
we set a special viewer for vectors from DMDACreateGlobalVector()
which puts it in lexicographic order on output.


> I'm also not sure about what the dof and stencil width arguments do.
>

'dof' is how many unknowns lie at each vertex. 'sw' is the width of the
ghost region for local vectors.


> I'm still unsure about the usage of MatCreateShell and
> MatShellSetOperation, since it seems that MyMatMult should still have 3
> inputs just like MatMult (the matrix and two vectors).
>

MatShell is a type where you provide your own function implementations,
rather than using those for a particular storage format, like AIJ.


> Since I'm not forming the matrix, does that mean the matrix input is
> meaningless but still needs to exist for the sake of this format?
>

It means you calculate the output yourself, using the input.


> After I create such a shell matrix, can I use it like a regular matrix in
> KSP and utilize preconditioners?
>

Many preconditioners want access to individual elements of the matrix,
which usually will not work with shell matrices, since the
user just wants to provide the multiply routine.

  Thanks,

    Matt


> Thanks!
> Yuyun
> ------------------------------
> *From:* petsc-users <petsc-users-bounces at mcs.anl.gov> on behalf of Yuyun
> Yang <yyang85 at stanford.edu>
> *Sent:* Sunday, February 16, 2020 3:12 AM
> *To:* Smith, Barry F. <bsmith at mcs.anl.gov>
> *Cc:* petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
> *Subject:* Re: [petsc-users] Matrix-free method in PETSc
>
> Thank you, that is very helpful information indeed! I will try it and send
> you my code when it works.
>
> Best regards,
> Yuyun
> ------------------------------
> *From:* Smith, Barry F. <bsmith at mcs.anl.gov>
> *Sent:* Saturday, February 15, 2020 10:02 PM
> *To:* Yuyun Yang <yyang85 at stanford.edu>
> *Cc:* petsc-users at mcs.anl.gov <petsc-users at mcs.anl.gov>
> *Subject:* Re: [petsc-users] Matrix-free method in PETSc
>
>   Yuyun,
>
>     If you are speaking about using a finite difference stencil on a
> structured grid where you provide the Jacobian vector products yourself by
> looping over the grid doing the stencil operation we unfortunately do not
> have exactly that kind of example.
>
>     But it is actually not difficult. I suggest starting with
> src/ts/examples/tests/ex22.c It computes the sparse matrix explicitly with
> FormIJacobian()
>
>     What you need to do is instead in main() use MatCreateShell() and
> MatShellSetOperation(,MATOP_MULT,(void (*)(void))MyMatMult) then provide
> the routine MyMatMult() to do your stencil based matrix free product; note
> that you can create this new routine by taking the structure of IFunction()
> and reorganizing it to do the Jacobian product instead. You will need to
> get the information about the shell matrix size on each process by calling
> DMDAGetCorners().
>
>     You will then remove the explicit computation of the Jacobian, and
> also remove the Event stuff since you don't need it.
>
>      Extending to 2 and 3d is straight forward.
>
>      Any questions let us know.
>
>    Barry
>
>    If you like this would make a great merge request with your code to
> improve our examples.
>
>
> > On Feb 15, 2020, at 9:42 PM, Yuyun Yang <yyang85 at stanford.edu> wrote:
> >
> > Hello team,
> >
> > I wanted to apply the Krylov subspace method to a matrix-free
> implementation of a stencil, such that the iterative method acts on the
> operation without ever constructing the matrix explicitly (for example,
> when doing backward Euler).
> >
> > I'm not sure whether there is already an example for that somewhere. If
> so, could you point me to a relevant example?
> >
> > Thank you!
> >
> > Best regards,
> > Yuyun
>
>

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

https://www.cse.buffalo.edu/~knepley/ <http://www.cse.buffalo.edu/~knepley/>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20200217/1c986654/attachment.html>