dolfin team mailing list archive

[Johan Hake <johan.hake@xxxxxxxxx>]

On Thursday 28 January 2010 00:30:06 Garth N. Wells wrote:
> Johan Hake wrote:
> > On Monday 25 January 2010 10:01:18 Garth N. Wells wrote:
> >> Johan Hake wrote:
> >>> On Monday 25 January 2010 03:10:46 Garth N. Wells wrote:
> >>>> Johan Hake wrote:
> >>>>> On Sunday 24 January 2010 16:39:14 Garth N. Wells wrote:
> >>>>>> Johan Hake wrote:
> >>>>>>> On Sunday 24 January 2010 16:12:37 Garth N. Wells wrote:
> >>>>>>>> Johan Hake wrote:
> >>>>>>>>> On Sunday 24 January 2010 00:03:41 Garth N. Wells wrote:
> >>>>>>>>>> Johan Hake wrote:
> >>>>>>>>>>> On Saturday 23 January 2010 14:55:08 Garth N. Wells wrote:
> >>>>>>>>>>>> Johan Hake wrote:
> >>>>>>>>>>>>> On Saturday 23 January 2010 08:42:14 Garth N. Wells wrote:
> >>>>>>>>>>>>>> Is it correct that behind the scenes that
> >>>>>>>>>>>>>>
> >>>>>>>>>>>>>>   U0 = Function(V)
> >>>>>>>>>>>>>>   U  = Function(V)
> >>>>>>>>>>>>>>   U0.vector()[:] = U.vector()[:]
> >>>>>>>>>>>>>>
> >>>>>>>>>>>>>> involves a GenericVector::get(..) call and a
> >>>>>>>>>>>>>> GenericVector::set(..) call? If so, it isn't ideal since it
> >>>>>>>>>>>>>> introduces unnecessary new/delete operations and unnecessary
> >>>>>>>>>>>>>> copying of data.
> >>>>>>>>>>>>>
> >>>>>>>>>>>>> None of GenericVector::get(..) or GenericVector::set(..) are
> >>>>>>>>>>>>> invoked, see __getslice__ and __setslice__ in la_post.i.
> >>>>>>>>>>>>>
> >>>>>>>>>>>>>    U0.vector()[:]
> >>>>>>>>>>>>>
> >>>>>>>>>>>>> involves
> >>>>>>>>>>>>>
> >>>>>>>>>>>>>   GenericVector::operator =(..)
> >>>>>>>>>>>>>
> >>>>>>>>>>>>> and
> >>>>>>>>>>>>>
> >>>>>>>>>>>>>   U.vector()[:]
> >>>>>>>>>>>>>
> >>>>>>>>>>>>> involves
> >>>>>>>>>>>>>
> >>>>>>>>>>>>>   GenericVector::copy()
> >>>>>>>>>>>>>
> >>>>>>>>>>>>> However the latter is unnecessary as you instead can do:
> >>>>>>>>>>>>>
> >>>>>>>>>>>>>   U0.vector()[:] = U.vector()
> >>>>>>>>>>>>>
> >>>>>>>>>>>>> invoking the assignment operator of U0's vector with U's
> >>>>>>>>>>>>> vector.
> >>>>>>>>>>>>
> >>>>>>>>>>>> What happens if I do
> >>>>>>>>>>>>
> >>>>>>>>>>>>   x = U.vector()[:]
> >>>>>>>>>>>
> >>>>>>>>>>> It just triggers the copy method of GenericVector, which is the
> >>>>>>>>>>> same behavior as for other itterable Python types.
> >>>>>>>>>>>
> >>>>>>>>>>>> ? Is x a numpy array?
> >>>>>>>>>>>
> >>>>>>>>>>> No you need to call array() to accomplish that.
> >>>>>>>>>>
> >>>>>>>>>> OK. What I'm trying to do is
> >>>>>>>>>>
> >>>>>>>>>>     # Get vectors
> >>>>>>>>>>     u_vec  = u.vector()[:]
> >>>>>>>>>>     u0_vec = u0.vector()[:]
> >>>>>>>>>>     v0_vec = v0.vector()[:]
> >>>>>>>>>>     a0_vec = a0.vector()[:]
> >>>>>>>>>
> >>>>>>>>> You should not need to make a copy of the vectors here.
> >>>>>>>>
> >>>>>>>> How can I avoid it?
> >>>>>>>
> >>>>>>> a_vec and v_vec are new vectors. None of the four vectors below get
> >>>>>>> modified by the a_vec and v_vec expressions so no need of copying,
> >>>>>>> and the v0 and a0 assignment should work with GenericVectors too.
> >>>>>>
> >>>>>> Do you mean that just
> >>>>>>
> >>>>>>   a_vec = 1.0/(2.0*beta)*((u - u0 - v0*dt)/(0.5*dt*dt) \
> >>>>>>            - (1.0-2.0*beta)*a0 )
> >>>>>>
> >>>>>> where u and u0 are GenericVectors should work?
> >>>>>
> >>>>> Have you tried?
> >>>>
> >>>> Can I somehow get a 'reference' to the vector so I don't have to use
> >>>> u.vector()[:] in the expressions?
> >>>
> >>> That is what u.vector() gives you.
> >>
> >> I figured it out eventually. There's nothing quite like a '&' symbol . .
> >> . . .
> >
> > Good!
> >
> > In Python every name in a namespace is a reference(pointer) to an
> > instance. Johan
> >
> >> Garth
> >>
> >>> Johan
> >>>
> >>>> Garth
> >>>>
> >>>>> As long as the rest (besides a0, which I assume also is a
> >>>>> GenericVector) are scalars everything should just work. The Python LA
> >>>>> interface (at least for GenericVector) should work more or less as
> >>>>> the NumPy interface which I think is nice :)
> 
> I have things working nicely now when working with GenericVector. I have:
> 
>   u_vec  = u.vector()
>   u0_vec = u0.vector()
>   v0_vec = v0.vector()
>   a0_vec = a0.vector()
> 
>   a_vec  = (1.0/(2.0*beta))*( (u_vec - u0_vec - v0_vec*dt)/(0.5*dt*dt) \
>           - (1.0-2.0*beta)*a0_vec )
>   v_vec = dt*((1.0-gamma)*a0_vec + gamma*a_vec) + v0_vec

