August 2012 - petsc-dev - mailman3.cels.anl.gov

Patch review
by Jed Brown 24 Aug '12

24 Aug '12

http://petsc.cs.iit.edu/petsc/petsc-dev/rev/067c5dbc352e Just a couple requests regarding this patch: 1. Please split independent functionality into separate patches. This does three different things. 2. Please cut out extra whitespace like this and follow conventions about not including args in the declarations. @@ -273,6 +274,12 @@ PETSC_EXTERN PetscErrorCode KSPMonitorDefault(KSP,PetscInt,PetscReal,void *); PETSC_EXTERN PetscErrorCode KSPLSQRMonitorDefault(KSP,PetscInt,PetscReal,void *); PETSC_EXTERN PetscErrorCode KSPMonitorRange(KSP,PetscInt,PetscReal,void *); + + +PETSC_EXTERN PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy); +PETSC_EXTERN PetscErrorCode KSPMonitorDynamicToleranceDestroy(void **dummy); + + 3. This doesn't look like an improvement: @@ -163,7 +163,7 @@ /* scalar updates */ omega = xi2 / xi3; beta = - xi4 / sigma; - rho = PetscSqrtReal(PetscAbsScalar(xi1 - omega * xi2)); /* residual norm */ + rho = sqrt(fabs(xi1 - omega * xi2)); /* residual norm */ /* vector updates */ 4. Better to raise a PETSC_ERR_SUP than to just leave this fragment in a comment + This code is only correct for the real case... Need to modify for complex numbers... 5. If this thing works (which putting it in the public interface usually indicates) then it needs a man page. --- a/src/ksp/ksp/interface/iterativ.c +++ b/src/ksp/ksp/interface/iterativ.c @@ -330,6 +330,42 @@ PetscFunctionReturn(0); } +#undef __FUNCT__ +#define __FUNCT__ "KSPMonitorDynamicTolerance" +/* + A hack to using dynamic tolerence in preconditioner + */ +PetscErrorCode KSPMonitorDynamicTolerance(KSP ksp,PetscInt its,PetscReal fnorm,void *dummy) { 6. Two different new features and no tests! I'm a little concerned that the current review process is all ex post facto. It makes more work for reviewers and the history ends up looking confusing. I think we should work out a review process that can happen before things are pushed to petsc-dev, so that the authors of the patches can clean them up in response to comments. When done right, this ends up being less work for both the reviewers and the authors.

1 0

Re: [petsc-dev] Issues with DMComplexSetCone
by Matthew Knepley 24 Aug '12

24 Aug '12

On Fri, Aug 24, 2012 at 3:48 PM, Blaise Bourdin <bourdin(a)lsu.edu> wrote: > I am attaching a debugged version. > Checked in as DMComplex ex1f90.F Matt > Blaise > > > > On Aug 24, 2012, at 3:30 PM, Matthew Knepley <knepley(a)gmail.com> wrote: > > On Fri, Aug 24, 2012 at 2:16 PM, Chris Eldred <chris.eldred(a)gmail.com>wrote: > >> The following code snippet does not appear to be working: >> > > 1) I hate Fortran. A lot. It is not debuggable. > > 2) Here is my code that should work, but the do loop somehow does not work: > > program main > implicit none > ! > #include <finclude/petsc.h90> > #include <finclude/petscdmcomplex.h90> > ! > DM dm > PetscInt, target, dimension(4) :: EC > PetscInt, pointer :: pEC(:) > PetscInt c, firstCell, numCells > PetscErrorCode ierr > > call PetscInitialize(PETSC_NULL_CHARACTER,ierr) > call DMComplexCreate(PETSC_COMM_WORLD, dm, ierr) > firstCell = 0 > numCells = 5 > call DMComplexSetChart(dm, 0, numCells, ierr) > do c=firstCell,numCells > write(*,*) 'c',c > call DMComplexSetConeSize(dm, c, 4) > end do > call DMSetUp(dm, ierr) > > EC(1) = 1 > EC(2) = 2 > EC(3) = 3 > EC(4) = 4 > pEC => EC > write(*,*) 'cell',c,pEC > c = 0 > call DMComplexSetCone(dm, c , pEC, ierr) > CHKERRQ(ierr) > call DMComplexGetCone(dm, c , pEC, ierr) > CHKERRQ(ierr) > write(*,*) 'cell',c,pEC > > call DMDestroy(dm,ierr) > call PetscFinalize(ierr) > end > > > Matt > > >> PetscInt, dimension(4) :: EC >> >> EC(1) = 1 >> EC(2) = 2 >> EC(3) = 3 >> EC(4) = 4 >> write(*,*) 'cell',EC >> call DMComplexSetCone(model_mesh, lcell , EC, ierr) >> CHKERRQ(ierr) >> call DMComplexGetCone(model_mesh, lcell , EC, ierr) >> CHKERRQ(ierr) >> write(*,*) 'cell', EC >> >> I get the following as output: >> >> cell 1 2 3 4 >> cell 18673792 0 -1 -1 >> >> DMComplexGetConeSize returns 4 as expected. I am using the latest >> version of petsc-dev. >> >> Later on there is a segmentation fault violation, probably related to >> this issue. >> >> Any ideas? >> >> -- >> Chris Eldred >> DOE Computational Science Graduate Fellow >> Graduate Student, Atmospheric Science, Colorado State University >> B.S. Applied Computational Physics, Carnegie Mellon University, 2009 >> chris.eldred(a)gmail.com >> > > > > -- > 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 > > > -- > Department of Mathematics and Center for Computation & Technology > Louisiana State University, Baton Rouge, LA 70803, USA > Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 > http://www.math.lsu.edu/~bourdin > > > > > > > > > -- 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] Issues with DMComplexSetCone
by Blaise Bourdin 24 Aug '12

