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