analyze preconditioned operator?

[Matthew Knepley]

On Tue, Nov 4, 2008 at 2:33 PM, Matt Funk <mafunk at nmsu.edu> wrote:
> 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.

Simple to explain. You get the ORIGINAL_MATRIX from a global location. If you
can't or won't use a global variable, then you need the context.

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

It will raise a PetscError, which SLEPc should propagate back to you.

  Matt

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



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