[Scilab-users] {EXT} need a more efficient and faster code: suggestions welcome

[Stéphane Mottelet]

Hi all,

I hope that the following will stop the (useless) competition : if speed 
is a real bottleneck for this particular computation, using plain C is 
the only option. The self-contained script (n=20000) given at the end of 
this message gives the following timings on various hardware (only one 
processor core is used) :

2,8 GHz Quad-Core Intel Xeon : 1.08 s (OSX)
2.2 GHz Xeon(R) CPU E5-2660 v2 @  : 1.84 s (Linux)

Don't use plain Scilab if speed is a crucial issue....

//  start of code
source=['#include <math.h>'
'#define SQR(x) ((x)*(x))'
'#define MIN(x, y) (((x) < (y)) ? (x) : (y))'
'void mindist(double d[],double x[],int *n) {'
'  int i,j,k;'
'  double dist;'
'  for (i=0;i<*n;i++) { '
'    d[i]=INFINITY;'
'    for (j=0;j<*n;j++) if (i!=j) { '
'      dist=0;'
'      for (k=0;k<3;k++) dist+=SQR(x[j*3+k]-x[i*3+k]);'
'      d[i]=MIN(d[i],dist);'
'    }'
'  }'
'}']
mputl(source,'mindist.c')
ilib_for_link('mindist','mindist.c',[],"c");
exec loader.sce

n=20000;
r=23;
radius = r*grand(n,1,'def').^(1/3);
phi = 2*%pi*grand(n,1,'def');
costheta = 1 - 2*grand(n,1,'def');
radsintheta = radius.*sin(acos(costheta));

X=[radsintheta.*cos(phi) radsintheta.*sin(phi) radius.*costheta];

n=size(X,1);
tic;
dist=sqrt(call("mindist",X',2,"d",n,3,"i","out",[1,n],1,"d"));
disp(toc())
// end of code

Le 31/01/2018 à 23:31, Rafael Guerra a écrit :
> Hi Heinz,
>
> Find herein a vectorised memory hog solution for Scilab 6, which is fast (~30 s for 20K points on Win7 laptop):
>
> // START OF CODE (Scilab 6)
> Clear
> n=20000;
> r=23;
> radius = r*grand(n,1,'def').^(1/3);
> phi = 2*%pi*grand(n,1, 'def');
> costheta = 1 - 2*grand(n,1, 'def');
> radsintheta = radius.*sin(acos(costheta));
> X = [radsintheta.*cos(phi), radsintheta.*sin(phi), radius.*costheta];
> tic;
> X2 = sum(X.*X, 'c');
> X3 = X2.*.ones(n,1)';
> X3 = X3 + X3';
> D = abs(X3 - 2*(X*X'))';
> D(1:n+1:$) = [];  // remove diagonal 0's
> D = matrix(D,n-1,n); // reshape D to (n-1) x n size
> MinDist1 = sqrt(min(D,'r'));
> t1=toc();
> disp(t1)
> jj=-0.05:0.1:3.05;
> H1 = histc(jj,MinDist1);
> disp([0.05+jj(1:$-1)' H1']);
> // END OF CODE
>
> Regards,
> Rafael
>
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> http://lists.scilab.org/mailman/listinfo/users


-- 
Stéphane Mottelet
Ingénieur de recherche
EA 4297 Transformations Intégrées de la Matière Renouvelable
Département Génie des Procédés Industriels
Sorbonne Universités - Université de Technologie de Compiègne
CS 60319, 60203 Compiègne cedex
Tel : +33(0)344234688
http://www.utc.fr/~mottelet