Comment 8 for bug 1088175

On Wed, Dec 12, 2012 at 08:58:26AM -0000, Garth Wells wrote:
> On Tue, Dec 11, 2012 at 11:21 PM, Anders Logg <email address hidden> wrote:
> > On Sun, Dec 09, 2012 at 04:02:03PM +0100, Anders Logg wrote:
> >> On Sun, Dec 09, 2012 at 02:32:19PM -0000, Garth Wells wrote:
> >> > On Sun, Dec 9, 2012 at 2:01 PM, Anders Logg <email address hidden> wrote:
> >> > > On Sun, Dec 09, 2012 at 01:27:50PM -0000, Garth Wells wrote:
> >> > >> On Sun, Dec 9, 2012 at 1:20 PM, Anders Logg <email address hidden> wrote:
> >> > >> > I can take a look but object to the description. It is indeed written
> >> > >> > to work in parallel. It has also been tested to work in parallel
> >> > >> > before.
> >> > >> >
> >> > >>
> >> > >> Tested where?
> >> > >
> >> > > I thought the unit test was running fine on the buildbots but I see
> >> > > now it was never added to the list of tests (in the main test.py).
> >> > >
> >> >
> >> > The unit test is a special case for which it may run on 3 processes
> >> > because the number of dofs is divisible by 3. In general, the parallel
> >> > layouts will not match.
> >>
> >> It doesn't even work for that case. So what happened is that I thought
> >> the test was run on the buildbot and therefore worked in parallel. I
> >> see now why it doesn't work but it should be relatively easy to fix.
> >
> > I've made an attempt to fix this but it still crashes.
>
> If you're using test/unit/la/python/LinearOperator.py to test, the line:
>
> b = Vector(V.dim())
>
> will certainly cause a problem. We have the function
> GenericMatrix::.resize(GenericVector&x, std::size_t dim) to re-size a
> vector consistently with a matrix operator,

Yes, that line was obviously wrong. I've replaced it with a proper
right-hand side now in the unit test but it still breaks.

It looks like I need to initialize the local rows and columns separately:

http://www.mcs.anl.gov/petsc/petsc-current/docs/manualpages/Mat/MatCreateShell.html

The local rows are determined by the result vector y and the local
columns by the vector x in y = Ax. But in a Krylov solve, the result
vector y will again be used to multiply the matrix, so it looks to me
that it needs to be the same?

--
Anders