%%**********************************************************
%% cheby0:
%%
%% minimize || p(d) ||_infty
%% p = polynomial of degree <= m such that p(0) = 1.
%%
%% Here d = n-vector
%%----------------------------------------------------------
%% [blk,Avec,C,b,X0,y0,Z0,objval,p] = cheby0(d,m,feas,solve);
%%
%% d = a vector.
%% m = degree of polynomial.
%% feas = 1 if want feasible starting point
%% = 0 if otherwise.
%% solve = 0 if just want initialization
%% = 1 if want to solve the problem
%%
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%**********************************************************

function [blk,Avec,C,b,X0,y0,Z0,objval,p] = cheby0(d,m,feas,solve)

if nargin <= 3; solve = 0; end;
if nargin <= 2; feas = 0; end;

if ~isreal(d); error(' d must be a real vector '); end;
if (size(d,1) < size(d,2)); d = d.'; end;
time=cputime;

n = length(d);
blk{1,1} = 's'; blk{1,2} = 2*ones(1,n);

e = ones(n,1);
V(:,1) = e/norm(e); R(1,1) = 1/norm(e);
for i =1:m
v = d.*V(:,i);
for j = 1:i %% Arnoldi iterations:
H(j,i) = (V(:,j))'*v; %% constructing upper-Hessenberg matrix.
v = v - H(j,i)*V(:,j); %% orthonormaliztion of Krylov basis.
end;
H(i+1,i) = norm(v);
V(:,i+1) = v/H(i+1,i);
R(1:i+1,i+1) = (1/H(i+1,i))*([0; R(1:i,i)] - [R(1:i,1:i)*H(1:i,i); 0]);
end;

v1(1:2:2*n,1) = ones(n,1);
v2 = v1;
d1 = [0; v1]; d2 = [v2; 0];
C = spdiags([d2 d1],[-1 1],2*n,2*n);
b = [zeros(2*m,1); 1];

iV = sqrt(-1)*V; v1 = zeros(2*n,1);
A = cell(1,2*m+1);
for k = 1:m
v1(1:2:2*n) = d.*V(:,k);
v2 = conj(v1);
d1 = [0; v1]; d2 = [v2; 0];
A{k} = spdiags([d2 d1],[-1 1],2*n,2*n);
end;
for k = 1:m
v1(1:2:2*n) = d.*iV(:,k);
v2 = conj(v1);
d1 = [0; v1]; d2 = [v2; 0];
A{k+m} = spdiags([d2 d1],[-1 1],2*n,2*n);
end;
A{2*m+1} = speye(2*n);

Avec = svec(blk,A,ones(size(blk,1),1));
if (feas == 1);
X0 = speye(2*n)/(2*n);
y0 = -1.1*normest(C)*[zeros(2*m,1); 1];
Z0 = ops(C,'-',Asum(blk,A,y0,1));
elseif (feas == 0);
[X0,y0,Z0] = infeaspt(blk,Avec,C,b);
end;
disp(sprintf('time taken to form A = %3.1f',cputime-time));

if solve;
[obj,X,y,Z] = sqlp(blk,Avec,C,b,X0,y0,Z0);
y = y(1:m) + sqrt(-1)*y(m+1:2*m);
x1 = R(1:m,1:m)*y(1:m);
p = [-x1(m:-1:1); 1];
objval = -mean(obj);
end;
%%==========================================================