<div dir="auto"><div>I've worked on a few different codes doing matrix assembly on GPU independently of petsc. In all instances to plug into petsc all I need are the device CSR pointers and some guarantee they don't move around (my first try without setpreallocation on CPU I saw the value array pointer move after the first solve). It would also be nice to have a guarantee there aren't any unnecessary copies since memory constraints are always a concern.</div><div dir="auto"><br></div><div dir="auto">Here I call</div><div dir="auto">MatCreateSeqAIJCUSPARSE</div><div dir="auto">MatSeqAIJSetPreallocationCSR (filled using a preexisting CSR on host using the correct index arrays and zeros for values)</div><div dir="auto">MatSeqAIJGetCSRAndMemType (grab the allocated device CSR pointers and use those directly)</div><div dir="auto"><br></div><div dir="auto">Then in the Jacobian evaluation routine I fill that CSR directly with no calls to MatSetValues, just </div><div dir="auto"><br></div><div dir="auto">MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);</div><div dir="auto"> MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);</div><div dir="auto"><br></div><div dir="auto">after to put it in the correct state.</div><div dir="auto"><br></div><div dir="auto">In this code to fill the CSR coefficients, each GPU thread gets one row and fills it. No race conditions to contend with. Technically I'm duplicating some computations (a given dof could fill its own row and column) but this is much faster than the linear solver anyway.</div><div dir="auto"><br></div><div dir="auto">Other mesh based codes did GPU assembly using either coloring or mutexes, but still just need the CSR value array to fill.</div><div dir="auto"><br><br><div class="gmail_quote" dir="auto"><div dir="ltr" class="gmail_attr">On Fri, Jan 6, 2023, 9:44 PM Junchao Zhang <<a href="mailto:junchao.zhang@gmail.com">junchao.zhang@gmail.com</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div dir="ltr"><div dir="ltr"><br><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Jan 6, 2023 at 7:35 PM Mark Lohry <<a href="mailto:mlohry@gmail.com" target="_blank" rel="noreferrer">mlohry@gmail.com</a>> wrote:<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>Well, I think it's a moderately crazy idea unless it's less painful to implement than I'm thinking. Is there a use case for a mixed device system where one petsc executable might be addressing both a HIP and CUDA device beyond some frankenstein test system somebody cooked up? In all my code I implicitly assume I have either have one host with one device or one host with zero devices. I guess you can support these weird scenarios, but why? Life is hard enough supporting one device compiler with one host compiler.<br></div><div><br></div><div>Many thanks Junchao -- with combinations of SetPreallocation I was able to grab allocated pointers out of petsc. Now I have all the jacobian construction on device with no copies.<br></div></div></blockquote><div>Hi, Mark, could you say a few words about how you assemble matrices on GPUs?  We ported MatSetValues like routines to GPUs but did not continue this approach since we have to resolve data races between GPU threads.</div><div> </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> </div></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Fri, Jan 6, 2023 at 12:27 AM Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank" rel="noreferrer">bsmith@petsc.dev</a>> wrote:<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><div><br></div>  So Jed's "everyone" now consists of "no one" and Jed can stop complaining that "everyone" thinks it is a bad idea.<div><br></div><div><br><div><br><blockquote type="cite"><div>On Jan 5, 2023, at 11:50 PM, Junchao Zhang <<a href="mailto:junchao.zhang@gmail.com" target="_blank" rel="noreferrer">junchao.zhang@gmail.com</a>> wrote:</div><br><div><div dir="ltr"><div dir="ltr"><br><br></div><br><div class="gmail_quote"><div dir="ltr" class="gmail_attr">On Thu, Jan 5, 2023 at 10:32 PM Barry Smith <<a href="mailto:bsmith@petsc.dev" target="_blank" rel="noreferrer">bsmith@petsc.dev</a>> wrote:<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>
<br>
> On Jan 5, 2023, at 3:42 PM, Jed Brown <<a href="mailto:jed@jedbrown.org" target="_blank" rel="noreferrer">jed@jedbrown.org</a>> wrote:<br>
> <br>
> Mark Adams <<a href="mailto:mfadams@lbl.gov" target="_blank" rel="noreferrer">mfadams@lbl.gov</a>> writes:<br>
> <br>
>> Support of HIP and CUDA hardware together would be crazy, <br>
> <br>
> I don't think it's remotely crazy. libCEED supports both together and it's very convenient when testing on a development machine that has one of each brand GPU and simplifies binary distribution for us and every package that uses us. Every day I wish PETSc could build with both simultaneously, but everyone tells me it's silly.<br>
<br>
  Not everyone at all; just a subset of everyone. Junchao is really the hold-out :-)<br></blockquote><div>I am not, instead I think we should try (I fully agree it can ease binary distribution).  But satish needs to install such a machine first :)</div><div>There are issues out of our control if we want to mix GPUs in execution.  For example, how to do VexAXPY on a cuda vector and a hip vector? Shall we do it on the host? Also, there are no gpu-aware MPI implementations supporting messages between cuda memory and hip memory.</div><blockquote class="gmail_quote" style="margin:0px 0px 0px 0.8ex;border-left:1px solid rgb(204,204,204);padding-left:1ex">
<br>
  I just don't care about "binary packages" :-); I think they are an archaic and bad way of thinking about code distribution (but yes the alternatives need lots of work to make them flawless, but I think that is where the work should go in the packaging world.)<br>
<br>
   I go further and think one should be able to automatically use a CUDA vector on a HIP device as well, it is not hard in theory but requires thinking about how we handle classes and subclasses a little to make it straightforward; or perhaps Jacob has fixed that also?</blockquote><div> </div></div></div>
</div></blockquote></div><br></div></div></blockquote></div>
</blockquote></div></div>
</blockquote></div></div></div>