Re: [petsc-dev] plans for preconditioners for SeqSELL

[Richard Tran Mills]

All,

I agree with Karl that we need an easier way to make use of SELL for the 
specific operations it helps accelerate, while facilitating the many 
more types of operations that AIJ supports. I'm not sure what the best 
approach is, or if there is more than one approach it might be 
appropriate to support.

Karl suggests making a subclass of AIJ that uses the SELL format 
selectively. I did something similar with CSRPERM (now AIJPERM) once: a 
student and I added an option to store each group of rows with the same 
number of nonzero elements in ELLPACK format. A SeqAIJPERM matrix always 
keeps everything around that a SeqAIJ one does, so this meant keeping 
two versions of the matrix stored. (For reasons including laziness, I 
never bothered to put this code back into PETSc -- the current AIJPERM 
code only stores the matrix entries in AIJ format.)

Karl, are you thinking of a matrix subclass that has everything that an 
AIJ matrix does, but also keeps a SELL copy around for operations like 
MatMult()? Would it make sense to just keep another Mat inside (like how 
MPIAIJ keeps multiple Mat instances) that *is* of type MATSELL, that 
gets built/updated as needed? Would this involve carrying too much 
baggage around, associated with a complete Mat instance? I like the idea 
of having a MATSELL type available that is lean (insofar as not having 
storing an AIJ matrix) for those cases when a user really knows that the 
AIJ stuff will not be needed. But maybe it makes sense to be able to use 
that inside another matrix class. Perhaps we could have something, 
called, say, MATAIJMUTABLE that uses AIJ but might also create copies in 
SELL (or other formats, potentially) when appropriate -- perhaps based 
on a some performance model indicating which format is fastest for 
MatMult() or whatever.

Having a class that's an AIJ but can also use SELL is more convenient 
than adding a fallback to AIJ format inside MATSELL. I wonder if the 
latter option might be preferable in some circumstances, however, 
because it can avoid the extra memory footprint of also keeping the 
matrix in AIJ -- maybe AIJ operations are rarely needed and the AIJ 
conversion can just happen on a lazy basis.

--Richard

On 3/4/18 2:58 AM, Karl Rupp wrote:
Hi all,

I'm getting increasingly concerned about SELL not being a subclass of 
AIJ. As such, we have to deal with all these fallback operations now, 
whereas as a subclass of AIJ we could just selectively make use of the 
SELL format where we really benefit from it. "Use AIJ by default 
unless we have something optimized for SELL" is just much more 
appropriate for the few use cases of SELL than the current "SELL has 
to implement everything and usually this means to manually convert 
back to AIJ".

If there are no objections I'd like to clean this up. (Subclassing AIJ 
was unfortunately not available at the time Hong started his great 
work on SELL)

Best regards,
Karli



On 03/03/2018 07:52 AM, Richard Tran Mills wrote:
Resurrecting a few weeks old thread:

