[petsc-dev] bug in MatZeroRowsColumns for MPIAIJ

[Stephan Kramer]

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).

Cheers
Stephan
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test.c
Type: text/x-csrc
Size: 6431 bytes
Desc: not available
URL: <http://lists.mcs.anl.gov/pipermail/petsc-dev/attachments/20131130/34f6da20/attachment.bin>