[Scilab-users] vectorization difficulty

classic Classic list List threaded Threaded
6 messages Options
Carrico, Paul Carrico, Paul
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] vectorization difficulty

Dear All

 

May I ask you advices in order to use vectorization instead of the loop? All the  trials I did have failed – the problem I’m confront to comes from the S matrix

 

Thanks for your feedback

 

Paul

 

########################################################

 

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);

for i = 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);

for i = 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)

 


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: vectorization difficulty

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);

for i = 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);

for i = 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
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
mottelet mottelet
Reply | Threaded
Open this post in threaded view
|

Re: vectorization difficulty

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);

for i = 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);

for i = 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
[hidden email]
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

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
fmiyara fmiyara
Reply | Threaded
Open this post in threaded view
|

[Scilab-users] compact layout for console


Dear all,

Is there any way to get a more compact presentation of results in the
console?

I mean to avoid the unnecessary extra blank lines between output lines.

Matlab allows so from the preferences, but I couldn't find a similar
option in Scilab. Maybe there is some undocumented feature somehere...?

Regards,

Federico Miyara

---
El software de antivirus Avast ha analizado este correo electrónico en busca de virus.
https://www.avast.com/antivirus

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
Samuel GOUGEON Samuel GOUGEON
Reply | Threaded
Open this post in threaded view
|

Re: compact layout for console

Hello Federico,

The most compact you can get is after having used
--> mode(1)

https://help.scilab.org/docs/6.0.1/en_US/mode.html

Regards
Samuel

Le 20/01/2019 à 09:46, Federico Miyara a écrit :

>
> Dear all,
>
> Is there any way to get a more compact presentation of results in the
> console?
>
> I mean to avoid the unnecessary extra blank lines between output lines.
>
> Matlab allows so from the preferences, but I couldn't find a similar
> option in Scilab. Maybe there is some undocumented feature somehere...?
>
> Regards,
>
> Federico Miyara
>
> ---
> El software de antivirus Avast ha analizado este correo electrónico en
> busca de virus.
> https://www.avast.com/antivirus
>
> _______________________________________________
> users mailing list
> [hidden email]
> http://lists.scilab.org/mailman/listinfo/users
>

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users
fmiyara fmiyara
Reply | Threaded
Open this post in threaded view
|

Re: compact layout for console

Samuel,

Thank you!

Regards,

Federico


On 20/01/2019 07:09, Samuel Gougeon wrote:
Hello Federico,

The most compact you can get is after having used
--> mode(1)

https://help.scilab.org/docs/6.0.1/en_US/mode.html

Regards
Samuel

Le 20/01/2019 à 09:46, Federico Miyara a écrit :

Dear all,

Is there any way to get a more compact presentation of results in the console?

I mean to avoid the unnecessary extra blank lines between output lines.

Matlab allows so from the preferences, but I couldn't find a similar option in Scilab. Maybe there is some undocumented feature somehere...?

Regards,

Federico Miyara

---
El software de antivirus Avast ha analizado este correo electrónico en busca de virus.
https://www.avast.com/antivirus

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users


_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users




Libre de virus. www.avast.com

_______________________________________________
users mailing list
[hidden email]
http://lists.scilab.org/mailman/listinfo/users