[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: for i=1:bignum, x(a(i)) += y(i); end
From: |
Miroslaw Kwasniak |
Subject: |
Re: for i=1:bignum, x(a(i)) += y(i); end |
Date: |
Thu, 19 May 2005 02:17:35 +0200 |
User-agent: |
Mutt/1.5.9i |
On Wed, May 18, 2005 at 02:39:14PM -0400, Etienne Grossmann wrote:
>
>
> Hi All,
>
> has anyone written an oct-file for this?
>
> for i=1:bignum, x(a(i)) += y(i); end
Not oct, but m-function with about 15x speedup (and 2x slowup for small
bignum:)
[xa,a_unique] = sum_by_index(y,a);
x(a_unique)+=xa;
Timing results (t0-loop, t1-sum_by_index):
octaveX:29> bignum=100;n_classes=10; sum_x
t0 = 0.012109 t1 = 0.021936 Delta_max = 0
octaveX:30> bignum=100;n_classes=100; sum_x
t0 = 0.012310 t1 = 0.022559 Delta_max = 0
octaveX:31> bignum=10000;n_classes=100; sum_x
t0 = 0.75447 t1 = 0.059115 Delta_max = 0
octaveX:32> bignum=10000;n_classes=10; sum_x
t0 = 0.75039 t1 = 0.058044 Delta_max = 0
octaveX:33> bignum=10000;n_classes=1000; sum_x
t0 = 0.74890 t1 = 0.064619 Delta_max = 0
octaveX:34> bignum=100000;n_classes=100; sum_x
t0 = 7.4814 t1 = 0.47814 Delta_max = 0
Mirek
%%%%%%%%%%%%%%%%
function [y,a] = sum_by_index(x,a)
[ a i ] = sort(a);
x = x(i);
% Find location of unique values of a
i = [ -2*abs(a(1))-1, a(:)' ];
i = find(diff(i));
a = a(i);
lx = length(x);
li = length(i);
% row indices
r = ones(1,lx);
r(i(2:end)) = 1-diff(i);
r = cumsum(r);
% column indices
c = zeros(1,lx);
c(i)= 1:li;
c(i(2:end)) -= 1:li-1;
c = cumsum(c);
% Place x columnwise & sum
R = max(r);
C = length(a);
y = zeros(1,R*C);
y(sub2ind([R C],r,c)) = x;
y = sum(reshape(y,R,C));
endfunction
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Test script
%bignum=100000
%n_classes=500
x=rand(1,bignum);
a=round(n_classes*x)+1;
x=a+(1:length(a))/100;
y1=y=zeros(1,max(a));
tic
for i=1:length(a), y(a(i)) += x(i); end
t0=toc
tic
[ya,a_unique] = sum_by_index(x,a);
y1(a_unique)+=ya;
t1=toc
Delta_max=max(abs(y-y1))
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-------------------------------------------------------------
Octave is freely available under the terms of the GNU GPL.
Octave's home on the web: http://www.octave.org
How to fund new projects: http://www.octave.org/funding.html
Subscription information: http://www.octave.org/archive.html
-------------------------------------------------------------