<div><br><div class="gmail_quote"><div dir="auto">On Sun, 6 May 2018 at 17:52, Robert Speck <<a href="mailto:r.speck@fz-juelich.de">r.speck@fz-juelich.de</a>> wrote:<br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex">OK, thanks, I see. This is basically what happens in the poisson2d.py<br>
example, too, right?<br>
<br>
I tried it with the shell matrix (?) used in the poisson2d example and<br>
it works right away, but then I fail to see how to make use of the<br>
preconditioners for KSP (see my original message)..</blockquote><div dir="auto"><br></div><div dir="auto">If you want assemble something within the operator created by DMCreateMatrix(), take a look at this function</div><div dir="auto"><br></div><div dir="auto"><a href="http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetValuesStencil.html">http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatSetValuesStencil.html</a><br></div><div dir="auto"><br></div><div dir="auto">Study some of the examples listed at the bottom of the page. </div><div dir="auto"><br></div><div dir="auto">(But you will have to concern yourself with the spatial decomposition, just like in the poisson2d and Bratu example.)</div><div dir="auto"><br></div><div dir="auto">Thanks,</div><div dir="auto"> Dave</div><div dir="auto"><br></div><div dir="auto"><br></div><div dir="auto"><br></div><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"><br>
<br>
Thanks again!<br>
-Robert-<br>
<br>
<br>
On 06.05.18 16:52, Jed Brown wrote:<br>
> Robert Speck <<a href="mailto:r.speck@fz-juelich.de" target="_blank">r.speck@fz-juelich.de</a>> writes:<br>
> <br>
>> Thanks for your reply and help. Yes, this is going to be a PDE solver<br>
>> for structured grids. The first goal would be IDC (or Crank-Nicholson)<br>
>> for the heat equation, which would require both solving a linear system<br>
>> and application of the matrix.<br>
>><br>
>> The code I wrote for testing parallel matrix-vector multiplication can<br>
>> be found here:<br>
>> <a href="https://github.com/Parallel-in-Time/pySDC/blob/petsc/pySDC/playgrounds/PETSc/playground_matmult.py" rel="noreferrer" target="_blank">https://github.com/Parallel-in-Time/pySDC/blob/petsc/pySDC/playgrounds/PETSc/playground_matmult.py</a><br>
>><br>
>> Both vectors and matrix come from the DMDA, but I guess filling them is<br>
>> done in a wrong way? Or do I need to convert global/natural vectors to<br>
>> local ones somewhere?<br>
> <br>
> Global and Natural are not the same (see user's manual for details).<br>
> The matrix acts on a Global vector. See<br>
> petsc4py/demo/bratu3d/bratu3d.py for examples of efficiently setting<br>
> values (and computing residuals) using Global vectors. It should be<br>
> simpler/cleaner code than you currently have.<br>
> <br>
>><br>
>> Best<br>
>> -Robert-<br>
>><br>
>><br>
>><br>
>> On 06.05.18 14:44, Dave May wrote:<br>
>>> On Sun, 6 May 2018 at 10:40, Robert Speck <<a href="mailto:r.speck@fz-juelich.de" target="_blank">r.speck@fz-juelich.de</a><br>
>>> <mailto:<a href="mailto:r.speck@fz-juelich.de" target="_blank">r.speck@fz-juelich.de</a>>> wrote:<br>
>>><br>
>>> Hi!<br>
>>><br>
>>> I would like to do a matrix-vector multiplication (besides using linear<br>
>>> solvers and so on) with petsc4py. I took the matrix from this example<br>
>>><br>
>>> (<a href="https://bitbucket.org/petsc/petsc4py/src/master/demo/kspsolve/petsc-mat.py" rel="noreferrer" target="_blank">https://bitbucket.org/petsc/petsc4py/src/master/demo/kspsolve/petsc-mat.py</a>)<br>
>>><br>
>>><br>
>>> This example only produces a matrix. And from the code the matrix<br>
>>> produced is identical in serial or parallel. <br>
>>><br>
>>><br>
>>><br>
>>> and applied it to a PETSc Vector. All works well in serial, but in<br>
>>> parallel (in particular if ordering becomes relevant) the resulting<br>
>>> vector looks very different. <br>
>>><br>
>>><br>
>>> Given this, the way you defined the x vector in y = A x must be<br>
>>> different when run on 1 versus N mpi ranks. <br>
>>><br>
>>><br>
>>> Using the shell matrix of this example<br>
>>> (<a href="https://bitbucket.org/petsc/petsc4py/src/master/demo/poisson2d/poisson2d.py" rel="noreferrer" target="_blank">https://bitbucket.org/petsc/petsc4py/src/master/demo/poisson2d/poisson2d.py</a>)<br>
>>> helps, but then I cannot use matrix-based preconditioners for KSP<br>
>>> directly (right?). I also tried using DMDA for creating vectors and<br>
>>> matrix and for taking care of their ordering (which seems to be my<br>
>>> problem here), but that did not help either.<br>
>>><br>
>>> So, my question is this: How do I do easy parallel matrix-vector<br>
>>> multiplication with petsc4py in a way that allows me to use parallel<br>
>>> linear solvers etc. later on? I want to deal with spatial decomposition<br>
>>> as little as possible. <br>
>>><br>
>>><br>
>>> What's the context - are you solving a PDE?<br>
>>><br>
>>> Assuming you are using your own grid object (e.g. as you might have if<br>
>>> solving a PDE), and assuming you are not solving a 1D problem, you<br>
>>> actually have to "deal" with the spatial decomposition otherwise<br>
>>> performance could be quite terrible - even for something simple like a 5<br>
>>> point Laplacian on a structured grid in 2D<br>
>>><br>
>>> What data structures should I use? DMDA or<br>
>>> PETSc.Vec() and PETSc.Mat() or something else?<br>
>>><br>
>>><br>
>>> The mat vec product is not causing you a problem. Your issue appears to<br>
>>> be that you do not have a way to label entries in a vector in a<br>
>>> consistent manner.<br>
>>><br>
>>> What's the objective? Are you solving a PDE? If yes, structured grid? If<br>
>>> yes again, use the DMDA. It takes care of all the local-to-global and<br>
>>> global-to-local mapping you need.<br>
>>><br>
>>> Thanks,<br>
>>> Dave<br>
>>><br>
>>><br>
>>><br>
>>> Thanks!<br>
>>> -Robert-<br>
>>><br>
>>> --<br>
>>> Dr. Robert Speck<br>
>>> Juelich Supercomputing Centre<br>
>>> Institute for Advanced Simulation<br>
>>> Forschungszentrum Juelich GmbH<br>
>>> 52425 Juelich, Germany<br>
>>><br>
>>> Tel: +49 2461 61 1644<br>
>>> Fax: +49 2461 61 6656<br>
>>><br>
>>> Email: <a href="mailto:r.speck@fz-juelich.de" target="_blank">r.speck@fz-juelich.de</a> <mailto:<a href="mailto:r.speck@fz-juelich.de" target="_blank">r.speck@fz-juelich.de</a>><br>
>>> Website: <a href="http://www.fz-juelich.de/ias/jsc/speck_r" rel="noreferrer" target="_blank">http://www.fz-juelich.de/ias/jsc/speck_r</a><br>
>>> PinT: <a href="http://www.fz-juelich.de/ias/jsc/pint" rel="noreferrer" target="_blank">http://www.fz-juelich.de/ias/jsc/pint</a><br>
>>><br>
>>><br>
>>><br>
>>> ------------------------------------------------------------------------------------------------<br>
>>> ------------------------------------------------------------------------------------------------<br>
>>> Forschungszentrum Juelich GmbH<br>
>>> 52425 Juelich<br>
>>> Sitz der Gesellschaft: Juelich<br>
>>> Eingetragen im Handelsregister des Amtsgerichts Dueren Nr. HR B 3498<br>
>>> Vorsitzender des Aufsichtsrats: MinDir Dr. Karl Eugen Huthmacher<br>
>>> Geschaeftsfuehrung: Prof. Dr.-Ing. Wolfgang Marquardt (Vorsitzender),<br>
>>> Karsten Beneke (stellv. Vorsitzender), Prof. Dr.-Ing. Harald Bolt,<br>
>>> Prof. Dr. Sebastian M. Schmidt<br>
>>> ------------------------------------------------------------------------------------------------<br>
>>> ------------------------------------------------------------------------------------------------<br>
>>><br>
>><br>
>> -- <br>
>> Dr. Robert Speck<br>
>> Juelich Supercomputing Centre<br>
>> Institute for Advanced Simulation<br>
>> Forschungszentrum Juelich GmbH<br>
>> 52425 Juelich, Germany<br>
>><br>
>> Tel: +49 2461 61 1644<br>
>> Fax: +49 2461 61 6656<br>
>><br>
>> Email: <a href="mailto:r.speck@fz-juelich.de" target="_blank">r.speck@fz-juelich.de</a><br>
>> Website: <a href="http://www.fz-juelich.de/ias/jsc/speck_r" rel="noreferrer" target="_blank">http://www.fz-juelich.de/ias/jsc/speck_r</a><br>
>> PinT: <a href="http://www.fz-juelich.de/ias/jsc/pint" rel="noreferrer" target="_blank">http://www.fz-juelich.de/ias/jsc/pint</a><br>
<br>
-- <br>
Dr. Robert Speck<br>
Juelich Supercomputing Centre<br>
Institute for Advanced Simulation<br>
Forschungszentrum Juelich GmbH<br>
52425 Juelich, Germany<br>
<br>
Tel: +49 2461 61 1644<br>
Fax: +49 2461 61 6656<br>
<br>
Email: <a href="mailto:r.speck@fz-juelich.de" target="_blank">r.speck@fz-juelich.de</a><br>
Website: <a href="http://www.fz-juelich.de/ias/jsc/speck_r" rel="noreferrer" target="_blank">http://www.fz-juelich.de/ias/jsc/speck_r</a><br>
PinT: <a href="http://www.fz-juelich.de/ias/jsc/pint" rel="noreferrer" target="_blank">http://www.fz-juelich.de/ias/jsc/pint</a><br>
<br>
</blockquote></div></div>