Merge lp:~jobh/dolfin/fast-array into lp:~fenics-core/dolfin/trunk

[Joachim Haga]

Changes Array from shared (using boost::shared_array) to a simpler representation. This change makes Poisson assembly 35% faster.

Creating a shared_array pointer requires heap allocations, and a lot of the time in assemble was spent just wrapping plain pointers in arrays. Furthermore, the shared semantics weren't used -- and if required the Array can be assigned to a shared_ptr.

Since this changes the binary interface to compiled modules, I added an interface version number (apart from the dolfin version number, which has more to do with the release schedule). Therefore instant modules will be recompiled.

[Garth Wells (garth-wells)]

Where in the code were the most Arrays being created?

[Joachim Haga (jobh)]

In GenericFunction::evaluate, via the chain

Assembler::assemble_cells
UFC::update
Expression::restrict
GenericFunction::restrict_as_ufc
GenericFunction::evaluate <===
ffc_form_xxx

-j.

On 24 February 2012 13:34, Garth Wells <email address hidden> wrote:

> Where in the code were the most Arrays being created?
> --
> https://code.launchpad.net/~jobh/dolfin/fast-array/+merge/94467
> You are the owner of lp:~jobh/dolfin/fast-array.
>

[Joachim Haga (jobh)]

On 24 February 2012 13:34, Garth Wells <email address hidden> wrote:

> Where in the code were the most Arrays being created?
>

I wrote:

> In GenericFunction::evaluate, via the chain
>
> Assembler::assemble_cells
> UFC::update
> Expression::restrict
>

... and indeed it's the rhs that involves an expression that's faster:
f = Expression("...")
L = f*v*dx
b = assemble(L)

goes from 8.0s on 64x64x64 to 2.7s. This accounts for most (all?) of the
35% speed-up for A and b combined.

-j.

[Garth Wells (garth-wells)]

I'm not keen on ownership flags - they were a real mess before and it's been a lot better since we go rid of them.

Would be it be possible to (1) re-use the Array so we avoid dynamic allocation or (2) have an Array-like structure that always owns the data? The class Array was introduced to better interface with Python, but I recall that we later became aware that the Numpy interface has has a flag to prevent resizing, which would make typemaps more robust.

[Joachim Haga (jobh)]

I'd opine that the easiest would be to make array never own the data,
rather than always own the data. In that way, it can wrap anything.

* TimeSeries returns (owned) Arrays in two places. These could return
references to its own data (const vector<double>&), that would save a copy
(they are currently copied twice, once into an owned Array, then again in
swig into a numpy array). But maybe the amount of data returned is small,
in which case a vector<double> will be fine.

* Can't support resize() without ownership, so everywhere resize is used on
a reference argument should pass a vector<double>& instead. I think this is
a few places only.

At that point, Array will be a plain wrapper of others' data.

As for your (1), it might well be possible (keep a scratch array in the
GenericFunction class perhaps). But the thing is, the shared version was no
more safe. You can't keep a reference to an Array because you don't know if
it's shared or just wrapping data that will be deleted from under you. It
was "shared or borrowed" rather than "owned or borrowed", only without the
flag to tell the cases apart.

-j.

On 24 February 2012 19:37, Garth Wells <email address hidden> wrote:

> I'm not keen on ownership flags - they were a real mess before and it's
> been a lot better since we go rid of them.
>
> Would be it be possible to (1) re-use the Array so we avoid dynamic
> allocation or (2) have an Array-like structure that always owns the data?
> The class Array was introduced to better interface with Python, but I
> recall that we later became aware that the Numpy interface has has a flag
> to prevent resizing, which would make typemaps more robust.
> --
> https://code.launchpad.net/~jobh/dolfin/fast-array/+merge/94467
> You are the owner of lp:~jobh/dolfin/fast-array.
>

[lp:~jobh/dolfin/fast-array]

6592. By Joachim Haga on 2012-02-26

Compilation fix (initialisation order, with -Wall -Werror)

[Anders Logg (logg)]

On Fri, Feb 24, 2012 at 07:01:57PM -0000, Joachim Haga wrote:
> I'd opine that the easiest would be to make array never own the data,
> rather than always own the data. In that way, it can wrap anything.
>
> * TimeSeries returns (owned) Arrays in two places. These could
> return references to its own data (const vector<double>&), that
> would save a copy (they are currently copied twice, once into an
> owned Array, then again in swig into a numpy array). But maybe the
> amount of data returned is small, in which case a vector<double>
> will be fine.

The return of arrays from TimeSeries is typically not performance
critical so a copy is fine (which it already is).

> * Can't support resize() without ownership, so everywhere resize is
> used on a reference argument should pass a vector<double>&
> instead. I think this is a few places only.
>
> At that point, Array will be a plain wrapper of others' data.
>
> As for your (1), it might well be possible (keep a scratch array in
> the GenericFunction class perhaps). But the thing is, the shared
> version was no more safe. You can't keep a reference to an Array
> because you don't know if it's shared or just wrapping data that
> will be deleted from under you. It was "shared or borrowed" rather
> than "owned or borrowed", only without the flag to tell the cases
> apart.

I'm not very happy with our current use of Array. We use std::vector
in some places and Array in other places. I don't remember what the
original intention was with Array. Is the purpose just to be able to
pass shared data as arrays back and forth to Python? Or should it be a
general array class used throughout DOLFIN C++?

If it is just for Python-wrapping, I think we should make it never own
data as Joachim suggests.

--
Anders

> -j.
>
> On 24 February 2012 19:37, Garth Wells <email address hidden> wrote:
>
> > I'm not keen on ownership flags - they were a real mess before and it's
> > been a lot better since we go rid of them.
> >
> > Would be it be possible to (1) re-use the Array so we avoid dynamic
> > allocation or (2) have an Array-like structure that always owns the data?
> > The class Array was introduced to better interface with Python, but I
> > recall that we later became aware that the Numpy interface has has a flag
> > to prevent resizing, which would make typemaps more robust.
> >
>

[Garth Wells (garth-wells)]

On 27 February 2012 13:40, Anders Logg <email address hidden> wrote:
> On Fri, Feb 24, 2012 at 07:01:57PM -0000, Joachim Haga wrote:
>> I'd opine that the easiest would be to make array never own the data,
>> rather than always own the data. In that way, it can wrap anything.
>>
>> * TimeSeries returns (owned) Arrays in two places. These could
>> return references to its own data (const vector<double>&), that
>> would save a copy (they are currently copied twice, once into an
>> owned Array, then again in swig into a numpy array). But maybe the
>> amount of data returned is small, in which case a vector<double>
>> will be fine.
>
> The return of arrays from TimeSeries is typically not performance
> critical so a copy is fine (which it already is).
>
>> * Can't support resize() without ownership, so everywhere resize is
>> used on a reference argument should pass a vector<double>&
>> instead. I think this is a few places only.
>>
>> At that point, Array will be a plain wrapper of others' data.
>>
>> As for your (1), it might well be possible (keep a scratch array in
>> the GenericFunction class perhaps). But the thing is, the shared
>> version was no more safe. You can't keep a reference to an Array
>> because you don't know if it's shared or just wrapping data that
>> will be deleted from under you. It was "shared or borrowed" rather
>> than "owned or borrowed", only without the flag to tell the cases
>> apart.
>
> I'm not very happy with our current use of Array. We use std::vector
> in some places and Array in other places. I don't remember what the
> original intention was with Array. Is the purpose just to be able to
> pass shared data as arrays back and forth to Python?

Yes.

> Or should it be a
> general array class used throughout DOLFIN C++?
>

No - there is no point in re-inventing std::vector. We should use
std::vector where possible. It's more flexible and can be used with
STL algorithms.

> If it is just for Python-wrapping, I think we should make it never own
> data as Joachim suggests.
>

Memory has to be created somewhere in order to be used, and it has to
be cleaned up.

There are two issues here:

1) Dynamic allocation in order to interface with UFC. Using plain
pointers rather than smart pointers is likely faster, but I think that
it's still not a solution. We should try to avoid these allocations in
loops.

2) We need to re-visit the NumPy interface. If we can create NumPy
arrays of fixed length, then we can just use std::vector, and let the
NumPy array wrap the pointer &x[0]. Maybe Johan Hake can comment on
this?

Garth

> --
> Anders
>
>
>> -j.
>>
>> On 24 February 2012 19:37, Garth Wells <email address hidden> wrote:
>>
>> > I'm not keen on ownership flags - they were a real mess before and it's
>> > been a lot better since we go rid of them.
>> >
>> > Would be it be possible to (1) re-use the Array so we avoid dynamic
>> > allocation or (2) have an Array-like structure that always owns the data?
>> > The class Array was introduced to better interface with Python, but I
>> > recall that we later became aware that the Numpy interface has has a flag
>> > to prevent resizing, which would make typemaps more robust.
>> >
>>
>
> --
> ht...

[Joachim Haga (jobh)]

>> I'm not very happy with our current use of Array. We use std::vector
>> in some places and Array in other places. I don't remember what the
>> original intention was with Array. Is the purpose just to be able to
>> pass shared data as arrays back and forth to Python?
>
> Yes.

The current usage is:
- As input (Python->C++) it wraps the numpy data.
- As output (C++->Python) the data is copied into a new numpy array.

I don't know off-hand if the input case is a requirement (that C++
modifies the data) or just an optimisation.

>> Or should it be a
>> general array class used throughout DOLFIN C++?
>
> No - there is no point in re-inventing std::vector. We should use
> std::vector where possible. It's more flexible and can be used with
> STL algorithms.

I agree in general, but note that there is no way to make a
std::vector wrap numpy data like in the input case above. An Array --
even if it's just a typedef to std::pair<T*,uint> -- may be nice to
standardise how wrapped data is passed.

> 1) Dynamic allocation in order to interface with UFC. Using plain
> pointers rather than smart pointers is likely faster, but I think that
> it's still not a solution. We should try to avoid these allocations in
> loops.

There is no dynamic allocation of actual data, not before and not now.
There was, however, dynamic allocation of the shared_array reference
count which is heap allocated. That is what this merge request fixes.

-j.

[Johan Hake (johan-hake)]