24 Aug '12

I am attaching a debugged version. Blaise On Aug 24, 2012, at 3:30 PM, Matthew Knepley <knepley(a)gmail.com> wrote: > On Fri, Aug 24, 2012 at 2:16 PM, Chris Eldred <chris.eldred(a)gmail.com> wrote: > The following code snippet does not appear to be working: > > 1) I hate Fortran. A lot. It is not debuggable. > > 2) Here is my code that should work, but the do loop somehow does not work: > > program main > implicit none > ! > #include <finclude/petsc.h90> > #include <finclude/petscdmcomplex.h90> > ! > DM dm > PetscInt, target, dimension(4) :: EC > PetscInt, pointer :: pEC(:) > PetscInt c, firstCell, numCells > PetscErrorCode ierr > > call PetscInitialize(PETSC_NULL_CHARACTER,ierr) > call DMComplexCreate(PETSC_COMM_WORLD, dm, ierr) > firstCell = 0 > numCells = 5 > call DMComplexSetChart(dm, 0, numCells, ierr) > do c=firstCell,numCells > write(*,*) 'c',c > call DMComplexSetConeSize(dm, c, 4) > end do > call DMSetUp(dm, ierr) > > EC(1) = 1 > EC(2) = 2 > EC(3) = 3 > EC(4) = 4 > pEC => EC > write(*,*) 'cell',c,pEC > c = 0 > call DMComplexSetCone(dm, c , pEC, ierr) > CHKERRQ(ierr) > call DMComplexGetCone(dm, c , pEC, ierr) > CHKERRQ(ierr) > write(*,*) 'cell',c,pEC > > call DMDestroy(dm,ierr) > call PetscFinalize(ierr) > end > > > Matt > > PetscInt, dimension(4) :: EC > > EC(1) = 1 > EC(2) = 2 > EC(3) = 3 > EC(4) = 4 > write(*,*) 'cell',EC > call DMComplexSetCone(model_mesh, lcell , EC, ierr) > CHKERRQ(ierr) > call DMComplexGetCone(model_mesh, lcell , EC, ierr) > CHKERRQ(ierr) > write(*,*) 'cell', EC > > I get the following as output: > > cell 1 2 3 4 > cell 18673792 0 -1 -1 > > DMComplexGetConeSize returns 4 as expected. I am using the latest > version of petsc-dev. > > Later on there is a segmentation fault violation, probably related to > this issue. > > Any ideas? > > -- > Chris Eldred > DOE Computational Science Graduate Fellow > Graduate Student, Atmospheric Science, Colorado State University > B.S. Applied Computational Physics, Carnegie Mellon University, 2009 > chris.eldred(a)gmail.com > > > > -- > 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 -- Department of Mathematics and Center for Computation & Technology Louisiana State University, Baton Rouge, LA 70803, USA Tel. +1 (225) 578 1612, Fax +1 (225) 578 4276 http://www.math.lsu.edu/~bourdin

