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

[Dmitry Karpeyev]

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/f53eb5ca/attachment.html>