%                             [dense,Adotdden] = getdense(At,Ablkjc,K,pars)
% GETDENSE  Creates dense.{l,cols,q}.
%   Try to find small proportion of the cone primitives that appear
%   in a large proportion of the primal constraints.
%
% SEE ALSO sedumi
% ******************** INTERNAL FUNCTION OF SEDUMI ********************

function [dense,Adotdden] = getdense(At,Ablkjc,K,pars)

 %  
 %   This file is part of SeDuMi 1.05
 %   Copyright (C) 2001 Jos F. Sturm
 %     Dept. Econometrics & O.R., Tilburg University, the Netherlands.
 %     Supported by the Netherlands Organization for Scientific Research (NWO).
 %   Affiliation SeDuMi 1.03 and 1.04Beta (2000):
 %     Dept. Quantitative Economics, Maastricht University, the Netherlands.
 %   Affiliations up to SeDuMi 1.02 (AUG1998):
 %     CRL, McMaster University, Canada.
 %     Supported by the Netherlands Organization for Scientific Research (NWO).
 % 
 %   This program is free software; you can redistribute it and/or modify
 %   it under the terms of the GNU General Public License as published by
 %   the Free Software Foundation; either version 2 of the License, or
 %   (at your option) any later version.
 % 
 %   This program is distributed in the hope that it will be useful,
 %   but WITHOUT ANY WARRANTY; without even the implied warranty of
 %   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 %   GNU General Public License for more details.
 % 
 %   You should have received a copy of the GNU General Public License
 %   along with this program; if not, write to the Free Software
 %   Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 %

% ------------------------------------------------------------
% Initialize: colnz = #nz-constraints per LP/LOR variable,
%  h = MAX(5, max{#nz-constraints | PSD-blocks})
% "5" is taken as a number that should not be beyond the denq-quantile.
% ------------------------------------------------------------
 NORMDEN = 5;
 colnz = sum(spones(extractA(At,Ablkjc,0,3,1,K.lq+1))')';
 h = max([NORMDEN, sum(findblks(At,Ablkjc,3,[],K.sblkstart)')]);
 [N,m] = size(At);
% ------------------------------------------------------------
% Replace colnz-entries for the Lorentz-trace ("x1") variables by
% nnz-constraints for that Lorentz block. Namely, these variables
% must be removed if the Lorentz block (d[k]*d[k]'*a[k]) is removed.
% ------------------------------------------------------------
 i1 = K.mainblks(1); i2 = K.mainblks(2);       % Bounds Lorentz trace.
 Ablkq = extractA(At,Ablkjc,1,2,i1,i2);
 if i1 < i2          % MATLAB 5.0, 5.1 cannot handle spones(empty)
   Ablkq2 = findblks(At,Ablkjc,2,3,K.qblkstart); % Tag if any nonzero in block
   Ablkq = spones(spones(Ablkq) + Ablkq2);
   colnz(i1:i2-1) = sum(Ablkq');             % And store at Lorentz trace pos.
 end
% ------------------------------------------------------------
% FIND THE denq-quantile for DENSE COLUMNS:
% If e.g. pars.denq = 0.75 and pars.denf = 20:
%    Find 75% quantile spquant. anything denser than 20*spquant is
%    tagged as dense.
% NB: spquant is chosen such that all columns with nnz <= h are to left of it.
%   This is because we don't allow PSD-block removal.
% ------------------------------------------------------------
 bigcolnz = colnz(find(colnz > h));              % h must be left anyhow
 denqN = ceil(pars.denq * length(colnz)) - (N-length(bigcolnz));
 if denqN < 1                % denqN is quantile on 1:length(bigcolnz) basis
   spquant = h;              % Take h as threshold
 else
   ordnzs = sort(bigcolnz);
   spquant = ordnzs(denqN);  % Take spquant > h as threshold
 end
% ------------------------------------------------------------
% Find LP cols, LORENTZ cols, and LORENTZ blocks with >denf*spquant nzs
% ------------------------------------------------------------
 dense.cols = find(colnz > pars.denf * spquant);     %dense LP, Q-blk, Q-cols
 dense.q = find(colnz(i1:i2-1) > pars.denf * spquant);        %dense Q-blks
 dense.l = length(find(dense.cols < i1));            %#dense LP cols.
% ----------------------------------------
% The number of dense columns should be small relative to
% the number of constraints - otherwise we better do 1 full Chol.
% ----------------------------------------
 if length(dense.cols) > m / 2;
   dense.l = 0; dense.cols = []; dense.q = [];
 end
 if isempty(dense.q)
   Adotdden = sparse(m,0);         % = Ablkq(dense.q,:)' block struct.
 else
% --------------------------------------------------
% HANDLE DENSE LORENTZ BLOCKS
% Let Adotdden lists the nz-constraints for each dense q-block.
% --------------------------------------------------
   Adotdden = Ablkq(dense.q,:)';
 end
