[scilab-Users] code which does not work
Ivan Maximov
vaborg at gmail.com
Wed Mar 18 17:52:47 CET 2009
Hi Francois,
> - how you run it
by ctrl+l in editor
> - in what Scilab version
4.2.1 it gives error message
5.1 it kills Scilab at all
> - on what platform
Fedore 10, x86
> - where you set the breakpoint and how you do it (from the debugger or
> directly using setbpt?)
from debugger
breakpoint in string 7 function krotov_csa.sci
below is code
__________________________________________________________________________________________________________
main script
__________________________________________________________________________________________________________
// INEPT and CSA
clearglobal;
funcprot(0);
global path dbg
path='/home/vaborg/scilab/csa/';
cd(path);
// Debug info
dbg=1;
exec(path+'system.sce');
exec(path+'H1.sci');
exec(path+'penalty.sci');
exec(path+'krotov_csa.sci');
disp('START GENERAL');
tic();
// ITERATION LOOP
w1=rf*rand(1,N); w2=rf*rand(1,N);
w3=rf*rand(1,N); w4=rf*rand(1,N);
[eff,r1,r2,r3,r4]=krotov_csa(w1,w2,w3,w4);
ti=toc();
disp('STOP GENERAL');
printf("Time is: %i\n",ti);
__________________________________________________________________________________________________________
system.sce
__________________________________________________________________________________________________________
global dt N dp csI csS csIN csSN l1 csN Hij
err=exec(path+'pauli.sce','errcatch');
if err~=0,
disp('ERROR: Can not find pauli.sce!');
abort();
end
ix=x.*.u; iy=y.*.u; iz=z.*.u;
sx=u.*.x; sy=u.*.y; sz=u.*.z;
uu=u.*.u;
dp=2*%pi;
init=sx;
tini=ix;
//------------------------------
// couplings time and shifts
//------------------------------
J=140;
T=1/J;
N=200; NIT=20;
l1=1e-5;
delta=1e-2; eta=delta;
dt=T/N; rf=1e3;
// shift arrays
csI=6:0.1:9;
csS=115:2.5:125;
wi=-500e6*dp;
ws=50.684e6*dp;
csIN=length(csI);
csSN=length(csS);
csN=csIN*csSN;
// Hamiltonians
Hamilt=dp*J*iz*sz;
Hij(:,:,csN)=uu;
U(:,:,csN)=uu;
B(:,:,csN)=uu;
for i=1:csIN,
for j=1:csSN,
Hij(:,:,(i-1)*csSN+j)=(-csI(i)*1e-6)*wi*iz+(-csS(j)*1e-6)*ws*sz;
U(:,:,(i-1)*csSN+j)=uu;
B(:,:,(i-1)*csSN+j)=uu;
end
end
__________________________________________________________________________________________________________
krotov_csa.sci
__________________________________________________________________________________________________________
function [ef,u1,u2,u3,u4]=krotov_csa(x1,x2,x3,x4)
u1=x1; u2=x2; u3=x3; u4=x4;
t1=x1; t2=x2; t3=x3; t4=x4;
PROP=uu;
TE=uu;
pause;
for i=1:N,
for j=1:csN,
PROP=expm(-%i*dt*(Hij(:,:,j)+Hamilt+H1(u1(i),u2(i),u3(i),u4(i))));
TE=PROP*U(:,:,j);
U(:,:,j)=TE;
end
end
eff=0; pen=0; fu=0;
for i=1:csN,
eff=eff+real(trace(tini*U(:,:,i)*init*U(:,:,i)'))/csN;
end
pen=penalty(u1,u2,u3,u4);
fu=1+eff-pen;
if dbg==1,
printf("EFF: %f \tF: %f \tP: %f\n",eff,fu,pen);
end
for i=1:csN,
B(:,:,i)=U(:,:,i)'+init*U(:,:,i)'*tini;
end
for i=N:(-1):1,
for j=1:csN,
PROP=expm(-%i*dt*(Hij(:,:,j)+Hamilt+H1(t1(i),t2(i),t3(i),t4(i))));
B(:,:,j)=B(:,:,j)*PROP;
end
end
//-----------------------------------------------
// ITERATION
//-----------------------------------------------
for iter=1:NIT,
for i=1:csN,
U(:,:,i)=uu;
end
PROP=zeros(4,4);
for i=1:csN,
PROP=PROP+B(:,:,i);
end
u1(1)=(1-delta)*t1(1)+delta/l1*imag(trace(ix*PROP));
u2(1)=(1-delta)*t2(1)+delta/l1*imag(trace(iy*PROP));
u3(1)=(1-delta)*t3(1)+delta/l1*imag(trace(sx*PROP));
u4(1)=(1-delta)*t4(1)+delta/l1*imag(trace(sy*PROP));
for i=2:N,
for j=1:csN,
PROP=expm(-%i*dt*(Hij(:,:,j)+Hamilt+H1(u1(i-1),u2(i-1),u3(i-1),u4
(i-1)))); U(:,:,j)=PROP*U(:,:,j);
PROP=expm(%i*dt*(Hij(:,:,j)+Hamilt+H1(t1(i-1),t2(i-1),t3(i-1),t4
(i-1)))); B(:,:,j)=B(:,:,j)*PROP;
end
PROP=zeros(4,4);
for j=1:csN,
PROP=PROP+U(:,:,j)*B(:,:,j);
end
u1(i)=(1-delta)*t1(i)+delta/l1*imag(trace(ix*PROP));
u2(i)=(1-delta)*t2(i)+delta/l1*imag(trace(iy*PROP));
u3(i)=(1-delta)*t3(i)+delta/l1*imag(trace(sx*PROP));
u4(i)=(1-delta)*t4(i)+delta/l1*imag(trace(sy*PROP));
end
for j=1:csN,
PROP=expm(-%i*dt*(Hij(:,:,j)+Hamilt+H1(u1(N),u2(N),u3(N),u4(N))));
U(:,:,j)=PROP*U(:,:,j);
end
temp=0;
for j=1:csN,
temp=temp+real(trace(tini*U(:,:,j)*init*U(:,:,j)'))/csN;
end
eff=cat(1,eff,temp);
pen=cat(1,pen,penalty(u1,u2,u3,u4));
fu=cat(1,fu,1+eff(iter+1)-pen(iter+1));
if dbg==1,
printf("EFF: %f \tF: %f \tP: %f\n",eff(iter+1),fu(iter+1),pen(iter+1));
end
if (fu(iter+1)-fu(iter))<0,
disp('ERROR: DISCONVERGENCY');
break;
elseif (fu(iter+1)-fu(iter))>0 & (fu(iter+1)-fu(iter))<1e-4,
disp('CONVERGENT!');
break;
end
for j=1:csN,
B(:,:,j)=U(:,:,j)'+init*U(:,:,j)'*tini;
end
PROP=zeros(4,4);
for j=1:csN,
PROP=PROP+U(:,:,j)*B(:,:,j);
end
t1(N)=(1-eta)*u1(N)+eta/l1*imag(trace(ix*PROP));
t2(N)=(1-eta)*u2(N)+eta/l1*imag(trace(iy*PROP));
t3(N)=(1-eta)*u3(N)+eta/l1*imag(trace(sx*PROP));
t4(N)=(1-eta)*u4(N)+eta/l1*imag(trace(sy*PROP));
for i=(N-1):(-1):1,
for j=1:csN,
PROP=expm(%i*dt*(Hij(:,:,j)+Hamilt+H1(u1(i+1),u2(i+1),u3(i+1),u4(i
+1)))); U(:,:,j)=PROP*U(:,:,j);
PROP=expm(-%i*dt*(Hij(:,:,j)+Hamilt+H1(t1(i+1),t2(i+1),t3(i+1),t4(i
+1)))); B(:,:,j)=B(:,:,j)*PROP;
end
PROP=zeros(4,4);
for j=1:csN,
PROP=PROP+U(:,:,j)*B(:,:,j);
end
t1(i)=(1-eta)*u1(i)+eta/l1*imag(trace(ix*PROP));
t2(i)=(1-eta)*u2(i)+eta/l1*imag(trace(iy*PROP));
t3(i)=(1-eta)*u3(i)+eta/l1*imag(trace(sx*PROP));
t4(i)=(1-eta)*u4(i)+eta/l1*imag(trace(sy*PROP));
end
for j=1:csN,
PROP=expm(%i*dt*(Hij(:,:,j)+Hamilt+H1(u1(1),u2(1),u3(1),u4(1))));
U(:,:,j)=PROP*U(:,:,j);
PROP=expm(-%i*dt*(Hij(:,:,j)+Hamilt+H1(t1(1),t2(1),t3(1),t4(1))));
B(:,:,j)=B(:,:,j)*PROP;
end
// disp(U);
end
ef=eff(length(eff));
endfunction
__________________________________________________________________________________________________________
H1.sci
__________________________________________________________________________________________________________
function [res]=H1(x1,x2,x3,x4)
res=x1*ix+x2*iy+x3*sx+x4*sy;
endfunction
__________________________________________________________________________________________________________
penalty.sci
__________________________________________________________________________________________________________
function [res]=penalty(x1,x2,x3,x4)
res=dt*l1*sum(x1.^2+x2.^2+x3.^2+x4.^2);
endfunction
__________________________________________________________________________________________________________
pauli.sce
__________________________________________________________________________________________________________
// Pauli matrix
x=[0,0.5;0.5,0];
y=[0,-0.5*%i;0.5*%i,0];
z=[0.5,0;0,-0.5];
u=eye(2,2);
More information about the users
mailing list