1 0

Re: [petsc-dev] Issues with DMComplexSetCone
by Matthew Knepley 24 Aug '12

24 Aug '12

On Fri, Aug 24, 2012 at 2:16 PM, Chris Eldred <chris.eldred(a)gmail.com>wrote: > The following code snippet does not appear to be working: > 1) I hate Fortran. A lot. It is not debuggable. 2) Here is my code that should work, but the do loop somehow does not work: program main implicit none ! #include <finclude/petsc.h90> #include <finclude/petscdmcomplex.h90> ! DM dm PetscInt, target, dimension(4) :: EC PetscInt, pointer :: pEC(:) PetscInt c, firstCell, numCells PetscErrorCode ierr call PetscInitialize(PETSC_NULL_CHARACTER,ierr) call DMComplexCreate(PETSC_COMM_WORLD, dm, ierr) firstCell = 0 numCells = 5 call DMComplexSetChart(dm, 0, numCells, ierr) do c=firstCell,numCells write(*,*) 'c',c call DMComplexSetConeSize(dm, c, 4) end do call DMSetUp(dm, ierr) EC(1) = 1 EC(2) = 2 EC(3) = 3 EC(4) = 4 pEC => EC write(*,*) 'cell',c,pEC c = 0 call DMComplexSetCone(dm, c , pEC, ierr) CHKERRQ(ierr) call DMComplexGetCone(dm, c , pEC, ierr) CHKERRQ(ierr) write(*,*) 'cell',c,pEC call DMDestroy(dm,ierr) call PetscFinalize(ierr) end Matt > PetscInt, dimension(4) :: EC > > EC(1) = 1 > EC(2) = 2 > EC(3) = 3 > EC(4) = 4 > write(*,*) 'cell',EC > call DMComplexSetCone(model_mesh, lcell , EC, ierr) > CHKERRQ(ierr) > call DMComplexGetCone(model_mesh, lcell , EC, ierr) > CHKERRQ(ierr) > write(*,*) 'cell', EC > > I get the following as output: > > cell 1 2 3 4 > cell 18673792 0 -1 -1 > > DMComplexGetConeSize returns 4 as expected. I am using the latest > version of petsc-dev. > > Later on there is a segmentation fault violation, probably related to > this issue. > > Any ideas? > > -- > Chris Eldred > DOE Computational Science Graduate Fellow > Graduate Student, Atmospheric Science, Colorado State University > B.S. Applied Computational Physics, Carnegie Mellon University, 2009 > chris.eldred(a)gmail.com > -- 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

Issues with DMComplexSetCone
by Chris Eldred 24 Aug '12

24 Aug '12

The following code snippet does not appear to be working: PetscInt, dimension(4) :: EC EC(1) = 1 EC(2) = 2 EC(3) = 3 EC(4) = 4 write(*,*) 'cell',EC call DMComplexSetCone(model_mesh, lcell , EC, ierr) CHKERRQ(ierr) call DMComplexGetCone(model_mesh, lcell , EC, ierr) CHKERRQ(ierr) write(*,*) 'cell', EC I get the following as output: cell 1 2 3 4 cell 18673792 0 -1 -1 DMComplexGetConeSize returns 4 as expected. I am using the latest version of petsc-dev. Later on there is a segmentation fault violation, probably related to this issue. Any ideas? -- Chris Eldred DOE Computational Science Graduate Fellow Graduate Student, Atmospheric Science, Colorado State University B.S. Applied Computational Physics, Carnegie Mellon University, 2009 chris.eldred(a)gmail.com

1 0

