function [A, b, C, X, y, Z, status] = genHardSDP(m, n, gap, NoP, NoD, dualHasInteriorFlag)
%This file generates general hard instance (Strict complementary
%slackness conditions fails) for SDPs
%   min trace CX s.t. Aop X = b  X psd
% input: 
%  m, n, the size of b and C.
%  gap: the desired strictly complementary gap; can be 0;
%  NoP: the maximal rank of the primal optimum
%  NoD: the maximal rank of the dual optimum;
%  dualhasInteriorFlag(optional): set to 1 means that the dual is
%    guranteed to to have interior; set to 0 means the dual might not
%    have interior.   
%Output: 
%  A hard SDP instance. with data A, b, C,  optimal solution
%  X, y, Z
%  where A{i} is matrix A_i.
%  status == 1, successful.  
%  status ==0, error.
  
% Author: Hua Wei, Henry Wolkowicz
%      University of Waterloo
%  Feb 1, 2006  
  
status =0;
  
%default parameters
if (nargin < 6) dualHasInteriorFlag =0; end

if (NoP+ NoD + gap ~= n)
    NoD= n-NoP - gap;
    %error(['error on the parameters: NoP + NoD + noofGap is not equal to n']);
    disp(['error on the parameters: NoP + NoD + noofGap is not equal to n']);
    disp(['NoD is reset:  NoD= n-NoP - gap']);
end

if (NoP<=0 | NoD <=0 | gap <0) 
  error(['NoP, NoP should be positive and gap should be non-negative.' ...
	 ' Check the input']);
end

% 1. construct Optimal X and Z.
Q = zeros(n,n);
while rank(Q) < n
  Q = rand(n,n); 
end
Q = orth(Q);


P = 1:NoP; G = NoP+1: NoP + gap ; D = NoP + gap +1 : n;

Dx = diag(rand(NoP,1) * 100 ) + 0.1;
Dz = diag(rand(NoD,1) * 100 ) + 0.1;

X = Q(:, P) * Dx * Q(:, P)';
Z = Q(:, D) * Dz * Q(:, D)';



% 2. generate random A_i
  A=cell(m,1);
  lo= -10000; up= 10000;
  Y1 = generateY1(gap,lo,up);
  Y2 = rand(NoD, NoP)*(up-lo)+lo; 
  %zeroIdx = find(sum(abs(Y2))==0);
  while ( Q(:,D)*Y2 == 0 )  
	'repeat here line 48 in genHardSDP'
	Y2 = rand(NoD, NoP)*(up-lo)+lo;
	zeroIdx = find(sum(abs(Y2))==0);
  end
  Y3 = rand(NoD, gap) * (up-lo) + lo;
  Y4 = rand(NoD, NoD) * (up-lo) + lo;
  Y4 = (Y4+Y4')/2;
  YY = [ zeros(NoP, NoP) zeros(NoP, gap) Y2' 
		 zeros(gap, NoP) Y1     Y3'
		 Y2  Y3               Y4];
  A{1} = Q * YY * Q';
  
  AQ = zeros(n*NoP, m);
  tmp= A{1}*Q(:,P);
  AQ(:,1) = tmp(:);

  i =2;
  while (i<= m)
	if dualHasInteriorFlag & i==2    
	  %setting this to 1 will not generate positive matrix
	  newMatrix = rand(n,n)*2*sqrt(up)+ (-sqrt(abs(lo)));
	  newMatrix = (newMatrix * newMatrix' + eye(n))/n;
	else	  
	  newMatrix = rand(n,n)*(up-lo)+lo;
	  newMatrix = (newMatrix + newMatrix')/2;
	end
    tmp = newMatrix*Q(:,P);
    AQ(:,i) = tmp(:);
    if (rank(AQ(:,1:i)) == i)
      A{i} = newMatrix;
      i = i+1;
    end
  end
  
  
b = zeros(m,1);  
for i=1:m
  b(i) = trace(A{i}* X);
end

C = zeros(n,n); 
y = rand(m,1)*(up-lo)+lo;
Tmp = zeros(n,n);
for i=1:m
  Tmp = Tmp + A{i} * y(i);
end
C = Tmp + Z;


status =1;  
  
return;  
  
function Y1 = generateY1(gap, lo, up )
  Y1 = rand(gap,gap)*(up-lo) + lo; %genearate Y1 in the range of (up-lo)
  Y1= (Y1+Y1')/2;  
  Y1 = Y1 + rand * ( 10 * (up-lo) )* eye(gap);
  tol = 0.1;
  while min(eig(Y1))< tol
	Y1 = rand(gap,gap)*(up-lo) + lo;
	Y1= (Y1+Y1')/2;  
	Y1 = Y1 + rand * (10* (up-lo) ) * eye(gap);
  end

  return;