[Scilab-users] vectorization difficulty

[Stéphane Mottelet]

Like this (I have simplified your script a little bit) for n=1000 there 
is a x60 speedup with Cardan's formulas:

clear

n  =  1000;
S11  =  rand(n,1);
S22  =  rand(n,1);
S33  =  rand(n,1);
S12  =  rand(n,1);
S23  =  rand(n,1);
S13  =  rand(n,1);
princ  =  zeros(n,3);

// with a loop

tic();
princ1  =  zeros(n,3);
for  i  =  1  :  n
     S=[S11(i)  S12(i)  S13(i)
        S12(i)  S22(i)  S23(i)
        S13(i)  S23(i)  S33(i)];
     princ1(i,:)  =  gsort(spec(S))';  
end

duration1  =  toc();  printf("Duration 1 = %g\n",duration1);

// using Cardan formulas

tic();
// characteristic polynomial is poly([d c b a],"x","coeff")
S13sq  =  S13.*S13;
S12sq  =  S12.*S12;
S23sq  =  S23.*S23;
//a=1; (not used) b=-S11-S22-S33;
c=S11.*S22+S11.*S33+S22.*S33-S23sq-S13sq-S12sq;
d=S11.*S23sq+S22.*S13sq+S12sq.*S33  -  S11.*S22.*S33-2*S13.*S12.*S23;

b2=b.*b;
p=-b2/3+c;
q=b.*(2*b2-9*c)/27+d;

//delta=-(4*p.*p.*p+27*q.*q) (not used since matrix is symetric =>real 
eigenvalues)

princ2=zeros(n,3);
theta=acos(1.5*q./p.*sqrt(-3./p))/3;
for  k=0:2
     princ2(:,k+1)=2*sqrt(-p/3).*cos(theta+2*k*%pi/3)-b/3;
end
princ2=gsort(princ2,"c")

duration2  =  toc();  printf("Duration 2 = %g\n",duration2);

disp(norm(princ1-princ2,%inf))
disp(duration1/duration2)

S.

Le 18/01/2019 à 15:21, Stéphane Mottelet a écrit :
> Hello Paul,
>
> If you stick to 3x3, you can vectorize the Cardan formulas applied to 
> the characteristic polynomial of each individual matrix.
>
> S.
>
>
>
> Le 15/01/2019 à 09:56, Carrico, Paul a écrit :
>>
>> clc
>>
>> mode(0)
>>
>> clear
>>
>> function*V*=_eigen_val_(*S11*, *S12*, *S13*, *S22*, *S23*, *S33*)
>>
>> S_u = [0 *S12* *S13* ; 0 0 *S23* ; 0. 0. 0.];
>>
>> S_d = [*S11* ; *S22* ; *S33*];
>>
>> S = S_u + S_u' + diag(S_d);
>>
>> *V* = gsort(spec(S),'lr','d')'; clear S; clear S_u; clear S_d;
>>
>> endfunction
>>
>> n= 10; m = 1;
>>
>> S11= rand(n,m);
>>
>> S22= rand(n,m);
>>
>> S33= rand(n,m);
>>
>> S12= rand(n,m);
>>
>> S23= rand(n,m);
>>
>> S13= rand(n,m);
>>
>> princ= zeros(n,3);
>>
>> /// with a loop/
>>
>> tic();
>>
>> princ1= zeros(n,3*m);
>>
>> fori = 1 : n
>>
>> S_u = [0 S12(i,m) S13(i,m) ; 0 0 S23(i,m) ; 0. 0. 0.];
>>
>> S_d = [S11(i,m) ; S22(i,1) ; S33(i,m)];
>>
>> S = S_u + S_u' + diag(S_d);
>>
>> princ1(i,:) = gsort(spec(S),'lr','d')'; clear S; clear S_u; clear S_d;
>>
>> end
>>
>> duration1= toc(); printf("Duration 1 = %g\n",duration1);
>>
>> /// using a function/
>>
>> tic();
>>
>> princ2= zeros(n,3*m);
>>
>> fori = 1 : n
>>
>> princ2(i,:) = 
>> _eigen_val_(S11(i,m),S12(i,m),S13(i,m),S22(i,m),S23(i,m),S33(i,m));
>>
>> end
>>
>> duration2= toc(); printf("Duration 2 = %g\n",duration2);
>>
>> isequal(princ1,princ2)
>>
>> /// using vectorization (in combination with the function ?)/
>>
>> i= (1:n)';
>>
>> S= zeros(3,3,n)
>
>
> -- 
> 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
>
> _______________________________________________
> users mailing list
> users at lists.scilab.org
> https://antispam.utc.fr/proxy/1/c3RlcGhhbmUubW90dGVsZXRAdXRjLmZy/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

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.scilab.org/pipermail/users/attachments/20190118/b6e0b040/attachment.htm>