%%*******************************************************************
%%  Read in a problem in SDPA sparse format.
%%
%%  [blk,Avec,C,b] = read_sdpa(fname)
%%
%%  Input: fname = name of the file containing SDP data in
%%                 SDPA foramt. 
%%
%% 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] = read_sdpa(fname); 

%%
%%  Open the file for input
%%
   compressed = 0; 
   if exist(fname)
      fid = fopen(fname,'r');
   elseif exist([fname,'.Z']); 
      compressed = 1; 
      unix(['uncompress ',fname,'.Z']);
      fid = fopen(fname,'r');
   elseif exist([fname,'.gz']); 
      compressed = 2; 
      unix(['gunzip ',fname,'.gz']);
      fid = fopen(fname,'r');
   else
      fprintf('** Problem not found, please specify the correct path. \n');
      return;
   end;
%%
%%  Clean up special characters and comments from the file 
%%
   [tmpr,count] = fscanf(fid,'%c');
   linefeeds = findstr(tmpr,char(10));
   comment_chars = '*"=';
   cumidx = [];
   for i=1:length(comment_chars)
      idx = findstr(tmpr,comment_chars(i));
      cumidx = [cumidx,idx];
   end;
   for j=length(cumidx):-1:1
      if (cumidx(j)==1) | (strcmp(tmpr(cumidx(j)-1),char(10)))
         tmpr(cumidx(j):linefeeds(min(find(cumidx(j)<linefeeds))))='';
      else
         tmpr(cumidx(j):linefeeds(min(find(cumidx(j)<linefeeds)))-1)='';
      end;
   end;
   special_chars = ',{}()';
   cumidx=[];
   for i=1:length(special_chars)
      idx = findstr(tmpr,special_chars(i));
      cumidx = [cumidx,idx];
   end;
   tmpr(cumidx) = blanks(length(cumidx));
   nztotal = length(tmpr);
%%
%% Close the file 
%%
   fclose('all');
   if compressed==1; unix(['compress ',fname]); end;
   if compressed==2; unix(['gzip ',fname]); end; 
%% 
%%  Next, read in basic problem size parameters.
%%
   datavec = sscanf(tmpr,'%f'); clear tmpr;
   if size(datavec,1) < size(datavec,2); datavec = datavec'; end; 
   m = datavec(1); 
   numblk  = datavec(2);
   blksize = datavec(2+[1:numblk]); 
   if size(blksize,1) > size(blksize,2); blksize = blksize'; end
%%
%% Get input  b.
%%
   idxstrt = 2+numblk; 
   b = datavec(idxstrt+[1:m]);    
   idxstrt = idxstrt+m; 
   b = -b;
%%
%% Construct blk
%%
   spblkidx = find( (blksize>0) & (blksize <= 20) ); 
   deblkidx = find( (blksize<0) | (blksize > 20) ); 
   denumblk = length(deblkidx); 
   spnumblk = length(spblkidx); 
   if ~isempty(spblkidx) 
      blk = cell(denumblk+1,2);  
      C = cell(denumblk+1,1); A = cell(denumblk+1,m); 
      blk{denumblk+1,1} = 's'; blk{denumblk+1,2} = blksize(spblkidx); 
   else
      blk = cell(denumblk,2);  
      C = cell(denumblk,1); A = cell(denumblk,m); 
   end
   for p = 1:length(deblkidx)
      n = blksize(deblkidx(p)); 
      if (n > 0); 
         blk{p,1} = 's'; blk{p,2} = n;
      else
         blk{p,1} = 'l'; blk{p,2} = abs(n); 
      end  
   end
%%
%% Construct single blocks of A,C
%%
   len = length(datavec);    
   Y = reshape(datavec(idxstrt+1:len),5,(len-idxstrt)/5)';   
   clear datavec;    
   Y = sortrows(Y,[1 2]); 
   matidx = [0; find(diff(Y(:,1)) ~= 0); size(Y,1)];
%%
   for k = 1:length(matidx)-1
      idx = [matidx(k)+1 : matidx(k+1)];       
      Ytmp  = Y(idx,1:5); 
      matno = Ytmp(1,1); 
      for p = 1:denumblk 
         n = blksize(deblkidx(p));   
         idx = find( Ytmp(:,2)==deblkidx(p) ); 
         Yp = Ytmp(idx,3:5); 
         len = length(idx); 
         if (matno == 0) 
            if (n > 0)
               tmp = -sparse(Yp(:,1),Yp(:,2),Yp(:,3),n,n); 
               C{p} = tmp + triu(tmp,1)'; 
            else
               tmp = -sparse(Yp(:,1),ones(len,1),Yp(:,3),abs(n),1); 
               C{p} = tmp; 
            end
         else
            if (n > 0)
               tmp = -sparse(Yp(:,1),Yp(:,2),Yp(:,3),n,n); 
               A{p,matno} = tmp + triu(tmp,1)'; 
            else
               tmp = -sparse(Yp(:,1),ones(len,1),Yp(:,3),abs(n),1);
               A{p,matno} = tmp;
            end
         end
      end
   end 
%%
%% Construct big sparse block of A,C 
%%
   if ~isempty(spblkidx) 
      cumspblksize = [0 cumsum(blksize(spblkidx))]; 
      n = sum(blksize(spblkidx)); 
      for k = 0:m
         idx = find( Y(:,1)==k ); 
         Ytmp = Y(idx,1:5); 
         Yp = []; 
         for p = 1:spnumblk 
            idx = find( Ytmp(:,2)==spblkidx(p) ); 
            offset = cumspblksize(p); 
            Yp = [Yp; Ytmp(idx,3)+offset, Ytmp(idx,4)+offset, Ytmp(idx,5)]; 
         end
         if (k == 0) 
            tmp = -sparse(Yp(:,1),Yp(:,2),Yp(:,3),n,n); 
            C{denumblk+1} = tmp + triu(tmp,1)';              
         else
            tmp = -sparse(Yp(:,1),Yp(:,2),Yp(:,3),n,n); 
            A{denumblk+1,k} = tmp + triu(tmp,1)'; 
         end
      end
   end 
%%  
%% 
   Avec = svec(blk,A,ones(size(blk,1),1)); 
%% 
%%******************************************************************








