[petsc-dev] bug in MatZeroRowsColumns for MPIAIJ

[Matthew Knepley]

On Sat, Nov 30, 2013 at 12:13 PM, Stephan Kramer <s.kramer at imperial.ac.uk>wrote:

> Dear petsc-devs
>
> Trying to write a test case for MatZeroRowColumns (as requested by Matt)
> in order to extend it to MPIBAIJ, I came across a bug in the current
> implementation for MPIAIJ. The attached test code assembles a simple 2D
> Laplacian (with the standard 4,-1,-1,-1 stencil) on the following grid:
>
> 0 1 2  3
> 4 5 6  7
> 8 9 10 11
>
> with boundary conditions applied to nodes 0,1,2,3,4 and 8 with values
> (stored in Vec x) 0., 1., 2., 3., 4. and 8. For simplicity I've kept a zero
> rhs (before the call to MatZeroRowsColumns)  After the call to
> MatZeroRowColumns I therefore expect the following rhs:
>
> 0. 1. 2. 3. 4. 5. 2. 3. 8. 8. 0. 0.
>
> The entries at 0,1,2,3,4 and 8 correspond to the boundary values
> (multiplied by a diagonal value of 1.0) and the values at 5,6,7 and 9 to
> the the lifted values: 1.+4.=5.,2.,3. and 8. However, if you run the
> attached example that's not what you get, the lifted values at 6 and 7 are
> wrong. This works fine in serial (you have to change nlocal=4 to get the
> same grid).
>
> I think this is caused by the way that entries in x are communicated with
> other processes that have off-diagonal entries referencing those. In line
> 944 of mpiaij.c it uses a VecScatter with ADD_VALUES to do this, but it
> uses the local work vector l->lvec without zeroing it first. This means it
> will still have values in it of whatever it was used for before. In this
> case that would be the MatMult of line 125 of test.c (in fact you can "fix"
> the example by adding the line 136-139 where it does a spurious MatMult
> with a zero vector).
>
> The attached code could be easily extended to test block matrices (simply
> increasing bs). The commented out lines in the assembly were basically to
> make the problem a bit more challenging (have an asymmetric stencil and the
> setting of off-process boundary nodes). The expected results for MPIBAIJ
> can be obtained by comparing with the MPIAIJ results (once that works
> correctly).
>

Okay, I have put this in knepley/feature-mat-zerorowscols-baij and pushed
to next. I also fixed the bug that you found. It should
not take much over the holidays to get the BAIJ functionality you want.

  Thanks,

     Matt


> Cheers
> Stephan
>



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