[petsc-users] Using a Matrix Shell with SLEPc

[Miguel Arriaga]

Hi,
Thank you so much for your reply.

1
One note that I forgot to mention is that I'm using FORTRAN (sorry!).
Is it possible to set up the context to pass M and the ksp for M? I also
try to avoid common blocks but in this case it seemed inevitable. Since M
and A are already global in the code that I'm using maybe I could just use
the kspM in the context.

The other side of my question was more on whether doing MatCreateShell(...)
and then MatShellSetOperation(...,MATOP_CREATE,...) would work because it
seems like a chicken-vs-egg kind of situation. Is the CREATE operation
executed after this? I couldn't find an example that uses this operation
and if I google MATOP_CREATE the only thing I get is the source code of the
file petscmat.h.

2
Thanks!
I'm sorry if this was already included in your answer but I just wanted to
make sure. For a moderately large sparse problem (100K equations)
considering that I'm using two KSPSolve in each call to MatMult, and I'm
going to use this Shell in an eigensolver (that will probably call MatMult
several times), wouldn't it end up being faster to just do a Sparse Direct
to get the LU decomposition of M and then just do the forward and back
substitution?

3
My colleague suggested that I do a first run without shift for
EPS_LARGEST_MAGNITUDE and then, if this value is negative, I would use the
eigenvalue as a shift and compute EPS_LARGEST_MAGNITUDE again. This seems
like a fine solution to me but I was wondering if just doing
EPS_LARGEST_REAL without shift would already do something like this but
more efficiently since, for example, you can detect very soon if your
eigenvalue is converging to a negative number, allowing with just a few
iterations to get the magnitude of the necessary shift. Should I just add a
very large shift to begin with? Or is EPS_LARGEST_REAL fine?

Thank for your time,
Miguel




On Fri, Nov 7, 2014 at 5:56 PM, Jed Brown <jed at jedbrown.org> wrote:

> Miguel Arriaga <miguelarriagaecunha at gmail.com> writes:
>
> > Hi there,
> > I have a matrix shell where the Mult operation looks like this:
> > y=1/2(K+K*) x
> > with K=M^-1 A and K* its transpose.
> >
> > 1- For setting-up the operation MatMult, since I need to do a KSPSolve on
> > M, should I set-up the KSP on the operation MATOP_CREATE and destroy it
> on
> > MATOP_DESTROY and store it as a global variable?
>
> Store it in the context.  That's why callbacks always have a context.
> Global variables are terrible.
>
> > 2- For "inverting" M, considering that M is symmetric, Should I (a) use
> LU,
> > or should I (b) use an iterative procedure like CG with ILU, since M is
> > symmetric positive definite? Is there a faster iterative method for a SPD
> > matrix?
>
> CG/Jacobi is usually fine for the mass matrix.
>
> > 3- I want to use this matrix shell in SLEPc to compute the largest
> positive
> > eigenvalue, which will be smaller in magnitude than the smallest negative
> > real (possibly very close to zero). Since it is a matrix shell I assume
> > that I can't really use Shift-and-invert without choosing an iterative
> > procedure that only requires MatMult.
>
> Start with a simple shift (not shift-and-invert).
>
> > What would be the best Method+Preconditioner for this case? Note that
> > this is a Real Symmetric matrix (but not Positive definite). Should I
> > set all of these as true?
> >
> > MATOP_IS_HERMITIAN
> > MATOP_IS_STRUCTURALLY_SYMMETRIC
> > MATOP_IS_SYMMETRIC
> >
> > Thank you,
> > Miguel
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141107/4a77772d/attachment.html>