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

Re: [petsc-dev] who broke petsc-dev
by Matthew Knepley 13 Nov '11

13 Nov '11

On Sun, Nov 13, 2011 at 10:44 PM, Jed Brown <jedbrown(a)mcs.anl.gov> wrote: > On Sun, Nov 13, 2011 at 16:36, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > >> 1) It is legitimate to call KSPSetOperators with a null Amat so this is a >> bug that crashes large parts of PETSc >> > > This logic should be moved much later (like to KSPSetUp()). > I disagree, since SetOperators() is where we know that the matrix is changing. Matt > > >> >> 2) Do you really want to do this? I am not sure I like it. >> > > I think that attaching the null space to the Mat should be the preferred > way to communicate it, especially in the case of a deeply nested hierarchy. > -- 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] developer hint
by Barry Smith 13 Nov '11

13 Nov '11

On Nov 13, 2011, at 4:42 PM, Jed Brown wrote: > On Sun, Nov 13, 2011 at 16:29, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > I have PETSC_OPTIONS set in my .bashrc file to -malloc -malloc_debug -malloc_dump by default. > > This means when I introduce memory leaks I find them immediately. Like someone just introduced in src/dm/examples/test/ex30.c but did not detect because they did not have these options set. > > This is my fault. In DMDestroy(), there is somewhat strange logic to count the number of circular references in the "gotten" vectors. > > DMSliced clones new Vecs off of it's own private Vec, which has yet another circular reference. Should we just artificially decrement the reference count on the DM, and then re-reference it in DMDestroy_Sliced before destroying slice->globalvector? I wonder if there is a more elegant way to handle these circular references. Why do it anyway different than DMDA? Better consistent and not a great way then 10 different inconsistent and slightly different and better and worse ways.

1 0

Re: [petsc-dev] who broke petsc-dev
by Matthew Knepley 13 Nov '11

13 Nov '11

On Sun, Nov 13, 2011 at 10:36 PM, Barry Smith <bsmith(a)mcs.anl.gov> wrote: 0) Me PetscErrorCode KSPSetOperators(KSP ksp,Mat Amat,Mat Pmat,MatStructure flag) > { > MatNullSpace nullsp; > PetscErrorCode ierr; > > PetscFunctionBegin; > PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); > if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); > if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); > if (Amat) PetscCheckSameComm(ksp,1,Amat,2); > if (Pmat) PetscCheckSameComm(ksp,1,Pmat,3); > if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);} > ierr = PCSetOperators(ksp->pc,Amat,Pmat,flag);CHKERRQ(ierr); > if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = > KSP_SETUP_NEWMATRIX; /* so that next solve call will call PCSetUp() on new > matrix */ > if (ksp->guess) { > ierr = KSPFischerGuessReset(ksp->guess);CHKERRQ(ierr); > } > ierr = MatGetNullSpace(Amat, &nullsp);CHKERRQ(ierr); > if (nullsp) { > ierr = KSPSetNullSpace(ksp, nullsp);CHKERRQ(ierr); > } > PetscFunctionReturn(0); > } > > > 1) It is legitimate to call KSPSetOperators with a null Amat so this is a > bug that crashes large parts of PETSc > What does a NULL A even mean here. I can understand a NULL P, which might mean A == P, but what would it mean to have no A? > 2) Do you really want to do this? I am not sure I like it. > Yes, I definitely do. How else will the solver get information about the null space of an operator. 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

Re: [petsc-dev] who broke petsc-dev
by Jed Brown 13 Nov '11

13 Nov '11

On Sun, Nov 13, 2011 at 16:36, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > 1) It is legitimate to call KSPSetOperators with a null Amat so this is a > bug that crashes large parts of PETSc > This logic should be moved much later (like to KSPSetUp()). > > 2) Do you really want to do this? I am not sure I like it. > I think that attaching the null space to the Mat should be the preferred way to communicate it, especially in the case of a deeply nested hierarchy.

1 0

Re: [petsc-dev] developer hint
by Jed Brown 13 Nov '11

13 Nov '11

On Sun, Nov 13, 2011 at 16:29, Barry Smith <bsmith(a)mcs.anl.gov> wrote: > I have PETSC_OPTIONS set in my .bashrc file to -malloc -malloc_debug > -malloc_dump by default. > > This means when I introduce memory leaks I find them immediately. Like > someone just introduced in src/dm/examples/test/ex30.c but did not detect > because they did not have these options set. > This is my fault. In DMDestroy(), there is somewhat strange logic to count the number of circular references in the "gotten" vectors. DMSliced clones new Vecs off of it's own private Vec, which has yet another circular reference. Should we just artificially decrement the reference count on the DM, and then re-reference it in DMDestroy_Sliced before destroying slice->globalvector? I wonder if there is a more elegant way to handle these circular references.

