[petsc-dev] bug in MatZeroRowsColumns for MPIAIJ

[Stephan Kramer]

On 19/12/13 21:03, Matthew Knepley wrote:
> On Sat, Nov 30, 2013 at 12:13 PM, Stephan Kramer <s.kramer at imperial.ac.uk <mailto: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

Great. Let me know if there's anything else I can do to help. Just a small heads up: the new src/mat/examples/tests/ex18.c has an uninitialized bool nonlocalBC
Cheers
Stephan