On 02/27/2012 03:01 PM, Garth Wells wrote:
> On 27 February 2012 13:40, Anders Logg<email address hidden> wrote:
>> On Fri, Feb 24, 2012 at 07:01:57PM -0000, Joachim Haga wrote:
>>> I'd opine that the easiest would be to make array never own the data,
>>> rather than always own the data. In that way, it can wrap anything.
>>>
>>> * TimeSeries returns (owned) Arrays in two places. These could
>>> return references to its own data (const vector<double>&), that
>>> would save a copy (they are currently copied twice, once into an
>>> owned Array, then again in swig into a numpy array). But maybe the
>>> amount of data returned is small, in which case a vector<double>
>>> will be fine.
>>
>> The return of arrays from TimeSeries is typically not performance
>> critical so a copy is fine (which it already is).
>>
>>> * Can't support resize() without ownership, so everywhere resize is
>>> used on a reference argument should pass a vector<double>&
>>> instead. I think this is a few places only.
>>>
>>> At that point, Array will be a plain wrapper of others' data.
>>>
>>> As for your (1), it might well be possible (keep a scratch array in
>>> the GenericFunction class perhaps). But the thing is, the shared
>>> version was no more safe. You can't keep a reference to an Array
>>> because you don't know if it's shared or just wrapping data that
>>> will be deleted from under you. It was "shared or borrowed" rather
>>> than "owned or borrowed", only without the flag to tell the cases
>>> apart.
>>
>> I'm not very happy with our current use of Array. We use std::vector
>> in some places and Array in other places. I don't remember what the
>> original intention was with Array. Is the purpose just to be able to
>> pass shared data as arrays back and forth to Python?
>
> Yes.
>
>> Or should it be a
>> general array class used throughout DOLFIN C++?
>>
>
> No - there is no point in re-inventing std::vector. We should use
> std::vector where possible. It's more flexible and can be used with
> STL algorithms.

Agree.

>> If it is just for Python-wrapping, I think we should make it never own
>> data as Joachim suggests.

Agree. But we need to go over the interface and check that it makes
sense everywhere. TimeSeries could probably just return a std::vector
which is copied in the Python interface.

> Memory has to be created somewhere in order to be used, and it has to
> be cleaned up.
>
> There are two issues here:
>
> 1) Dynamic allocation in order to interface with UFC. Using plain
> pointers rather than smart pointers is likely faster, but I think that
> it's still not a solution. We should try to avoid these allocations in
> loops.
>
> 2) We need to re-visit the NumPy interface. If we can create NumPy
> arrays of fixed length, then we can just use std::vector, and let the
> NumPy array wrap the pointer&x[0]. Maybe Johan Hake can comment on
> this?

The problem is to pass a NumPy array to the C++ interface when it wants
a std::vector&. Then we need to copy the values. A light weight Array
class fixes that issue.

Wrapping a std::vector& into a Numpy array is possible without copying
by just passing the pointer. This would be similar to what we do today
with th...

[Anders Logg (logg)]

On Mon, Feb 27, 2012 at 03:32:04PM -0000, Johan Hake wrote:
> On 02/27/2012 03:01 PM, Garth Wells wrote:
> > On 27 February 2012 13:40, Anders Logg<email address hidden> wrote:
> >> On Fri, Feb 24, 2012 at 07:01:57PM -0000, Joachim Haga wrote:
> >>> I'd opine that the easiest would be to make array never own the data,
> >>> rather than always own the data. In that way, it can wrap anything.
> >>>
> >>> * TimeSeries returns (owned) Arrays in two places. These could
> >>> return references to its own data (const vector<double>&), that
> >>> would save a copy (they are currently copied twice, once into an
> >>> owned Array, then again in swig into a numpy array). But maybe the
> >>> amount of data returned is small, in which case a vector<double>
> >>> will be fine.
> >>
> >> The return of arrays from TimeSeries is typically not performance
> >> critical so a copy is fine (which it already is).
> >>
> >>> * Can't support resize() without ownership, so everywhere resize is
> >>> used on a reference argument should pass a vector<double>&
> >>> instead. I think this is a few places only.
> >>>
> >>> At that point, Array will be a plain wrapper of others' data.
> >>>
> >>> As for your (1), it might well be possible (keep a scratch array in
> >>> the GenericFunction class perhaps). But the thing is, the shared
> >>> version was no more safe. You can't keep a reference to an Array
> >>> because you don't know if it's shared or just wrapping data that
> >>> will be deleted from under you. It was "shared or borrowed" rather
> >>> than "owned or borrowed", only without the flag to tell the cases
> >>> apart.
> >>
> >> I'm not very happy with our current use of Array. We use std::vector
> >> in some places and Array in other places. I don't remember what the
> >> original intention was with Array. Is the purpose just to be able to
> >> pass shared data as arrays back and forth to Python?
> >
> > Yes.
> >
> >> Or should it be a
> >> general array class used throughout DOLFIN C++?
> >>
> >
> > No - there is no point in re-inventing std::vector. We should use
> > std::vector where possible. It's more flexible and can be used with
> > STL algorithms.
>
> Agree.
>
> >> If it is just for Python-wrapping, I think we should make it never own
> >> data as Joachim suggests.
>
> Agree. But we need to go over the interface and check that it makes
> sense everywhere. TimeSeries could probably just return a
> std::vector which is copied in the Python interface.

Yes, definitely for TimeSeries. Perhaps we could settle for that
everytime we want to pass output data to Python?

Then Array input would essentially be limited to a few particular
callbacks that need to be made efficient (GenericFunction::eval).

--
Anders

> > Memory has to be created somewhere in order to be used, and it has
to
> > be cleaned up.
> >
> > There are two issues here:
> >
> > 1) Dynamic allocation in order to interface with UFC. Using plain
> > pointers rather than smart pointers is likely faster, but I think that
> > it's still not a solution. We should try to avoid these allocations in
> > loops.
> >
> > 2) We need to re-visit the NumPy interface. If we can create NumPy
> > arrays of fixed length,...

[Garth Wells (garth-wells)]

On 27 February 2012 17:01, Anders Logg <email address hidden> wrote:
> On Mon, Feb 27, 2012 at 03:32:04PM -0000, Johan Hake wrote:
>> On 02/27/2012 03:01 PM, Garth Wells wrote:
>> > On 27 February 2012 13:40, Anders Logg<email address hidden>  wrote:
>> >> On Fri, Feb 24, 2012 at 07:01:57PM -0000, Joachim Haga wrote:
>> >>> I'd opine that the easiest would be to make array never own the data,
>> >>> rather than always own the data. In that way, it can wrap anything.
>> >>>
>> >>> * TimeSeries returns (owned) Arrays in two places. These could
>> >>> return references to its own data (const vector<double>&), that
>> >>> would save a copy (they are currently copied twice, once into an
>> >>> owned Array, then again in swig into a numpy array). But maybe the
>> >>> amount of data returned is small, in which case a vector<double>
>> >>> will be fine.
>> >>
>> >> The return of arrays from TimeSeries is typically not performance
>> >> critical so a copy is fine (which it already is).
>> >>
>> >>> * Can't support resize() without ownership, so everywhere resize is
>> >>> used on a reference argument should pass a vector<double>&
>> >>> instead. I think this is a few places only.
>> >>>
>> >>> At that point, Array will be a plain wrapper of others' data.
>> >>>
>> >>> As for your (1), it might well be possible (keep a scratch array in
>> >>> the GenericFunction class perhaps). But the thing is, the shared
>> >>> version was no more safe. You can't keep a reference to an Array
>> >>> because you don't know if it's shared or just wrapping data that
>> >>> will be deleted from under you. It was "shared or borrowed" rather
>> >>> than "owned or borrowed", only without the flag to tell the cases
>> >>> apart.
>> >>
>> >> I'm not very happy with our current use of Array. We use std::vector
>> >> in some places and Array in other places. I don't remember what the
>> >> original intention was with Array. Is the purpose just to be able to
>> >> pass shared data as arrays back and forth to Python?
>> >
>> > Yes.
>> >
>> >> Or should it be a
>> >> general array class used throughout DOLFIN C++?
>> >>
>> >
>> > No - there is no point in re-inventing std::vector. We should use
>> > std::vector where possible. It's more flexible and can be used with
>> > STL algorithms.
>>
>> Agree.
>>
>> >> If it is just for Python-wrapping, I think we should make it never own
>> >> data as Joachim suggests.
>>
>> Agree. But we need to go over the interface and check that it makes
>> sense everywhere. TimeSeries could probably just return a
>> std::vector which is copied in the Python interface.
>
> Yes, definitely for TimeSeries. Perhaps we could settle for that
> everytime we want to pass output data to Python?
>
> Then Array input would essentially be limited to a few particular
> callbacks that need to be made efficient (GenericFunction::eval).
>

In the case of GenericFunction::eval we can use std::vector if, in the
wrapping process, we can prevent resizing on the Python side.

Maybe we can remove Array and support (a) std::vectors via copy and
(b) references to a std::vector with resizing prevented.

Garth

> --
> Anders
>
>
>> > Memory has to be created somewhere in order to be used, ...

[Anders Logg (logg)]

