I couldn't find the definition of sum right now, but I guess it is calling
BLAS or something like that, which might explain the difference?
Consider for example
function trapz3{T<:Number}(x::Vector{T}, y::Vector{T})
local n = length(x)
if length(y) != n
error("Vectors 'x', 'y' must be of same length.")
end
if n == 1; return 0.0; end
int = (x[2:end]-x[1:end-1]) .* (y[2:end]+y[1:end-1])
r = 0.0
for i in 1:length(int)
r += int[i]
end
r / 2.0, BLAS.asum(int)/2.0, sum(int)/2.0
end
which gives me
julia> trapz3(x, y)
(1.9999999999984162,1.9999999999983549,1.999999999998355)
Best,
Alex.
PS: Just a minor remark: one should probably use zero(T) instead of 0.0 (or
even zero(zero(T)*zero(T))).
On Thursday, 24 April 2014 10:28:20 UTC+2, Hans W Borchers wrote:
>
> Just being curious, I compared these two version on accuracy and execution
> time, having grown up in Matlab as well. Let's call the first, vectorized
> version by Tomas as 'trapz1', the second version by Cameron as 'trapz2',
> see below. Then I get the following results:
>
> julia> x = linspace(0, pi, 1000000);
>
> julia> y = sin(x);
>
> julia> trapz1(x, y)
> 1.999999999998355
>
> julia> trapz2(x, y)
> 1.9999999999984162
>
> julia> @time for i in 1:100; trapz1(x, y); end
> elapsed time: 4.858258575 seconds (5600105600 bytes allocated)
>
> julia> @time for i in 1:100; trapz2(x, y); end
> elapsed time: 0.41866061 seconds (1600 bytes allocated)
>
>
> and in general the unvectorized version is more than 10 times faster. But
> why are the two results different (though only by 6.1e-14)? It looks like
> they are performing exactly the same computations. Does the sequence of
> elementary operations play a role?
> Is it to be expected that a version with a localized length is faster? I
> didn't see significantly different running times.
>
>
> PS: The function definitions I used here are:
>
> function trapz1{T<:Number}(x::Vector{T}, y::Vector{T})
> local n = length(x)
> if length(y) != n
> error("Vectors 'x', 'y' must be of same length.")
> end
> if n == 1; return 0.0; end
> sum((x[2:end]-x[1:end-1]) .* (y[2:end]+y[1:end-1])) / 2.0
> end
>
>
> function trapz2{T<:Number}(x::Vector{T}, y::Vector{T})
> local n = length(x)
> if (length(y) != n)
> error("Vectors 'x', 'y' must be of same length")
> end
> if n == 1; return 0.0; end
> r = 0.0
> for i in 2:n
> r += (x[i] - x[i-1]) * (y[i] + y[i-1])
> end
> r / 2.0
> end
>
>
>
> On Wednesday, April 23, 2014 11:20:50 PM UTC+2, Tomas Lycken wrote:
>>
>>
>> On Wednesday, April 23, 2014 11:10:15 PM UTC+2, Cameron McBride wrote:
>>>
>>> Or you can use the non-vectorized version and save the overhead of the
>>> temporary arrays being created by the addition and multiplication steps.
>>>
>>
>> There's really no way I can hide that I learnt scientific computing in
>> Matlab, is there? :P
>>
>>>
>>>