Re: [petsc-dev] [petsc-maint #129715] __float128 complex
by Alexander Grayver 24 Aug '12

24 Aug '12

On 24.08.2012 17:53, Jose E. Roman wrote: > I guess PETSC_SCALAR should be correctly defined for the case complex+__float128 in petscsys.h. Maybe a new PetscDataType PETSC___COMPLEX128? This would require changing PetscMPIDataTypeToPetscDataType and then PetscBinaryRead... > > It is better if someone else can be in charge of this change. It is not necessary, I can write a kind of convertor myself. Otherwise complex128 seems to be working. -- Regards, Alexander

1 0

Re: [petsc-dev] [petsc-maint #129715] __float128 complex
by Jose E. Roman 24 Aug '12

24 Aug '12

El 24/08/2012, a las 13:57, Alexander Grayver escribió: > Jose, > > Thanks! > > One issue. I saw there is special treatment for quad precision in PetscBinaryRead: > > #if defined(PETSC_USE_REAL___FLOAT128) > double *ppp; > #endif > ... > #if defined(PETSC_USE_REAL___FLOAT128) > /* If using __float128 precision we still read in doubles from file */ > if (type == PETSC_SCALAR) { > m = m/2; > ierr = PetscMalloc(n*sizeof(double),&ppp);CHKERRQ(ierr); > pp = (char*)ppp; > } > #endif > > As far as I can see it doesn't take into account complex numbers currently? I guess PETSC_SCALAR should be correctly defined for the case complex+__float128 in petscsys.h. Maybe a new PetscDataType PETSC___COMPLEX128? This would require changing PetscMPIDataTypeToPetscDataType and then PetscBinaryRead... It is better if someone else can be in charge of this change. Jose

1 0

Re: [petsc-dev] [petsc-maint #129715] __float128 complex
by Alexander Grayver 24 Aug '12

24 Aug '12

Jose, Thanks! One issue. I saw there is special treatment for quad precision in PetscBinaryRead: #if defined(PETSC_USE_REAL___FLOAT128) double *ppp; #endif ... #if defined(PETSC_USE_REAL___FLOAT128) /* If using __float128 precision we still read in doubles from file */ if (type == PETSC_SCALAR) { m = m/2; ierr = PetscMalloc(n*sizeof(double),&ppp);CHKERRQ(ierr); pp = (char*)ppp; } #endif As far as I can see it doesn't take into account complex numbers currently? On 23.08.2012 18:02, Jose E. Roman wrote: > I have pushed a change to petsc-dev that allows to configure --with-precision=__float128 --with-scalar-type=complex > > I have tested this with src/ksp/ksp/examples/tutorials/ex12.c with 1 and 2 processes. It seems to work. > > Please review the change, since I had to modify pinit.c > > A couple of comments: > > 1) The f2cblaslapack interface to quad precision subroutines (q* and w*) is defined only in the case of "uscore" mangling. Maybe a test in configure should assure this. > > 2) I could not make it work --with-clanguage=c++ - I get compilation errors like this. Any idea? -- Regards, Alexander

1 0

Re: [petsc-dev] Adjacency relations in Sieve and stencils
by Jed Brown 23 Aug '12

23 Aug '12

On Thu, Aug 23, 2012 at 6:16 PM, Matthew Knepley <knepley(a)gmail.com> wrote: > On Thu, Aug 23, 2012 at 5:59 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > >> On Thu, Aug 23, 2012 at 5:50 PM, Chris Eldred <chris.eldred(a)gmail.com>wrote: >> >>> For the extended operators, I was going to write: >>> >>> 1) DMComplexGetLabelClosure- returns the points from a transitive >>> closure with a given value for a given label (this will implement both >>> star(vertex) \ support(vertex) and closure(cell) \ cone(cell); along >>> with many other operators) >>> 2) DMComplexGetNeighbors- returns U cone(support()) OR U >>> support(cone()) by setting a Boolean flag >>> >>> These two operators should cover "nearest neighbors" of all types. >>> >> >> What about having >> >> typedef enum {DMCOMPLEX_VERTEX = 0x1,DMCOMPLEX_EDGE = 0x2,DMCOMPLEX_FACE >> = 0x4,DMCOMPLEX_CELL = 0x8} DMComplexPointType; >> >> DMComplexGetClosure(dm,point_type,label,array,basepoint,&npoints,&values) >> >> You can always use point_type = my_type-1 to get values from all points >> in the closure, but this allows you to jump directly to particular levels. >> Similar for Support. Then you could get rid of cone and star because nobody >> can remember those anyway. You just have Closure and Support with results >> restricted to points of a given type (topological dimension). >> > > Never. > Note that you can do bitwise or and that mytype-1 gives the current Closure behavior while mytype>>1 gives Cone. (Wrap the shifts in convenience functions if you like.)

1 0

Re: [petsc-dev] Adjacency relations in Sieve and stencils
by Matthew Knepley 23 Aug '12

23 Aug '12

On Thu, Aug 23, 2012 at 5:59 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Thu, Aug 23, 2012 at 5:50 PM, Chris Eldred <chris.eldred(a)gmail.com>wrote: > >> For the extended operators, I was going to write: >> >> 1) DMComplexGetLabelClosure- returns the points from a transitive >> closure with a given value for a given label (this will implement both >> star(vertex) \ support(vertex) and closure(cell) \ cone(cell); along >> with many other operators) >> 2) DMComplexGetNeighbors- returns U cone(support()) OR U >> support(cone()) by setting a Boolean flag >> >> These two operators should cover "nearest neighbors" of all types. >> > > What about having > > typedef enum {DMCOMPLEX_VERTEX = 0x1,DMCOMPLEX_EDGE = 0x2,DMCOMPLEX_FACE = > 0x4,DMCOMPLEX_CELL = 0x8} DMComplexPointType; > > DMComplexGetClosure(dm,point_type,label,array,basepoint,&npoints,&values) > > You can always use point_type = my_type-1 to get values from all points in > the closure, but this allows you to jump directly to particular levels. > Similar for Support. Then you could get rid of cone and star because nobody > can remember those anyway. You just have Closure and Support with results > restricted to points of a given type (topological dimension). > Never. Matt > >> On Thu, Aug 23, 2012 at 12:39 PM, Chris Eldred <chris.eldred(a)gmail.com> >> wrote: >> > Sounds good- I will write the extended topological operators ( >> > star(vertex) \ support(vertex) ; closure(cell) \ cone(cell) and U >> > cone(support(edge)) ) and let you know if I have questions/run into >> > issues! >> > >> > -Chris >> > >> > On Thu, Aug 23, 2012 at 12:22 PM, Matthew Knepley <knepley(a)gmail.com> >> wrote: >> >> On Thu, Aug 23, 2012 at 12:19 PM, Chris Eldred <chris.eldred(a)gmail.com >> > >> >> wrote: >> >>> >> >>> I am working with a 2D unstructured mesh using sieve and I wanted to >> >>> get some advice on how to determine adjacency relations using the >> >>> topological operators provided by Sieve (cone, support, closure, >> >>> star). I have figured out the following: >> >>> >> >>> Vertices for a given edge: cone(edge) >> >>> Cells for a given edge: support(edge) >> >>> >> >>> Cells for a given vertex: star(vertex) \ support(vertex) >> >>> Edges for a given vertex: support(vertex) >> >>> >> >>> Edges for a given cell: cone(cell) >> >>> Vertices for a given cell: closure(cell) \ cone(cell) >> >>> >> >>> The set of edges associated with the cells for a given edge: U >> >>> cone(support(edge)) - ie the union of the cones of the support of the >> >>> edge >> >>> >> >>> Are there other topological operators that more naturally express >> >>> these relations (especially the ones involving unions or complements)? >> >> >> >> >> >> The only other builtin operations are meet and join. You are correct >> that >> >> you >> >> want to use union there. We could add routines like >> >> >> >> DMComplexGetConeUnion(dm, [p], &unionSize, &unionArray); >> >> >> >> but it was not clear that would do anything other than bloat the >> interface. >> >> I >> >> think the best way to procede is to write these convenience routines >> using >> >> the lower level primitives, and move them into PETSc if they turn out >> to be >> >> useful. I can help you write them. >> >> >> >> Matt >> >> >> >>> >> >>> >> >>> -- >> >>> Chris Eldred >> >>> DOE Computational Science Graduate Fellow >> >>> Graduate Student, Atmospheric Science, Colorado State University >> >>> B.S. Applied Computational Physics, Carnegie Mellon University, 2009 >> >>> chris.eldred(a)gmail.com >> >> >> >> >> >> >> >> >> >> -- >> >> 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 >> > >> > >> > >> > -- >> > Chris Eldred >> > DOE Computational Science Graduate Fellow >> > Graduate Student, Atmospheric Science, Colorado State University >> > B.S. Applied Computational Physics, Carnegie Mellon University, 2009 >> > chris.eldred(a)gmail.com >> >> >> >> -- >> Chris Eldred >> DOE Computational Science Graduate Fellow >> Graduate Student, Atmospheric Science, Colorado State University >> B.S. Applied Computational Physics, Carnegie Mellon University, 2009 >> chris.eldred(a)gmail.com >> > > -- 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