%%*******************************************************************
%%  Read in a problem in SeDuMi format.
%%
%%  [blk,A,C,b] = read_sedumi(fname)
%%
%%  Input: fname.mat = name of the file containing SDP data in
%%                     SeDuMi format.
%%
%% SDPT3: version 3.0
%% Copyright (c) 2000 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%******************************************************************

  function [blk,Avec,C,b]=read_sedumi(fname);
%%
%%  First, load the matlab file containing At, c, b, and K
%%
  K.q = [];
  K.l = [];
  A = 0;
  At = 0;
  compressed = 0; 
  if exist([fname,'.mat.gz']); 
     compressed = 1; 
     unix(['gunzip ', fname,'.mat.gz']);
  elseif exist([fname,'.gz']); 
     compressed = 2; 
     unix(['gunzip ', fname,'.gz']);
  elseif exist([fname,'.mat.Z']); 
     compressed = 3; 
     unix(['uncompress ', fname,'.mat.Z']);
  elseif exist([fname,'.Z']); 
     compressed = 4; 
     unix(['uncompress ', fname,'.Z']);
  end
  if exist([fname,'.mat']) | exist(fname) 
     eval(['load ', fname]);
  else
     fprintf('** Problem not found, please specify the correct path. \n');
     return;
  end
  if (compressed == 1)
     unix(['gzip ', fname,'.mat']);
  elseif (compressed == 2)
     unix(['gzip ', fname]);
  elseif (compressed == 3)
     unix(['compress ', fname,'.mat']);
  elseif (compressed == 4)
     unix(['compress ', fname]);
  end

%%
  if (size(c,1) == 1), c = c'; end;
  if (size(b,1) == 1), b = b'; end;
  if (At == 0), At = A'; end;  clear A;
  [nn,mm] = size(At); if (max(size(c)) == 1); c = c*ones(nn,1); end; 

  if ~isfield(K,'l'); K.l = 0; end   
  if ~isfield(K,'q'); K.q = 0; end
  if ~isfield(K,'s'); K.s = 0; end
  if K.l == 0 | isempty(K.l); K.l = 0; end;
  if sum(K.q) == 0 | isempty(K.q); K.q = 0; end
  if sum(K.s) == 0 | isempty(K.s); K.s = 0; end 
%%
%%
%%
   m = length(b);
   rowidx = 0;  idxblk = 0;  
   if ~(K.l == 0) 
      len = K.l;   
      idxblk = idxblk + 1; 
      blk{idxblk,1} = 'l'; blk{idxblk,2} = K.l; 
      Avec{idxblk,1} = At(rowidx+[1:len],:); 
      C{idxblk,1} = c(rowidx+[1:len]); 
      rowidx = rowidx + len; 
   end
   if ~(K.q == 0) 
      len = sum(K.q); 
      idxblk = idxblk + 1; 
      blk{idxblk,1} = 'q'; blk{idxblk,2} = K.q; 
      Avec{idxblk,1} = At(rowidx+[1:len],:); 
      C{idxblk,1} = c(rowidx+[1:len]); 
      rowidx = rowidx + len;
   end
   if ~(K.s == 0) 
      smblkdim = 30; 
      blksize = K.s;  
      if size(blksize,2) == 1; blksize = blksize'; end
      blknnz = [0 cumsum(blksize.*blksize)];   
      deblkidx = find(blksize > smblkdim); 
      if ~isempty(deblkidx)
         for p = 1:length(deblkidx)
             idxblk = idxblk + 1; 
             n = blksize(deblkidx(p)); 
             pblk{1,1} = 's'; pblk{1,2} = n;
             blk(idxblk,:) = pblk; 
             Atmp = At(rowidx+blknnz(deblkidx(p))+[1:n*n],:); 
             Avec{idxblk,1} = sparse(n*(n+1)/2,m);
             warning_yes = 1;  
             for k = 1:m 
                 Ak = mexmat(pblk,Atmp(:,k),1);
                 if (norm(Ak-Ak','fro') > 1e-13);
                    if warning_yes;
                       fprintf(' matrix not symmetric, use 0.5*(Ak+Ak'')\n');
                       warning_yes = 0; 
                    end
                    Ak = 0.5*(Ak+Ak'); 
                 end
                 Avec{idxblk,1}(:,k) = svec(pblk,Ak,1); 
             end
             Ctmp = c(rowidx+blknnz(deblkidx(p))+[1:n*n]);
             Ctmp = mexmat(pblk,Ctmp,1);
             C{idxblk,1} = 0.5*(Ctmp+Ctmp'); 
          end 
      end
      spblkidx = find(blksize <= smblkdim);
      if ~isempty(spblkidx)
         pos = []; len = 0;
         for p = 1:length(spblkidx)
             n = blksize(spblkidx(p)); 
             len = len + n*(n+1)/2; 
             pos = [pos, rowidx+blknnz(spblkidx(p))+[1:n*n]];
         end 
         idxblk = idxblk + 1; 
         blk{idxblk,1} = 's';  blk{idxblk,2} = blksize(spblkidx); 
         Avec{idxblk,1} = sparse(len,m); 
         Atmp = At(pos,:); 
         warning_yes = 1; 
         for k = 1:m
             Ak = mexmat(blk(idxblk,:),Atmp(:,k),1); 
             if (norm(Ak-Ak','fro') > 1e-13);
                if warning_yes;
                   fprintf(' matrix not symmetric, use 0.5*(Ak+Ak'')\n');
                   warning_yes = 0; 
                end
                Ak = 0.5*(Ak+Ak'); 
             end
             Avec{idxblk,1}(:,k) = svec(blk(idxblk,:),Ak,1); 
         end
         Ctmp = c(pos); 
         Ctmp = mexmat(blk(idxblk,:),Ctmp,1); 
         C{idxblk,1} = 0.5*(Ctmp+Ctmp');                 
      end
   end   
%%