On Mon, Feb 27, 2012 at 07:52:27PM -0000, Garth Wells wrote:
> On 27 February 2012 17:01, Anders Logg <email address hidden> wrote:
> > On Mon, Feb 27, 2012 at 03:32:04PM -0000, Johan Hake wrote:
> >> On 02/27/2012 03:01 PM, Garth Wells wrote:
> >> > On 27 February 2012 13:40, Anders Logg<email address hidden>  wrote:
> >> >> On Fri, Feb 24, 2012 at 07:01:57PM -0000, Joachim Haga wrote:
> >> >>> I'd opine that the easiest would be to make array never own the data,
> >> >>> rather than always own the data. In that way, it can wrap anything.
> >> >>>
> >> >>> * TimeSeries returns (owned) Arrays in two places. These could
> >> >>> return references to its own data (const vector<double>&), that
> >> >>> would save a copy (they are currently copied twice, once into an
> >> >>> owned Array, then again in swig into a numpy array). But maybe the
> >> >>> amount of data returned is small, in which case a vector<double>
> >> >>> will be fine.
> >> >>
> >> >> The return of arrays from TimeSeries is typically not performance
> >> >> critical so a copy is fine (which it already is).
> >> >>
> >> >>> * Can't support resize() without ownership, so everywhere resize is
> >> >>> used on a reference argument should pass a vector<double>&
> >> >>> instead. I think this is a few places only.
> >> >>>
> >> >>> At that point, Array will be a plain wrapper of others' data.
> >> >>>
> >> >>> As for your (1), it might well be possible (keep a scratch array in
> >> >>> the GenericFunction class perhaps). But the thing is, the shared
> >> >>> version was no more safe. You can't keep a reference to an Array
> >> >>> because you don't know if it's shared or just wrapping data that
> >> >>> will be deleted from under you. It was "shared or borrowed" rather
> >> >>> than "owned or borrowed", only without the flag to tell the cases
> >> >>> apart.
> >> >>
> >> >> I'm not very happy with our current use of Array. We use std::vector
> >> >> in some places and Array in other places. I don't remember what the
> >> >> original intention was with Array. Is the purpose just to be able to
> >> >> pass shared data as arrays back and forth to Python?
> >> >
> >> > Yes.
> >> >
> >> >> Or should it be a
> >> >> general array class used throughout DOLFIN C++?
> >> >>
> >> >
> >> > No - there is no point in re-inventing std::vector. We should use
> >> > std::vector where possible. It's more flexible and can be used with
> >> > STL algorithms.
> >>
> >> Agree.
> >>
> >> >> If it is just for Python-wrapping, I think we should make it never own
> >> >> data as Joachim suggests.
> >>
> >> Agree. But we need to go over the interface and check that it makes
> >> sense everywhere. TimeSeries could probably just return a
> >> std::vector which is copied in the Python interface.
> >
> > Yes, definitely for TimeSeries. Perhaps we could settle for that
> > everytime we want to pass output data to Python?
> >
> > Then Array input would essentially be limited to a few particular
> > callbacks that need to be made efficient (GenericFunction::eval).
> >
>
> In the case of GenericFunction::eval we can use std::vector if, in
> the wrapping process, we can prevent resizing on the Python side.

Don't know if that's po...

[Johan Hake (johan-hake)]

On 02/27/2012 09:16 PM, Anders Logg wrote:
> On Mon, Feb 27, 2012 at 07:52:27PM -0000, Garth Wells wrote:
>> On 27 February 2012 17:01, Anders Logg<email address hidden> wrote:
>>> On Mon, Feb 27, 2012 at 03:32:04PM -0000, Johan Hake wrote:
>>>> On 02/27/2012 03:01 PM, Garth Wells wrote:
>>>>> On 27 February 2012 13:40, Anders Logg<email address hidden> wrote:
>>>>>> On Fri, Feb 24, 2012 at 07:01:57PM -0000, Joachim Haga wrote:
>>>>>>> I'd opine that the easiest would be to make array never own the data,
>>>>>>> rather than always own the data. In that way, it can wrap anything.
>>>>>>>
>>>>>>> * TimeSeries returns (owned) Arrays in two places. These could
>>>>>>> return references to its own data (const vector<double>&), that
>>>>>>> would save a copy (they are currently copied twice, once into an
>>>>>>> owned Array, then again in swig into a numpy array). But maybe the
>>>>>>> amount of data returned is small, in which case a vector<double>
>>>>>>> will be fine.
>>>>>>
>>>>>> The return of arrays from TimeSeries is typically not performance
>>>>>> critical so a copy is fine (which it already is).
>>>>>>
>>>>>>> * Can't support resize() without ownership, so everywhere resize is
>>>>>>> used on a reference argument should pass a vector<double>&
>>>>>>> instead. I think this is a few places only.
>>>>>>>
>>>>>>> At that point, Array will be a plain wrapper of others' data.
>>>>>>>
>>>>>>> As for your (1), it might well be possible (keep a scratch array in
>>>>>>> the GenericFunction class perhaps). But the thing is, the shared
>>>>>>> version was no more safe. You can't keep a reference to an Array
>>>>>>> because you don't know if it's shared or just wrapping data that
>>>>>>> will be deleted from under you. It was "shared or borrowed" rather
>>>>>>> than "owned or borrowed", only without the flag to tell the cases
>>>>>>> apart.
>>>>>>
>>>>>> I'm not very happy with our current use of Array. We use std::vector
>>>>>> in some places and Array in other places. I don't remember what the
>>>>>> original intention was with Array. Is the purpose just to be able to
>>>>>> pass shared data as arrays back and forth to Python?
>>>>>
>>>>> Yes.
>>>>>
>>>>>> Or should it be a
>>>>>> general array class used throughout DOLFIN C++?
>>>>>>
>>>>>
>>>>> No - there is no point in re-inventing std::vector. We should use
>>>>> std::vector where possible. It's more flexible and can be used with
>>>>> STL algorithms.
>>>>
>>>> Agree.
>>>>
>>>>>> If it is just for Python-wrapping, I think we should make it never own
>>>>>> data as Joachim suggests.
>>>>
>>>> Agree. But we need to go over the interface and check that it makes
>>>> sense everywhere. TimeSeries could probably just return a
>>>> std::vector which is copied in the Python interface.
>>>
>>> Yes, definitely for TimeSeries. Perhaps we could settle for that
>>> everytime we want to pass output data to Python?
>>>
>>> Then Array input would essentially be limited to a few particular
>>> callbacks that need to be made efficient (GenericFunction::eval).
>>>
>>
>> In the case of GenericFunction::eval we can use std::vector if, in
>> the wrapping process, we can prevent resizing on the Python side.
>
> Don't know if ...

[Anders Logg (logg)]

On Tue, Feb 28, 2012 at 05:08:11AM +0100, Johan Hake wrote:
> On 02/27/2012 09:16 PM, Anders Logg wrote:
> >On Mon, Feb 27, 2012 at 07:52:27PM -0000, Garth Wells wrote:
> >>On 27 February 2012 17:01, Anders Logg<email address hidden> wrote:
> >>>On Mon, Feb 27, 2012 at 03:32:04PM -0000, Johan Hake wrote:
> >>>>On 02/27/2012 03:01 PM, Garth Wells wrote:
> >>>>>On 27 February 2012 13:40, Anders Logg<email address hidden> wrote:
> >>>>>>On Fri, Feb 24, 2012 at 07:01:57PM -0000, Joachim Haga wrote:
> >>>>>>>I'd opine that the easiest would be to make array never own the data,
> >>>>>>>rather than always own the data. In that way, it can wrap anything.
> >>>>>>>
> >>>>>>>* TimeSeries returns (owned) Arrays in two places. These could
> >>>>>>>return references to its own data (const vector<double>&), that
> >>>>>>>would save a copy (they are currently copied twice, once into an
> >>>>>>>owned Array, then again in swig into a numpy array). But maybe the
> >>>>>>>amount of data returned is small, in which case a vector<double>
> >>>>>>>will be fine.
> >>>>>>
> >>>>>>The return of arrays from TimeSeries is typically not performance
> >>>>>>critical so a copy is fine (which it already is).
> >>>>>>
> >>>>>>>* Can't support resize() without ownership, so everywhere resize is
> >>>>>>>used on a reference argument should pass a vector<double>&
> >>>>>>>instead. I think this is a few places only.
> >>>>>>>
> >>>>>>>At that point, Array will be a plain wrapper of others' data.
> >>>>>>>
> >>>>>>>As for your (1), it might well be possible (keep a scratch array in
> >>>>>>>the GenericFunction class perhaps). But the thing is, the shared
> >>>>>>>version was no more safe. You can't keep a reference to an Array
> >>>>>>>because you don't know if it's shared or just wrapping data that
> >>>>>>>will be deleted from under you. It was "shared or borrowed" rather
> >>>>>>>than "owned or borrowed", only without the flag to tell the cases
> >>>>>>>apart.
> >>>>>>
> >>>>>>I'm not very happy with our current use of Array. We use std::vector
> >>>>>>in some places and Array in other places. I don't remember what the
> >>>>>>original intention was with Array. Is the purpose just to be able to
> >>>>>>pass shared data as arrays back and forth to Python?
> >>>>>
> >>>>>Yes.
> >>>>>
> >>>>>>Or should it be a
> >>>>>>general array class used throughout DOLFIN C++?
> >>>>>>
> >>>>>
> >>>>>No - there is no point in re-inventing std::vector. We should use
> >>>>>std::vector where possible. It's more flexible and can be used with
> >>>>>STL algorithms.
> >>>>
> >>>>Agree.
> >>>>
> >>>>>>If it is just for Python-wrapping, I think we should make it never own
> >>>>>>data as Joachim suggests.
> >>>>
> >>>>Agree. But we need to go over the interface and check that it makes
> >>>>sense everywhere. TimeSeries could probably just return a
> >>>>std::vector which is copied in the Python interface.
> >>>
> >>>Yes, definitely for TimeSeries. Perhaps we could settle for that
> >>>everytime we want to pass output data to Python?
> >>>
> >>>Then Array input would essentially be limited to a few particular
> >>>callbacks that need to be made efficient (GenericFunction::eval).
> >>>
> >>
> >>In the cas...

[lp:~jobh/dolfin/fast-array]

6593. By Joachim Haga on 2012-03-01

Fix missing return statement

6594. By Joachim Haga on 2012-03-01

Merge with trunk

6595. By Joachim Haga on 2012-03-01

Change TimeSeries interface to return std::vector rather than Array

[Joachim Haga (jobh)]

This seems to have stalled without a consensus. Still, I think the improvements are too dramatic to just drop. Assembling a "v*f*dx" RHS on a 64^3 cube (twice):

6.95 s before
2.30 s after

With Garth's suggestion of reusing an Array: 2.54s, but not thread safe without some extra complexity (likely compiler-dependent).

I'd like to repeat my argument that my proposed change isn't any less safe than what's there already. The old Array used reference_to_no_delete_ptr in several constructors, thus the memory was managed by the caller -- exactly the same as it is here, if the Array is not the owner. Arrays were never shared [*], so the "shared" variant is exactly the same as the new "owned" variant.