1 0

who broke petsc-dev
by Barry Smith 13 Nov '11

13 Nov '11

PetscErrorCode KSPSetOperators(KSP ksp,Mat Amat,Mat Pmat,MatStructure flag) { MatNullSpace nullsp; PetscErrorCode ierr; PetscFunctionBegin; PetscValidHeaderSpecific(ksp,KSP_CLASSID,1); if (Amat) PetscValidHeaderSpecific(Amat,MAT_CLASSID,2); if (Pmat) PetscValidHeaderSpecific(Pmat,MAT_CLASSID,3); if (Amat) PetscCheckSameComm(ksp,1,Amat,2); if (Pmat) PetscCheckSameComm(ksp,1,Pmat,3); if (!ksp->pc) {ierr = KSPGetPC(ksp,&ksp->pc);CHKERRQ(ierr);} ierr = PCSetOperators(ksp->pc,Amat,Pmat,flag);CHKERRQ(ierr); if (ksp->setupstage == KSP_SETUP_NEWRHS) ksp->setupstage = KSP_SETUP_NEWMATRIX; /* so that next solve call will call PCSetUp() on new matrix */ if (ksp->guess) { ierr = KSPFischerGuessReset(ksp->guess);CHKERRQ(ierr); } ierr = MatGetNullSpace(Amat, &nullsp);CHKERRQ(ierr); if (nullsp) { ierr = KSPSetNullSpace(ksp, nullsp);CHKERRQ(ierr); } PetscFunctionReturn(0); } 1) It is legitimate to call KSPSetOperators with a null Amat so this is a bug that crashes large parts of PETSc 2) Do you really want to do this? I am not sure I like it.

1 0

developer hint
by Barry Smith 13 Nov '11

13 Nov '11

I have PETSC_OPTIONS set in my .bashrc file to -malloc -malloc_debug -malloc_dump by default. This means when I introduce memory leaks I find them immediately. Like someone just introduced in src/dm/examples/test/ex30.c but did not detect because they did not have these options set. Please fix your own messes before other people slip on them. [ 0]80 bytes DMCreate_Sliced() line 277 in /Users/barrysmith/Src/petsc-dev/src/dm/impls/sliced/sliced.c [0] DMSetType() line 1377 in /Users/barrysmith/Src/petsc-dev/src/dm/interface/dm.c [0] DMSlicedCreate() line 310 in /Users/barrysmith/Src/petsc-dev/src/dm/impls/sliced/sliced.c [ 0]16 bytes PetscCommDuplicate() line 152 in /Users/barrysmith/Src/petsc-dev/src/sys/objects/tagm.c [0] PetscHeaderCreate_Private() line 30 in /Users/barrysmith/Src/petsc-dev/src/sys/objects/inherit.c [0] DMCreate() line 32 in /Users/barrysmith/Src/petsc-dev/src/dm/interface/dm.c [0] DMSlicedCreate() line 310 in /Users/barrysmith/Src/petsc-dev/src/dm/impls/sliced/sliced.c [ 0]224 bytes DMCreate() line 39 in /Users/barrysmith/Src/petsc-dev/src/dm/interface/dm.c [0] DMSlicedCreate() line 310 in /Users/barrysmith/Src/petsc-dev/src/dm/impls/sliced/sliced.c [ 0]64 bytes DMCreate() line 39 in /Users/barrysmith/Src/petsc-dev/src/dm/interface/dm.c [0] DMSlicedCreate() line 310 in /Users/barrysmith/Src/petsc-dev/src/dm/impls/sliced/sliced.c [ 0]3776 bytes DMCreate() line 39 in /Users/barrysmith/Src/petsc-dev/src/dm/interface/dm.c [0] DMSlicedCreate() line 310 in /Users/barrysmith/Src/petsc-dev/src/dm/impls/sliced/sliced.c Barry

1 0

Re: [petsc-dev] [petsc-users] Any suggestion for this kinds of matrix?
by Jed Brown 11 Nov '11

11 Nov '11

