Ok, thanks for the info (and @inbounds does improve it a bit). I usually
follow your advice and fuse the operations together when I need the speed,
but since I do all manner of combinations of vectorized operations
throughout my module I tend to prefer using .*=, ./=, etc unless I need it.
Is there a way to define function .*!=, ./!= and call it in the form A .*!=
B? If I try the following it gives me an error:
function .*!=(A,B)
@inbounds for ix in eachindex(A,B)
A[ix] = A[ix] * B[ix]
end
end
On Friday, December 18, 2015 at 10:20:31 AM UTC-8, Yichao Yu wrote:
On Fri, Dec 18, 2015 at 1:09 PM, Ethan Anderes <[email protected]
> <javascript:>> wrote:
> > Hi everyone. I have a pretty basic question. Are there inplace
> vectorized
> > operation for things like A .*= B, A ./= B when A and B are
> Array{Float64,
> > d} for general d? I’m ok with writing these myself (and I love the fact
> that
> > Julia allows me to do this and have it be fast) but I have the feeling
> like
> > it must be somewhere in Base and I’m missing it. Also, with regard to
> > writing my own version, I would like to avoid the tiny bit of additional
> > technical debt incurred by defining these myself (and the fact that A
> .*= B,
> > A ./= B, are just so easy to read).
>
> I usually just write out the loop myself. It has the addition benefit
> of being able to add other operations at the same time which I often
> end up doing (e.g. do A = A .* B * c + d instead)
> Ref https://github.com/JuliaLang/julia/issues/249
>
> P.S. you can add `@inbounds` to your `myownscale!` which might speed
> it up by a factor of a few (due to automatic simd)
>
> >
> > Note: I get an error if I do A[:] .*= B for Array{Float64, d}, d>1. I
> could
> > do A[:,...,:] .*= B or A[:] .*= B[:] but the former needs to work for
> > generic dimension and I’m worried the later has issues with regard to
> > looping, linear indexing and still creats temporaries. Also, scale!(A,B)
> > doesn’t work for me when both A and B are of the same dimension.
> >
> > This snippit gives the timings.
> >
> > function myownscale!(A,B)
> > for ix in eachindex(A,B)
> > A[ix] = A[ix] * B[ix]
> > end
> > end
> >
> > function test1!(A,B)
> > A[:] .*= B # <-- gives error
> > end
> >
> > function test2!(A,B)
> > A[:,:] .*= B
> > end
> >
> > function test3!(A,B)
> > A[:] .*= B[:]
> > end
> >
> > A, B = rand(1_000, 1_000), rand(1_000, 1_000);
> >
> > test1!(A,B); #<-- error
> >
> > # warmup
> > @time test2!(A,B);
> > @time test3!(A,B);
> > @time myownscale!(A,B);
> >
> > @time test2!(A,B);
> > @time test3!(A,B);
> > @time myownscale!(A,B);
> >
> > The last three commands result in:
> >
> > julia> @time test2!(A,B);
> > 0.011629 seconds (22 allocations: 15.260 MB, 34.11% gc time)
> >
> > julia> @time test3!(A,B);
> > 0.007802 seconds (27 allocations: 22.889 MB, 19.09% gc time)
> >
> > julia> @time myownscale!(A,B);
> > 0.001341 seconds (4 allocations: 160 bytes)
> >
> > So, am I missing something in Base or is idiomatic Julia just defining
> these
> > myself?
>
