On 19/12/13 21:03, Matthew Knepley wrote:
On Sat, Nov 30, 2013 at 12:13 PM, Stephan Kramer <[email protected] <mailto:[email protected]>> 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