Stefano, did you get around to coding something up to do an automatic 
conversion to SeqAIJ for operations unsupported by SELL format? I did 
some hacking the other day to try to get PCGAMG to use SELL inside 
the smoothers and this turns out to be way more complicated than I'd 
like and very bug prone (I haven't found all of mine, anyway). I 
think it may be preferable to be able to pass a SELL matrix to PCGAMG 
and have an internal conversion happen in the SELL matrix to AIJ 
format for doing the MatPtAP and LU solves. Support for this would 
certainly make it easier for users in a lot other cases as well, and 
might make the use of SELL much more likely. If no one has already 
done some work on this, I'll take a stab at it.

--Richard

On Mon, Feb 12, 2018 at 10:04 AM, Richard Tran Mills <rtmi...@anl.gov 
<mailto:rtmi...@anl.gov>> wrote:

    On Mon, Feb 12, 2018 at 8:47 AM, Smith, Barry F. <bsm...@mcs.anl.gov
    <mailto:bsm...@mcs.anl.gov>> wrote:



        > On Feb 12, 2018, at 10:25 AM, Stefano Zampini 
<stefano.zamp...@gmail.com <mailto:stefano.zamp...@gmail.com>> wrote:
        >
        > Barry,
        >
        > for sure Amat,Pmat is the right approach; however, with 
complicated user codes, we are not always in control of having a 
different Jacobian matrix.
        > Since Mat*SELL does not currently support any 
preconditioning except PCSOR and PCJACOBI, we ask the user to put 
codes like
        >
        > if (type is SELL)
        >  create two matrices (and maybe modify the code in many 
other parts)
        > else
        >   ok with the previous code

            I don't disagree with what you are saying and am not opposed
        to the proposed work.

            Perhaps we need to do a better job with making the mat,pmat
        approach simpler or better documented so more people use it
        naturally in their applications.


    I wrote some code like that in some of the Jacobian/function
    routines in PFLOTRAN to experiment with MATSELL, and it works, but
    looks and feels pretty hacky. And if I wanted to support it for all
    of the different systems that PFLOTRAN can model, then I'd have to
    reproduce that it in many different Jacobian and function evaluation
    routines. I also don't like that it makes it awkward to play with
    the many combinations of matrix types and preconditioners that PETSc
    allows: The above pseudocode should really say "if (type is SELL)
    and (preconditioner is not PCSOR or PCJACOBI)". I do think that
    Amat,Pmat is a good approach in many situations, but it's easy to
    construct scenarios in which it falls short.

    In some situations, what I'd like to have happen is what Stefano is
    talking about, with an automatic conversion to AIJ happening if SELL
    doesn't support an operation. But, ideally, I think this sort of
    implicit format conversion shouldn't be something hard-coded into
    the workings of SELL. Instead, there should be some general
    mechanism by which PETSc recognizes that a particular operation is
    unsupported for a given matrix format, and then it can (optionally)
    copy/convert to a different matrix type (probably default to AIJ,
    but it shouldn't have to be AIJ) that supports the operation. This
    sort of implicit data rearrangement game may actually become more
    important if future computer architectures strongly prefer different
    data layouts different types of operations (though let's not get
    ahead of ourselves).

    --Richard


             Barry

         >
         > Just my two cents.
         >
         >
         > 2018-02-12 19:10 GMT+03:00 Smith, Barry F.
        <bsm...@mcs.anl.gov <mailto:bsm...@mcs.anl.gov>>:
         >
         >
         > > On Feb 12, 2018, at 9:59 AM, Stefano Zampini
        <stefano.zamp...@gmail.com <mailto:stefano.zamp...@gmail.com>>
        wrote:
         > >
         > > FYI, I just checked and MatSOR_*SELL does not use any
        vectorized instruction.
         > > Why just not converting to SeqAIJ, factor and then use the
        AIJ implementation for MatSolve for the moment?
         >
         >   Why not use the mat, pmat feature of the solvers to pass in
        both matrices and have the solvers handle using two formats
        simultaneously instead of burdening the MatSELL code with tons
        of special code for automatically converting to AIJ for 
solvers etc?
         >
         >
         > >
         > > 2018-02-12 18:06 GMT+03:00 Stefano Zampini
        <stefano.zamp...@gmail.com <mailto:stefano.zamp...@gmail.com>>:
         > >
         > >
         > > 2018-02-12 17:36 GMT+03:00 Jed Brown <j...@jedbrown.org
        <mailto:j...@jedbrown.org>>:
         > > Karl Rupp <r...@iue.tuwien.ac.at
        <mailto:r...@iue.tuwien.ac.at>> writes:
         > >
         > > > Hi Stefano,
         > > >
         > > >> Is there any plan to write code for native ILU/ICC etc
        for SeqSELL, at least to have BJACOBI in parallel?
         > > >
         > > > (imho) ILU/ICC is a pain to do with SeqSELL. Point-Jacobi
        should be
         > > > possible, yes. SELL is really just tailored to MatMults
        and a pain for
         > > > anything that is not very similar to a MatMult...
         > >
         > > There is already MatSOR_*SELL. MatSolve_SeqSELL wouldn't
        be any harder.
         > > I think it would be acceptable to convert to SeqAIJ,
        factor, and convert
         > > the factors back to SELL.
         > >
         > > Yes, this was my idea. Today I have started coding
        something. I'll push the branch whenever I have anything working
         > >
         > >
         > >
         > > --
         > > Stefano
         > >
         > >
         > >
         > > --
         > > Stefano
         >
         >
         >
         >
         > --
         > Stefano