ex5.cxx
Description: Binary data
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> >>> >>> >> >> >
