%%*****************************************************************************
%% checkdeprows: compute AAt to determine if the 
%%               matrices Ak are linearly independent. 
%%              
%% SDPT3: version 3.0 
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last Modified: 7 May 01
%%*****************************************************************************

   function  [Avec,b,y,indeprows,feasible] = checkdeprows(blk,Avec,b,y,rmdepconstr);
   
   global spdensity chol_options
  
   m = length(b);  
   AAt = sparse(m,m);  numdencol = 0; UU = []; 
   for p = 1:size(blk,1)
      pblk = blk(p,:); 
      if strcmp(pblk{1},'s'); 
         AAt = AAt + Avec{p}'*Avec{p};  
      else 
         [decolidx,spcolidx] = checkdense(Avec{p}');          
         if ~isempty(decolidx); 
            n2 = size(Avec{p},1); 
            dd = ones(n2,1); 
            len= length(decolidx); 
            dd(decolidx) = zeros(len,1);  
            AAt = AAt + (Avec{p}' *spdiags(dd,0,n2,n2)) *Avec{p};
            tmp = Avec{p}(decolidx,:)'; 
            UU = [UU tmp]; 
            numdencol = numdencol + len; 
         else
            AAt = AAt + Avec{p}'*Avec{p};  
         end
      end
   end
   if (numdencol > 0); fprintf('\n  number of dense column = %d.',numdencol); end;
%%
%%
%%
   feasible = 1;
   indeprows = [];
   if (nnz(AAt) > spdensity*m^2); 
      if issparse(AAt); AAt = full(AAt); end; 
      chol_options = 'matlab';
   else
      if ~issparse(AAt), AAt = sparse(AAt); end;
      if (norm(AAt-spdiags(diag(AAt),0,m,m),inf) == 0)
         chol_options = 'matlab'; 
      else
         chol_options = 'lipsol'; 
      end
   end
   if chol_options == 'lipsol';  
      symbfct(AAt);  
      btmp = randn(m,1); btmp = btmp/norm(btmp); 
      sol = linsysolve(AAt,UU,UU,btmp,1); 
      if any(abs(sol) > 1e10) | any(sol == nan); 
         if isempty(UU);    
            if (rmdepconstr); 
               fprintf('\n  ... removing dependent constraints. ');    
               [Avec,b,y,indeprows,feasible] = rmdeprows(blk,Avec,b,y);
            else
               fprintf('\n warning: constraints may be nearly linearly dependent.');
               fprintf('\n To remove these constraints,');
               fprintf(' re-run sqlp.m with OPTIONS.rmdepconstr = 1.');
            end
         else
            fprintf('\n warning: the sparse part of AAt may be nearly singular'); 
         end
      end
   elseif chol_options == 'matlab'; 
      if ~isempty(UU); AAt = AAt + UU*UU'; end; 
      [UA,indef] = chol(AAt);
      if indef | (min(abs(diag(UA))) < 1e-6); 
         if (rmdepconstr); 
            fprintf('\n  ... removing dependent constraints. '); 
            [Avec,b,y,indeprows,feasible] = rmdeprows(blk,Avec,b,y);
         else
            fprintf('\n warning: constraints may be nearly linearly dependent.');
            fprintf('\n To remove these constraints,');
            fprintf(' re-run sqlp.m with OPTIONS.rmdepconstr = 1.');
         end
      end
   end
%%
%%*****************************************************************************
%% rmdeprows: eliminate dependent matrices
%%
%% [Avec,b,y0,indeprows,feasible] = rmdeprows(blk,Avec,b,y0);
%%
%% Input:  A,C,b: data for the SDP instance.
%% Output: A,C,b: data for the SDP instance with dependent constraints 
%%                eliminated
%%         feasible : feasibility flag
%% 
%% SDPT3: version 3.0
%% Copyright (c) 1997 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: 2 Feb 01
%%*****************************************************************************

   function [Avec,b,y0,indeprows,feasible] = rmdeprows(blk,Avec,b,y0);

   feasible=1;
%%
%%-----------------------------------------
%% delete zero matrices
%%-----------------------------------------
%%
   m = length(b); 
   normA = zeros(1,m);
   for p = 1:size(blk,1)
       normA = normA + sqrt(sum(Avec{p}.*Avec{p})); 
   end
   if any(normA == 0)
      zerorows=find(normA == 0);
      if any(b(zerorows) ~= 0)
         fprintf('\n rmdeprows: problem infeasible\n');
         fprintf(' zero row(s) with nonzero rhs\n');
         fprintf(1,'%c',7);  % ring a bell
         feasible = 0; return;
      end;
      nonzerorows=find(normA > 0);
      for p = 1:size(blk,1)
          Avec{p} = Avec{p}(:,nonzerorows);
      end
      b  = b(nonzerorows);
      y0 = y0(nonzerorows); 
      fprintf('\n rmdeprows: %i zero-row deleted\n',length(zerorows));
   else
      nonzerorows = [1:m];
   end;
%%
%%-----------------------------------------
%% eliminate dependent matrices 
%%-----------------------------------------
%%  
    vecA = []; 
    for p = 1:size(blk,1); 
        vecA = [vecA Avec{p}']; 
    end
%%
%% row reduce [ vecA | b ] to remove dependent rows/detect infeasibility
%%
   [indeprows,feasible]=rowreduce(vecA,b);
   if feasible
      for p = 1:size(blk,1) 
         Avec{p} = Avec{p}(:,indeprows);
      end
      b  = b(indeprows);
      y0 = y0(indeprows); 
      m  = length(b);
      indeprows = nonzerorows(indeprows);  
   else
      indeprows = nonzerorows(indeprows);  
      fprintf('\n rmdeprows: problem infeasible\n');
      fprintf('Inconsistent constraints\n');
      fprintf(1,'%c',7);  % ring a bell
      return;
   end;
%%
%%*****************************************************************************
%% rowreduce: eliminates dependent rows/ detects inconsistencies
%%            adapted from the MATLAB m file rref.m
%%
%% SDPT3: version 2.1 
%% Copyright (c) 1998 by
%% K.C. Toh, M.J. Todd, R.H. Tutuncu
%% Last modified: July 1998
%%*****************************************************************************

   function [indeprows,feasible] = rowreduce(A,b,tol);
%%
   [m,n] = size(A);
   mn=min(m,n);
%% Compute the default tolerance if none was provided.
   if (nargin < 3), tol = max(m,n)*eps*norm(A,'inf'); end

   PsLs=speye(m); % matrix of permutations and lower tri. elementary operations
   P=1:m;         % vector of permutations
   indeprows=P; 
%%
%% Loop over the m x min(m,n) (sub)matrix.
%%
   i = 1;
   j = 1;
   jb = [];
   while (i <= m) & (j <= mn)  
      %% Find value and index of largest element in the remainder of column j.
      [p,k] = max(abs(A(i:m,j))); k = k+i-1;
      if (p <= tol)
         % The column is negligible, zero it out.
         A(i:m,j) = zeros(m-i+1,1);
         j = j + 1;
      else
         jb = [jb j];
         %% Swap i-th and k-th rows.
         %%
         PsLs([i k],:) = PsLs([k i],:);
         P([i k]) = P([k i]);
         A([i k],j:mn) = A([k i],j:mn);
      %%
      %% Subtract multiples of the first m columns of 
      %% pivot row from all the higher indexed rows
      %%
         nzidxA=find(A(i,j:mn))+j-1;
         nzidxP=find(PsLs(i,:));
         for k = [i+1:m];
            if abs(A(k,j)) > tol
               PsLs(k,nzidxP) = PsLs(k,nzidxP)-(A(k,j)/A(i,j))*PsLs(i,nzidxP);
               A(k,nzidxA) = A(k,nzidxA) - (A(k,j)/A(i,j))*A(i,nzidxA);
            end;      
         end;
         i = i + 1;
         j = j + 1;
      end;
   end;
%%
%% now the rest of the columns (if any and if necessary)
%%
   while (j<=n) & length(jb)<mn
      % first calculate the lower part of "reduced" column
        A(i:m,j)=PsLs(i:m,:)*A(:,j);
        [p,k] = max(abs(A(i:m,j))); k = k+i-1;
        if (p <= tol)
      	   % The column is negligible, move on.
           j = j + 1;
      	else
	   % Remember column index
  	   jb = [jb j];
	   % Swap i-th and k-th rows.
	   % 
           if k~= i
              A([i k],j) = A([k i],j);
	      PsLs([i k],:) = PsLs([k i],:);
              P([i k]) = P([k i]);
           end;
           nzidxP=find(PsLs(i,:));
           for k = [i+1:m]
	      if abs(A(k,j)) > tol
	         PsLs(k,nzidxP)= ...
                 PsLs(k,nzidxP)-(A(k,j)/A(i,j))*PsLs(i,nzidxP);
	      end      
	   end
	   i = i + 1;
	   j = j + 1;
	end;
    end;
%%
   feasible=1;
   if i<m+1
      b=PsLs*b;
      if any(abs(b(i:m))>tol)
         feasible=0;
      else
         fprintf('\n rmdeprows: %d dependent constraints removed.\n',m-i+1);
         indeprows=sort(P(1:i-1));
      end;
   end;
   return; 
%%=============================================================================




