On Monday, April 27, 2015 at 2:52:26 PM UTC-5, Alex wrote: > > On Monday, 27 April 2015 21:05:33 UTC+2, François Fayard wrote: > > If one implements y = y + delta_t * dy_dt in the code with a Vector, you > create a lot of heap arrays that has a huge cost (It seems that ODE.jl is > doing that by the way). So you are forced to loop over the vector elements, > which does not look really good. How about implementing a > add_linear_combination!(y, delta_t, dy_dt) which would make things easier > to extend for any kind of "Vector". This is exactly what people from odeint > have been doing and it looks like a nice idea. If you want to go crazy you > can event implement a MPI vector whose memory is scattered around different > nodes. Then you just implement the add_linear_combination for your array, > and you are done. > > The problem is that `y = y + delta_t * dy_dt` currently makes two copies, > as you have observed. There is an open discussion about in-place assignment > operators [1], which would help. There is also Devectorize.jl [2]. However, > I agree that having something like add!(y, x, beta) would be nice and I > think this also came up in some discussion. Funnily enough we just recently > discussed this strategy for ODE.jl [3] ...
There is http://julia.readthedocs.org/en/latest/stdlib/linalg/#Base.LinAlg.BLAS.axpy! if you're using BLAS-compatible array types.
