On Sat, Dec 28, 2013 at 9:01 AM, Stephan Kramer <[email protected]>wrote:
> 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 > Okay, I have checked in preliminary code. It passes the tests for Mat ex18. I have merged it to next. Let me know how it works for you. 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
