Hi Barry, 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.
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. Thanks again for your help, Gianluca 2015-11-12 12:25 GMT-08:00 Barry Smith <[email protected]>: > I forgot the attached code > > > On Nov 12, 2015, at 2:24 PM, Barry Smith <[email protected]> wrote: > > > > > > This is kind of tedious to get right. Plus it unfortunately will not > work for all possible partitionings > > > > Issues: > > > > 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() > > > > 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. > > > > 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. > > > > 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 > > > > > > $ petscmpiexec -n 1 ./ex5 > > Vec Object: 1 MPI processes > > type: seq > > Vec Object:Vec_0x84000000_0 1 MPI processes > > type: mpi > > Process [0] > > 0. > > 0. > > 0. > > 0. > > 0. > > 0. > > 1. > > 1. > > 1. > > 0. > > 0. > > 1. > > 1. > > 1. > > 0. > > 0. > > 0. > > 0. > > 0. > > 0. > > Vec Object: 1 MPI processes > > type: seq > > Vec Object:Vec_0x84000000_1 1 MPI processes > > type: mpi > > Process [0] > > 1. > > 1. > > 1. > > 1. > > 1. > > 1. > > Mat Object: 1 MPI processes > > type: seqaij > > row 0: > > row 1: > > row 2: > > row 3: > > row 4: > > row 5: > > row 6: (0, 1.) > > row 7: (1, 1.) > > row 8: (2, 1.) > > row 9: > > row 10: > > row 11: (3, 1.) > > row 12: (4, 1.) > > row 13: (5, 1.) > > row 14: > > row 15: > > row 16: > > row 17: > > row 18: > > row 19: > > ~/Src/petsc/test-dir (master=) arch-double > > $ petscmpiexec -n 2 ./ex5 > > Vec Object: 2 MPI processes > > type: mpi > > Vec Object:Vec_0x84000004_0 2 MPI processes > > type: mpi > > Process [0] > > 0. > > 0. > > 0. > > 0. > > 0. > > 0. > > 1. > > 1. > > 2. > > 0. > > 0. > > 1. > > Process [1] > > 1. > > 2. > > 0. > > 0. > > 0. > > 0. > > 0. > > 0. > > Vec Object: 2 MPI processes > > type: mpi > > Vec Object:Vec_0x84000004_1 2 MPI processes > > type: mpi > > Process [0] > > 1. > > 1. > > 2. > > 1. > > Process [1] > > 1. > > 2. > > Mat Object: 2 MPI processes > > type: mpiaij > > row 0: > > row 1: > > row 2: > > row 3: > > row 4: (0, 1.) > > row 5: (1, 1.) > > row 6: > > row 7: (2, 1.) > > row 8: (3, 1.) > > row 9: > > row 10: > > row 11: > > row 12: > > row 13: > > row 14: (4, 2.) > > row 15: > > row 16: (5, 2.) > > row 17: > > row 18: > > row 19: > > ~/Src/petsc/test-dir (master=) arch-double > > $ petscmpiexec -n 4 ./ex5 > > Vec Object: 4 MPI processes > > type: mpi > > Vec Object:Vec_0x84000004_0 4 MPI processes > > type: mpi > > Process [0] > > 0. > > 0. > > 0. > > 0. > > 0. > > 0. > > Process [1] > > 1. > > 1. > > 2. > > 0. > > Process [2] > > 0. > > 3. > > 3. > > 4. > > 0. > > 0. > > Process [3] > > 0. > > 0. > > 0. > > 0. > > Vec Object: 4 MPI processes > > type: mpi > > Vec Object:Vec_0x84000004_1 4 MPI processes > > type: mpi > > Process [0] > > 1. > > 1. > > Process [1] > > 2. > > Process [2] > > 3. > > 3. > > Process [3] > > 4. > > Mat Object: 4 MPI processes > > type: mpiaij > > row 0: > > row 1: > > row 2: > > row 3: > > row 4: (0, 1.) > > row 5: (1, 1.) > > row 6: > > row 7: > > row 8: (2, 2.) > > row 9: > > row 10: > > row 11: (3, 3.) > > row 12: (4, 3.) > > row 13: > > row 14: > > row 15: > > row 16: (5, 4.) > > row 17: > > row 18: > > row 19: > > > > > > It does NOT run correctly with 3 processes > > > > $ petscmpiexec -n 3 ./ex5 > > [1]PETSC ERROR: --------------------- Error Message > -------------------------------------------------------------- > > [1]PETSC ERROR: Argument out of range > > [1]PETSC ERROR: Local index 2 too large 1 (max) at 0 > > [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html > for trouble shooting. > > [1]PETSC ERROR: Petsc Development GIT revision: v3.6.2-1539-gf82f7f7 > GIT Date: 2015-11-09 20:26:06 -0600 > > [1]PETSC ERROR: ./ex5 on a arch-double named > visitor098-088.wl.anl-external.org by barrysmith Thu Nov 12 14:18:56 2015 > > [1]PETSC ERROR: Configure options --download-hwloc --download-hypre > --download-mpich --with-afterimage PETSC_ARCH=arch-double > > [1]PETSC ERROR: #1 ISLocalToGlobalMappingApply() line 423 in > /Users/barrysmith/Src/PETSc/src/vec/is/utils/isltog.c > > [1]PETSC ERROR: #2 MatSetValuesLocal() line 2020 in > /Users/barrysmith/Src/PETSc/src/mat/interface/matrix.c > > [1]PETSC ERROR: --------------------- Error Message > -------------------------------------------------------------- > > [1]PETSC ERROR: Argument out of range > > [1]PETSC ERROR: Local index 2 too large 1 (max) at 0 > > [1]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html > for trouble shooting. > > [1]PETSC ERROR: Petsc Development GIT revision: v3.6.2-1539-gf82f7f7 > GIT Date: 2015-11-09 20:26:06 -0600 > > [1]PETSC ERROR: ./ex5 on a arch-double named > visitor098-088.wl.anl-external.org by barrysmith Thu Nov 12 14:18:56 2015 > > [1]PETSC ERROR: Configure options --download-hwloc --download-hypre > --download-mpich --with-afterimage PETSC_ARCH=arch-double > > [1]PETSC ERROR: #3 ISLocalToGlobalMappingApply() line 423 in > /Users/barrysmith/Src/PETSc/src/vec/is/utils/isltog.c > > [1]PETSC ERROR: #4 VecSetValuesLocal() line 1058 in > /Users/barrysmith/Src/PETSc/src/vec/vec/interface/rvector.c > > Vec Object: 3 MPI processes > > type: mpi > > Vec Object:Vec_0x84000004_0 3 MPI processes > > type: mpi > > Process [0] > > 0. > > 0. > > 0. > > 0. > > 0. > > 0. > > 1. > > 2. > > Process [1] > > 2. > > 0. > > 0. > > 1. > > 2. > > 2. > > 0. > > 0. > > Process [2] > > 0. > > 0. > > 0. > > 0. > > Vec Object: 3 MPI processes > > type: mpi > > Vec Object:Vec_0x84000004_1 3 MPI processes > > type: mpi > > Process [0] > > 1. > > 2. > > Process [1] > > 0. > > 1. > > Process [2] > > 4. > > 0. > > Mat Object: 3 MPI processes > > type: mpiaij > > row 0: > > row 1: > > row 2: > > row 3: (0, 1.) > > row 4: > > row 5: (1, 1.) > > row 6: > > row 7: > > row 8: > > row 9: > > row 10: (2, 2.) > > row 11: (3, 2.) > > row 12: (3, 2.) > > row 13: > > row 14: > > row 15: > > row 16: > > row 17: > > row 18: > > row 19: > > > > 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 > > 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. > > > > Barry > > > > > >> On Nov 12, 2015, at 12:23 PM, Gianluca Meneghello <[email protected]> > wrote: > >> > >> Hi Barry, > >> > >> 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). > >> > >> Am I mistaken? Is there any example I could start from? > >> > >> Thanks again, > >> > >> Gianluca > >> > >> On Wed, Nov 11, 2015 at 10:19 PM, Barry Smith <[email protected]> > wrote: > >> DMDAGetOwnershipRanges > >> > >>> On Nov 11, 2015, at 10:47 PM, Gianluca Meneghello <[email protected]> > wrote: > >>> > >>> Hi, > >>> > >>> thanks for the very quick reply. > >>> > >>> 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. > >>> > >>> Thanks again for your help > >>> > >>> Gianluca > >>> > >>> On Wed, Nov 11, 2015 at 8:12 PM, Barry Smith <[email protected]> > wrote: > >>> > >>> 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. > >>> > >>> 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. > >>> > >>> The 2 DM don't have any magically way of knowing that you created > another DMDA and want it to be compatible automatically. > >>> > >>> Barry > >>> > >>> DMDACreate2d(PETSC_COMM_WORLD , DM_BOUNDARY_GHOSTED , > DM_BOUNDARY_GHOSTED , DMDA_STENCIL_BOX , > >>> Mx , Nx , PETSC_DECIDE , PETSC_DECIDE , 1 , 0 , 0 , 0 , > &dax); > >>> DMDACreate2d(PETSC_COMM_WORLD , DM_BOUNDARY_NONE , > DM_BOUNDARY_NONE , DMDA_STENCIL_BOX , > >>> Mx-2*bs , Nx-2*bs , PETSC_DECIDE , PETSC_DECIDE , 1 , 0 , 0 , 0 , > &daf); > >>> > >>>> On Nov 11, 2015, at 10:05 PM, Gianluca Meneghello <[email protected]> > wrote: > >>>> > >>>> Hi, > >>>> > >>>> I am trying to do something apparently really simple but with no > success. > >>>> > >>>> 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. > >>>> > >>>> 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. > >>>> > >>>> 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). > >>>> > >>>> I attach a short example where B is intended to map f to the interior > of x. > >>>> > >>>> mpirun -n 1 ./test -draw_pause -1 works fine while > >>>> mpirun -n 2 ./test -draw_pause -1 shows the problem > >>>> > >>>> I have not found any example with non square matrices in the src > folder, any help is very welcome. > >>>> > >>>> Thanks for your time, > >>>> > >>>> Gianluca > >>>> <test.cpp> > >>> > >>> > >> > >> > > > >
