<div dir="ltr"><div>Dave convinced me that it's so relatively easy to use MatPreallocator that we might as well do it here even if we're going to do something more general later. Here's a draft of how it looks in 1D - most of the changes in the MR are just putting the guts of the matrix assembly in a separate function so we can call it twice. The part that actually uses the MatPreallocator API is very contained and so doesn't seem like it would be difficult to refactor.<br></div><div> <a href="https://gitlab.com/petsc/petsc/-/merge_requests/4769">https://gitlab.com/petsc/petsc/-/merge_requests/4769</a></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Am Mi., 12. Jan. 2022 um 14:53 Uhr schrieb Patrick Sanan <<a href="mailto:patrick.sanan@gmail.com">patrick.sanan@gmail.com</a>>:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><div dir="ltr"><div>Thanks a lot, all! So given that there's still some debate about whether we should even use MATPREALLOCATOR or  a better integration of that hash logic , as in Issue 852, I'll proceed with simply aping what DMDA does (with apologies for all this code duplication).<br></div><div><br></div><div>One thing I had missed, which I just added, is respect for DMSetMatrixPreallocation() / -dm_preallocate_only . Now, when that's activated, the previous behavior should be recovered, but by default it'll assemble a matrix which can be used for coloring, as is the main objective of this thread.<br></div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">Am Mi., 12. Jan. 2022 um 08:53 Uhr schrieb Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>>:<br></div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex"><br>
  Actually given the Subject of this email "Finite difference approximation of Jacobian" what I suggested is A complete story anyways for Patrick since the user cannot provide numerical values (of course Patrick could write a general new hash matrix type and then use it with zeros but that is asking a lot of him).<br>
<br>
  Regarding Fande's needs. One could rig things so that later one could "flip" back the matrix to again use the hasher for setting values when the contacts change. <br>