Note [*]: It was used implicitly to transfer ownership in the TimeSeries return values. Since there's agreement than an extra copy does no harm there, I've changed those to return std::vector instead.

Then it can be deprecated or redesigned or whatever later, but in the meantime dolfin will be much faster :)

PS. The "interface version" thing seems to not work for all instant files, so instant-clean is still required. How should this be handled?

[Johan Hake (johan-hake)]

I am in favor of applying the patch.

On Tue, Mar 6, 2012 at 12:20 PM, Joachim Haga <email address hidden> wrote:
> This seems to have stalled without a consensus. Still, I think the improvements are too dramatic to just drop. Assembling a "v*f*dx" RHS on a 64^3 cube (twice):
>
> 6.95 s before
> 2.30 s after
>
> With Garth's suggestion of reusing an Array: 2.54s, but not thread safe without some extra complexity (likely compiler-dependent).
>
> I'd like to repeat my argument that my proposed change isn't any less safe than what's there already. The old Array used reference_to_no_delete_ptr in several constructors, thus the memory was managed by the caller -- exactly the same as it is here, if the Array is not the owner. Arrays were never shared [*], so the "shared" variant is exactly the same as the new "owned" variant.
>
> Note [*]: It was used implicitly to transfer ownership in the TimeSeries return values. Since there's agreement than an extra copy does no harm there, I've changed those to return std::vector instead.
>
> Then it can be deprecated or redesigned or whatever later, but in the meantime dolfin will be much faster :)
>
> PS. The "interface version" thing seems to not work for all instant files, so instant-clean is still required. How should this be handled?

Not sure what is not working, but it should be enough to change the
DOLFIN_VERSION_MICRO, as this is eventualy included in the md5sum
which determines if a modules should be recompiled or not.

Johan

> --
> https://code.launchpad.net/~jobh/dolfin/fast-array/+merge/94467
> Your team DOLFIN Core Team is requested to review the proposed merge of lp:~jobh/dolfin/fast-array into lp:dolfin.

[Anders Logg (logg)]

I'm also very much in favor.

--
Anders

On Tue, Mar 06, 2012 at 11:53:31AM -0000, Johan Hake wrote:
> I am in favor of applying the patch.
>
> On Tue, Mar 6, 2012 at 12:20 PM, Joachim Haga <email address hidden> wrote:
> > This seems to have stalled without a consensus. Still, I think the improvements are too dramatic to just drop. Assembling a "v*f*dx" RHS on a 64^3 cube (twice):
> >
> > 6.95 s before
> > 2.30 s after
> >
> > With Garth's suggestion of reusing an Array: 2.54s, but not thread safe without some extra complexity (likely compiler-dependent).
> >
> > I'd like to repeat my argument that my proposed change isn't any less safe than what's there already. The old Array used reference_to_no_delete_ptr in several constructors, thus the memory was managed by the caller -- exactly the same as it is here, if the Array is not the owner. Arrays were never shared [*], so the "shared" variant is exactly the same as the new "owned" variant.
> >
> > Note [*]: It was used implicitly to transfer ownership in the TimeSeries return values. Since there's agreement than an extra copy does no harm there, I've changed those to return std::vector instead.
> >
> > Then it can be deprecated or redesigned or whatever later, but in the meantime dolfin will be much faster :)
> >
> > PS. The "interface version" thing seems to not work for all instant files, so instant-clean is still required. How should this be handled?
>
> Not sure what is not working, but it should be enough to change the
> DOLFIN_VERSION_MICRO, as this is eventualy included in the md5sum
> which determines if a modules should be recompiled or not.
>
> Johan
>
>

[Garth Wells (garth-wells)]

On 6 March 2012 12:11, Anders Logg <email address hidden> wrote:
> I'm also very much in favor.
>

I'm happy in principle, but want some clarity on Array. Can we make it
that Array is always a view (i.e., never owns that data), and
therefore cannot be resized? This would clean up some const hacks.

Garth

> --
> Anders
>
> On Tue, Mar 06, 2012 at 11:53:31AM -0000, Johan Hake wrote:
>> I am in favor of applying the patch.
>>
>> On Tue, Mar 6, 2012 at 12:20 PM, Joachim Haga <email address hidden> wrote:
>> > This seems to have stalled without a consensus. Still, I think the improvements are too dramatic to just drop. Assembling a "v*f*dx" RHS on a 64^3 cube (twice):
>> >
>> > 6.95 s before
>> > 2.30 s after
>> >
>> > With Garth's suggestion of reusing an Array: 2.54s, but not thread safe without some extra complexity (likely compiler-dependent).
>> >
>> > I'd like to repeat my argument that my proposed change isn't any less safe than what's there already. The old Array used reference_to_no_delete_ptr in several constructors, thus the memory was managed by the caller -- exactly the same as it is here, if the Array is not the owner. Arrays were never shared [*], so the "shared" variant is exactly the same as the new "owned" variant.
>> >
>> > Note [*]: It was used implicitly to transfer ownership in the TimeSeries return values. Since there's agreement than an extra copy does no harm there, I've changed those to return std::vector instead.
>> >
>> > Then it can be deprecated or redesigned or whatever later, but in the meantime dolfin will be much faster :)
>> >
>> > PS. The "interface version" thing seems to not work for all instant files, so instant-clean is still required. How should this be handled?
>>
>> Not sure what is not working, but it should be enough to change the
>> DOLFIN_VERSION_MICRO, as this is eventualy included in the md5sum
>> which determines if a modules should be recompiled or not.
>>
>> Johan
>>
>>
>
> --
> https://code.launchpad.net/~jobh/dolfin/fast-array/+merge/94467
> Your team DOLFIN Core Team is requested to review the proposed merge of lp:~jobh/dolfin/fast-array into lp:dolfin.

[Anders Logg (logg)]

On Tue, Mar 06, 2012 at 12:28:20PM -0000, Garth Wells wrote:
> On 6 March 2012 12:11, Anders Logg <email address hidden> wrote:
> > I'm also very much in favor.
> >
>
> I'm happy in principle, but want some clarity on Array. Can we make it
> that Array is always a view (i.e., never owns that data), and
> therefore cannot be resized? This would clean up some const hacks.

That's fine with me. And in addition we should limit the use of Array
to just the few places where we need it for passing data through
the SWIG layer without copying.

--
Anders

> Garth
>
> >
> > On Tue, Mar 06, 2012 at 11:53:31AM -0000, Johan Hake wrote:
> >> I am in favor of applying the patch.
> >>
> >> On Tue, Mar 6, 2012 at 12:20 PM, Joachim Haga <email address hidden> wrote:
> >> > This seems to have stalled without a consensus. Still, I think the improvements are too dramatic to just drop. Assembling a "v*f*dx" RHS on a 64^3 cube (twice):
> >> >
> >> > 6.95 s before
> >> > 2.30 s after
> >> >
> >> > With Garth's suggestion of reusing an Array: 2.54s, but not thread safe without some extra complexity (likely compiler-dependent).
> >> >
> >> > I'd like to repeat my argument that my proposed change isn't any less safe than what's there already. The old Array used reference_to_no_delete_ptr in several constructors, thus the memory was managed by the caller -- exactly the same as it is here, if the Array is not the owner. Arrays were never shared [*], so the "shared" variant is exactly the same as the new "owned" variant.
> >> >
> >> > Note [*]: It was used implicitly to transfer ownership in the TimeSeries return values. Since there's agreement than an extra copy does no harm there, I've changed those to return std::vector instead.
> >> >
> >> > Then it can be deprecated or redesigned or whatever later, but in the meantime dolfin will be much faster :)
> >> >
> >> > PS. The "interface version" thing seems to not work for all instant files, so instant-clean is still required. How should this be handled?
> >>
> >> Not sure what is not working, but it should be enough to change the
> >> DOLFIN_VERSION_MICRO, as this is eventualy included in the md5sum
> >> which determines if a modules should be recompiled or not.
> >>
> >> Johan
> >>
> >>
> >
>

[Johan Hake (johan-hake)]

On Tue, Mar 6, 2012 at 1:28 PM, Garth Wells <email address hidden> wrote:
> On 6 March 2012 12:11, Anders Logg <email address hidden> wrote:
>> I'm also very much in favor.
>>
>
> I'm happy in principle, but want some clarity on Array. Can we make it
> that Array is always a view (i.e., never owns that data), and
> therefore cannot be resized? This would clean up some const hacks.

I agree. For all argments needing to be resized should use
std::vector. Not sure where the resize functionality is/was used.

Johan

[Joachim Haga (jobh)]

>
> > I'm happy in principle, but want some clarity on Array. Can we make it
> > that Array is always a view (i.e., never owns that data), and
> > therefore cannot be resized? This would clean up some const hacks.
>
> I agree. For all argments needing to be resized should use
> std::vector. Not sure where the resize functionality is/was used.

It is possible (of course). There are quite a few user-exposed interfaces
that need changing (f.x. all the GeneralVector::set/get methods, since it
would be strange to just change the get methods).

And it may be a lot of grunt work. But I'll have a go at it if you say
that's the way to go (even if it requires change to user code).

-j.

[Joachim Haga (jobh)]

After looking at la_get_set_items.i: _get_vector_values(), shown below, I
think I'll punt on this. It's just too much work for an uncertain outcome
(it is quite possible that a lot of complexity has to be added back to deal
with cases like this, making this pointless).

// Returns the values from a Vector.
dolfin::Array<double>& _get_vector_values( dolfin::GenericVector* self)
{
  dolfin::Array<double>* values;
  try
  {
    // Try accessing the value pointer directly [jobh: return borrowed data]
    double* data_values = self->data();
    values = new dolfin::Array<double>(self->size(), data_values);
  }
  catch (std::runtime_error e)
  {
    // We couldn't access the values directly [jobh: return owned data]
    values = new dolfin::Array<double>(self->size());
    self->get_local(*values);
  }
  return *values;
}

On 6 March 2012 15:01, Joachim Berdal Haga <email address hidden> wrote:

