Logo Search packages:      
Sourcecode: semidef-oct version File versions  Download package

etp.m

function [d,Q] = etp(A)
% [d,Q] = etp(A);
%
% solves the eductional testing problem:
%
% (primal problem) maximize     e'*d 
%                  subject to   A - diag(d) >= 0
%                               d >= 0
%
% (dual problem)   minimize     Tr A*Q
%                  subject to   Q >= 0
%                               diag(Q) >= e
% 
% (e is the vector with all components one)
% 
% Input arguments:
% - A:             nxn, positive definite matrix.
%
% Output arguments:
% - d:             n-vector, solution for primal problem.
% - Q:             nxn, solution for dual problem.
%
% The relative duality gap is less than 0.1 percent.

n = size(A,1);  
if (size(A,2) ~= n) error('A must be square.'); end;
if (norm(A - A') > 1e-10*norm(A)) error('A must be symmetric.'); end;
lambdamin = min(eig(A)); 
if lambdamin < 1e-8, error('A must be positive definite.');  end;

% c = -e
c = -ones(n,1);    

% F = [A -E1 ... -En]  (Ei: zero except for E_ii = 1) 
F = [A zeros(n,n*n)];
F(:,n+[1:n+1:n*n]) = -eye(n);

#printf("F=");
#F
#prompt
   
% F = [ A(:) -E1(:) ... -En(:); 0 e1 ... en ] (ei: i-th unit vector)
F = [reshape(F,n*n,n+1); zeros(n,1) eye(n)];
blck_szs = [n; ones(n,1)];

#printf("reshaped F=");
#F
#prompt

for mnum = 1:columns(F)
  Fi = 0;
  i0 = 0;
  j0 = 0;
  for bnum = 1:length(blck_szs)
    bs = blck_szs(bnum);
    bs2 = bs*bs;
    fmidx = i0 + (1:bs);
    Fi(fmidx,fmidx) = reshape(F(j0+(1:bs2),mnum),bs,bs);
    i0 = i0 + bs;
    j0 = j0 +bs2;
  endfor
  #printf("F_%d=\n",mnum);
  #Fi
  #prompt
endfor
   
% primal initial point
x0 = 0.9*lambdamin*ones(n,1); 

% dual initial point
Z0 = ;

#printf("initial parameters:")
#x0
#Z0
#prompt

nu=50;  abstol=1e-8;  reltol=1e-3;  tv=0.0;  maxiters=100; 
[d,Z,ul,info,time] = sp(F, blck_szs, c, x0, Z0, nu, abstol, reltol,...
                        tv, maxiters);
#d
#Z
#ul
#info
#time
#prompt

if(time(3) == maxiters)
  error("maximum number of iterations exceeded");
endif
Q = reshape(Z(1:n*n),n,n);

Generated by  Doxygen 1.6.0   Back to index