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

[Dmitry Karpeyev]

On Fri, Mar 6, 2015 at 9:02 PM Barry Smith <bsmith at mcs.anl.gov> wrote:

>
>   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.
>
> Sure.  What about the MatSetType();MatSetType() route to a potentially
ill-defined state?  That sequence will effectively destroy and recreate the
inner data structures, instead of recursively resetting them.

Simpler still, MatSetType(mat,other) will destroy one implementation,
create another (e.g., BAIJ instead of AIJ), and we are free to create an
identical nonzerostate in a matrix of a different type!  The PC will only
know that the Mat's object state got incremented and will be very
confused.  Does that mean we should disallow the resetting of matrix types?

>
>   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.)
> > >
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20150307/64ae82ff/attachment.html>