> > I'm happy in principle, but want some clarity on Array. Can we make it
>> > that Array is always a view (i.e., never owns that data), and
>> > therefore cannot be resized? This would clean up some const hacks.
>>
>> I agree. For all argments needing to be resized should use
>> std::vector. Not sure where the resize functionality is/was used.
>
>
> It is possible (of course). There are quite a few user-exposed interfaces
> that need changing (f.x. all the GeneralVector::set/get methods, since it
> would be strange to just change the get methods).
>
> And it may be a lot of grunt work. But I'll have a go at it if you say
> that's the way to go (even if it requires change to user code).
>
> -j.
>
>

[Garth Wells (garth-wells)]

On 6 March 2012 15:13, Joachim Haga <email address hidden> wrote:
> After looking at la_get_set_items.i: _get_vector_values(), shown below, I
> think I'll punt on this. It's just too much work for an uncertain outcome
> (it is quite possible that a lot of complexity has to be added back to deal
> with cases like this, making this pointless).
>
> // Returns the values from a Vector.
> dolfin::Array<double>& _get_vector_values( dolfin::GenericVector* self)
> {
>  dolfin::Array<double>* values;
>  try
>  {
>    // Try accessing the value pointer directly [jobh: return borrowed data]
>    double* data_values = self->data();
>    values = new dolfin::Array<double>(self->size(), data_values);
>  }
>  catch (std::runtime_error e)
>  {
>    // We couldn't access the values directly [jobh: return owned data]
>    values = new dolfin::Array<double>(self->size());
>    self->get_local(*values);
>  }
>  return *values;
> }
>

We should just be able to use std::vector in the C++ interface for this.

Garth

>
> On 6 March 2012 15:01, Joachim Berdal Haga <email address hidden> wrote:
>
>> > I'm happy in principle, but want some clarity on Array. Can we make it
>>> > that Array is always a view (i.e., never owns that data), and
>>> > therefore cannot be resized? This would clean up some const hacks.
>>>
>>> I agree. For all argments needing to be resized should use
>>> std::vector. Not sure where the resize functionality is/was used.
>>
>>
>> It is possible (of course). There are quite a few user-exposed interfaces
>> that need changing (f.x. all the GeneralVector::set/get methods, since it
>> would be strange to just change the get methods).
>>
>> And it may be a lot of grunt work. But I'll have a go at it if you say
>> that's the way to go (even if it requires change to user code).
>>
>> -j.
>>
>>
>
> --
> https://code.launchpad.net/~jobh/dolfin/fast-array/+merge/94467
> Your team DOLFIN Core Team is requested to review the proposed merge of lp:~jobh/dolfin/fast-array into lp:dolfin.

[Anders Logg (logg)]

On Tue, Mar 06, 2012 at 02:10:28PM -0000, Joachim Haga wrote:
> >
> > > I'm happy in principle, but want some clarity on Array. Can we make it
> > > that Array is always a view (i.e., never owns that data), and
> > > therefore cannot be resized? This would clean up some const hacks.
> >
> > I agree. For all argments needing to be resized should use
> > std::vector. Not sure where the resize functionality is/was used.
>
>
> It is possible (of course). There are quite a few user-exposed interfaces
> that need changing (f.x. all the GeneralVector::set/get methods, since it
> would be strange to just change the get methods).
>
> And it may be a lot of grunt work. But I'll have a go at it if you say
> that's the way to go (even if it requires change to user code).

I suggest we merge this now, and think about cleaning up the use of
Array later.

Just to check: do we break any user interfaces with this patch? I
assume not.

--
Anders

[Johan Hake (johan-hake)]

On Tue, Mar 6, 2012 at 3:10 PM, Joachim Haga <email address hidden> wrote:
>>
>> > I'm happy in principle, but want some clarity on Array. Can we make it
>> > that Array is always a view (i.e., never owns that data), and
>> > therefore cannot be resized? This would clean up some const hacks.
>>
>> I agree. For all argments needing to be resized should use
>> std::vector. Not sure where the resize functionality is/was used.
>
> It is possible (of course). There are quite a few user-exposed interfaces
> that need changing (f.x. all the GeneralVector::set/get methods, since it
> would be strange to just change the get methods).

I did not anticipate using std::vector for these methods as these
methods assume the user provides arrays with the correct size. It
would cost too much to resize std::vector for these operations. Maybe
we can use Array<double> for these but for now they just work with
NumPy arrays directly.

Johan

> And it may be a lot of grunt work. But I'll have a go at it if you say
> that's the way to go (even if it requires change to user code).

> -j.
>
> --
> https://code.launchpad.net/~jobh/dolfin/fast-array/+merge/94467
> Your team DOLFIN Core Team is requested to review the proposed merge of lp:~jobh/dolfin/fast-array into lp:dolfin.

[Johan Hake (johan-hake)]

On Tue, Mar 6, 2012 at 8:19 PM, Anders Logg <email address hidden> wrote:
> On Tue, Mar 06, 2012 at 02:10:28PM -0000, Joachim Haga wrote:
>> >
>> > > I'm happy in principle, but want some clarity on Array. Can we make it
>> > > that Array is always a view (i.e., never owns that data), and
>> > > therefore cannot be resized? This would clean up some const hacks.
>> >
>> > I agree. For all argments needing to be resized should use
>> > std::vector. Not sure where the resize functionality is/was used.
>>
>>
>> It is possible (of course). There are quite a few user-exposed interfaces
>> that need changing (f.x. all the GeneralVector::set/get methods, since it
>> would be strange to just change the get methods).
>>
>> And it may be a lot of grunt work. But I'll have a go at it if you say
>> that's the way to go (even if it requires change to user code).
>
> I suggest we merge this now, and think about cleaning up the use of
> Array later.

I agree. I do not have a working DOLFIN right now so I will
unfortunately not be able to do it.

Johan

> Just to check: do we break any user interfaces with this patch? I
> assume not.
>
> --
> Anders
>
> --
> https://code.launchpad.net/~jobh/dolfin/fast-array/+merge/94467
> Your team DOLFIN Core Team is requested to review the proposed merge of lp:~jobh/dolfin/fast-array into lp:dolfin.

[Anders Logg (logg)]

On Tue, Mar 06, 2012 at 08:16:47PM -0000, Johan Hake wrote:
> On Tue, Mar 6, 2012 at 8:19 PM, Anders Logg <email address hidden> wrote:
> > On Tue, Mar 06, 2012 at 02:10:28PM -0000, Joachim Haga wrote:
> >> >
> >> > > I'm happy in principle, but want some clarity on Array. Can we make it
> >> > > that Array is always a view (i.e., never owns that data), and
> >> > > therefore cannot be resized? This would clean up some const hacks.
> >> >
> >> > I agree. For all argments needing to be resized should use
> >> > std::vector. Not sure where the resize functionality is/was used.
> >>
> >>
> >> It is possible (of course). There are quite a few user-exposed interfaces
> >> that need changing (f.x. all the GeneralVector::set/get methods, since it
> >> would be strange to just change the get methods).
> >>
> >> And it may be a lot of grunt work. But I'll have a go at it if you say
> >> that's the way to go (even if it requires change to user code).
> >
> > I suggest we merge this now, and think about cleaning up the use of
> > Array later.
>
> I agree. I do not have a working DOLFIN right now so I will
> unfortunately not be able to do it.

I can merge if there are no objections.

--
Anders

> Johan
>
> > Just to check: do we break any user interfaces with this patch? I
> > assume not.
> >
> >
>

[Joachim Haga (jobh)]

>
> >> I agree. For all argments needing to be resized should use
> >> std::vector. Not sure where the resize functionality is/was used.
> >
> > It is possible (of course). There are quite a few user-exposed interfaces
> > that need changing (f.x. all the GeneralVector::set/get methods, since it
> > would be strange to just change the get methods).
>
> I did not anticipate using std::vector for these methods as these
> methods assume the user provides arrays with the correct size. It
> would cost too much to resize std::vector for these operations. Maybe
>

Not sure if we're talking about the same thing? PETScVector::get_local for
example uses resize() on the Array& that is passed as parameter.

[Joachim Haga (jobh)]

On 6 March 2012 21:31, Anders Logg <email address hidden> wrote:

> > > I suggest we merge this now, and think about cleaning up the use of
> > > Array later.
> >
> > I agree. I do not have a working DOLFIN right now so I will
> > unfortunately not be able to do it.
>
> I can merge if there are no objections.
>

Certainly not from me! But hold on a bit, tomorrow morning I will check the
instant issue, and see if it's possible to disable the copy constructor (to
ensure it's not used in unanticipated ways).

-j.

[Johan Hake (johan-hake)]

On Tue, Mar 6, 2012 at 10:04 PM, Joachim Haga <email address hidden> wrote:
>>
>> >> I agree. For all argments needing to be resized should use
>> >> std::vector. Not sure where the resize functionality is/was used.
>> >
>> > It is possible (of course). There are quite a few user-exposed interfaces
>> > that need changing (f.x. all the GeneralVector::set/get methods, since it
>> > would be strange to just change the get methods).
>>
>> I did not anticipate using std::vector for these methods as these
>> methods assume the user provides arrays with the correct size. It
>> would cost too much to resize std::vector for these operations. Maybe
>>
>
> Not sure if we're talking about the same thing? PETScVector::get_local for
> example uses resize() on the Array& that is passed as parameter.

You are right. We do not talk about the same thing.

Here I think we should just check that the user pass an Array of the
correct size avoiding resizing.

Johan

[Joachim Haga (jobh)]

On 6 March 2012 17:55, Garth Wells <email address hidden> wrote:

> On 6 March 2012 15:13, Joachim Haga <email address hidden> wrote:
> > After looking at la_get_set_items.i: _get_vector_values(), shown below, I
> > think I'll punt on this. It's just too much work for an uncertain outcome
> > (it is quite possible that a lot of complexity has to be added back to
> deal
> > with cases like this, making this pointless).
> >
> > // Returns the values from a Vector.
> > dolfin::Array<double>& _get_vector_values( dolfin::GenericVector* self)
>
> We should just be able to use std::vector in the C++ interface for this.
>

Hm, yes, I guess an extra copy won't matter that much here. I'll see how it
goes.

-j.

