%
 M = [M  q;-q' 0];
 [m,n] = size(M);
 if exist('toler') ~= 1 
   toler=1.e-6; 
 end
 if exist('gamma') ~= 1 
   gamma=n/(1+n);    
 end
 if exist('alpha') ~= 1 
   alpha=0.7;    
 end
 x = ones(n,1);
 s = ones(n,1);
 mu= 1;
 gap = 1;
 iter =0;
 while gap >= toler,
%
   iter = iter + 1;
   XX=sparse(1:n,1:n,x);
   SS=sparse(1:n,1:n,s);
   dx=[XX*M+SS]\[gamma*mu*ones(n,1)-(1-gamma)*x.*s];
   ds=M*dx;
%
%  choose step-size
%
   nora = min([dx;ds]);
   if nora/mu > -toler,
     if dx(n) > toler/n,
       s=max(ds(1:n-1),zeros(n-1,1))/abs(dx(n));
       x=max(dx(1:n-1),zeros(n-1,1))/abs(dx(n));
       M=M(1:n-1,1:n-1);
       disp('Find an approximate maximal feasible solution');
       return;
     elseif dx(1:n-1)'*ds(1:n-1) < toler,
       x=max(dx(1:n-1),zeros(n-1,1))/abs(ds(n));
       s=max(ds(1:n-1),zeros(n-1,1))/abs(ds(n));
       M=M(1:n-1,1:n-1);
       disp('The problem is infeasible');
       return;
     end;
   end;
   nora = min([dx./x;ds./s]);
   nora = abs(alpha/nora);
   x = x + nora*dx;
   s = s + nora*ds;
   mu = x'*s/n;
   gap = 1/mu,
 end;
%
%  output solution
%
 M=M(1:n-1,1:n-1);
 if x(n)/s(n) > toler,
   s=s(1:n-1)/x(n);
   x=x(1:n-1)/x(n);
   disp('Find an approximate complemetarity or maximal-feasible solution');
 else
   x=x(1:n-1)/s(n);
   s=s(1:n-1)/s(n);
   disp('The problem is infeasible');
 end;
 return
%
%  This program solves the LCP in standard form:
%
%     s = Mx + q, x^Ts=0, (x,s) >= 0.
%
%  Input 
%      M: Sparse n-dimensional monotone matrix.
%      q: n-dimensional column vector
%      toler: relative stopping tolerance: the objective value close to 
%             the local optimal one in the range of tolerance. 
%             Default value: 1.e-6.
%      alpha: step size: 0 < alpha < 1. Default value: .7.
%      gamma: weight parameter: 0< gamma <=1. Default value: n/(n+1).
%  Output
%     x,s: a maximal feasible or complementarity solution.
%
%  Homogeneous self-dual for LCP
%
%  Comment: Use eta =0 and gamma > 1.
%