2018-03-03 9:52 GMT+03:00 Richard Tran Mills <[email protected]>: > 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? >
No, I just added a faster conversion routine to SeqAIJ and support for sequential MUMPS (which requires COO storage) and SUPERLU (the PETSc interface already has a spot in the data structure for a copy of the matrix) here https://bitbucket.org/petsc/petsc/branch/stefano_zampini/feature-factor-sell > 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 <[email protected]> > wrote: > >> On Mon, Feb 12, 2018 at 8:47 AM, Smith, Barry F. <[email protected]> >> wrote: >> >>> >>> >>> > On Feb 12, 2018, at 10:25 AM, Stefano Zampini < >>> [email protected]> 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. <[email protected]>: >>> > >>> > >>> > > On Feb 12, 2018, at 9:59 AM, Stefano Zampini < >>> [email protected]> 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 < >>> [email protected]>: >>> > > >>> > > >>> > > 2018-02-12 17:36 GMT+03:00 Jed Brown <[email protected]>: >>> > > Karl Rupp <[email protected]> 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 >>> >>> >> > -- Stefano
