% genlp
% --- Clovis Gonzaga, 3-22-2003
% Generates an lp problem with data read from lpdata.m, which must be edited.
% The results are saved in lp.mat and the memory is cleared.
% The command load lp loads the problem.
% Problem: min c'x subject to Ax=b, x>=0, where
% A is sparse and has a given density,
% the dimensions of the primal and dual faces are given
% The parameter 'option' (from lpdata.m) determines how the initial point is
% chosen:
%    option=0: no initial feasible point.
%    option=1: (x0,s0)=(e,e).
%    option=2: random initial feasible point satisfying
%              x0's0=n and x0'sopt+s0'xopt=n.

% The program reads the data from the file lpdata.m, which
% should be edited. The problem generated is saved in the
% file lp.mat. The command  'load lp' loads the following data:
%    A, b, c, m, n, option, pdim, ddim, xopt sopt p rankAB x0 s0,
% where pdim is the dimension of the primal optimal face,
%       ddim is the dimension of the dual optimal face,
%       xopt is an optimal (interior) primal solution,
%       sopt is an optimal (interior) dual solution,
%       x0, s0 is an interior initial primal-dual solution (if option>0).
%       p and rankAB are defined in the paper, and may be useful.

%read problem data
   lpdata;

%Number of large primal variables:
   if pdim>n-m-1,
     disp('Dimension of primal face too large')
     pdim=n-m-1;
   end
   if ddim>m-1,  
     disp('Dimension of dual face too large')
     ddim=m-1;
   end
   p=m+pdim-ddim;
   rankAB=p-pdim;
%The submatrix A_B will have rankAB linearly independent rows.
%optimal primal solution xopt'=[ >0 | 0 ]
%                                   p   n

   temp=xopt_l + (xopt_u-xopt_l)*rand(p,1);
   xopt=[temp; zeros(n-p,1)];

%optimal dual solution  sopt'= [ 0  |  >0 ]
%                                   p     n

   temp=sopt_l + (sopt_u-sopt_l)*rand(n-p,1);
   sopt=[zeros(p,1); temp];

% Generate the initial point if option>0.
   if option>0
      if option==1, x0=ones(n,1); s0=x0; 
      else 
         x0=xopt_l + (xopt_u-xopt_l)*rand(n,1);
         s0=sopt_l + (sopt_u-sopt_l)*rand(n,1);
         temp=sqrt(x0'*s0/n);
         x0=x0/temp; s0=s0/temp;
      end;
   %Now x0's0=n. Enforce x0'sopt+s0'xopt=n.
      temp=(x0'*sopt+s0'*xopt)/n;
      xopt=xopt/temp; sopt=sopt/temp;
   end   

%  Matrix generation: first generate a random entry in each column
%  vectors: col=1:n, line(i)=position in line i, s random vector.
%  rankAB<m is possible.
   col=1:n;
   line=zeros(1,n);
   for i=1:rankAB,
      line(i)=i;
   end
   for i=1:m-rankAB,
      line(p+i)=rankAB+i;
   end   
% choose random line positions for remaining columns
   for i=[rankAB+1:p],
      line(i)= floor(1+rand*m);
   end
   for i=[p+m-rankAB+1:n],
      line(i)= floor(1+rand*m);
   end
% Construct a random vector s and the sparse matrix
   s=randn(1,n);
   A=sparse(line,col,s);
   A=A+sprandn(m,n,a_dens);
% Enforce the rank condition
A11=A(1:rankAB, 1:p);
A21=sparse(m-rankAB,p);
A2=A(1:m,p+1:n);
A=[[A11;A21],A2 ];
%Now the random matrix entries and the vector b and c:
   b=A*xopt;
   c=sopt;
   
   if option>0
   %Change last column of A so that Ax0=b.
      A(:,n)=(b-A(:,1:n-1)*x0(1:n-1))/x0(n) ;
   %Change first row of A so that A'[-1|..0..]'+s0=c.
      A(1,:)=(s0-c)';
      b(1)=A(1,:)*xopt; 
   end   
if make_figure,
   clf;
   spy(A,'k')
   ax1=axis;
   hold on
   ax=[-0.1*n 1.1*n -0.1*m 1.1*m];
   axis(ax);
   plot([0 n n 0 0],[0 0 m m 0],'-k');
   plot([m m],[0,m],'--r','linewidth',1);
   plot([p p],[0,m],'--r','linewidth',1);
   plot([rankAB rankAB],[0,m],'--r','linewidth',1);
   text(m,1.05*m,'m');
   text(p,1.05*m,'p');
   text(rankAB,1.05*m,'r');
   text(n,1.05*m,'n');
   text(1.05*n,m,'m');
   axis off;
   axis(ax1);

   spy(A)
   axis(ax);
   plot([0 p p 0 0],[0 0 m m 0],'-k','linewidth',2);
   title('Matrix A')
end;

%Save the results and clear memory.
   if option==0
   save lp A b c m n option pdim ddim xopt sopt p rankAB;
   else
   save lp A b c m n option pdim ddim xopt sopt p rankAB x0 s0;
   end
   clear
