Comment 9 for bug 1088175

On Wed, Dec 12, 2012 at 10:39:37AM +0100, Anders Logg wrote:
> 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?

Works now. The problem was this...

  m_local = local_range.first;
  n_local = local_range.second;

:-)

Now replaced by this:

  m_local = local_range_y.second - local_range_y.first;
  n_local = local_range_x.second - local_range_x.first;

A LinearOperator must now be initialized with vectors x and y matching
the product y = Ax. In most cases, it will work fine to use the same
vector (vector b or u.vector()) since they will have the same parallel
layout.

--
Anders