analyze preconditioned operator?

Tue Nov 4 13:33:11 CST 2008

Hi,

first of all: sorry to be picking this back up so late.
(If you guys cannot answer because of that i understand. I was busy with other 
stuff and am now just returning to this ... )

As was described earlier (see below) i create a Shell Matrix and registered 
the mulitplying fcn with it via MatSetShellOperation(...).

Regarding the myApply(...) fcn (which provides the custom multiply)
I looked at the Fortran example for using MattSetShellOperations at:
http://www-unix.mcs.anl.gov/petsc/petsc-as/snapshots/petsc-current/src/ksp/ksp/examples/tutorials/ex14f.F.html

Question #1:

In it (=>the fortran code) the MatGetContext fcn (which was included in the 
myApply fcn (see below)) is never used. So i took it out as well. Things 
appear to work (at least nothing crashes and i get an answer), so am i still 
missing something? 

(I simply do this:
1) ORIGINAL_MATRIX*VECTOR_X => WORKSPACE
2) PCAPPLY(PC, WORKSPACE, VECTOR_Y)

where ORIGINAL_MATRIX is a pointer to the original matrix. This pointer is 
stored in a global structure that can be accessed in the myApply fcn)

I guess i never understood why it (=>MatGetContext fcn) was used anyway. The 
main point of the myApply() fcn below is to return the vector y, correct? So 
what is the purpose of the MatGetContext fcn? 
I'd like to understand it, but if it is too complicated to explain, then don't 
worry about it.

One sanity check i did to see whether things make sense was to replace:

...
MatMult(ctx->M, x, ctx->work);
PCApply(ctx->pc, ctx->work, y);
...

in the myApply fcn with simply:

...
MatMult(ctx->M, x, y);
...

which should then return the same answer as if i had just used my original 
matrix in the first place: This checks out. 

So things look ok, however, just out of curiosity i have another question:

Question #2 :

I registered the MULT operation with the shell matrix. However, since there 
are many other operations defined in petscmat.h and i don't know whether 
SLEPC uses those or not, what happens if were to to use another operation?

Will it fall back to the original operator since it is not defined for the new 
(shell) matrix? Will it crash? 
What is the behavior going to be?

thanks for all the help
matt

On Friday 10 October 2008, Matthew Knepley wrote:
> On Fri, Oct 10, 2008 at 3:11 PM, Matt Funk <mafunk at nmsu.edu> wrote:
> > Ok,
> >
> > thanks, i think i finally got the main idea (the emphasize here being on
> > 'i think' which might not mean much ...)
> >
> > One other question though:
> > What do i declare ctx to be: PetscObject * ctx ?
>
> No, you declare a struct which holds
>
>   a) The original matrix
>
>   b) The preconditioner
>
>   c) A work vector
>
> or you would not be able to code the apply method.
>
> > Another question related to KSPComputeExplicitOperator:
> >
> > I was trying to use KSPComputeExplicitOperator and i am having issues
> > with it as well. My matrix is a sparse matrix of length ((2*73)^3) with a
> > maximum of ~30 entries per row.
>
> This is way way way too big to form the operator explicitly.
>
> > I believe that in order to use KSPComputeExplicitOperator i need to set
> > up a matrix, allocate the memory for it and then pass it to the routine.
> > From what i gather on the reference page this matrix needs to be of type
> > MatDense? (if so then i cannot use KSPComputeExplicitOperator due to the
> > size of the resulting matrix)
> > However, it said on the website that when multiple procs are used it uses
> > a sparse format.
> > So why cannot i not use a sparse format in serial?
>
> 1) This operator is in general not sparse at all
>
> 2) We only use a sparse format in parllel because matMPIDense was
> having problems
>
> It is instructive to look at what the code is doing. The function is
> very short and you
> can get to it right from the link on the manpage.
>
>   Matt
>
> > thanks
> > matt
> >
> > On Thursday 09 October 2008, Matthew Knepley wrote:
> >> On Thu, Oct 9, 2008 at 5:58 PM, Matt Funk <mafunk at nmsu.edu> wrote:
> >> > So then,
> >> >
> >> > is it then correct to say that this code "registers" this extra
> >> > operation (i.e. applying the preconditioner) with the matrix context
> >> > such that whenever a matrix operation involving this matrix is invoked
> >> > (like MatMult for example) the pc-applying fcn (i.e. myApply() ) is
> >> > called first?
> >>
> >> The idea here is to replace a given matrix A, which you are passing to
> >> SLEPc, with another matrix M, which is a shell matrix. When MatMult() is
> >> called on M, we call MatMult on A, and then PCApply on the result.
> >>
> >>   Matt
> >>
> >> > matt
> >> >
> >> > ps: sorry for being a little slow on the uptake here ...
> >> >
> >> > On Thursday 09 October 2008, Matthew Knepley wrote:
> >> >> On Thu, Oct 9, 2008 at 5:25 PM, Matt Funk <mafunk at nmsu.edu> wrote:
> >> >> > Hi Matt,
> >> >> >
> >> >> > so, the basic idea with this code is to apply the pc to each column
> >> >> > vector of the matrix? Is that right?
> >> >>
> >> >> No, that is what KSPGetExplicitOperator() does. This applies the
> >> >> operator in a matrix-free way.
> >> >>
> >> >> > Also, in your example: when is myApply actually invoked? I also
> >> >> > looked at the example listed under the MatShellSetOperation
> >> >> > reference page.
> >> >>
> >> >> When MatMult(A) is called.
> >> >>
> >> >>   Matt
> >> >>
> >> >> > Is the function then actually internally called when
> >> >> > MatShellSetOperation is called, or when KSPSetOperators is called
> >> >> > or KSPSolve?
> >> >> >
> >> >> > The reason i am asking is that if it is called when KSPSolve called
> >> >> > then there is a problem because for the analysis i never call
> >> >> > KSPSolve directly.
> >> >> >
> >> >> > thanks
> >> >> > matt
> >> >> >
> >> >> > On Thursday 09 October 2008, Matthew Knepley wrote:
> >> >> >> On Thu, Oct 9, 2008 at 4:32 PM, Matt Funk <mafunk at nmsu.edu> wrote:
> >> >> >> > mmhh,
> >> >> >> >
> >> >> >> > i think i am missing something. Doesn't PCApply() apply the
> >> >> >> > preconditioner to a vector? So how would that work (easily) with
> >> >> >> > a matrix?
> >> >> >>
> >> >> >> You do not apply it to the matrix. Here is a skeleton (maybe has
> >> >> >> mistakes)
> >> >> >>
> >> >> >> void myApply(Mat A, Vec x, Vec y) {
> >> >> >>   MatShellGetContext(A, &ctx);
> >> >> >>   MatMult(ctx->M, x, ctx->work);
> >> >> >>   PCApply(ctx->pc, ctx->work, y);
> >> >> >> }
> >> >> >>
> >> >> >> MatShellSetOperation(A, MATOP_MULT, myApply)
> >> >> >>
> >> >> >>   Matt
> >> >> >>
> >> >> >> > matt
> >> >> >> >
> >> >> >> > On Thursday 09 October 2008, Matthew Knepley wrote:
> >> >> >> >> CApply().