The usage of construct like zero(zero(T)*zero(T)) is not really easy to
understand. I assume you propose this due to the upcasting arithmetic
integer behavior of non-native integer types (int16, ...), right?
Am Donnerstag, 24. April 2014 11:46:37 UTC+2 schrieb Alex:
>
> 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
>>>
>>>>
>>>>
