<div dir="ltr"><div class="gmail_quote"><div dir="ltr">On Fri, Jul 27, 2018 at 4:25 AM Mark Olesen <<a href="mailto:Mark.Olesen@esi-group.com">Mark.Olesen@esi-group.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">Hi Barry,<br>
<br>
Thanks for the information. One of the major takeaways from talking to <br>
various people at the workshop in June was about the imperative for <br>
proper pre-allocation sizing.<br>
As I started looking into how the matrices are being assembled, it is <br>
quickly apparently that I need to go off and determine the number of <br>
non-zeroes everywhere. After which I would use this information for <br>
petsc for provide good preallocation values and use the matrix set <br>
routines to set/copy the values. The hurdle, from there, to just sending <br>
off the appropriate CSR chunks for on-diagonal and off-diagonal doesn't <br>
seem so much higher. Unfortunately it appears that I'll need to approach <br>
it more indirectly.<br></blockquote><div><br></div><div>While its true that proper preallocation is essential, I want to point out that</div><div>computing the nonzero pattern usually has a negligible cost. If you just run through</div><div>your assembly routine twice, the first time leaving out the value computation, its</div><div>the easiest route.</div><div><br></div><div>Also, Lisandro has now given us really nice hash table support. I use this to calculate the nonzero</div><div>pattern in several places in PETSc.</div><div><br></div><div> Thanks,</div><div><br></div><div> Matt</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
Thanks,<br>
/mark<br>
<br>
<br>
<br>
On 07/24/18 18:02, Smith, Barry F. wrote:<br>
> <br>
> Mark,<br>
> <br>
> I think you are over-optimizing your matrix assembly leading to complicated, fragile code. Better just to create the matrix and use MatSetValues() to set values into the matrix and not to work directly with the various sparse matrix data structures. If you wish to work directly with the sparse matrix data structures then you are largely on your own and need to figure out yourself how to use them. Plus you will only get a small benefit time wise in going the much more complicated route.<br>
> <br>
> Barry<br>
> <br>
> <br>
>> On Jul 24, 2018, at 6:52 AM, Mark Olesen <<a href="mailto:Mark.Olesen@esi-group.com" target="_blank">Mark.Olesen@esi-group.com</a>> wrote:<br>
>><br>
>> I'm still at the beginning phase of looking at PETSc and accordingly have some beginner questions. My apologies if these are FAQs, but I didn't find much to address these specific questions.<br>
>><br>
>> My simulation matrices are sparse, and will generally (but not always) be generated in parallel. There is currently no conventional internal storage format (something like a COO variant), but lets just assume that I have CSR format for the moment.<br>
>><br>
>> I would like the usual combination of convenience and high efficiency, but efficiency (speed, memory) is the main criterion.<br>
>><br>
>> For serial, MatCreateSeqAIJWithArrays() looks like the thing to be using. It would provide a very thin wrapper around my CSR matrix without much additional allocation. The only extra allocation appears to be a precompute column range size per row (ilen) instead of doing it on-the-fly. If my matrix is actually to be considered symmetric, then use MatCreateSeqSBAIJWithArrays() instead.<br>
>> This all seems pretty clear.<br>
>><br>
>><br>
>> For parallel, MatCreateMPIAIJWithSplitArrays() appears to be the equivalent for efficiency, but I also read the note discouraging its use, which I fully appreciate. It also leads neatly into my question. I obviously will have fairly ready access to my on-processor portions of the matrix, but collecting the information for the off-processor portions is required. What would a normal or recommended approach look like?<br>
>><br>
>> For example,<br>
>> ====<br>
>> Mat A = MatCreateSeqAIJWithArrays() to wrap the local CSR.<br>
>><br>
>> Mat B = MatCreateSeqAIJ(). Do some preallocation for num non-zeroes, use MatSetValues() to fill in. Need extra garray[] as linear lookup for the global column numbers of B.<br>
>><br>
>> Or as an alternative, calculate the off-diagonal as a CSR by hand and use Mat B = MatCreateSeqAIJWithArrays() to wrap it.<br>
>><br>
>> Finally,<br>
>> Use MatCreateMPIAIJWithSeqAIJ() to produce the full matrix.<br>
>><br>
>> Assuming that I used MatCreateSeqAIJWithArrays() to create both the A and B matrices, then they both hold a shallow copy of my own storage.<br>
>> In MatCreateSeqAIJWithArrays(), I can't really tell what happens to the A matrix. For the B matrix, it appears that its column entries are changed to be those of the global columns and its data values are handed off to another MatCreateSeqAIJ() as the off-diagonal bits. The original B matrix is tagged up to avoid any deletion, and the shallow copied part is tagged to be deleted. If I follow this properly, this implies that if I was managing the storage of the original B matrix myself, I now have double deletion?<br>
>><br>
>> I would have expected something like this instead (around line 3431 of mpiaij.c in master):<br>
>><br>
>> /* Retain original memory management */<br>
>> bnew->singlemalloc = b->singlemalloc;<br>
>> bnew->free_a = b->free_a;<br>
>> bnew->free_ij = b->free_ij;<br>
>><br>
>> /* B arrays are shared by Bnew */<br>
>> b->singlemalloc = PETSC_FALSE;<br>
>> b->free_a = PETSC_FALSE;<br>
>> b->free_ij = PETSC_FALSE;<br>
>> ierr = MatDestroy(&B);CHKERRQ(ierr);<br>
>><br>
>><br>
>> Have I gone off in completely the wrong direction here?<br>
>> Is there a better method of approaching this?<br>
>><br>
>> Cheers,<br>
>> /mark<br>
> <br>
<br>
-- <br>
Dr Mark OLESEN<br>
Principal Engineer, ESI-OpenCFD<br>
ESI GmbH | Einsteinring 24 | 85609 Munich | GERMANY<br>
Mob. +49 171 9710 149<br>
<a href="http://www.openfoam.com" rel="noreferrer" target="_blank">www.openfoam.com</a> | <a href="http://www.esi-group.com" rel="noreferrer" target="_blank">www.esi-group.com</a> | <a href="mailto:mark.olesen@esi-group.com" target="_blank">mark.olesen@esi-group.com</a><br>
</blockquote></div><br clear="all"><div><br></div>-- <br><div dir="ltr" class="gmail_signature" data-smartmail="gmail_signature"><div dir="ltr"><div><div dir="ltr"><div>What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead.<br>-- Norbert Wiener</div><div><br></div><div><a href="http://www.caam.rice.edu/~mk51/" target="_blank">https://www.cse.buffalo.edu/~knepley/</a><br></div></div></div></div></div></div>