[petsc-users] Using a Matrix Shell with SLEPc

[Jed Brown]

Miguel Arriaga <miguelarriagaecunha at gmail.com> writes:

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

Fortran 90 has "derived types" that you can use for this.  If you insist
on using Fortran 77, bad software design is inevitable (you can use
magic arrays or globals).

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

MATOP_CREATE was removed (it existed for something different anyway and
is unused now).  It should have been deleted so that people like you
wouldn't stumble upon it.

I recommend doing that setup up-front and set the context when you set
MATOP_MULT.

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

M is pretty well conditioned and direct solves are can be mighty
expensive in 3D.  So try it.  For a 2D 5-point stencil, direct might
win.  For a 3D 13-point stencil (second neighbors, a class of fourth
order methods), iterative is probably vastly better.

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

Start simple, EPS_LARGEST_REAL.
-------------- next part --------------
A non-text attachment was scrubbed...
Name: signature.asc
Type: application/pgp-signature
Size: 818 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-users/attachments/20141108/39864aef/attachment.pgp>