[Johan Hake (johan-hake)]

On Tue, Mar 6, 2012 at 10:50 PM, Joachim Haga <email address hidden> wrote:
> On 6 March 2012 17:55, Garth Wells <email address hidden> wrote:
>
>> On 6 March 2012 15:13, Joachim Haga <email address hidden> wrote:
>> > After looking at la_get_set_items.i: _get_vector_values(), shown below, I
>> > think I'll punt on this. It's just too much work for an uncertain outcome
>> > (it is quite possible that a lot of complexity has to be added back to
>> deal
>> > with cases like this, making this pointless).
>> >
>> > // Returns the values from a Vector.
>> > dolfin::Array<double>& _get_vector_values( dolfin::GenericVector* self)
>>
>> We should just be able to use std::vector in the C++ interface for this.
>>
>
> Hm, yes, I guess an extra copy won't matter that much here. I'll see how it
> goes.

We can change this syntax to:

  void _get_vector_values(dolfin::GenericVector* self,
dolfin::Array<double>& values)

where values now is just a view of a numpy array. In the Python layer
we then create the NumPy array of correct size, call the function and
populate the values.

Johan

[lp:~jobh/dolfin/fast-array]

6596. By Joachim Haga on 2012-03-07

Disable Array copy constructor, remove utility methods that are not used

6597. By Joachim Haga on 2012-03-07

Merge with trunk

[Joachim Haga (jobh)]

I've now pushed a version which is merged with trunk (the merge wasn't clean so I thought it wise) and has been fully tested. The recent changes are:

- Disable the Array copy constructor, to make it safer (no unanticipated sharing or copying).
- Allow resize() to same size even for un-owned Arrays, it turned out that the ale code used this.
- Remove unused utility methods (Array::zero, Array::min, etc)

It passes all tests after merge. (KrylovSolver and one CGAL demo fail, but that's not new.)

Johan: The interface version I mentioned is something I added to compilemodules.py which changes the hash without changing the visible version. If you think it better to change DOLFIN_VERSION_MICRO, I'll do that instead (but if so, what comes after "+"?)

[Anders Logg (logg)]

I'll merge this later if there are now objections.

--
Anders

On Wed, Mar 07, 2012 at 01:28:34PM -0000, Joachim Haga wrote:
> I've now pushed a version which is merged with trunk (the merge wasn't clean so I thought it wise) and has been fully tested. The recent changes are:
>
> - Disable the Array copy constructor, to make it safer (no unanticipated sharing or copying).
> - Allow resize() to same size even for un-owned Arrays, it turned out that the ale code used this.
> - Remove unused utility methods (Array::zero, Array::min, etc)
>
> It passes all tests after merge. (KrylovSolver and one CGAL demo fail, but that's not new.)
>
> Johan: The interface version I mentioned is something I added to compilemodules.py which changes the hash without changing the visible version. If you think it better to change DOLFIN_VERSION_MICRO, I'll do that instead (but if so, what comes after "+"?)

[Anders Logg (logg)]

Merged.

Joachim, could you send a patch for ChangeLog to include one-line
summaries of all your changes. Thanks.

--
Anders

On Wed, Mar 07, 2012 at 01:28:34PM -0000, Joachim Haga wrote:
> I've now pushed a version which is merged with trunk (the merge wasn't clean so I thought it wise) and has been fully tested. The recent changes are:
>
> - Disable the Array copy constructor, to make it safer (no unanticipated sharing or copying).
> - Allow resize() to same size even for un-owned Arrays, it turned out that the ale code used this.
> - Remove unused utility methods (Array::zero, Array::min, etc)
>
> It passes all tests after merge. (KrylovSolver and one CGAL demo fail, but that's not new.)
>
> Johan: The interface version I mentioned is something I added to compilemodules.py which changes the hash without changing the visible version. If you think it better to change DOLFIN_VERSION_MICRO, I'll do that instead (but if so, what comes after "+"?)

[Joachim Haga (jobh)]

=== modified file 'ChangeLog'
--- ChangeLog 2012-02-13 21:49:10 +0000
+++ ChangeLog 2012-03-08 12:12:20 +0000
@@ -1,3 +1,8 @@
+ - MPI functionality for distributing values between neighbours
+ - SystemAssembler now works in parallel with topological/geometric
boundary search
+ - New symmetric assembler with ability for stand-alone RHS assemble
+ - Major speed-up of DirichletBC computation and mesh marking
+ - Major speed-up of assembly of functions and expressions
  - Major speed-up of mesh topology computation

On 8 March 2012 11:46, Anders Logg <email address hidden> wrote:
> Merged.
>
> Joachim, could you send a patch for ChangeLog to include one-line
> summaries of all your changes. Thanks.
>
> --
> Anders
>
>
> On Wed, Mar 07, 2012 at 01:28:34PM -0000, Joachim Haga wrote:
>> I've now pushed a version which is merged with trunk (the merge wasn't clean so I thought it wise) and has been fully tested. The recent changes are:
>>
>> - Disable the Array copy constructor, to make it safer (no unanticipated sharing or copying).
>> - Allow resize() to same size even for un-owned Arrays, it turned out that the ale code used this.
>> - Remove unused utility methods (Array::zero, Array::min, etc)
>>
>> It passes all tests after merge. (KrylovSolver and one CGAL demo fail, but that's not new.)
>>
>> Johan: The interface version I mentioned is something I added to compilemodules.py which changes the hash without changing the visible version. If you think it better to change DOLFIN_VERSION_MICRO, I'll do that instead (but if so, what comes after "+"?)

=== modified file 'dolfin/adaptivity/TimeSeries.cpp'
--- dolfin/adaptivity/TimeSeries.cpp	2012-03-01 12:01:06 +0000
+++ dolfin/adaptivity/TimeSeries.cpp	2012-03-07 13:15:44 +0000
@@ -23,6 +23,7 @@
 #include <sstream>
 #include <boost/scoped_ptr.hpp>
 
+#include <dolfin/common/constants.h>
 #include <dolfin/io/File.h>
 #include <dolfin/io/BinaryFile.h>
 #include <dolfin/la/GenericVector.h>
@@ -207,20 +208,14 @@
   file >> mesh;
 }
 //-----------------------------------------------------------------------------
-Array<double> TimeSeries::vector_times() const
+std::vector<double> TimeSeries::vector_times() const
 {
-  Array<double> times(_vector_times.size());
-  for (uint i = 0; i < _vector_times.size(); i++)
-    times[i] = _vector_times[i];
-  return times;
+  return _vector_times;
 }
 //-----------------------------------------------------------------------------
-Array<double> TimeSeries::mesh_times() const
+std::vector<double> TimeSeries::mesh_times() const
 {
-  Array<double> times(_mesh_times.size());
-  for (uint i = 0; i < _mesh_times.size(); i++)
-    times[i] = _mesh_times[i];
-  return times;
+  return _mesh_times;
 }
 //-----------------------------------------------------------------------------
 void TimeSeries::clear()
 
=== modified file 'dolfin/adaptivity/TimeSeries.h'
--- dolfin/adaptivity/TimeSeries.h	2011-10-24 04:11:41 +0000
+++ dolfin/adaptivity/TimeSeries.h	2012-03-07 13:15:44 +0000
@@ -22,7 +22,7 @@
 #define __TIME_SERIES_H
 
 #include <string>
-#include <dolfin/common/Array.h>
+#include <vector>
 #include <dolfin/common/Variable.h>
 
 namespace dolfin
@@ -104,16 +104,16 @@
     /// Return array of sample times for vectors
     ///
     /// *Returns*
-    ///     _Array_ <double>
+    ///     std::vector<double>
     ///         The times.
-    Array<double> vector_times() const;
+    std::vector<double> vector_times() const;
 
     /// Return array of sample times for meshes
     ///
     /// *Returns*
-    ///     _Array_ <double>
+    ///     std::vector<double>
     ///         The times.
-    Array<double> mesh_times() const;
+    std::vector<double> mesh_times() const;
 
     /// Clear time series
     void clear();
 
=== modified file 'dolfin/ale/HarmonicSmoothing.cpp'
--- dolfin/ale/HarmonicSmoothing.cpp	2012-02-01 07:31:55 +0000
+++ dolfin/ale/HarmonicSmoothing.cpp	2012-03-07 13:15:44 +0000
@@ -108,7 +108,7 @@
     solve(A, x, b, "gmres", prec);
 
     // Get new coordinates
-    Array<double> _new_coordinates(N, new_coordinates.data().get() + dim*N);
+    Array<double> _new_coordinates(N, new_coordinates.data() + dim*N);
     x.get_local(_new_coordinates);
   }
 
 
=== modified file 'dolfin/common/Array.h'
--- dolfin/common/Array.h	2012-03-01 12:01:06 +0000
+++ dolfin/common/Array.h	2012-03-07 13:15:44 +0000
@@ -16,9 +16,10 @@
 // along with DOLFIN. If not, see <http://www.gnu.org/licenses/>.
 //
 // Modified by Anders Logg 2010-2011
+// Modified by Joachim B Haga 2012
 //
 // First added:  2009-12-06
-// Last changed: 2011-11-13
+// Last changed: 2012-02-23
 
 #ifndef __DOLFIN_ARRAY_H
 #define __DOLFIN_ARRAY_H
@@ -26,8 +27,7 @@
 #include <sstream>
 #include <string>
 #include <utility>
-#include <boost/shared_array.hpp>
-
+#include <vector>
 
 #include <dolfin/common/constants.h>
 #include <dolfin/common/types.h>