<br>
> On Jan 12, 2022, at 2:48 AM, Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> wrote:<br>
> <br>
> <br>
> <br>
>> On Jan 12, 2022, at 2:22 AM, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> wrote:<br>
>> <br>
>> Because if a user jumps right into MatSetValues() and then wants to use the result, it better have the values. It's true that DM preallocation normally just "inserts zeros" and thus matches what MATPREALLOCATOR does.<br>
> <br>
>  I am not saying the user jumps right into MatSetValues(). I am saying they do a "value" free assembly to have the preallocation set up and then they do the first real assembly; hence very much JUST a refactorization of MATPREALLOCATE. So the user is "preallocation-aware" but does not have to do anything difficult to determine their pattern analytically.<br>
> <br>
>  Of course handling the values also is fine and simplifies user code and the conceptual ideas (since preallocation ceases to exist as a concept for users). I have no problem with skipping the "JUST a refactorization of MATPREALLOCATE code" but for Patrick's current pressing need I think a refactorization of MATPREALLOCATE would be fastest to develop hence I suggested it.<br>
> <br>
>  I did not look at Jed's branch but one way to eliminate the preallocation part completely would be to have something like MATHASH and then when the first assembly is complete do a MatConvert to AIJ or whatever.  The conversion could be hidden inside the assemblyend of the hash so the user need not be aware that something different was done the first time through.  Doing it this way would easily allow multiple MATHASH* implementations (written by different people) that would compete on speed and memory usage.  So there could be MatSetType() and a new MatSetInitialType() where the initial type would be automatically set to some MATHASH if preallocation info is not provided.<br>
> <br>
>> <br>
>> Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> writes:<br>
>> <br>
>>> Why does it need to handle values? <br>
>>> <br>
>>>> On Jan 12, 2022, at 12:43 AM, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank">jed@jedbrown.org</a>> wrote:<br>
>>>> <br>
>>>> I agree with this and even started a branch jed/mat-hash in (yikes!) 2017. I think it should be default if no preallocation functions are called. But it isn't exactly the MATPREALLOCATOR code because it needs to handle values too. Should not be a lot of code and will essentially remove this FAQ and one of the most irritating subtle aspects of new codes using PETSc matrices.<br>
>>>> <br>
>>>> <a href="https://petsc.org/release/faq/#assembling-large-sparse-matrices-takes-a-long-time-what-can-i-do-to-make-this-process-faster-or-matsetvalues-is-so-slow-what-can-i-do-to-speed-it-up" rel="noreferrer" target="_blank">https://petsc.org/release/faq/#assembling-large-sparse-matrices-takes-a-long-time-what-can-i-do-to-make-this-process-faster-or-matsetvalues-is-so-slow-what-can-i-do-to-speed-it-up</a><br>
>>>> <br>
>>>> Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank">bsmith@petsc.dev</a>> writes:<br>
>>>> <br>
>>>>> I think the MATPREALLOCATOR as a MatType is a cumbersome strange thing and would prefer it was just functionality that Mat provided directly; for example MatSetOption(mat, preallocator_mode,true); matsetvalues,... MatAssemblyBegin/End. Now the matrix has its proper nonzero structure of whatever type the user set initially, aij, baij, sbaij, .... And the user can now use it efficiently.<br>
>>>>> <br>
>>>>> Barry<br>
>>>>> <br>
>>>>> So turning on the option just swaps out temporarily the operations for MatSetValues and AssemblyBegin/End to be essentially those in MATPREALLOCATOR. The refactorization should take almost no time and would be faster than trying to rig dmstag to use MATPREALLOCATOR as is.<br>
>>>>> <br>
>>>>> <br>
>>>>>> On Jan 11, 2022, at 9:43 PM, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>> wrote:<br>
>>>>>> <br>
>>>>>> On Tue, Jan 11, 2022 at 12:09 PM Patrick Sanan <<a href="mailto:patrick.sanan@gmail.com" target="_blank">patrick.sanan@gmail.com</a> <mailto:<a href="mailto:patrick.sanan@gmail.com" target="_blank">patrick.sanan@gmail.com</a>>> wrote:<br>
>>>>>> Working on doing this incrementally, in progress here: <a href="https://gitlab.com/petsc/petsc/-/merge_requests/4712" rel="noreferrer" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/4712</a> <<a href="https://gitlab.com/petsc/petsc/-/merge_requests/4712" rel="noreferrer" target="_blank">https://gitlab.com/petsc/petsc/-/merge_requests/4712</a>><br>
>>>>>> <br>
>>>>>> This works in 1D for AIJ matrices, assembling a matrix with a maximal number of zero entries as dictated by the stencil width (which is intended to be very very close to what DMDA would do if you <br>
>>>>>> associated all the unknowns with a particular grid point, which is the way DMStag largely works under the hood).<br>
>>>>>> <br>
>>>>>> Dave, before I get into it, am I correct in my understanding that MATPREALLOCATOR would be better here because you would avoid superfluous zeros in the sparsity pattern,<br>
>>>>>> because this routine wouldn't have to assemble the Mat returned by DMCreateMatrix()?<br>
>>>>>> <br>
>>>>>> Yes, here is how it works. You throw in all the nonzeros you come across. Preallocator is a hash table that can check for duplicates. At the end, it returns the sparsity pattern.<br>
>>>>>> <br>
>>>>>> Thanks,<br>
>>>>>> <br>
>>>>>>   Matt<br>
>>>>>> <br>
>>>>>> If this seems like a sane way to go, I will continue to add some more tests (in particular periodic BCs not tested yet) and add the code for 2D and 3D. <br>
>>>>>> <br>
>>>>>> <br>
>>>>>> <br>
>>>>>> Am Mo., 13. Dez. 2021 um 20:17 Uhr schrieb Dave May <<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a> <mailto:<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a>>>:<br>
>>>>>> <br>
>>>>>> <br>
>>>>>> On Mon, 13 Dec 2021 at 20:13, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a> <mailto:<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>>> wrote:<br>
>>>>>> On Mon, Dec 13, 2021 at 1:52 PM Dave May <<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a> <mailto:<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a>>> wrote:<br>
>>>>>> On Mon, 13 Dec 2021 at 19:29, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a> <mailto:<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>>> wrote:<br>
>>>>>> On Mon, Dec 13, 2021 at 1:16 PM Dave May <<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a> <mailto:<a href="mailto:dave.mayhem23@gmail.com" target="_blank">dave.mayhem23@gmail.com</a>>> wrote:<br>
>>>>>> <br>
>>>>>> <br>
>>>>>> On Sat 11. Dec 2021 at 22:28, Matthew Knepley <<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a> <mailto:<a href="mailto:knepley@gmail.com" target="_blank">knepley@gmail.com</a>>> wrote:<br>
>>>>>> On Sat, Dec 11, 2021 at 1:58 PM Tang, Qi <<a href="mailto:tangqi@msu.edu" target="_blank">tangqi@msu.edu</a> <mailto:<a href="mailto:tangqi@msu.edu" target="_blank">tangqi@msu.edu</a>>> wrote:<br>
>>>>>> Hi,<br>
>>>>>> Does anyone have comment on finite difference coloring with DMStag? We are using DMStag and TS to evolve some nonlinear equations implicitly. It would be helpful to have the coloring Jacobian option with that.<br>
>>>>>> <br>
>>>>>> Since DMStag produces the Jacobian connectivity,<br>
>>>>>> <br>
>>>>>> This is incorrect.<br>
>>>>>> The DMCreateMatrix implementation for DMSTAG only sets the number of nonzeros (very inaccurately). It does not insert any zero values and thus the nonzero structure is actually not defined. <br>
>>>>>> That is why coloring doesn’t work.<br>
>>>>>> <br>
>>>>>> Ah, thanks Dave.<br>
>>>>>> <br>
>>>>>> Okay, we should fix that.It is perfectly possible to compute the nonzero pattern from the DMStag information.<br>
>>>>>> <br>
>>>>>> Agreed. The API for DMSTAG is complete enough to enable one to<br>
>>>>>> loop over the cells, and for all quantities defined on the cell (centre, face, vertex), <br>
>>>>>> insert values into the appropriate slot in the matrix. <br>
>>>>>> Combined with MATPREALLOCATOR, I believe a compact and readable<br>
>>>>>> code should be possible to write for the preallocation (cf DMDA).<br>
>>>>>> <br>
>>>>>> I think the only caveat with the approach of using all quantities defined on the cell is <br>
>>>>>> It may slightly over allocate depending on how the user wishes to impose the boundary condition,<br>
>>>>>> or slightly over allocate for says Stokes where there is no pressure-pressure coupling term.<br>
>>>>>> <br>
>>>>>> Yes, and would not handle higher order stencils.I think the overallocating is livable for the first imeplementation.<br>
>>>>>> <br>
>>>>>> <br>
>>>>>> Sure, but neither does DMDA. <br>
>>>>>> <br>
>>>>>> The user always has to know what they are doing and set the stencil width accordingly.<br>
>>>>>> I actually had this point listed in my initial email (and the stencil growth issue when using FD for nonlinear problems),<br>
>>>>>> however I deleted it as all the same issue exist in DMDA and no one complains (at least not loudly) :D<br>
>>>>>> <br>
>>>>>> <br>
>>>>>> <br>
>>>>>> <br>
>>>>>> Thanks,<br>
>>>>>> <br>
>>>>>>   Matt<br>
>>>>>> <br>
>>>>>> Thanks,<br>
>>>>>> Dave<br>
>>>>>> <br>
>>>>>> <br>
>>>>>> Paging Patrick :)<br>
>>>>>> <br>
>>>>>> Thanks,<br>
>>>>>> <br>
>>>>>>  Matt<br>
>>>>>> <br>
>>>>>> Thanks,<br>
>>>>>> Dave<br>
>>>>>> <br>
>>>>>> <br>
>>>>>> you can use -snes_fd_color_use_mat. It has many options. Here is an example of us using that:<br>
>>>>>> <br>
>>>>>> <a href="https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex19.c#L898" rel="noreferrer" target="_blank">https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex19.c#L898</a> <<a href="https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex19.c#L898" rel="noreferrer" target="_blank">https://gitlab.com/petsc/petsc/-/blob/main/src/snes/tutorials/ex19.c#L898</a>><br>
>>>>>> <br>
>>>>>> Thanks,<br>
>>>>>> <br>
>>>>>>   Matt<br>
>>>>>> <br>
>>>>>> Thanks,<br>
>>>>>> Qi<br>
>>>>>> <br>
>>>>>> <br>
>>>>>>> On Oct 15, 2021, at 3:07 PM, Jorti, Zakariae via petsc-users <<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a> <mailto:<a href="mailto:petsc-users@mcs.anl.gov" target="_blank">petsc-users@mcs.anl.gov</a>>> wrote:<br>
>>>>>>> <br>
>>>>>>> Hello,<br>
>>>>>>> <br>
>>>>>>> Does the Jacobian approximation using coloring and finite differencing of the function evaluation work in DMStag?   <br>
>>>>>>> Thank you.<br>
>>>>>>> Best regards,<br>
>>>>>>> <br>
>>>>>>> Zakariae <br>
>>>>>> <br>
>>>>>> <br>
>>>>>> <br>
>>>>>> -- <br>
>>>>>> 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<br>
>>>>>> <br>
>>>>>> <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a> <<a href="http://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">http://www.cse.buffalo.edu/~knepley/</a>><br>
>>>>>> <br>
>>>>>> <br>
>>>>>> -- <br>
>>>>>> 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<br>
>>>>>> <br>
>>>>>> <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a> <<a href="http://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">http://www.cse.buffalo.edu/~knepley/</a>><br>
>>>>>> <br>
>>>>>> <br>
>>>>>> -- <br>
>>>>>> 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<br>
>>>>>> <br>
>>>>>> <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a> <<a href="http://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">http://www.cse.buffalo.edu/~knepley/</a>><br>
>>>>>> <br>
>>>>>> <br>
>>>>>> -- <br>
>>>>>> 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<br>
>>>>>> <br>
>>>>>> <a href="https://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">https://www.cse.buffalo.edu/~knepley/</a> <<a href="http://www.cse.buffalo.edu/~knepley/" rel="noreferrer" target="_blank">http://www.cse.buffalo.edu/~knepley/</a>><br>
> <br>
<br>
</blockquote></div>
</blockquote></div>