<div dir="ltr">Hi Barry,<div><br></div><div>thanks for your help. I am going to have a look at the code right away, I am curious to understand better how PETSc works. </div><div><br></div><div>I have in the mean time tried a different approach, modifying DMCreateInterpolation_DA_2D_Q0 to provide the mapping (B) between two dms which are not one the refinement of the other. It looks like it is working. If that can be useful, I will clean it up and mail it to you.</div><div><br></div><div>Thanks again for your help,</div><div><br></div><div>Gianluca</div><div><br></div><div><br></div><div><br></div><div><br></div><div><br></div><div> </div></div><div class="gmail_extra"><br><div class="gmail_quote">2015-11-12 12:25 GMT-08:00 Barry Smith <span dir="ltr"><<a href="mailto:bsmith@mcs.anl.gov" target="_blank">bsmith@mcs.anl.gov</a>></span>:<br><blockquote class="gmail_quote" style="margin:0 0 0 .8ex;border-left:1px #ccc solid;padding-left:1ex"> I forgot the attached code<br>
<div class="HOEnZb"><div class="h5"><br>
> On Nov 12, 2015, at 2:24 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
><br>
><br>
> This is kind of tedious to get right. Plus it unfortunately will not work for all possible partitionings<br>
><br>
> Issues:<br>
><br>
> 1) The MatSetValues() works with global number but PETSc global numbering is per process so with DMDA does not normally match the "natural" ordering. You are trying to use global numbering in the natural ordering. To get it to work you need to use local numbering and MatSetValuesLocal()<br>
><br>
> 2) it is next to impossible to debug as written. So what I did to help debug is to put values using the local numbering into the two vectors f and x and print those vectors. This makes it easy to see when values are put in the wrong place.<br>
><br>
> 3) VecView() for DMDA prints out the vectors in the global natural ordering so it is easy to see if the vectors have the correct values in the correct locations. MatView() for DMDA however prints it only in the PETSc ordering by process so one needs to manual translate to make sure the matrix is correct.<br>
><br>
> Anyways I've attached your code with small Mx=5 and My=4 it runs correctly with 1,2 and 4 processes here is the output<br>
><br>
><br>
> $ petscmpiexec -n 1 ./ex5<br>
> Vec Object: 1 MPI processes<br>
> type: seq<br>
> Vec Object:Vec_0x84000000_0 1 MPI processes<br>
> type: mpi<br>
> Process [0]<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 1.<br>
> 1.<br>
> 1.<br>
> 0.<br>
> 0.<br>
> 1.<br>
> 1.<br>
> 1.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> Vec Object: 1 MPI processes<br>
> type: seq<br>
> Vec Object:Vec_0x84000000_1 1 MPI processes<br>
> type: mpi<br>
> Process [0]<br>
> 1.<br>
> 1.<br>
> 1.<br>
> 1.<br>
> 1.<br>
> 1.<br>
> Mat Object: 1 MPI processes<br>
> type: seqaij<br>
> row 0:<br>
> row 1:<br>
> row 2:<br>
> row 3:<br>
> row 4:<br>
> row 5:<br>
> row 6: (0, 1.)<br>
> row 7: (1, 1.)<br>
> row 8: (2, 1.)<br>
> row 9:<br>
> row 10:<br>
> row 11: (3, 1.)<br>
> row 12: (4, 1.)<br>
> row 13: (5, 1.)<br>
> row 14:<br>
> row 15:<br>
> row 16:<br>
> row 17:<br>
> row 18:<br>
> row 19:<br>
> ~/Src/petsc/test-dir (master=) arch-double<br>
> $ petscmpiexec -n 2 ./ex5<br>
> Vec Object: 2 MPI processes<br>
> type: mpi<br>
> Vec Object:Vec_0x84000004_0 2 MPI processes<br>
> type: mpi<br>
> Process [0]<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 1.<br>
> 1.<br>
> 2.<br>
> 0.<br>
> 0.<br>
> 1.<br>
> Process [1]<br>
> 1.<br>
> 2.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> Vec Object: 2 MPI processes<br>
> type: mpi<br>
> Vec Object:Vec_0x84000004_1 2 MPI processes<br>
> type: mpi<br>
> Process [0]<br>
> 1.<br>
> 1.<br>
> 2.<br>
> 1.<br>
> Process [1]<br>
> 1.<br>
> 2.<br>
> Mat Object: 2 MPI processes<br>
> type: mpiaij<br>
> row 0:<br>
> row 1:<br>
> row 2:<br>
> row 3:<br>
> row 4: (0, 1.)<br>
> row 5: (1, 1.)<br>
> row 6:<br>
> row 7: (2, 1.)<br>
> row 8: (3, 1.)<br>
> row 9:<br>
> row 10:<br>
> row 11:<br>
> row 12:<br>
> row 13:<br>
> row 14: (4, 2.)<br>
> row 15:<br>
> row 16: (5, 2.)<br>
> row 17:<br>
> row 18:<br>
> row 19:<br>
> ~/Src/petsc/test-dir (master=) arch-double<br>
> $ petscmpiexec -n 4 ./ex5<br>
> Vec Object: 4 MPI processes<br>
> type: mpi<br>
> Vec Object:Vec_0x84000004_0 4 MPI processes<br>
> type: mpi<br>
> Process [0]<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> Process [1]<br>
> 1.<br>
> 1.<br>
> 2.<br>
> 0.<br>
> Process [2]<br>
> 0.<br>
> 3.<br>
> 3.<br>
> 4.<br>
> 0.<br>
> 0.<br>
> Process [3]<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> Vec Object: 4 MPI processes<br>
> type: mpi<br>
> Vec Object:Vec_0x84000004_1 4 MPI processes<br>
> type: mpi<br>
> Process [0]<br>
> 1.<br>
> 1.<br>
> Process [1]<br>
> 2.<br>
> Process [2]<br>
> 3.<br>
> 3.<br>
> Process [3]<br>
> 4.<br>
> Mat Object: 4 MPI processes<br>
> type: mpiaij<br>
> row 0:<br>
> row 1:<br>
> row 2:<br>
> row 3:<br>
> row 4: (0, 1.)<br>
> row 5: (1, 1.)<br>
> row 6:<br>
> row 7:<br>
> row 8: (2, 2.)<br>
> row 9:<br>
> row 10:<br>
> row 11: (3, 3.)<br>
> row 12: (4, 3.)<br>
> row 13:<br>
> row 14:<br>
> row 15:<br>
> row 16: (5, 4.)<br>
> row 17:<br>
> row 18:<br>
> row 19:<br>
><br>
><br>
> It does NOT run correctly with 3 processes<br>
><br>
> $ petscmpiexec -n 3 ./ex5<br>
> [1]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>
> [1]PETSC ERROR: Argument out of range<br>
> [1]PETSC ERROR: Local index 2 too large 1 (max) at 0<br>
> [1]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>
> [1]PETSC ERROR: Petsc Development GIT revision: v3.6.2-1539-gf82f7f7 GIT Date: 2015-11-09 20:26:06 -0600<br>
> [1]PETSC ERROR: ./ex5 on a arch-double named <a href="http://visitor098-088.wl.anl-external.org" rel="noreferrer" target="_blank">visitor098-088.wl.anl-external.org</a> by barrysmith Thu Nov 12 14:18:56 2015<br>
> [1]PETSC ERROR: Configure options --download-hwloc --download-hypre --download-mpich --with-afterimage PETSC_ARCH=arch-double<br>
> [1]PETSC ERROR: #1 ISLocalToGlobalMappingApply() line 423 in /Users/barrysmith/Src/PETSc/src/vec/is/utils/isltog.c<br>
> [1]PETSC ERROR: #2 MatSetValuesLocal() line 2020 in /Users/barrysmith/Src/PETSc/src/mat/interface/matrix.c<br>
> [1]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------<br>
> [1]PETSC ERROR: Argument out of range<br>
> [1]PETSC ERROR: Local index 2 too large 1 (max) at 0<br>
> [1]PETSC ERROR: See <a href="http://www.mcs.anl.gov/petsc/documentation/faq.html" rel="noreferrer" target="_blank">http://www.mcs.anl.gov/petsc/documentation/faq.html</a> for trouble shooting.<br>
> [1]PETSC ERROR: Petsc Development GIT revision: v3.6.2-1539-gf82f7f7 GIT Date: 2015-11-09 20:26:06 -0600<br>
> [1]PETSC ERROR: ./ex5 on a arch-double named <a href="http://visitor098-088.wl.anl-external.org" rel="noreferrer" target="_blank">visitor098-088.wl.anl-external.org</a> by barrysmith Thu Nov 12 14:18:56 2015<br>
> [1]PETSC ERROR: Configure options --download-hwloc --download-hypre --download-mpich --with-afterimage PETSC_ARCH=arch-double<br>
> [1]PETSC ERROR: #3 ISLocalToGlobalMappingApply() line 423 in /Users/barrysmith/Src/PETSc/src/vec/is/utils/isltog.c<br>
> [1]PETSC ERROR: #4 VecSetValuesLocal() line 1058 in /Users/barrysmith/Src/PETSc/src/vec/vec/interface/rvector.c<br>
> Vec Object: 3 MPI processes<br>
> type: mpi<br>
> Vec Object:Vec_0x84000004_0 3 MPI processes<br>
> type: mpi<br>
> Process [0]<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 1.<br>
> 2.<br>
> Process [1]<br>
> 2.<br>
> 0.<br>
> 0.<br>
> 1.<br>
> 2.<br>
> 2.<br>
> 0.<br>
> 0.<br>
> Process [2]<br>
> 0.<br>
> 0.<br>
> 0.<br>
> 0.<br>
> Vec Object: 3 MPI processes<br>
> type: mpi<br>
> Vec Object:Vec_0x84000004_1 3 MPI processes<br>
> type: mpi<br>
> Process [0]<br>
> 1.<br>
> 2.<br>
> Process [1]<br>
> 0.<br>
> 1.<br>
> Process [2]<br>
> 4.<br>
> 0.<br>
> Mat Object: 3 MPI processes<br>
> type: mpiaij<br>
> row 0:<br>
> row 1:<br>
> row 2:<br>
> row 3: (0, 1.)<br>
> row 4:<br>
> row 5: (1, 1.)<br>
> row 6:<br>
> row 7:<br>
> row 8:<br>
> row 9:<br>
> row 10: (2, 2.)<br>
> row 11: (3, 2.)<br>
> row 12: (3, 2.)<br>
> row 13:<br>
> row 14:<br>
> row 15:<br>
> row 16:<br>
> row 17:<br>
> row 18:<br>
> row 19:<br>
><br>
> The reason is that in this case the DMDAs are decomposed into three strips, for x the strips are xi = 0,1 then xi = 2,3 then xi= 4<br>
> for f the strips are fi=0, fi=1, fi=2 so there is no way to get a consistent local numbering for both x and f at the same time with 3 strips. So like the DMDA interpolation routines your code can only work for certain decompositions. Forget what I said about lx and ly I don't think that is relevant for what you are trying to do.<br>
><br>
> Barry<br>
><br>
><br>
>> On Nov 12, 2015, at 12:23 PM, Gianluca Meneghello <<a href="mailto:gianmail@gmail.com">gianmail@gmail.com</a>> wrote:<br>
>><br>
>> Hi Barry,<br>
>><br>
>> sorry, but I still cannot make it. I guess what I need is something similar to MatRestrict/MatInterpolate (and B is something similar to what is created from DMCreateInterpolation, except for the fact that the nonzero entries are distributed differently).<br>
>><br>
>> Am I mistaken? Is there any example I could start from?<br>
>><br>
>> Thanks again,<br>
>><br>
>> Gianluca<br>
>><br>
>> On Wed, Nov 11, 2015 at 10:19 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
>> DMDAGetOwnershipRanges<br>
>><br>
>>> On Nov 11, 2015, at 10:47 PM, Gianluca Meneghello <<a href="mailto:gianmail@gmail.com">gianmail@gmail.com</a>> wrote:<br>
>>><br>
>>> Hi,<br>
>>><br>
>>> thanks for the very quick reply.<br>
>>><br>
>>> One more question: is there a way to get the lx and ly from the first dm and use them (modified) for the second dm? DMDAGetInfo does not seem to provide this information.<br>
>>><br>
>>> Thanks again for your help<br>
>>><br>
>>> Gianluca<br>
>>><br>
>>> On Wed, Nov 11, 2015 at 8:12 PM, Barry Smith <<a href="mailto:bsmith@mcs.anl.gov">bsmith@mcs.anl.gov</a>> wrote:<br>
>>><br>
>>> When you create the 2 DM you must be set the lx, ly arguments (the ones you set to 0) in your code carefully to insure that the vectors for the 2 DM you create have compatible layout to do the matrix vector product.<br>
>>><br>
>>> You can run a very small problem with 2 processors and printing out the vectors to see the layout to make sure you get it correct.<br>
>>><br>
>>> The 2 DM don't have any magically way of knowing that you created another DMDA and want it to be compatible automatically.<br>
>>><br>
>>> Barry<br>
>>><br>
>>> DMDACreate2d(PETSC_COMM_WORLD , DM_BOUNDARY_GHOSTED , DM_BOUNDARY_GHOSTED , DMDA_STENCIL_BOX ,<br>
>>> Mx , Nx , PETSC_DECIDE , PETSC_DECIDE , 1 , 0 , 0 , 0 , &dax);<br>
>>> DMDACreate2d(PETSC_COMM_WORLD , DM_BOUNDARY_NONE , DM_BOUNDARY_NONE , DMDA_STENCIL_BOX ,<br>
>>> Mx-2*bs , Nx-2*bs , PETSC_DECIDE , PETSC_DECIDE , 1 , 0 , 0 , 0 , &daf);<br>
>>><br>
>>>> On Nov 11, 2015, at 10:05 PM, Gianluca Meneghello <<a href="mailto:gianmail@gmail.com">gianmail@gmail.com</a>> wrote:<br>
>>>><br>
>>>> Hi,<br>
>>>><br>
>>>> I am trying to do something apparently really simple but with no success.<br>
>>>><br>
>>>> I need to perform a matrix-vector multiplication x = B f , where the length of x is bigger than the length of f (or viceversa). Thus, B cannot be created using DMCreateMatrix.<br>
>>>><br>
>>>> Both x and f are obtained from different DMs, the smaller covering only a subdomain of the larger. The application is to apply a control f to a system, e.g. \dot{x} = A x + B f.<br>
>>>><br>
>>>> The problem is, when running on more than one core, the vector x is not organized as I would expect (everything works fine on a single core).<br>
>>>><br>
>>>> I attach a short example where B is intended to map f to the interior of x.<br>
>>>><br>
>>>> mpirun -n 1 ./test -draw_pause -1 works fine while<br>
>>>> mpirun -n 2 ./test -draw_pause -1 shows the problem<br>
>>>><br>
>>>> I have not found any example with non square matrices in the src folder, any help is very welcome.<br>
>>>><br>
>>>> Thanks for your time,<br>
>>>><br>
>>>> Gianluca<br>
>>>> <test.cpp><br>
>>><br>
>>><br>
>><br>
>><br>
><br>
<br>
</div></div></blockquote></div><br></div>