@@ -44,55 +44,46 @@
 
   template <typename T> class Array
   {
+
   public:
 
-    /// Create empty array
-    Array() : _size(0), x(0) {}
-
-    /// Create array of size N
-    explicit Array(uint N) : _size(N), x(new T[N]) {}
-
-    /// Copy constructor (arg name need to have a different name that 'x')
-    Array(const Array& other) : _size(0), x(0)
-    { *this = other; }
-
-    /// Construct array from a shared pointer
-    Array(uint N, boost::shared_array<T> x) : _size(N), x(x) {}
-
-    /// Construct array from a pointer. Array will not take ownership.
-    Array(uint N, T* x) : _size(N), x(boost::shared_array<T>(x, NoDeleter())) {}
-
-    /// Assignment operator
-    const Array& operator= (const Array& x)
+    /// Create an empty array, with ownership. Must be resized or assigned to to be of any use.
+    explicit Array() : _size(0), _x(NULL), _owner(true) {}
+
+    /// Create array of size N. Array has ownership.
+    explicit Array(uint N) : _size(N), _x(new T[N]), _owner(true) {}
+
+    /// Create an array wrapping a std::vector. Array does not take ownership.
+    explicit Array(std::vector<T>& vec) : _size(vec.size()), _x(&vec[0]), _owner(false) {}
+
+    /// Construct array from a pointer. Array does not take ownership.
+    Array(uint N, T* x) : _size(N), _x(x), _owner(false) {}
+
+    /// Destructor.
+    ~Array()
     {
-      // Resize if necessary
-      if (x.empty() && !x.x)
-      {
-        this->x.reset();
-        this->_size = 0;
-      }
-      else if (this->_size != x.size())
-      {
-        this->x.reset(new T[x.size()]);
-        this->_size = x.size();
-      }
-
-      // Copy data
-      if (_size > 0)
-      {
-        dolfin_assert(this->x);
-        dolfin_assert(x.x);
-        std::copy(&x.x[0], &x.x[_size], &this->x[0]);
-      }
-
-      return *this;
+      if (_owner && _x)
+        delete[] _x;
     }
 
-    /// Construct array from a pointer. Array will not take ownership.
-    void update(uint N, T* _x)
+    /// Resize the array. If the new size if different from the old size, the
+    /// Array must be owner; the old contents are then lost.
+    void resize(uint N)
     {
+      if (_size == N)
+        return;
+      if (!_owner)
+        dolfin_error("Array.h", "resize", "Only owned arrays can be resized");
+      if (_x)
+        delete[] _x;
+      _x = (N == 0 ? NULL : new T[N]);
       _size = N;
-      x.reset(_x, NoDeleter());
+    }
+
+    /// Clear the array (resize to 0).
+    void clear()
+    {
+      resize(0);
     }
 
     /// Return informal string representation (pretty-print).
@@ -116,101 +107,64 @@
       return s.str();
     }
 
-    /// Clear array
-    void clear()
-    {
-      this->x.reset();
-      this->_size = 0;
-    }
-
-    /// Resize array to size N. If size changes, contents will be destroyed.
-    void resize(uint N)
-    {
-      // Special case
-      if (N == _size)
-        return;
-
-      // Special case
-      if (N == 0)
-      {
-        clear();
-        return;
-      }
-
-      // FIXME: Do we want to allow resizing of shared data?
-      if (x.unique())
-      {
-        _size = N;
-        x.reset(new T[N]);
-      }
-      else
-      {
-        dolfin_error("Array.h",
-                     "resize Array",
-                     "Data is shared");
-      }
-    }
-
-    /// Return true if empty
-    bool empty() const
-    { return _size == 0; }
-
     /// Return size of array
     uint size() const
     { return _size; }
 
-    /// Zero array
-    void zero()
-    { dolfin_assert(x); std::fill(&x[0], &x[_size], 0.0); }
-
     /// Set entries which meet (abs(x[i]) < eps) to zero
     void zero_eps(double eps=DOLFIN_EPS);
 
-    /// Return minimum value of array
-    T min() const
-    { dolfin_assert(x); return *std::min_element(&x[0], &x[_size]); }
-
-    /// Return maximum value of array
-    T max() const
-    { dolfin_assert(x); return *std::max_element(&x[0], &x[_size]); }
-
     /// Access value of given entry (const version)
     const T& operator[] (uint i) const
-    { dolfin_assert(x); dolfin_assert(i < _size); return x[i]; }
+    { dolfin_assert(i < _size); return _x[i]; }
 
     /// Access value of given entry (non-const version)
     T& operator[] (uint i)
+    { dolfin_assert(i < _size); return _x[i]; }
+
+    /// Assignment operator. If resize is required, the Array must be owner.
+    Array& operator= (const Array& other)
     {
-      dolfin_assert(x);
-      dolfin_assert(i < _size);
-      return x[i];
+      resize(other._size);
+      if (_size > 0)
+        std::copy(&other._x[0], &other._x[_size], &_x[0]);
+      return *this;
     }
 
-    /// Assign value to all entries
-    const Array<T>& operator= (T& x)
+    /// Assignment operator from std::vector. If resize is required, the Array must be owner.
+    Array& operator= (const std::vector<T>& other)
     {
-      dolfin_assert(this->x);
-      for (uint i = 0; i < _size; ++i)
-        this->x[i] = x;
+      resize(other.size());
+      if (_size > 0)
+        std::copy(&other[0], &other[_size], &_x[0]);
       return *this;
     }
 
     /// Return pointer to data (const version)
-    const boost::shared_array<T> data() const
-    { return x; }
+    const T* data() const
+    { return _x; }
 
     /// Return pointer to data (non-const version)
-    boost::shared_array<T> data()
-    { return x; }
+    T* data()
+    { return _x; }
+
+   private:
+
+    /// Disable copy construction, to avoid unanticipated sharing or
+    /// copying of data. This means that an Array must always be passed as
+    /// reference, or as a (possibly shared) pointer.
+    Array(const Array& other) /* leave body undefined */;
 
   private:
 
-    // Length of array
+    /// Length of array
     uint _size;
 
-    // Array data
-    boost::shared_array<T> x;
+    /// Array data
+    T* _x;
 
+    /// True if instance is owner of data
+    bool _owner;
   };
 
   template <typename T>
@@ -224,11 +178,10 @@
   template <>
   inline void Array<double>::zero_eps(double eps)
   {
-    dolfin_assert(x);
     for (uint i = 0; i < _size; ++i)
     {
-      if (std::abs(x[i]) < eps)
-        x[i] = 0.0;
+      if (std::abs(_x[i]) < eps)
+        _x[i] = 0.0;
     }
   }
 
 
=== modified file 'dolfin/function/Function.cpp'
--- dolfin/function/Function.cpp	2012-03-01 12:01:06 +0000
+++ dolfin/function/Function.cpp	2012-03-07 13:15:44 +0000
@@ -315,7 +315,7 @@
   const Mesh& mesh = *_function_space->mesh();
 
   // Find the cell that contains x
-  const double* _x = x.data().get();
+  const double* _x = x.data();
   const Point point(mesh.geometry().dim(), _x);
   int id = mesh.intersected_cell(point);
 
@@ -440,7 +440,7 @@
   dolfin_assert(_function_space->mesh());
   const Mesh& mesh = *_function_space->mesh();
 
-  const double* _x = x.data().get();
+  const double* _x = x.data();
   const uint dim = mesh.geometry().dim();
   const Point point(dim, _x);
 
 
=== modified file 'dolfin/graph/ZoltanInterface.cpp'
--- dolfin/graph/ZoltanInterface.cpp	2011-11-18 12:14:50 +0000
+++ dolfin/graph/ZoltanInterface.cpp	2012-03-07 13:15:44 +0000
@@ -72,7 +72,7 @@
   // Call Zoltan function to compute coloring
   int num_id = 1;
   int rc = zoltan.Color(num_id, graph.size(), &global_ids[0],
-                        reinterpret_cast<int*>(colors.data().get()));
+                        reinterpret_cast<int*>(colors.data()));
   if (rc != ZOLTAN_OK)
   {
     dolfin_error("ZoltanInterface.cpp",
 
=== modified file 'dolfin/io/BinaryFile.cpp'
--- dolfin/io/BinaryFile.cpp	2012-03-01 12:01:06 +0000
+++ dolfin/io/BinaryFile.cpp	2012-03-07 13:15:44 +0000
@@ -65,7 +65,7 @@
 
   const uint n = read_uint();
   Array<double> values(n);
-  read_array(n, values.data().get());
+  read_array(n, values.data());
 
   vector.resize(n);
   vector.set_local(values);
@@ -148,7 +148,7 @@
   Array<double> values(n);
   vector.get_local(values);
   write_uint(n);
-  write_array(n, values.data().get());
+  write_array(n, values.data());
 
   close_write();
 }
 
=== modified file 'dolfin/io/XMLFunctionData.cpp'
--- dolfin/io/XMLFunctionData.cpp	2011-11-18 12:14:50 +0000
+++ dolfin/io/XMLFunctionData.cpp	2012-03-07 13:15:44 +0000
@@ -51,8 +51,8 @@
   const FunctionSpace& V = *u.function_space();
 
   std::vector<std::pair<uint, uint> > global_to_cell_dof;
-  Array<double> x;
-  Array<uint> indices;
+  std::vector<double> x;
+  std::vector<uint> indices;
 
   const uint num_dofs = V.dim();
 
@@ -115,7 +115,7 @@
       indices[i] = new_index;
     }
 
-    vector.set(x.data().get(), x.size(), indices.data().get());
+    vector.set(&x[0], x.size(), &indices[0]);
   }
 
   // Finalise vector
 
=== modified file 'dolfin/io/XMLVector.cpp'
--- dolfin/io/XMLVector.cpp	2011-11-18 12:14:50 +0000
+++ dolfin/io/XMLVector.cpp	2012-03-07 13:15:44 +0000
@@ -44,7 +44,7 @@
   read(data, indices, xml_dolfin);
 
   // Set data (GenericVector::apply will be called by calling function)
