<div dir="ltr"><div>Hi Everyone (especially Herr Doktor Hong Zhang),<br></div><div><br></div><div>It sure took me a while to get around to it, but I coded up (in branch rmills/feature-add-mataijsell) the skeleton of what I'm currently calling the 'AIJSELL' matrix type, which is a "meta" matrix type that inherits from AIJ but maintains a "shadow" MATSEQSELL copy of the matrix for use when that is appropriate. This works similarly to how AIJMKL works, insofar as the PetscObjectState is tracked and used to determine when the SELL "shadow" copy needs to be updated before an operation like MatMult() is applied. (Note that what I have so far only does MatMult, but it should be easy to add other operations once I've verified that my implementation is working.) I decided on the approach of using a MATSEQSELL instance inside an AIJ-derived matrix type because I still want to be able to support standalone SELL matrices when there really isn't a need for any of the AIJ-only operations.</div><div><br></div><div>My AIJSELL implementation seems to work fine in many cases, but I've spent several hours chasing some bugs that crop up in certain cases. I now think that the bugs are not a problem with my AIJSELL implementation (though I'd welcome another pair of eyes to help me verify this), but are lurking somewhere in the SELL code and haven't been hit before because AIJSELL makes it easier to exercise the SELL class on more problems.</div><div><br></div><div>I've been using the well-loved SNES ex19 example (that is used for 'make check') for some of my debugging. I've found that everything works great when I size the problem grid such that the resulting matrices have a row count that is evenly divisible by 8 (when using double precision scalars). So</div><div><br></div><div><span style="font-family:monospace,monospace">./ex19 -ksp_type gmres -pc_type none -snes_monitor -ksp_monitor -da_grid_x 4 -da_grid_y 5 -mat_seqaij_type seqaijsell -snes_max_it 1 -ksp_monitor</span></div><div><br></div><div>which results in an 80x80 Jacobian (since there are 4 degrees of freedom per grid point) runs fine, but</div><div><br></div><div><span style="font-family:monospace,monospace">./ex19 -ksp_type gmres -pc_type none -snes_monitor -ksp_monitor
-da_grid_x 3 -da_grid_y 5 -mat_seqaij_type seqaijsell -snes_max_it 1
-ksp_monitor</span></div><div><br></div><div>results in either a segfault or the GMRES(30) solve terminating at 10,000 iterations because the residual norm just keeps bouncing up and down.</div><div><br></div><div>Sticking some MatView() calls in, it appears that the MatMult() results inside the KSP solve are incorrect because the SELL copy of the matrix generated by my MatConvert() call is incorrect when the row count isn't evenly divisible by eight. In the case of SNES ex19, MatConvert_Basic() is currently used because the block size is greater than unity. So I think something is going wrong either in the MatSetValues() or (my guess) MatAssemblyBegin()/End().</div><div><br></div><div>Poking around with a debugger, I see that the column index array ("colidx") that is part of the SEQSELLHEADER ends up with nonsense values, and the segfaults I see tend to occur when indexing into the x vector with this.</div><div><br></div><div>I'm not familiar enough with the internals of SELL to look further into debugging this without a lot of effort. Hong, can you help me look at this?</div><div><br></div><div>Best regards,</div><div>Richard<br></div></div><div class="gmail_extra"><br><div class="gmail_quote">On Tue, Mar 6, 2018 at 3:59 AM, Karl Rupp <span dir="ltr"><<a href="mailto:rupp@iue.tuwien.ac.at" target="_blank">rupp@iue.tuwien.ac.at</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><span class=""><br>
> Karl, are you thinking of a matrix subclass that has everything that an <br>
> AIJ matrix does, but also keeps a SELL copy around for operations like <br>
> MatMult()? Would it make sense to just keep another Mat inside (like how <br>
> MPIAIJ keeps multiple Mat instances) that *is* of type MATSELL, that <br>
> gets built/updated as needed? Would this involve carrying too much <br>
> baggage around, associated with a complete Mat instance?<br>
<br>
</span>What I have in mind is to put the SELL datastructures into a A->spptr, <br>
very much like you did for AIJMKL.<br>
<span class=""><br>
<br>
> I like the idea <br>
> of having a MATSELL type available that is lean (insofar as not having <br>
> storing an AIJ matrix) for those cases when a user really knows that the <br>
> AIJ stuff will not be needed. But maybe it makes sense to be able to use <br>
> that inside another matrix class. Perhaps we could have something, <br>
> called, say, MATAIJMUTABLE that uses AIJ but might also create copies in <br>
> SELL (or other formats, potentially) when appropriate -- perhaps based <br>
> on a some performance model indicating which format is fastest for <br>
> MatMult() or whatever.<br>
<br>
</span>The actual overhead of also storing a SELL datastructure in terms of <br>
memory footprint is at most 2x. When you keep in mind that extra memory <br>
during the matrix assembly is also needed in the stash, then the impact <br>
on overall memory consumption is about 50 percent extra. Given that SELL <br>
is a performance optimization for SpMV (and hence the SELL datastructure <br>
is only populated if you call MatMult on the particular matrix), I'm not <br>
too worried about the increased memory footprint at this time.<br>
<span class=""><br>
<br>
> Having a class that's an AIJ but can also use SELL is more convenient <br>
> than adding a fallback to AIJ format inside MATSELL. I wonder if the <br>
> latter option might be preferable in some circumstances, however, <br>
> because it can avoid the extra memory footprint of also keeping the <br>
> matrix in AIJ -- maybe AIJ operations are rarely needed and the AIJ <br>
> conversion can just happen on a lazy basis.<br>
<br>
</span>I advocate an 'optimize as needed' approach here. Let's first make SELL <br>
compatible with the full range of AIJ operations and preconditioners.<br>
<br>
Best regards,<br>
Karli<br>
<span class=""><br>
<br>
<br>
> <br>
> --Richard<br>
> <br>
> On 3/4/18 2:58 AM, Karl Rupp wrote:<br>
>> Hi all,<br>
>><br>
>> I'm getting increasingly concerned about SELL not being a subclass of <br>
>> AIJ. As such, we have to deal with all these fallback operations now, <br>
>> whereas as a subclass of AIJ we could just selectively make use of the <br>
>> SELL format where we really benefit from it. "Use AIJ by default <br>
>> unless we have something optimized for SELL" is just much more <br>
>> appropriate for the few use cases of SELL than the current "SELL has <br>
>> to implement everything and usually this means to manually convert <br>
>> back to AIJ".<br>
>><br>
>> If there are no objections I'd like to clean this up. (Subclassing AIJ <br>
>> was unfortunately not available at the time Hong started his great <br>
>> work on SELL)<br>
>><br>
>> Best regards,<br>
>> Karli<br>
>><br>
>><br>
>><br>
>> On 03/03/2018 07:52 AM, Richard Tran Mills wrote:<br>
</span><span class="">>>> Resurrecting a few weeks old thread:<br>
>>><br>
>>> Stefano, did you get around to coding something up to do an automatic <br>
>>> conversion to SeqAIJ for operations unsupported by SELL format? I did <br>
>>> some hacking the other day to try to get PCGAMG to use SELL inside <br>
>>> the smoothers and this turns out to be way more complicated than I'd <br>
>>> like and very bug prone (I haven't found all of mine, anyway). I <br>
>>> think it may be preferable to be able to pass a SELL matrix to PCGAMG <br>
>>> and have an internal conversion happen in the SELL matrix to AIJ <br>
>>> format for doing the MatPtAP and LU solves. Support for this would <br>
>>> certainly make it easier for users in a lot other cases as well, and <br>
>>> might make the use of SELL much more likely. If no one has already <br>
>>> done some work on this, I'll take a stab at it.<br>
>>><br>
>>> --Richard<br>
>>><br>
>>> On Mon, Feb 12, 2018 at 10:04 AM, Richard Tran Mills <<a href="mailto:rtmills@anl.gov">rtmills@anl.gov</a> <br>
</span><span class="">>>> <mailto:<a href="mailto:rtmills@anl.gov">rtmills@anl.gov</a>>> wrote:<br>
>>><br>
>>> On Mon, Feb 12, 2018 at 8:47 AM, Smith, Barry F. <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a><br>
</span><span class="">>>> <mailto:<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>>> wrote:<br>
>>><br>
>>><br>
>>><br>
>>> > On Feb 12, 2018, at 10:25 AM, Stefano Zampini <br>
</span><div><div class="h5">>>> <<a href="mailto:stefano.zampini@gmail.com">stefano.zampini@gmail.com</a> <mailto:<a href="mailto:stefano.zampini@gmail.com">stefano.zampini@gmail.<wbr>com</a>>> wrote:<br>
>>> ><br>
>>> > Barry,<br>
>>> ><br>
>>> > for sure Amat,Pmat is the right approach; however, with <br>
>>> complicated user codes, we are not always in control of having a <br>
>>> different Jacobian matrix.<br>
>>> > Since Mat*SELL does not currently support any <br>
>>> preconditioning except PCSOR and PCJACOBI, we ask the user to put <br>
>>> codes like<br>
>>> ><br>
>>> > if (type is SELL)<br>
>>> > create two matrices (and maybe modify the code in many <br>
>>> other parts)<br>
>>> > else<br>
>>> > ok with the previous code<br>
>>><br>
>>> I don't disagree with what you are saying and am not opposed<br>
>>> to the proposed work.<br>
>>><br>
>>> Perhaps we need to do a better job with making the mat,pmat<br>
>>> approach simpler or better documented so more people use it<br>
>>> naturally in their applications.<br>
>>><br>
>>><br>
>>> I wrote some code like that in some of the Jacobian/function<br>
>>> routines in PFLOTRAN to experiment with MATSELL, and it works, but<br>
>>> looks and feels pretty hacky. And if I wanted to support it for all<br>
>>> of the different systems that PFLOTRAN can model, then I'd have to<br>
>>> reproduce that it in many different Jacobian and function evaluation<br>
>>> routines. I also don't like that it makes it awkward to play with<br>
>>> the many combinations of matrix types and preconditioners that PETSc<br>
>>> allows: The above pseudocode should really say "if (type is SELL)<br>
>>> and (preconditioner is not PCSOR or PCJACOBI)". I do think that<br>
>>> Amat,Pmat is a good approach in many situations, but it's easy to<br>
>>> construct scenarios in which it falls short.<br>
>>><br>
>>> In some situations, what I'd like to have happen is what Stefano is<br>
>>> talking about, with an automatic conversion to AIJ happening if SELL<br>
>>> doesn't support an operation. But, ideally, I think this sort of<br>
>>> implicit format conversion shouldn't be something hard-coded into<br>
>>> the workings of SELL. Instead, there should be some general<br>
>>> mechanism by which PETSc recognizes that a particular operation is<br>
>>> unsupported for a given matrix format, and then it can (optionally)<br>
>>> copy/convert to a different matrix type (probably default to AIJ,<br>
>>> but it shouldn't have to be AIJ) that supports the operation. This<br>
>>> sort of implicit data rearrangement game may actually become more<br>
>>> important if future computer architectures strongly prefer different<br>
>>> data layouts different types of operations (though let's not get<br>
>>> ahead of ourselves).<br>
>>><br>
>>> --Richard<br>
>>><br>
>>><br>
>>> Barry<br>
>>><br>
>>> ><br>
>>> > Just my two cents.<br>
>>> ><br>
>>> ><br>
>>> > 2018-02-12 19:10 GMT+03:00 Smith, Barry F.<br>
</div></div>>>> <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a> <mailto:<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>>>:<br>
<span class="">>>> ><br>
>>> ><br>
>>> > > On Feb 12, 2018, at 9:59 AM, Stefano Zampini<br>
</span>>>> <<a href="mailto:stefano.zampini@gmail.com">stefano.zampini@gmail.com</a> <mailto:<a href="mailto:stefano.zampini@gmail.com">stefano.zampini@gmail.<wbr>com</a>>><br>
<span class="">>>> wrote:<br>
>>> > ><br>
>>> > > FYI, I just checked and MatSOR_*SELL does not use any<br>
>>> vectorized instruction.<br>
>>> > > Why just not converting to SeqAIJ, factor and then use the<br>
>>> AIJ implementation for MatSolve for the moment?<br>
>>> ><br>
>>> > Why not use the mat, pmat feature of the solvers to pass in<br>
>>> both matrices and have the solvers handle using two formats<br>
>>> simultaneously instead of burdening the MatSELL code with tons<br>
>>> of special code for automatically converting to AIJ for <br>
>>> solvers etc?<br>
>>> ><br>
>>> ><br>
>>> > ><br>
>>> > > 2018-02-12 18:06 GMT+03:00 Stefano Zampini<br>
</span>>>> <<a href="mailto:stefano.zampini@gmail.com">stefano.zampini@gmail.com</a> <mailto:<a href="mailto:stefano.zampini@gmail.com">stefano.zampini@gmail.<wbr>com</a>>>:<br>
<span class="im HOEnZb">>>> > ><br>
>>> > ><br>
>>> > > 2018-02-12 17:36 GMT+03:00 Jed Brown <<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a><br>
</span><span class="im HOEnZb">>>> <mailto:<a href="mailto:jed@jedbrown.org">jed@jedbrown.org</a>>>:<br>
>>> > > Karl Rupp <<a href="mailto:rupp@iue.tuwien.ac.at">rupp@iue.tuwien.ac.at</a><br>
</span><div class="HOEnZb"><div class="h5">>>> <mailto:<a href="mailto:rupp@iue.tuwien.ac.at">rupp@iue.tuwien.ac.at</a>><wbr>> writes:<br>
>>> > ><br>
>>> > > > Hi Stefano,<br>
>>> > > ><br>
>>> > > >> Is there any plan to write code for native ILU/ICC etc<br>
>>> for SeqSELL, at least to have BJACOBI in parallel?<br>
>>> > > ><br>
>>> > > > (imho) ILU/ICC is a pain to do with SeqSELL. Point-Jacobi<br>
>>> should be<br>
>>> > > > possible, yes. SELL is really just tailored to MatMults<br>
>>> and a pain for<br>
>>> > > > anything that is not very similar to a MatMult...<br>
>>> > ><br>
>>> > > There is already MatSOR_*SELL. MatSolve_SeqSELL wouldn't<br>
>>> be any harder.<br>
>>> > > I think it would be acceptable to convert to SeqAIJ,<br>
>>> factor, and convert<br>
>>> > > the factors back to SELL.<br>
>>> > ><br>
>>> > > Yes, this was my idea. Today I have started coding<br>
>>> something. I'll push the branch whenever I have anything working<br>
>>> > ><br>
>>> > ><br>
>>> > ><br>
>>> > > --<br>
>>> > > Stefano<br>
>>> > ><br>
>>> > ><br>
>>> > ><br>
>>> > > --<br>
>>> > > Stefano<br>
>>> ><br>
>>> ><br>
>>> ><br>
>>> ><br>
>>> > --<br>
>>> > Stefano<br>
>>><br>
>>><br>
>>><br>
> <br>
</div></div></blockquote></div><br></div>