[petsc-dev] What should MatReset() do?

[Barry Smith]

  Note that the other XXXReset() do have code for each type, for example

PetscErrorCode PCReset_BJacobi_Singleblock(PC pc)
{
  PC_BJacobi             *jac  = (PC_BJacobi*)pc->data;
  PC_BJacobi_Singleblock *bjac = (PC_BJacobi_Singleblock*)jac->data;
  PetscErrorCode         ierr;

  PetscFunctionBegin;
  ierr = KSPReset(jac->ksp[0]);CHKERRQ(ierr);
  ierr = VecDestroy(&bjac->x);CHKERRQ(ierr);
  ierr = VecDestroy(&bjac->y);CHKERRQ(ierr);
  PetscFunctionReturn(0);
}

  in fact the goal of XXReset() is to free up all the "data" but keep the various options that had been set. Your "simple minded" reset would for example lose any matrix options that had been set by the user on the inner matrices. So, no you should not just be doing a "simple minded" reset. For example what would your simple minded reset do for a nested matrix? Throw away the four levels of nesting the user carefully crafted? We don't want that.


  Barry



> On Mar 6, 2015, at 8:53 PM, Dmitry Karpeyev <karpeev at mcs.anl.gov> wrote:
> 
> 
> 
> On Fri, Mar 6, 2015 at 8:33 PM Barry Smith <bsmith at mcs.anl.gov> wrote:
> 
> > On Mar 6, 2015, at 8:22 PM, Dmitry Karpeyev <karpeev at mcs.anl.gov> wrote:
> >
> >
> >
> > On Fri, Mar 6, 2015 at 6:03 PM Barry Smith <bsmith at mcs.anl.gov> wrote:
> >
> > > On Mar 6, 2015, at 4:41 PM, Jed Brown <jed at jedbrown.org> wrote:
> > >
> > > Dmitry Karpeyev <karpeev at mcs.anl.gov> writes:
> > >> This is trickier than it might appear: nonzerostate effectively counts the
> > >> global number of nonzeros.
> >
> >    No it does not. Note in MatZeroRows_SeqAIJ() when entries are deleted from the matrix we still increase the nonzerostate.
> > From MatAssemblyEnd_MPIAIJ():
> >
> >   if ((!mat->was_assembled && mode == MAT_FINAL_ASSEMBLY) || !((Mat_SeqAIJ*)(aij->A->data))->nonew) {
> >     PetscObjectState state = aij->A->nonzerostate + aij->B->nonzerostate;
> 
>   This is ok. So long as the two sub matrix states are always increasing the state of the total matrix will increase.
> 
> >     ierr = MPI_Allreduce(&state,&mat->nonzerostate,1,MPIU_INT64,MPI_SUM,PetscObjectComm((PetscObject)mat));CHKERRQ(ierr);
> >   }
> >
> > Barring MatZeroRows MATSEQAIJ matrices aij->A and aij->B will simply increment their nonzerostates on each new nonzero insertion in MatSetVAlues_SeqAIJ(), so the containing MATMPIAIJ M ends up with the total number of nonzeros in its nonzerostate.
> >
> > Now I reset my M.  That will blow away aij->A and aij->B.
> 
>   They should not be "blown away". They should also be reset and in being reset their nonzerostate will never get smaller.
> Okay, this requires implementing MatReset_XXX() that will carefully clean things out for each matrix type, whereas I was using a simple-minded approach of just destroying and recreating the object.  Partly because this is more or less the equivalent of the only way to "reset" the matrix now: 
> MatSetType(M,othertype);MatSetType(M,oldtype) (*).
> 
> 
> >  I can now insert the same number of nonzeros, but in a different pattern.  M will end up with the same nonzerostate as before the reset and confuse the PC, no?
> 
>    Why are you reseting the nonzerostate to zero, just don't do that.
> I'm not resetting it to anything: nonzerostate will get recomputed from those of A&B using Allreduce.  They (A&B), however, will get rebuilt from scratch, if all of M's  guts are destroyed and recreated.  Even if we insist on writing complicated MatReset_XXX() variants that avoid that,   we still need to make sure that (*) above (which will destroy and recreate the inner A&B) doesn't confuse the preconditioner.  Or do we simply tell the user not to do that :-)?
> 
> >
> > This started out as a discussion about MatReset(), but I think this _may_ be a bug we are seeing in one of the elastic contact applications: PCASM tries to rebuild itself with MAT_REUSE_MATRIX when subdomain matrices actually have different numbers of nonzeros.  I have to say I haven't ascertained that an inconsistent nonzerostate cases the problem, yet -- reproducible test cases that trigger the problem are still too big to debug.
> 
>   It is possible that somewhere the state is not properly handled by being incremented.
> 
>   Barry
> 
> >
> > >> The PC will rebuild if its state is stale, but
> > >> it will reuse matrices (e.g., subdomain matrices in PCASM) if nonzerostate
> > >> is up to date.  This works if the sparsity pattern never drops nonzeros,
> > >> but that's no longer true if reset is allowed. I can reset a matrix,
> > >> preallocate and assemble it so that the global number of nonzeros will be
> > >> the same as before the resetting, but local sparsity patterns will change.
> > >> This could happen, for example, when I have moving particles or, less
> > >> exotically, when I have elastic contact and nodes move past each other.
> > >> That will break PCASM.
> >
> >   Just increase the nonzerostate flag by one on a reset (that is there is no reason to ever set it back to zero). Now nonzerostate is monotonically increasing.
> >
> >   Barry
> >
> >
> >   Barry
> >
> > >
> > > On pretty simple and reliable solution would be to take a cryptographic
> > > hash of the row/col arrays.  I assume BG is really atrocious at hashing,
> > > but is it so bad that this is not viable?  (There are several places
> > > where we use kinda fragile state counters or trust the user, but hashes
> > > would make rebuilding more reliable and transparent.)
> >