-  x.set(data.data().get(), data.size(), indices.data().get());
+  x.set(data.data(), data.size(), indices.data());
 }
 //-----------------------------------------------------------------------------
 void XMLVector::read(Array<double>& x, Array<uint>& indices,
 
=== modified file 'dolfin/la/EpetraVector.cpp'
--- dolfin/la/EpetraVector.cpp	2012-03-01 12:01:06 +0000
+++ dolfin/la/EpetraVector.cpp	2012-03-07 13:15:44 +0000
@@ -326,7 +326,7 @@
 
   values.resize(x->MyLength());
 
-  const int err = x->ExtractCopy(values.data().get(), 0);
+  const int err = x->ExtractCopy(values.data(), 0);
   if (err!= 0)
   {
     dolfin_error("EpetraVector.cpp",
@@ -470,7 +470,7 @@
   Epetra_SerialComm serial_comm = f.get_serial_comm();
 
   // Create map for y
-  const int* _indices = reinterpret_cast<const int*>(indices.data().get());
+  const int* _indices = reinterpret_cast<const int*>(indices.data());
   Epetra_BlockMap target_map(indices.size(), indices.size(), _indices, 1, 0, serial_comm);
 
   // Reset vector y
 
=== modified file 'dolfin/la/MUMPSLUSolver.cpp'
--- dolfin/la/MUMPSLUSolver.cpp	2011-11-02 16:03:33 +0000
+++ dolfin/la/MUMPSLUSolver.cpp	2012-03-07 13:15:44 +0000
@@ -159,7 +159,7 @@
   // Gather RHS on root process and attach
   Array<double> _b;
   b.gather_on_zero(_b);
-  data.rhs = _b.data().get();
+  data.rhs = _b.data();
 
   // Scaling strategy (77 is default)
   data.ICNTL(8) = 77;
@@ -171,8 +171,8 @@
 
   // Attach solution data to MUMPS object
   data.lsol_loc = local_x_size;
-  data.sol_loc  = x_local_vals.data().get();
-  data.isol_loc = reinterpret_cast<int*>(x_local_indices.data().get());
+  data.sol_loc  = x_local_vals.data();
+  data.isol_loc = reinterpret_cast<int*>(x_local_indices.data());
 
   // Solve problem
   data.job = 3;
@@ -185,8 +185,8 @@
     x_local_indices[i]--;
 
   // Set x values
-  x.set(x_local_vals.data().get(), x_local_indices.size(),
-        x_local_indices.data().get());
+  x.set(x_local_vals.data(), x_local_indices.size(),
+        x_local_indices.data());
   x.apply("insert");
 
   // Clean up
 
=== modified file 'dolfin/la/PETScVector.cpp'
--- dolfin/la/PETScVector.cpp	2012-03-01 12:01:06 +0000
+++ dolfin/la/PETScVector.cpp	2012-03-07 13:15:44 +0000
@@ -228,7 +228,7 @@
   for (uint i = 0; i < local_size; ++i)
     rows[i] = i + n0;
 
-  VecGetValues(*x, local_size, &rows[0], values.data().get());
+  VecGetValues(*x, local_size, &rows[0], values.data());
 }
 //-----------------------------------------------------------------------------
 void PETScVector::set_local(const Array<double>& values)
@@ -251,7 +251,7 @@
   for (uint i = 0; i < local_size; ++i)
     rows[i] = i + n0;
 
-  VecSetValues(*x, local_size, &rows[0], values.data().get(), INSERT_VALUES);
+  VecSetValues(*x, local_size, &rows[0], values.data(), INSERT_VALUES);
 }
 //-----------------------------------------------------------------------------
 void PETScVector::add_local(const Array<double>& values)
@@ -274,7 +274,7 @@
   for (uint i = 0; i < local_size; ++i)
     rows[i] = i + n0;
 
-  VecSetValues(*x, local_size, &rows[0], values.data().get(), ADD_VALUES);
+  VecSetValues(*x, local_size, &rows[0], values.data(), ADD_VALUES);
 }
 //-----------------------------------------------------------------------------
 void PETScVector::get_local(double* block, uint m, const uint* rows) const
@@ -681,7 +681,7 @@
   }
 
   // Prepare data for index sets (global indices)
-  const int* global_indices = reinterpret_cast<const int*>(indices.data().get());
+  const int* global_indices = reinterpret_cast<const int*>(indices.data());
 
   // Prepare data for index sets (local indices)
   const int n = indices.size();
 
=== modified file 'dolfin/mesh/MeshSmoothing.cpp'
--- dolfin/mesh/MeshSmoothing.cpp	2011-06-02 19:26:59 +0000
+++ dolfin/mesh/MeshSmoothing.cpp	2012-03-07 13:15:44 +0000
@@ -157,13 +157,12 @@
   BoundaryMesh boundary(mesh);
 
   const uint dim = mesh.geometry().dim();
-  Array<double> x;
 
   // Smooth boundary
   MeshGeometry& geometry = boundary.geometry();
   for (uint i = 0; i < boundary.num_vertices(); i++)
   {
-    x.update(dim, geometry.x(i));
+    Array<double> x(dim, geometry.x(i));
     sub_domain.snap(x);
   }
 
 
=== modified file 'dolfin/mesh/SubDomain.cpp'
--- dolfin/mesh/SubDomain.cpp	2012-03-06 13:55:20 +0000
+++ dolfin/mesh/SubDomain.cpp	2012-03-07 13:15:44 +0000
@@ -162,9 +162,6 @@
   // Extract exterior (non shared) facets markers
   const MeshFunction<bool>& exterior = mesh.parallel_data().exterior_facet();
 
-  // Array for vertex coordinate
-  Array<double> x;
-
   // Speed up the computation by only checking each vertex once (or twice if it
   // is on the boundary for some but not all facets).
   RangedIndexSet boundary_visited(mesh.num_vertices());
@@ -197,7 +194,7 @@
       {
         if (is_visited.insert(vertex->index()))
         {
-          x.update(_geometric_dimension, const_cast<double*>(vertex->x()));
+          Array<double> x(_geometric_dimension, const_cast<double*>(vertex->x()));
           is_inside[vertex->index()] = inside(x, on_boundary);
         }
 
@@ -212,8 +209,8 @@
     // Check midpoint (works also in the case when we have a single vertex)
     if (all_points_inside)
     {
-      x.update(_geometric_dimension,
-               const_cast<double*>(entity->midpoint().coordinates()));
+      Array<double> x(_geometric_dimension,
+                      const_cast<double*>(entity->midpoint().coordinates()));
       if (!inside(x, on_boundary))
         all_points_inside = false;
     }
 
=== modified file 'dolfin/swig/common/post.i'
--- dolfin/swig/common/post.i	2012-01-17 22:38:08 +0000
+++ dolfin/swig/common/post.i	2012-03-07 13:15:44 +0000
@@ -55,8 +55,6 @@
 // in any typemaps
 %feature("valuewrapper") dolfin::Array<TYPE>;
 
-%ignore dolfin::Array<TYPE>::Array(uint N, boost::shared_array<TYPE> x);
-
 // Cannot construct an Array from another Array.
 // Use NumPy Array instead
 %ignore dolfin::Array<TYPE>::Array(const Array& other);
@@ -72,7 +70,7 @@
   void __setitem__(unsigned int i, const TYPE& val) { (*self)[i] = val; }
 
   PyObject * array(){
-    return %make_numpy_array(1, TYPE_NAME)(self->size(), self->data().get(), true);
+    return %make_numpy_array(1, TYPE_NAME)(self->size(), self->data(), true);
   }
 
 }
@@ -83,7 +81,7 @@
 //-----------------------------------------------------------------------------
 CONST_ARRAY_IGNORES(double)
 ARRAY_EXTENSIONS(double, Double, double)
-ARRAY_EXTENSIONS(const double, ConstDouble, double)
+//ARRAY_EXTENSIONS(const double, ConstDouble, double)
 ARRAY_EXTENSIONS(unsigned int, UInt, uint)
 ARRAY_EXTENSIONS(int, Int, int)
 
 
=== modified file 'dolfin/swig/la/la_get_set_items.i'
--- dolfin/swig/la/la_get_set_items.i	2012-01-17 22:38:08 +0000
+++ dolfin/swig/la/la_get_set_items.i	2012-03-07 13:15:44 +0000
@@ -413,10 +413,10 @@
   dolfin::Array<double>* values = new dolfin::Array<double>(inds->size());
   if (row)
     // If returning a single row
-    self->get(values->data().get(), 1, &single, inds->size(), indices);
+    self->get(values->data(), 1, &single, inds->size(), indices);
   else
     // If returning a single column
-    self->get(values->data().get(), inds->size(), indices, 1, &single);
+    self->get(values->data(), inds->size(), indices, 1, &single);
 
   // Create the return vector and set the values
   boost::shared_ptr<dolfin::GenericVector> return_vec = self->factory().create_vector();
 
=== modified file 'dolfin/swig/typemaps/array.i'
--- dolfin/swig/typemaps/array.i	2012-01-30 22:44:17 +0000
+++ dolfin/swig/typemaps/array.i	2012-03-07 13:15:44 +0000
@@ -90,11 +90,11 @@
 // Director typemaps for dolfin::Array
 //-----------------------------------------------------------------------------
 %typemap(directorin) const dolfin::Array<double>& {
-  $input = %make_numpy_array(1, double)($1_name.size(), $1_name.data().get(), false);
+  $input = %make_numpy_array(1, double)($1_name.size(), $1_name.data(), false);
  }
 
 %typemap(directorin) dolfin::Array<double>& {
-  $input = %make_numpy_array(1, double)($1_name.size(), $1_name.data().get(), true);
+  $input = %make_numpy_array(1, double)($1_name.size(), $1_name.data(), true);
  }
 
 //-----------------------------------------------------------------------------
 
=== modified file 'site-packages/dolfin/compilemodules/compilemodule.py'
--- site-packages/dolfin/compilemodules/compilemodule.py	2012-01-30 21:20:53 +0000
+++ site-packages/dolfin/compilemodules/compilemodule.py	2012-03-07 13:15:44 +0000
@@ -42,6 +42,11 @@
            "expression_to_code_fragments",
            "math_header"]
 
+# Bump the interface version if anything changes that invalidates cached
+# modules (not required for change in generated code, swig version or dolfin
+# version)
+_interface_version = 0
+
 # A list of supported math builtins
 _math_builtins = [
     # cmath functions:
@@ -416,7 +421,8 @@
     if module_name is "":
         module_name = "dolfin_compile_code_%s" % hashlib.md5(repr(code) +\
                                             instant.get_swig_version() +\
-                                            dolfin.__version__).hexdigest()
+                                            dolfin.__version__ +\
+                                            str(_interface_version)).hexdigest()
     # Check the handed import files
     if dolfin_module_import:
         interface_files = []