<div class="gmail_quote">On Fri, Feb 10, 2012 at 19:33, Dmitry Karpeev <span dir="ltr"><<a href="mailto:karpeev@mcs.anl.gov">karpeev@mcs.anl.gov</a>></span> wrote:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">
<div>How can there be a single callback, when the two DMs could be of different types?</div><div>I guess it could be done using the double-dispatch mechanism: </div>

<div>DMComputeOperators(DM dmJ, DM dmP, Mat *J, Mat *P) would dispatch based on the types of the two DMs, </div><div>and, if no such method exists, defaulting internally to calling</div><div>DMComputeMatrix(DM dm, Mat *A) for both J and P.</div>
</blockquote><div><br></div><div>Or call the user's function which just assembles into the P matrix.</div><div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">


<div>I still think that DMComputeMatrix() should be different from DMCreateMatrix() -- the latter just preallocates --</div><div>eliminating the need for DMSetPreallocateOnly().</div></blockquote><div><br></div><div>How many times do we have to explain that DMCreateMatrix() *NEVER* fills in useful values? It has two capabilities: preallocate (just specify number of nonzeros per row) and also setting the column indices. It does both by default because that is better for performance profiling and better to catch setting values outside of the allocated sparsity pattern at the first incident instead of whenever a row overflows. Under no circumstances does it assemble physics and under no circumstances does DMComputeJacobian() allocate a matrix.</div>
<div> </div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><div><br></div><div>If the two DMs are of the same type and, say, share a mesh, then the dual dispatch can take advantage of</div>


<div>a corresponding implementation.</div></blockquote></div><br><div>The user should always be able to get the callback to do both.</div>