February 2011 - petsc-dev - mailman3.cels.anl.gov

Configure argument check
by Matthew Knepley 22 Feb '11

22 Feb '11

I believe Jed has suggested before that we multiplex this check over different architectures. This is also starting to bite me and I looked at MPI.py. I think just by changing args.petsc to args.<arch>.petsc we can achieve this. Of course, it will not work for the unfortunate poop that use make.inc, but we can clean that up later. Any objections to this change? Matt -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener

1 0

Overzealous guards in DMGetInjection_DA
by Jed Brown 21 Feb '11

21 Feb '11

ierr = DMDAGetInfo(dac,&dimc,&Mc,&Nc,&Pc,&mc,&nc,&pc,&dofc,&sc,&wrapc,&stc);CHKERRQ(ierr); ierr = DMDAGetInfo(daf,&dimf,&Mf,&Nf,&Pf,&mf,&nf,&pf,&doff,&sf,&wrapf,&stf);CHKERRQ(ierr); if (dimc != dimf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Dimensions of DMDA do not match %D %D",dimc,dimf);CHKERRQ(ierr); if (dofc != doff) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"DOF of DMDA do not match %D %D",dofc,doff);CHKERRQ(ierr); if (sc != sf) SETERRQ2(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil width of DMDA do not match %D %D",sc,sf);CHKERRQ(ierr); if (wrapc != wrapf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Periodic type different in two DMDAs");CHKERRQ(ierr); if (stc != stf) SETERRQ(((PetscObject)daf)->comm,PETSC_ERR_ARG_INCOMP,"Stencil type different in two DMDAs");CHKERRQ(ierr); if (Mc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in x direction"); if (dimc > 1 && Nc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in y direction"); if (dimc > 2 && Pc < 2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_OUTOFRANGE,"Coarse grid requires at least 2 points in z direction"); Several of these are not necessary. DMGetInjection_DA_3D: ierr = DMDAGetInfo(dac,0,&Mx,&My,&Mz,0,0,0,0,0,&pt,0);CHKERRQ(ierr); ierr = DMDAGetInfo(daf,0,&mx,&my,&mz,0,0,0,&dof,0,0,0);CHKERRQ(ierr); Seems silly to fail when not necessary. I'm inclined to move the guards into the specific implementation so that only the necessary amount of consistency needs to be checked. How about getting rid of DMGetInjection_DA entirely and just set the appropriate ops slot?

1 0

DA refine/coarsen with stencil_width=0 for nproc>1
by Lisandro Dalcin 21 Feb '11

21 Feb '11

Perhaps I'm missing something obvious, but DA refine/coarsen with stencil_width=0 for nproc>1 is failing, see the example below [dalcinl@nikola petsc4py-dev]$ cat test-da-rc.py import sys, petsc4py petsc4py.init(sys.argv) from petsc4py import PETSc OptDB = PETSc.Options() n = OptDB.getInt('n', 100) da = PETSc.DA().create([n], dof=1, stencil_width=0) rda = da.refine() cda = da.coarsen() [dalcinl@nikola petsc4py-dev]$ mpiexec -n 1 python test-da-rc.py [dalcinl@nikola petsc4py-dev]$ mpiexec -n 2 python test-da-rc.py Traceback (most recent call last): File "test-da-rc.py", line 10, in <module> rda = da.refine() Traceback (most recent call last): File "test-da-rc.py", line 10, in <module> rda = da.refine() File "DA.pyx", line 408, in petsc4py.PETSc.DA.refine (src/petsc4py.PETSc.c:115026) File "DA.pyx", line 408, in petsc4py.PETSc.DA.refine (src/petsc4py.PETSc.c:115026) petsc4py.PETSc.Errorpetsc4py.PETSc.Error: error code 60 [0] DMRefine() line 541 in /usr/local/petsc/dev/src/dm/interface/dm.c [0] DMRefine_DA() line 770 in /usr/local/petsc/dev/src/dm/impls/da/da.c [0] DMDARefineOwnershipRanges() line 668 in /usr/local/petsc/dev/src/dm/impls/da/da.c [0] Nonconforming object sizes [0] Could not find a compatible refined ownership range : error code 60 [1] DMRefine() line 541 in /usr/local/petsc/dev/src/dm/interface/dm.c [1] DMRefine_DA() line 770 in /usr/local/petsc/dev/src/dm/impls/da/da.c [1] DMDARefineOwnershipRanges() line 668 in /usr/local/petsc/dev/src/dm/impls/da/da.c [1] Nonconforming object sizes [1] Could not find a compatible refined ownership range -- Lisandro Dalcin --------------- CIMEC (INTEC/CONICET-UNL) Predio CONICET-Santa Fe Colectora RN 168 Km 472, Paraje El Pozo 3000 Santa Fe, Argentina Tel: +54-342-4511594 (ext 1011) Tel/Fax: +54-342-4511169

1 0

Re: [petsc-dev] vec/vec/examples/tests/ex37.c "Nest vector functionality. has bugs
by Jed Brown 20 Feb '11

20 Feb '11

On Sun, Feb 20, 2011 at 22:58, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > PetscObjectDecreaseReference() and have it > zero the pointer only if it got down to zero. > Using this across the board probably also requires making it beginner level. Part of the temptation of XXXDestroy() is that it is shorter (long names take longer to read so it helps for common things to have short names) and the explicit cast causes visual noise that obscures the actual argument of the function and prompts the reader to confirm that the cast is indeed safe. Removing the need for an explicit cast would reduce the noise, but I'm not sure it's a good idea. > > I really hate code that does XXXDestroy(x) then XXXXSomething(x); it is > fundamentally confusing > Very confusing without a comment.

1 0

Re: [petsc-dev] vec/vec/examples/tests/ex37.c "Nest vector functionality. has bugs
by Barry Smith 20 Feb '11

20 Feb '11

On Feb 20, 2011, at 2:07 PM, Jed Brown wrote: > On Sun, Feb 20, 2011 at 20:52, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > ierr = VecDestroy(a);CHKERRQ(ierr); ierr = VecDestroy(b);CHKERRQ(ierr); > > ierr = VecAssemblyBegin(X);CHKERRQ(ierr); > ierr = VecAssemblyEnd(X);CHKERRQ(ierr); > > ierr = VecMax( b, &index, &val );CHKERRQ(ierr); > PetscPrintf( PETSC_COMM_WORLD, "(max-b) = %f : index = %d \n", val, index ); > > b is destroyed and then used later. > > Well, sort of. The Vec is not actually destroyed because the VecNest X retains a reference. VecDestroy is called to relinquish ownership, not to actually free the memory. Although this is perhaps not very pretty, neither is calling PetscObjectDereference((PetscObject)b) for the same purpose. I think that waiting until the end of the function body to call VecDestroy is harder to reason about: when you read the function from the top down, ownership becomes ambiguous and you have to scan to the end of the function body to find out who keeps ownership. > > In this example, the max could easily be computed before creating X, but I guess we need a style guideline for other cases where we want to relinquish ownership, but still be able to inspect. The verbose choice is to VecDestroy and then immediately call VecNestGetSubVec (a low-level function that I would prefer people not need to call very often) or VecGetSubVector (the generic function, but you need index sets for that, not just PETSC_NULL) to get a non-ownership reference to the same thing you passed in. > > Or still use VecDestroy, but always include a comment about it just relinquishing ownership. I realize that usage conflicts with your plan to make XXDestroy() or XXDestroy_() or something zero the pointer (which by no means prevents stale references, but takes care of some of the obvious cases). I originally planned that PetscObjectDereference() was the thing to use in this case. From its manual page PetscObjectDereference - Indicates to any PetscObject that it is being referenced by one less PetscObject. This decreases the reference count for that object by one. but I now realize the name is not great. It implies that this particular reference is being removed. Could change it to PetscObjectDecreaseReference() and have it zero the pointer only if it got down to zero. I really hate code that does XXXDestroy(x) then XXXXSomething(x); it is fundamentally confusing Barry

1 0

Re: [petsc-dev] vec/vec/examples/tests/ex37.c "Nest vector functionality. has bugs
by Jed Brown 20 Feb '11

20 Feb '11

On Sun, Feb 20, 2011 at 20:52, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > ierr = VecDestroy(a);CHKERRQ(ierr); ierr = VecDestroy(b);CHKERRQ(ierr); > > ierr = VecAssemblyBegin(X);CHKERRQ(ierr); > ierr = VecAssemblyEnd(X);CHKERRQ(ierr); > > ierr = VecMax( b, &index, &val );CHKERRQ(ierr); > PetscPrintf( PETSC_COMM_WORLD, "(max-b) = %f : index = %d \n", val, index > ); > > b is destroyed and then used later. > Well, sort of. The Vec is not actually destroyed because the VecNest X retains a reference. VecDestroy is called to relinquish ownership, not to actually free the memory. Although this is perhaps not very pretty, neither is calling PetscObjectDereference((PetscObject)b) for the same purpose. I think that waiting until the end of the function body to call VecDestroy is harder to reason about: when you read the function from the top down, ownership becomes ambiguous and you have to scan to the end of the function body to find out who keeps ownership. In this example, the max could easily be computed before creating X, but I guess we need a style guideline for other cases where we want to relinquish ownership, but still be able to inspect. The verbose choice is to VecDestroy and then immediately call VecNestGetSubVec (a low-level function that I would prefer people not need to call very often) or VecGetSubVector (the generic function, but you need index sets for that, not just PETSC_NULL) to get a non-ownership reference to the same thing you passed in. Or still use VecDestroy, but always include a comment about it just relinquishing ownership. I realize that usage conflicts with your plan to make XXDestroy() or XXDestroy_() or something zero the pointer (which by no means prevents stale references, but takes care of some of the obvious cases).

1 0

vec/vec/examples/tests/ex37.c "Nest vector functionality. has bugs
by Barry Smith 20 Feb '11

20 Feb '11

ierr = VecDestroy(a);CHKERRQ(ierr); ierr = VecDestroy(b);CHKERRQ(ierr); ierr = VecAssemblyBegin(X);CHKERRQ(ierr); ierr = VecAssemblyEnd(X);CHKERRQ(ierr); ierr = VecMax( b, &index, &val );CHKERRQ(ierr); PetscPrintf( PETSC_COMM_WORLD, "(max-b) = %f : index = %d \n", val, index ); b is destroyed and then used later. Barry

1 0

Re: [petsc-dev] XXXDestroy() mistaken design in PETSc
by Barry Smith 19 Feb '11

19 Feb '11

On Feb 19, 2011, at 10:17 PM, Matthew Knepley wrote: > On Sat, Feb 19, 2011 at 10:08 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > What about > > extern PetscErrorCode VecDestroy_(Vec); > #define VecDestroy(a) (VecDestroy_(a) || (((a) = 0),0)) > > Not exactly PETSc style, but allows the switch without changing the API. > > Why not let VecDestroy_() take the pointer? Could do that. I think it is equivalent. > > Also, what about having a configure mode that goes back to the old interface (on by default), > but allow the new one to be configured. Might involve putting in #define stuff, or changing > files a little, but it seems cleaner. Won't that configure code also have to change all the calls to Destroy() to pass in the &? My goal was to not change the API (since you feared the universal breakage) hence not add the * to the API. And I hate to have two valid ways of doing the same thing VecFree(Vec*) (for new users) and VecDestroy(Vec) for old users. Barry > > Matt > > > Barry > > On Feb 15, 2011, at 4:47 PM, Barry Smith wrote: > > > > > In MPI one calls MPI_Comm_free(&comm) to allow the MPI implementation to set the pointer explicitly to 0 after the object is destroyed. > > > > In Petsc XXXDestroy() does not pass the pointer (because it seemed too unnatural to me in 1994) thus not allowing 0ing the pointer. > > > > Was this a bad design decision? Should it be revisited? > > > > Barry > > > > Two use cases > > > > 1) error detection when someone tries to reuse a freed object > > > > 2) when removing some objects from a data structure that will be used data one currently needs to do > > > > XXXXDestroy(mystruct->something);CHKERRQ(ierr); mystruct->something = 0; > > > > instead of the cleaner XXXDestroy(&mystruct->something);CHKERRQ(ierr); > > > > > -- > What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. > -- Norbert Wiener

1 0

Re: [petsc-dev] XXXDestroy() mistaken design in PETSc
by Matthew Knepley 19 Feb '11

19 Feb '11

On Sat, Feb 19, 2011 at 10:08 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > > What about > > extern PetscErrorCode VecDestroy_(Vec); > #define VecDestroy(a) (VecDestroy_(a) || (((a) = 0),0)) > > Not exactly PETSc style, but allows the switch without changing the API. Why not let VecDestroy_() take the pointer? Also, what about having a configure mode that goes back to the old interface (on by default), but allow the new one to be configured. Might involve putting in #define stuff, or changing files a little, but it seems cleaner. Matt > > Barry > > On Feb 15, 2011, at 4:47 PM, Barry Smith wrote: > > > > > In MPI one calls MPI_Comm_free(&comm) to allow the MPI implementation to > set the pointer explicitly to 0 after the object is destroyed. > > > > In Petsc XXXDestroy() does not pass the pointer (because it seemed too > unnatural to me in 1994) thus not allowing 0ing the pointer. > > > > Was this a bad design decision? Should it be revisited? > > > > Barry > > > > Two use cases > > > > 1) error detection when someone tries to reuse a freed object > > > > 2) when removing some objects from a data structure that will be used > data one currently needs to do > > > > XXXXDestroy(mystruct->something);CHKERRQ(ierr); mystruct->something = 0; > > > > instead of the cleaner XXXDestroy(&mystruct->something);CHKERRQ(ierr); > > -- What most experimenters take for granted before they begin their experiments is infinitely more interesting than any results to which their experiments lead. -- Norbert Wiener

1 0

Re: [petsc-dev] XXXDestroy() mistaken design in PETSc
by Barry Smith 19 Feb '11

19 Feb '11

What about extern PetscErrorCode VecDestroy_(Vec); #define VecDestroy(a) (VecDestroy_(a) || (((a) = 0),0)) Not exactly PETSc style, but allows the switch without changing the API. Barry On Feb 15, 2011, at 4:47 PM, Barry Smith wrote: > > In MPI one calls MPI_Comm_free(&comm) to allow the MPI implementation to set the pointer explicitly to 0 after the object is destroyed. > > In Petsc XXXDestroy() does not pass the pointer (because it seemed too unnatural to me in 1994) thus not allowing 0ing the pointer. > > Was this a bad design decision? Should it be revisited? > > Barry > > Two use cases > > 1) error detection when someone tries to reuse a freed object > > 2) when removing some objects from a data structure that will be used data one currently needs to do > > XXXXDestroy(mystruct->something);CHKERRQ(ierr); mystruct->something = 0; > > instead of the cleaner XXXDestroy(&mystruct->something);CHKERRQ(ierr);

1 0