%%********************************************************************* 
%% chebymat: compute the Chebyshev polynomial of 
%%           deg m of a matrix B.
%%        
%% [blk,Avec,C,b,X0,y0,Z0,objval,p] = chebymat(B,m,feas,solve);
%%
%% B = a square matrix.
%% m = degree of polynomial. 
%% feas  = 1 if want feasible starting point
%%       = 0 if otherwise.
%% solve = 0 just to initialize
%%       = 1 if want to solve using sdp.m
%%       = 2 if want to solve using sdphlf.m
%%
%% p = Chebyshev polynomial of B in Matlab format.
%%
%% 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] = chebymat(B,m,feas,solve);

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

   if ~isreal(B); error(' B must be a real matrix '); end; 
   N = length(B); 
   blk{1,1} = 's'; blk{1,2} = 2*N;
   if (m >= N); error('degree >= size of B'); end;

   A = cell(1,m+1); 
   [V,H,R] = orthbasis(B,m); 
   A{1} = V{m+1};
   for k = [2:m+1];  A{k} = V{k-1};  end;
   [blk,Avec,C,b,X0,y0,Z0] = norm_min(A,feas,0); 

   if solve; 
      [obj,X,y,Z] = sqlp(blk,Avec,C,b,X0,y0,Z0);
      if ~isreal(B); y = y(1:m) + sqrt(-1)*y(m+1:2*m); end;
      h  = R(1:m+1,1:m+1) \ [zeros(m,1); 1]; 
      x1 = R(1:m,1:m)*(h(m+1)*y(1:m) + h(1:m));
      p = [1; -x1(m:-1:1)]; 
      objval = norm(polyvalm(p,B)); 
   end; 
%%===================================================================== 




