function varargout = reducedLP(varargin);

% reducedLP  Solves linear programming problems in standard form.
%            The algorithm uses primal-dual methods where only the `M'
%            most nearly active constraints of the dual are taken into
%            account in the normal equation matrix.
%            Speedups in the time per iteration are thus expected
%            for linear programs where the number of (primal)
%            variables is much larger than the number of constraints.
%            Preliminary test results show that the number of iterations
%            to solution hardly increases, sometimes even decreases, for
%            M as small as 5 or 10 percent of the total number of constraints.
%            Accordingly, significant speedups in the time to solution
%            can be expected. 
%
%   [X,Y] = reducedLP(A,B,C) returns the primal optimal solution X
%   and the dual optimal solution Y to the linear programming problem
%   min C'X  subject to AX = B, X >= 0.
%
%   [X,Y,LOGS] = reducedLP(A,B,C) further returns 
%   a structure containing information collected during runtime.
%
%   reducedLP(A,B,C,OPTS) specifies options:
%     In general, the value 1 means 'active' and 0 means 'inactive'.
%     OPTS.M: number of constraints to keep [integer | {size(A,2)}]
%     OPTS.M_rel: M/size(A,2) [scalar | {OPTS.M/size(A,2)}]
%     OPTS.disp: verbose printing at each outer step [0 | {1} | 2]
%     OPTS.disp_fig: displays figure [{0} | 1]
%     OPTS.show_title: shows figure title [0 | {1}]
%     OPTS.tol: convergence tolerance [scalar | {1e-6}]
%     OPTS.tol_soft_constr: tolerance for soft constraint heuristic 
%          [scalar | {1e-14}]
%     OPTS.iter_ubound: maximal number of iterations 
%          [integer | {1e2}]
%     OPTS.beta: beta parameter [scalar | {.99}]
%     OPTS.xmax: xmax parameter [scalar | {1e15}]
%     OPTS.xunder: xunder [scalar | {1e-4}]
%     OPTS.x0: initial x (primal variables) [{automatically
%          selected}]
%     OTPS.y0: initial y (dual variables) [{automatically selected}]
%     OPTS.method: Primal-Dual Affine Scaling 
%          or Mehrotra Predictor Corrector  [PDAS | {MPC}]
%     OPTS.mu: way the centering parameter is computed [{Wright} | Mehrotra]
%     OPTS.SFSOP: use Mehrotra's SFSOP routine [{0} | 1]
%     OPTS.sparse: sparse equation solver [0 | {1}]     
%     OPTS.Chol: [0 for backslash operator | 1 for Choleski (CHOL)  
%          | {2} for CHOLINC(.,'inf')]
%     OPTS.loading: diagonal loading [scalar | {0}]
%     OPTS.nb_rep: number of repetitions (in order to stabilize
%          measured CPU times) [integer | {1}]
%     OPTS.log_scalar: return LOGS for scalar values [0 | {1}]
%     OPTS.log_dual: return LOGS for y [{0} | 1]
%     OPTS.log_primal: return LOGS for x and s [{0} | 1]
%
% Please report bugs to P.-A. Absil
% (http://www.montefiore.ulg.ac.be/~absil/)

% based on ipas35.m
% This is a copy of ipas35_fn03.m
% Copyright P.-A. Absil, 2005.


% ipas35_fn03.m - PA - Started Mon 06 Jun 2005, CSIT.
%   I have to test now.
% ipas35_fn02.m - PA - Started Mon 06 Jun 2005, kag1.
%   Go on.
% ipas35_fn01.m - PA - Started Sun 05 Jun 2005, CSIT.
%   Turn it into a function.


%  *********************************************************
%  *********************************************************
%  Process inputs and do error checking
%  *********************************************************
%  *********************************************************

if nargin<3
  error('Not enough input arguments.')
end

A = varargin{1};
b = varargin{2};
c = varargin{3};

if nargin==4
  opts = varargin{4};
else
  opts = [];
end

[m,n] = size(A);

% The variables set below determine the behaviour of the code. 
% Most of the variable names start with "sw_", which stands for "switch".
% Particular attention should be paid to the choices commented with
% "!!!".
% In general, the value '1' means "yes", "on", and zero means the contrary.

% Choice of M (or k here)

if isfield(opts,'M')
  k = opts.M;
elseif isfield(opts,'M_rel')
  k = round(n*opts.M_rel);
else
  k = n; % Default choice
end

% Display figure:
if isfield(opts,'disp_fig')
  disp_fig = opts.disp_fig;
else
  disp_fig = 0;
end

% Print info during run
if isfield(opts,'disp')
  sw_verbose = opts.disp
else
  sw_verbose = 0;  % !!! has to be checked
end

% Turn figure title on/off
if isfield(opts,'show_title')
  show_title = opts.show_title;
else
  show_title = 1;
end

% Convergence tolerance: (default in paper: 1e-8)
if isfield(opts,'tol')
  tol_conv = opts.tol;
else
  tol_conv = 1e-6;
end

% Tolerance for soft constraints.
if isfield(opts,'tol_soft_constr')
  tol_soft_constr = opts.tol_soft_constr;
else
  tol_soft_constr = 1e-14;  % Default: 1e-14
end

% Upper bound on number of iterations
if isfield(opts,'iter_ubound')
  iter_ubound = opts.iter_ubound;
else
  iter_ubound = 1e2;  % Default: 1e2 
end

% ** Parameters of the algorithm

if isfield(opts,'beta')
  beta = opts.beta;
else
  beta = .99;
end

if isfield(opts,'xmax')
  xmax = opts.xmax;
else 
  xmax = 1e15; 
end

if isfield(opts,'xunder')
  xunder = opts.xunder;
else
  xunder = 1e-4;
end


% Choice of x0
if isfield(opts,'x0')
  x0 = opts.x0;
  x0 = x0(:);
end

% Choice of y0
if isfield(opts,'y0')
  y0 = opts.y0;
  y0 = y0(:);
end

% Automatic choices of x0 and y0
% -1: x0, y0 given.
% 0: x0 = ones(m,1);
% 1: x0 min norm solution of Ax=b
% 2: x0 as in Mehrotra (SIOPT, 1992), but modified for dual-feasibility.
% 3: x0, y0, s0 as in Mehrotra. 
% Note that in this case y_0 is not always dual feasible 
% and s_0 is not the exact slack.

sw_xy0 = -1;
if ~exist('x0')
  sw_xy0 = 2;
  if ~exist('y0')
    sw_xy0 = 3;
  end
end

%  xxx Test feasibility if IPAS.

% Choice of method
% 0: Affine scaling (no mu parameter)
% 1: Mehrotra predictor-corrector, as described e.g.
%    in [NW99, Chap 14].

if isfield(opts,'method')
  if isequal(opts.method,'PDAS') | isequal(opts.method,'pdas')
    sw_MPC = 0;
  else
    sw_MPC = 1;
  end
else
  sw_MPC = 1;
end
  
% Choice of the way the mu parameter is selected:
% 0: as in Wright
% 1: as in Mehrotra.
if isfield(opts,'mu')
  if opts.mu == 'Wright'
    sw_mu_Meh = 0;
  elseif opts.mu == 'Mehrotra'
    sw_mu_Meh = 1;
  else 
    error(['OPTS.mu = ', opts.mu ,' not defined'])
  end
else
  sw_mu_Meh = 0;
end


% SFSOP is a subroutine defined in [Mehrotra 1992].
% 0: as in Wright (linear) - this is simpler than Mehrotra's SFSOP
% 1: as in [Mehrotra 1992] (quadratic expansion of central path)
% SFSOP stands for Step From Second-Order Polynomial
%
% Our experiments suggest that sw_mu_Meh and sw_Meh_SFSOP don't help.

if isfield(opts,'SFSOP')
  if opts.SFSOP == 1
    sw_Meh_SFSOP = 1;
  else
    sw_Meh_SFSOP = 0;
  end
else
  sw_Meh_SFSOP = 0;
end

sw_Meh_GTSF = 0;  % Default: 0  (1 not implemented)
% 0: "Classical" update.
% 1: With GTSF heuristic in [Meh92, p588]. "The step factor in the primal space is
% chosen so that the product of primal blocking variable (ls) and the
% corresponding dual slack is nearly equal to the value their product
% would take at the point on the central trajectory at which the
% duality gap is equal to (x-p_x)'(s-p_s)/(n*gamma_a)."


% Choose dual update rule
if sw_MPC == 0
  sw_d_update = 1;
else
  sw_d_update = 0;
end
% Default: 1 (if sw_MPC=0) / 0 (if sw_MPC=1)  !!!
% 0: Classical update rule (clipped below by beta tbar).
% 1: Update rule involving tbar-||Delta y||, for superlinear convergence.

% Choose primal update rule
if sw_MPC == 0
  sw_p_update = 2;
else
  sw_p_update = 0;
end
% Default: 2 (if sw_MPC=0) / 0 (if sw_MPC=1) !!!
% 0: Classical update rule along Delta x.
% 1: First update rule (not along Delta x, without x_-).
% 2: Update rule involving x_-.

sw_sc_dy = 0; % Default: 0.   
% Turns norm(dy) into norm(dy./dy0) in update rules.
% New, 10 May 2005.

sw_sc_Q = 0; % Default: 0.  
% Selects Q on the basis of (c-A'*y)./abs(c-A'*y0)
% New, 10 May 2005.

sw_force_decr_xs = 0;  % Default: 0
% 1: Scale the updates (new 17 Oct 2004)
psi_xs = .75;  % must be in (.5,1)
% Goal: get global convergence for affine scaling non feasible,
% by forcing x's to decrease (the p and d residuals will decrease,
% too, by linearity).
% Result: does not work well. We hit the constraints x>0, s>0 too
% soon, i.e., t primal and t dual are very small.
% This is because the Newton step tries to achieve c-A'y-s=0. But if
% c-A'y has negative components, then s will try to go negative...

% Choose dense/sparse solver
% Involved in the computation of the primal-dual matrix.
% 0: No sparse (full matrices).
% 1: Sparse.

if isfield(opts,'sparse') && opts.sparse == 0
  sw_sparse = 0;
else
  sw_sparse = 1;
end

% xxx Use && everywhere.

if isfield(opts,'Chol')
  sw_Chol = opts.Chol;
else
  sw_Chol = 2;
end
% Default: 1/2  !!!
% 0: Solve systems using Matlab's backslash operator.
% 1: Solve systems using chol.
% 2: Force Cholinc.

if isfield(opts,'loading')
  val_loading = opts.loading;
else
  val_loading = 0;
end
% !!! Default: 0 (or 1e-8)
% Diagonal loading.
% Purpose: help the approximate normal matrix be pd

sw_Q = 1;  % Default: 1
% Select heuristic for the choice of the index set Q.
% 0: Q = {1,...,n}
% 1: Q = indices of the k smallest entries of s.
% 2: Q = (k largest x/s) U (m smallest s).
% 3: Q = k largest x/s.
% 4: Q = (ind k smallest s) U (ind m lin ind col) 

sw_soft_constr = 1;  % Default: 1
% Soft constraints allow the constraints to slightly recede.
% 0: Hard constraints.
% 1: Soft constraints: s = max(1e-14,s)

sw_xtilde = 0;  % Default: 0
% 0: Compute full xtilde.
% 1: Compute xtildeQ and put the rest to zero.

% First feasibility check
d_feasible = 1; % Default: 1
% Assume that the problem is dual-feasible.
% Will be modified later if necessary.
% This variable has to be eventually set to 0 if y0 is not feasible.

if isfield(opts,'y0') && any(c-A'*y0<0)
    d_feasible = 0;
end

if ~sw_MPC & ~d_feasible
  error('PDAS with dual infeasible initial conditions')
end

% Choose number of repetitions
if isfield(opts,'nb_rep')
  nb_rep = opts.nb_rep;
else
  nb_rep = 1;
end
% This makes it possible to run problems several times and return the
% mean CPU over all runs.
%   xxx May be changed later?

% Determine what the LOG structure should return
if nargout >= 3
  if isfield(opts,'log_scalar')
    sw_log_scalar = opts.log_scalar;
  else
    sw_log_scalar = 1;
  end
  if isfield(opts,'log_dual')
    sw_log_dual = opts.log_dual;
  else
    sw_log_dual = 0;
  end
  if isfield(opts,'log_primal')
    sw_log_primal = opts.log_primal;
  else
    sw_log_primal = 0;
  end
else
  sw_log_scalar = 0;
  sw_log_dual = 0;
  sw_log_primal = 0;
end
   
    
%  *********************************************************
%  *********************************************************
%  Various stuff
%  *********************************************************
%  *********************************************************


% ******* If requested, find m lin ind cols of A:
%if sw_Q==4
%  load woodw_Q_li
%end
if sw_Q==4 % & 1==0
  first_zero = 1;
  bl_nr = 0;
  while first_zero == 1
    bl_nr = bl_nr+1;
    [Qt,Rt,Et] = qr(full(A(:,(bl_nr-1)*m+1:bl_nr*m)));
    ind_zero = find(abs(diag(Rt))<1e-5);
    first_zero = ind_zero(1)
  end
  Rt = Rt(:,1:first_zero-1);
  Q_li = [1:m]*Et;
  Q_li = Q_li(1:first_zero-1);
  Q_li = Q_li + (bl_nr-1)*m; 
  size_Q_li = length(Q_li);
  j = bl_nr*m;
  while size_Q_li < m
    j = j+1;
    if ~mod(j,100), j
    end
    [Qtt,Rtt] = qrinsert(Qt,Rt,size_Q_li+1,full(A(:,j)));
    width_Rtt = size(Rtt,2);
    if abs(Rtt(width_Rtt,width_Rtt)) > 1e-5;
      Qt = Qtt;
      Rt = Rtt;
      Q_li = [Q_li,j];
      size_Q_li = size_Q_li+1
      if rank(full(A(:,Q_li)))==size_Q_li
	disp('Full rank')
      else
	disp('Rank deficient!')
	pause
      end
    end
  end
end  % if sw_Q==4

% ******* clear all evol variables (quite tricky!)
to_del = char(who('*evol'));
for i=1:size(to_del,1)
  var_to_del = to_del(i,:);
  j = length(var_to_del);
  while var_to_del(j)==' '
    j=j-1;
  end
  var_to_del = var_to_del(1:j);
  clear(var_to_del)
end
% ******* clear all _v variables (quite tricky!)
to_del = char(who('*_v'));
for i=1:size(to_del,1)
  var_to_del = to_del(i,:);
  j = length(var_to_del);
  while var_to_del(j)==' '
    j=j-1;
  end
  var_to_del = var_to_del(1:j);
  clear(var_to_del)
end



if sw_verbose > 0
  k
end

CPU_tic = cputime;
CPU_normal_tot = 0;
CPU_solve_tot = 0;
CPU_step2_tot = 0;
CPU_x0_tot = 0;
CPU_Q_tot = 0;
CPU_stop_tot = 0;
iter_tot = 0;

%  *********************************************************
%  *********************************************************
%  LOOP ON REPETITIONS
%  *********************************************************
%  *********************************************************

for rep=1:nb_rep
  

% ** Data
CPUt_x0=cputime;
switch sw_xy0
 case 0
  x = ones(n,1);
 case 1
  [U,R] = qr(A',0);
  x = U*(R'\b);
  x = x + max(0,-1.5*min(x));
  x0 = x;
 case 2  % Mehrotra for x_0 only (not y_0 nor s_0)
  % [U,R] = qr(A',0);  % This is very slow !
  % Use this instead:
  if ~exist('s0')
    s0 = c-A'*y0;
  end
  Af = full(A);
  [U,R] = qr(Af',0);
  x = U*(R'\b);
  xpdx = x + max(0,-1.5*min(x));
  x = xpdx + .5*(xpdx'*s0)/sum(s0);
 case 3  % Mehrotra SIOPT 2(4):575-601, 1992, p589.
  % Beware! This modifies initial y and s, too.
  % [U,R] = qr(A',0);  % This is very slow !
  % Use this instead:
  Af = full(A);
  [U,R] = qr(Af',0);
  y0 = R\(U'*c);
  s_tilde = c-A'*y0;
  x_tilde = U*(R'\b);
  dx = max(0,-1.5*min(x_tilde));
  ds = max(0,-1.5*min(s_tilde));
  dx_tilde = dx + .5*(x_tilde+dx)'*(s_tilde+ds)/sum(s_tilde+ds);
  ds_tilde = ds + .5*(x_tilde+dx)'*(s_tilde+ds)/sum(x_tilde+dx);  
  % y = y_0;
  s = s_tilde + ds_tilde;
  x = x_tilde + dx_tilde;
  d_feasible = 0;
end
y = y0;
if d_feasible
  s = c-A'*y;  % Ensure dual feasibility (dual equality satisfied)
end
CPU_x0_tot = CPU_x0_tot + cputime - CPUt_x0;


% ** Prepare storage:
%D2AT = zeros(n,m);


test_conv = 0;
test_div = 0;
iter = 0;

% *** DISPLAY INFOS
if sw_verbose==2
  disp('***')
  iter
  %err
  %max(y)
  disp(['c-A^Ty>=0? ', num2str(all(c-A'*y>=0))])
  disp(['c-A^Ty-s=0? ', num2str(all(abs(c-A'*y-s)<=1e-6))])
  %that_p
  %that_d
  %beta_xs
  disp(['xTs = ', num2str(x'*s)])
  disp(['Primal resid = ', num2str(norm(b-A*x))])
  disp(['Dual resid = ', num2str(norm(c-A'*y-s))])
end

%  *********************************************************
%  LOOP ON ITERATES (y,x,s)
%  *********************************************************

% ************* Loop on iterates (y,x,s) *************

while ~test_conv & ~test_div

  if d_feasible
    cmAty = s;  % cmAty = c-A'*y
  else
    cmAty = c-A'*y;
  end

  if sw_sc_Q & iter==0
    cmAty0 = cmAty;
  end
  
  CPUtt_tic = cputime;
  % ** Choose Q:
  switch sw_Q
   case 0
    Q = [1:n];
   case 1
    if ~sw_sc_Q
      [ignore,Qfull] = sort(cmAty);
    else
      [ignore,Qfull] = sort(cmAty./abs(cmAty0));
    end
    Q = Qfull(1:k);
   case 2  % the k largest x/s and the m smallest s
    [ignore,ind_s_sort] = sort(cmAty);
    [ignore,ind_xovers_sort] = sort(x./cmAty);
    Q = union(ind_s_sort(1:m),ind_xovers_sort(end-k+1:end));
   case 3
    [ignore,ind_xovers_sort] = sort(x./cmAty);
    Q = ind_xovers_sort(end-k+1:end);
   case 4   % (ind k smallest s) U (ind m lin ind col)  
    [ignore,Qfull] = sort(cmAty);
    Qtmp = Qfull(1:k);
    Q = union(Qtmp,Q_li);
  end
  CPU_Q_tot = CPU_Q_tot + cputime - CPUtt_tic;
  
  % ** Soft constraints:
  if sw_soft_constr
    s = max(tol_soft_constr,s);
  end
    
  % *** STEP 1
  
  % ** Compute the approximate normal matrix:
  CPUtt_tic = cputime;
  lQ = length(Q);
  
  switch sw_sparse
   case 0
    %CPUt = cputime;
    D2ATQ = zeros(lQ,m);
    %cputime-CPUt
    %CPUt = cputime;
    for iQ = 1:lQ
      nloc = Q(iQ);
      D2ATQ(iQ,:) = x(nloc)/s(nloc)*transpose(A(:,nloc));
    end
    %cputime-CPUt
    %CPUt=cputime;
    AD2ATQ = A(:,Q)*D2ATQ;
    %cputime-CPUt
    %disp('**')
   case 1
    %CPUt = cputime;
    D2ATQ = spdiags(x(Q)./s(Q),0,lQ,lQ)*A(:,Q)';
    %cputime-CPUt
    %CPUt = cputime;
    AD2ATQ = A(:,Q)*spdiags(x(Q)./s(Q),0,lQ,lQ)*A(:,Q)';
    %cputime-CPUt
    %disp('**')
  end
  %cputime-CPUtt_tic
  CPU_normal_tot = CPU_normal_tot + cputime - CPUtt_tic;  

  % ** Diagonal loading
  if val_loading ~= 0
    AD2ATQ = AD2ATQ + val_loading * eye(m);
  end
  
  % ** Symmetrize:  !!
  %AD2ATQ = .5*(AD2ATQ+AD2ATQ');   % !!!!
    
  % ** Solve normal equations:
  CPUtt_tic = cputime;

  if d_feasible
    rhs = b;
  else
    dual_diff = cmAty-s;
    rhs = b+A*(1./s.*x.*dual_diff);
  end
  
  if sw_Chol==2 & ~issparse(AD2ATQ);  % sw_Chol=2: force Cholinc
    AD2ATQ = sparse(AD2ATQ);
  end
  switch sw_Chol
   case 0
    dy = AD2ATQ\rhs;
   case {1,2}  
    if issparse(AD2ATQ)
      R = cholinc(AD2ATQ,'inf');
      if any(any(isinf(R)))
	%disp('Singular PD matrix!!')
      end
    else
      [R,pp] = chol(AD2ATQ);
      if pp>0
	pp
	test_div = 1;
	%pause
      end
    end
    if test_div>0
      break
    end
    warning off
    dy = R\(R'\rhs);   % - p_pi1 in [Meh92]
    warning on
  end

  if d_feasible
    ds = -A'*dy;
  else
    ds = -A'*dy + dual_diff;  % -p_s1 in [Meh92]
  end
  
  switch sw_xtilde
   case 0
    xtilde = -x./s.*ds;
   case 1
    xtilde = zeros(n,1);
    xtilde(Q) = -x(Q)./s(Q).*ds(Q);
  end
  dx = xtilde-x;  % -p_x1 in [Meh92]
  CPU_solve_tot = CPU_solve_tot + cputime - CPUtt_tic;

  if sw_sc_dy & iter==0
    dy0 = dy;
  end

  ind_s_decr = find(ds<0);
  
  switch sw_MPC
   case 0  % Affine scaling
    % Nothing to do.
   case 1  % Predictor-corrector
    ind_x_decr = find(dx<0);
    ind_s_decr = find(ds<0);
    a_aff_p = min(min(-x(ind_x_decr)./dx(ind_x_decr)),1);
    if isempty(a_aff_p)
      a_aff_p = Inf;
    end
    a_aff_d = min(min(-s(ind_s_decr)./ds(ind_s_decr)),1);
    if isempty(a_aff_d)
      a_aff_d = Inf;
    end
    
    % Step MPC-2: compute mu, sigma
    mu = x'*s/n;
    mu_aff = (x+a_aff_p*dx)'*(s+a_aff_d*ds)/n;  % mdg in [Meh92, p 585]
    sigma = (mu_aff/mu)^3;  % sigma = Meh's (mdg/x's)^nu
                            % sigma*mu = Meh's mu
    dxds = dx.*ds;
    if sw_mu_Meh 
      % Do step 4-5 of Mehrotra (1992) p. 585, which is not in
      % Wright's book. 
      ef = (dx'*s./x.*dx + ds'*x./s.*ds)/(mu*n);  % ef in [Meh92, p585]
      if ef > 1.1
	sigma = sigma/min(a_aff_p,a_aff_d);
      end
    end
    
    %sigma = 0; 
    %dxds = zeros(n,1);
        
    % * Version 1: compute the cc's:
    % Step MPC-3:
    CPUtt_tic = cputime;
    rhs_cc = -2*(A*(s.^(-1).*(sigma*mu*ones(n,1)-dxds)));
    switch sw_Chol
     case 0
      dy_cc = AD2ATQ\rhs_cc;
     case {1,2}
      dy_cc = R\(R'\rhs_cc);
    end
      % -p_pi2 in [Meh92]
    ds_cc = -A'*dy_cc;  % (+)p_s2 in [Meh92]
    dx_cc = 2*s.^(-1).*(sigma*mu*ones(n,1)-dxds)-s.^(-1).*x.*ds_cc;
      % p_x2 in [Meh92] (in Step 3 of [Meh92], read p_x2)
    CPU_solve_tot = CPU_solve_tot + cputime - CPUtt_tic;

    switch sw_Meh_SFSOP
     case 0  % Do it simple, as in Wright's book
      dx_mpc = dx+.5*dx_cc;
      dy_mpc = dy+.5*dy_cc;
      ds_mpc = ds+.5*ds_cc;
     case 1  % Do it as in Mehrotra (1992)
      % Compute roots of quadratic polyns to find where the 
      % quadratic leaves the feasible set:
      % The equation is like x+ gamma_p dx + .5 gamma_d^2 dx_cc = 0
      % in the unknown gamma.
      gamma_p_plus = (-dx+sqrt(dx.^2-2*dx_cc.*dx))./dx_cc; 
      % !! divide by zero
      gamma_p_minus = (-dx-sqrt(dx.^2-2*dx_cc.*dx))./dx_cc;
      gamma_p_nonclipped = min([gamma_p_plus(isreal(gamma_p_plus)&(gamma_p_plus>0));...
			      gamma_p_minus(isreal(gamma_p_minus)&(gamma_p_minus>0))]);
      if ~isempty(gamma_p_nonclipped)
	gamma_p = min(gamma_p_nonclipped,1);
      else
	gamma_p = 1;
      end
      gamma_d_plus = (-ds+sqrt(ds.^2-2*ds_cc.*ds))./ds_cc; 
      % !! divide by zero
      gamma_d_minus = (-ds-sqrt(ds.^2-2*ds_cc.*ds))./ds_cc;
      gamma_d_nonclipped = min([gamma_d_plus(isreal(gamma_d_plus)&(gamma_d_plus>0));...
			      gamma_d_minus(isreal(gamma_d_minus)&(gamma_d_minus>0))]);
      if ~isempty(gamma_d_nonclipped)
	gamma_d = min(gamma_d_nonclipped,1);
      else
	gamma_d = 1;
      end
      dx_mpc = gamma_p*dx+.5*gamma_p^2*dx_cc;
      dy_mpc = gamma_d*dy+.5*gamma_d^2*dy_cc;   
      ds_mpc = gamma_d*ds+.5*gamma_d^2*ds_cc;
    end  % sw_Meh_SFSOP
    
    xtilde_mpc = x+dx_mpc;

    % * Version 2: Skip the cc and directly go for mpc.
    % !! Not ok for infeasible yet.
%    switch sw_Chol
%     case 0
%      dy_mpc = AD2ATQ\(b-A*(s.^(-1).*(sigma*mu*ones(n,1)-dxds)));
%     case 1
%      dy_mpc = R\(R'\(b-A*(s.^(-1).*(sigma*mu*ones(n,1)-dxds))));
%    end
%    ds_mpc = -A'*dy_mpc;
%    xtilde_mpc = s.^(-1).*(sigma*mu*ones(n,1)-dxds-x.*ds_mpc);
%    dx_mpc = xtilde_mpc-x;
    
    % Rename the updates:
    dy_aff = dy; ds_aff = ds; dx_aff = dx; xtilde_aff = xtilde;
    dy = dy_mpc; ds = ds_mpc; dx = dx_mpc; xtilde = xtilde_mpc;
    % Step MPC-4-5: apply the updates:
    %   Done later.
  end  % sw_MPC
  
  % *** STEP 2
  
  CPUtt_tic = cputime;
  if sw_sc_dy
    norm_dy = norm(dy./dy0);
  else
    norm_dy = norm(dy);
  end
  switch sw_Meh_GTSF
   case 0   % "Simple" version
    % ** (i) Compute largest feasible stepsize:
    ind_s_decr = find(ds<0);
    if isempty(ind_s_decr)
      tbar_d = Inf;
    else
      tbar_d = min(-s(ind_s_decr)./ds(ind_s_decr));
    end
    switch sw_d_update
     case 0
      that_d = min(beta*tbar_d,1);
     case 1  
      that_d = min(max(beta*tbar_d,tbar_d-norm_dy),1);
    end
    
    switch sw_p_update
     case 0  % Primal update along Delta x.
      ind_x_decr = find(dx<0);
      if isempty(ind_x_decr)
	tbar_p = Inf;
      else
	tbar_p = min(-x(ind_x_decr)./dx(ind_x_decr));
      end
      that_p = min(beta*tbar_p,1);  % !! Forgot a term... To do...
     case 1  % Former method (without x_-)
      % ** (ii) Compute xplus:
      xplus_tmp = min(max(min((norm_dy)^2+(norm(dx))^2,xunder),xtilde),xmax);
      dx = xplus_tmp-x;
      that_p = 1;
     case 2  % New update method (04 fev 2005), with x_-.
      xtilde_minus = min(xtilde,0);
      xplus_tmp = min(max(min((norm_dy)^2+(norm(xtilde_minus))^2,xunder),xtilde),xmax);
      dx = xplus_tmp-x;
      that_p = 1;
    end
    
     % ** (iii) Choose new Q:
     % See beginning of the loop.

    % ** Ensure decrease of x^Ts:
    if sw_force_decr_xs
      dxds = dx'*ds;
      if dxds > 0
	beta_xs = - psi_xs * (x'*that_d*ds+s'*that_p*dx)/... 
		  (that_p*that_d*dx'*ds);
      else
	beta_xs = Inf;
      end
      beta_xs = min(beta_xs,1);
    else
      beta_xs = 1;
    end
    if beta_xs < 0
      disp('beta_xs < 0 !')
      pause
    end
    
    % ** Updates: 
    xplus = x+beta_xs*that_p*dx;
    yplus = y+beta_xs*that_d*dy;
    splus = s+beta_xs*that_d*ds;
    
   case 1  % Update as GTSF in [Meh92, p 588]
    % ?? We don't understand how it could work (notably, why would f
    % be smaller than 1?)
  end
  
  CPU_step2_tot = CPU_step2_tot + cputime - CPUtt_tic;
  
  % *** STOPPING CRITERION
  
  CPUtt_tic = cputime;
  gap_Meh = (c'*xplus-b'*yplus)/(1+abs(b'*yplus));
  %resid = norm(b-A*xplus)+norm(xplus.*splus);
  err = norm(b-A*xplus)/(1+norm(xplus))+norm(c-A'*yplus-splus)/(1+norm(splus))+abs(gap_Meh);
  % !! Comput for line above to be optimized.
  %test_conv = (resid < tol_resid);
  %test_conv = gap_Meh < tol_gap;
  test_conv = err<tol_conv;
  if iter > iter_ubound
    test_div = 2;
  end
  CPU_stop_tot = CPU_stop_tot + cputime - CPUtt_tic;
  
  % *** DEBUG
  if sw_log_scalar
    err_evol(iter+1) = err;
    min_x_evol(iter+1) = min(x);
    min_s_evol(iter+1) = min(s);
    fn_primal_evol(iter+1) = c'*x; 
    fn_dual_evol(iter+1) = b'*y;
  end
  if sw_log_dual
    y_evol(:,iter+1) = y;
  end
  if sw_log_primal
    x_evol(:,iter+1) = x;
    s_evol(:,iter+1) = s;
  end
  
  if isnan(err)
    disp('NaN!')
    pause
  end
  
  
  % *** ACTIVATE UPDATES
  y = yplus;
  x = xplus;
  s = splus;
  iter = iter+1;

  % *** DISPLAY INFOS
  if sw_verbose==2
    disp('***')
    iter
    err
    max(y)
    disp(['c-A^Ty>=0? ', num2str(all(c-A'*y>=0))])
    disp(['c-A^Ty-s=0? ', num2str(all(abs(c-A'*y-s)<=1e-6))])
    that_p
    that_d
    beta_xs
    disp(['xTs = ', num2str(x'*s)])
    disp(['Primal resid = ', num2str(norm(b-A*x))])
    disp(['Dual resid = ', num2str(norm(c-A'*y-s))])
  end
  

end  % loop on iterates

%  *********************************************************
%  END LOOP ON ITERATES (y,x,s)
%  *********************************************************

% *** DEBUG

if sw_log_scalar
  err_evol(iter+1) = err;
  min_x_evol(iter+1) = min(x);
  min_s_evol(iter+1) = min(s);
  fn_primal_evol(iter+1) = c'*x; 
  fn_dual_evol(iter+1) = b'*y;
end
if sw_log_dual
  y_evol(:,iter+1) = y;
end
if sw_log_primal
  x_evol(:,iter+1) = x;
  s_evol(:,iter+1) = s;
end

iter_tot = iter_tot + iter;
if test_div 
  break 
end
end  % of for rep=1:nb_rep

%  *********************************************************
%  *********************************************************
%  END LOOP ON REPETITIONS
%  *********************************************************
%  *********************************************************

CPU_total = (cputime-CPU_tic)
disp(['Mean nb iter = ', num2str(iter_tot/nb_rep)])

%  *********************************************************
%  *********************************************************
%  Outputs
%  *********************************************************
%  *********************************************************

varargout{1} = x;
varargout{2} = y;

if nargout >= 3
%  varargout{3} = struct('nb_iter',[],'cpu_total',[],'cpu_normal',[], ...
%			'cpu_solve',[],'cpu_step2',[],'cpu_x0',[], ...
%			'cpu_Q',[],'cpu_stop',[],'err',[], ...
%			'fn_primal',[],'fn_dual',[],'min_x',[], ...
%			'min_y',[],'x',[],'y',[],'s',[]);
  if sw_log_scalar
    varargout{3}.nb_iter = iter_tot/nb_rep;
    varargout{3}.cpu_total = CPU_total/nb_rep;
    varargout{3}.cpu_normal = CPU_normal_tot/nb_rep;
    varargout{3}.cpu_solve = CPU_solve_tot/nb_rep;
    varargout{3}.cpu_step2 = CPU_step2_tot/nb_rep;
    varargout{3}.cpu_x0 = CPU_x0_tot/nb_rep;
    varargout{3}.cpu_Q = CPU_Q_tot/nb_rep;
    varargout{3}.cpu_stop = CPU_stop_tot/nb_rep;
    varargout{3}.err = err_evol;
    varargout{3}.fn_primal = fn_primal_evol;
    varargout{3}.fn_dual = fn_dual_evol;
    varargout{3}.min_x = min_x_evol;
    varargout{3}.min_s = min_s_evol;
  end
  
  if sw_log_dual
    varargout{3}.y = y_evol;
  end
  if sw_log_primal
    varargout{3}.x = x_evol;
    varargout{3}.s = s_evol;
  end
end






% Implementation of IPAS method proposed by Andre Tits.

% ipas35.m - PA - Started Mon 30 May 2005, CSIT.
%   Modify labels of figures.
%   Slightly modify vals_k.
% ipas34.m - PA - Started Tue 10 May 2005, CSIT.
%   Try scale invariance modifs.
% ipas33.m - PA - Started Tue 10 May 2005, CSIT.
%   From *32tt02.m.
% ipas32.m - PA - Started Sat 09 Apr 2005, CSIT.
%   Assess execution times in details, especially for sparse.
%   Count "cc" computation times in system solve timing.
%   Modify y0 for scsd1.
% ipas31.m - PA - Started Wed 09 Feb 2005, CSIT.
%   Add comments.
% ipas30.m - PA - Started Mon 07 Feb 2005, CSIT.
%   New version for safety.
% ipas29.m - PA - Started Sun 06 Feb 2005, CSIT.
%   Prepare for automatic runs.
% ipas28.m - PA - Started Sun 06 Feb 2005, CSIT.
%   More tests on Mehrotra heuristics, including non-dual-feasible
%   case.
%   Slighly modified definition of 'err' to fit better with [Meh92].
% ipas27.m - PA - Started Sat 05 Feb 2005, CSIT.
%   Restart working on Mehrotra as in [Mehrotra 1992].
% ipas26.m - PA - Started Thu 03 Feb 2005, CSIT.
%   New version...
%   Modify primal (x) update rule to use "x_-".
%   Force Cholinc.
%   Correct error in err=...
% ipas25.m - PA - Started Sun 17 Oct 2004, College Park.
%   Force decrease of x's.
% ipas24.m - PA - Started Sat 16 Oct 2004, College Park.
%   Test values of k in decreasing order and stop when failure.
% ipas23.m - PA - Started Thu 14 Oct 2004, College Park.
%   Infeasible version.
%   Note that Q is chosen from c-A'y, not from s.
%   Also, mu as in Mehrotra's paper.
% ipas22.m - PA - Started Fri 20 Aug 2004, CSIT.
%   Add diagonal loading heuristic.
%   Added also was initialization just as in Mehrotra.
%   Difficulty: y_0 is not necessarily feasible.
% ipas21.m - PA - Started Sun 01 Aug 2004, Houckaye.
%   Heuristic to ensure full-rank of A(:,Q).
% ipas20.m - PA - Started Sun 18 Jul 2004, Houckaye.
%   More tests on CPU times.
% ipas19.m - PA - Started Sat 17 Jul 2004, Houckaye.
%   Test problems FIT1D and FIT2D from Netlib.
% ipas18.m - PA - Started Sat 17 Jul 2004, Houckaye.
%   Test problem scsd6 from Netlib.
% ipas17.m - PA - Started Sat 26 Jun 2004, CSIT.
%   Test sw_problem = 2;
% ipas16.m - PA - Started Sat 19 Jun 2004, CSIT.
%   Other tests.
% ipas15.m - PA - Started Fri 18 Jun 2004, CSIT.
%   New version for security.
% ipas14.m - PA - Started Wed 16 Jun 2004, CSIT.
%   Alternate x0 (one of them like in Mehrotra).
% ipas13.m - PA - Started Wed 16 Jun 2004, CSIT.
%   Trying to debug the Mehrotra version.
% ipas12.m - PA - Started Wed 16 Jun 2004, CSIT.
%   Include Mehrotra predictor-corrector approach.
% ipas11.m - PA - Started Tue 15 Jun 2004, CSIT.
%   Completely erase previous version.
%   Implement according to [Woessner-Tits, 18 Jul 2003],
%   with choice between some variants.
%
% PA Absil, September 3, 1999.
% Version ipas08, September 15.
% What's new in ipas06: ok for n>2; delta_x <> dx; cost displayed.
% In ipas08: new data generator for "random-regular" polytope.
% Caution: this m-file contains a clear and a close('all')...

% [Meh 92] and [Mehrotra 1992] are "On the implementation of a
% primal-dual interior point method", Sanjay Mehrotra, SIOPT 2(4)
% 575-601, November 1992.