On Fri, Mar 6, 2015 at 8:33 PM Barry Smith <[email protected]> wrote:
> > > On Mar 6, 2015, at 8:22 PM, Dmitry Karpeyev <[email protected]> wrote: > > > > > > > > On Fri, Mar 6, 2015 at 6:03 PM Barry Smith <[email protected]> wrote: > > > > > On Mar 6, 2015, at 4:41 PM, Jed Brown <[email protected]> wrote: > > > > > > Dmitry Karpeyev <[email protected]> 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.) > > > >
