maxbisep.m

该算法是用于量子运算的matlab程序

% maxbisep   Maximum for biseparable multi-qubit states.
%   maxbisep(op,list) gives maximum value for an
%   operator for biseparable states; 
%   list determines the biparitioning.
%   Uses numerical search. Parameters can be set in
%   the beginning of maxbisep.m.

function m=maxbisep(op1,list)

% Copyright (C) 2005  Geza Toth    E.mail: [e-mail]
%
% This program is free software; you can redistribute it and/or
% modify it under the terms of the GNU General Public License
% as published by the Free Software Foundation; see gpl.txt
% of this subroutine package.
%
% This program is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with this program; if not, write to the Free Software
% Foundation, Inc., 51 Franklin Street, Fifth Floor, 
% Boston, MA  02110-1301, USA.

% Parameters for the minimization
Delta=0.01;
Nit1=5000;
Nit2=20000;

% Qubits
d=2;

[sx,sy]=size(op1);
N=log2(sx)/log2(d);
N=floor(N+0.5);

listneg=setdiff(N:-1:1,list);
listneg=-sort(-listneg);
Nr=length(listneg);
if isequal([listneg list],N:-1:1),
  op=op1;
else
  op=reorder(op1,[listneg list],d);
end %if

% Dimensions
d1=d^length(listneg);
d2=d^length(list);

rmax=-Inf;

for n=1:Nit1
     if mod(n,100)==0,  randn('state',sum(100*clock));  end %if
     f1=(rand(d1,1)-0.5)+i*(rand(d1,1)-0.5);
     f2=(rand(d2,1)-0.5)+i*(rand(d2,1)-0.5);
     f=kron(f1,f2);
     r=real((f'*op*f)/(f'*f));
     if r>rmax,
         rmax=r;
         f1max=f1;
         f2max=f2;        
     end %if
 end %for
 
 f10=f1max;
 f20=f2max;
 r0=rmax;
 inital_guess=r0;
 
 % Second phase of the search
 for n=1:Nit2
     if mod(n,100)==0,  randn('state',sum(100*clock));  end %if
     f1=(rand(d1,1)-0.5)+i*(rand(d1,1)-0.5);
     f2=(rand(d2,1)-0.5)+i*(rand(d2,1)-0.5);
     f1=f10+Delta*f1;
     f2=f20+Delta*f2;
     f=kron(f1,f2);
     r=real((f'*op*f)/(f'*f));
     if r>r0,
         r0=r;
         f10=f1;
         f20=f2;        
     end %if
end %for
 
m=r0;
% s=kron(f10,f20);