[petsc-dev] bug in MatZeroRowsColumns for MPIAIJ

[Matthew Knepley]

On Sat, Dec 28, 2013 at 9:01 AM, Stephan Kramer <s.kramer at imperial.ac.uk>wrote:

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

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