Re: [petsc-users] Matrix-free method in PETSc

[Smith, Barry F. via petsc-users]

> On Feb 17, 2020, at 7:56 AM, Yuyun Yang <yyan...@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?

   Yes, but you have to make sure the ones you create have the same sizes and 
parallel layouts. Generally best to get them from the DMDA or VecDuplicate() 
than the hassle of figuring out sizes.

> Or the vectors have to be formed by DMDACreateGlobalVector? I'm also not sure 
> about what the dof and stencil width arguments do.
> 
> 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). 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?

    Well the matrix input is your shell matrix so it likely has information you 
need to do your multiply routine. MatShellGetContext() (No you do not want to 
put your information about the matrix stencil inside global variables!)


> 
> After I create such a shell matrix, can I use it like a regular matrix in KSP 
> and utilize preconditioners?
> 
> Thanks!
> Yuyun
> From: petsc-users <petsc-users-boun...@mcs.anl.gov> on behalf of Yuyun Yang 
> <yyan...@stanford.edu>
> Sent: Sunday, February 16, 2020 3:12 AM
> To: Smith, Barry F. <bsm...@mcs.anl.gov>
> Cc: petsc-users@mcs.anl.gov <petsc-users@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. <bsm...@mcs.anl.gov>
> Sent: Saturday, February 15, 2020 10:02 PM
> To: Yuyun Yang <yyan...@stanford.edu>
> Cc: petsc-users@mcs.anl.gov <petsc-users@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 <yyan...@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