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 : [email protected] > Unsubscribe : https://launchpad.net/~dolfin > More help : https://help.launchpad.net/ListHelp > _______________________________________________ Mailing list: https://launchpad.net/~dolfin Post to : [email protected] Unsubscribe : https://launchpad.net/~dolfin More help : https://help.launchpad.net/ListHelp