On Fri, Nov 11, 2011 at 11:17, Mark F. Adams <mark.adams(a)columbia.edu>wrote: > All I recall was that I was confused because the residual dropped 10 > orders of magnitude, apparently, in the first iteration and Barry explained > the issue, I added 'unpreconditioned' residual and it was fine. My recollection with that problem was that the initial guess "almost" satisfied the PDE in the interior, but did not satisfy the boundary conditions. The preconditioned residual made the boundaries and interior of roughly equal importance, so the initial preconditioned residual looked large due to the big disagreement near boundaries. In unpreconditioned residuals, you never see the boundary equations in the norm, but they happen to be satisfied because the preconditioner tends to fix them. Or something like this, but with ksp_diagonal_scale, the residual drop was similar for preconditioned and unpreconditioned residuals. I want methods to be scale invariant, but I don't think it's practical to ask algorithms to behave the same way if you scale some equations very differently from others. That changes the conditioning of the operator and makes dot products with the operator meaningless. I think it's achievable to ask people to scale their equations properly if they want good performance. The questions is which causes less confusion when people mess up. I'm not sold either way regarding preconditioned or unpreconditioned residuals.

1 0

Re: [petsc-dev] [petsc-users] Any suggestion for this kinds of matrix?
by Mark F. Adams 11 Nov '11

11 Nov '11

On Nov 11, 2011, at 11:20 AM, Jed Brown wrote: > On Fri, Nov 11, 2011 at 10:08, Mark F. Adams <mark.adams(a)columbia.edu> wrote: > First, "my" matrices (actually XGC1 matrices from one of my projects and I don't even know exactly what they do) put 1.0 on the diagonal for BCs and god knows how the thing is scaled (1e13 apparently). But this does not matter because the RHS and initial guess are 0.0 so the BCs have been completely removed from the algebra, even if they are still in the data structures. > > So what happens when you use the same method with inhomogeneous Dirichlet conditions? Jed, I agreed with this point in my email. And I can now see your point that using the preconditioned residual is a good default for PETSc, or at least I'm willing to defer to your judgment. > This is especially bad if the interior is scaled by, say, 1e-13 instead of 1e+13, because the boundary conditions dominate the initial residual. > > > As I recall the preconditioned residual was confusing because since these are Laplacian matrices with a scale of 1e13, as you found, the residual dropped like 10 orders of magnitude in the first iteration, which was pretty confusing. > > The preconditioned residual fixes the scaling. If you evaluate the initial unpreconditioned residual, then you see the confusing scaling. But preconditioning, even with just Jacobi, fixes the scaling. All I recall was that I was confused because the residual dropped 10 orders of magnitude, apparently, in the first iteration and Barry explained the issue, I added 'unpreconditioned' residual and it was fine. So if I am recalling this correctly this default confused me but apparently, or what I infer from this discussion, more users get confused with using the unpreconditioned residual so you've chosen preconditioned residual as the default. Mark > > > So I don't see these XGC1 problems as being arguments for preconditioned residual, in fact they argue against it, right? > > The preconditioned residual fixes the scaling, assuming it is no worse than Jacobi.

1 0

Re: [petsc-dev] [petsc-users] Any suggestion for this kinds of matrix?
by Jed Brown 11 Nov '11

11 Nov '11

On Fri, Nov 11, 2011 at 10:08, Mark F. Adams <mark.adams(a)columbia.edu>wrote: > First, "my" matrices (actually XGC1 matrices from one of my projects and I > don't even know exactly what they do) put 1.0 on the diagonal for BCs and > god knows how the thing is scaled (1e13 apparently). But this does not > matter because the RHS and initial guess are 0.0 so the BCs have been > completely removed from the algebra, even if they are still in the data > structures. > So what happens when you use the same method with inhomogeneous Dirichlet conditions? This is especially bad if the interior is scaled by, say, 1e-13 instead of 1e+13, because the boundary conditions dominate the initial residual. > > As I recall the preconditioned residual was confusing because since these > are Laplacian matrices with a scale of 1e13, as you found, the residual > dropped like 10 orders of magnitude in the first iteration, which was > pretty confusing. > The preconditioned residual fixes the scaling. If you evaluate the initial unpreconditioned residual, then you see the confusing scaling. But preconditioning, even with just Jacobi, fixes the scaling. > > So I don't see these XGC1 problems as being arguments for preconditioned > residual, in fact they argue against it, right? > The preconditioned residual fixes the scaling, assuming it is no worse than Jacobi.

1 0