Good!

> where beta, gamma and dt are floats. I'd like now to know a little more
> about what's going on behind the scenes, especially with respect to the
> creation of temporaries. Do all the operations just call the relevant
> functions from the GenericVector interface? 

Yes, eventually we end up calling either operator *=() or operator +=().

> Will each operation create a temporary?

All binary operators creates a temporary, also the scalar multiplier. I think 
this would be the same as in NumPy. The way we do it is to create a copy of 
one of the vectors and then we apply the unary operator on that one, and 
return the result.

So the last expression would create:

  1 : (1.0-gamma)*a0_vec
  2 : gamma*a_vec
  3 : (1.0-gamma)*a0_vec + gamma*a_vec
  4 : dt*((1.0-gamma)*a0_vec + gamma*a_vec)
  5 : dt*((1.0-gamma)*a0_vec + gamma*a_vec) + v0_vec

temporaries (if I am not totally off). You should be able to speed this up by 
using axpy.

Johan


> Garth
> 
> >>>>> We cannot take 1./v, where v is a GenericVector.
> >>>>>
> >>>>>>> Do you get any error messages?
> >>>>>>
> >>>>>> No errors. What I have now seems to work fine.
> >>>>>
> >>>>> Ok, and that is because what you do obviously works for NumPy arrays.
> >>>>>
> >>>>> Johan
> >>>>>
> >>>>>> Garth
> 
> _______________________________________________
> Mailing list: https://launchpad.net/~dolfin
> Post to     : dolfin@xxxxxxxxxxxxxxxxxxx
> Unsubscribe : https://launchpad.net/~dolfin
> More help   : https://help.launchpad.net/ListHelp
>