%
 M = [M  q;-q' 0];
 [m,n] = size(M);
 if exist('toler') ~= 1 
   toler=1.e-6; 
 end
 if exist('gamma') ~= 1 
   gamma=1/sqrt(n);    
 end
 if exist('alpha') ~= 1 
   alpha=0.9;    
 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)-x.*s-(1-gamma)*XX*(M*x-s)];
   ds=M*dx+(1-gamma)*(M*x-s);
%
%  choose step-size
%
   nora = min([dx./x;ds./s]);
%
%  Check if feasibility can be reached early.
%
   if abs(nora) <= (1-gamma),
     nora = (1-gamma);
     s = s + ds/nora;
     x = x + dx/nora;
     M=M(1:n-1,1:n-1);
     if x(n) > toler^2,
       disp('Find a maximal feasible solution');
       s = s(1:n-1)/x(n);
       x = x(1:n-1)/x(n);
     else
       disp('The problem is infeasible');
       x = x(1:n-1);
       s = s(1:n-1);
     end;
     return;
   else
     nora = abs(alpha/nora);
     x = x + nora*dx;
     s = s + nora*ds;
     mu = x'*s/n;
     gap = mu,
   end;
 end;
%
% output solution
%
 n = n-1;
 M=M(1:n,1:n);
 if x(n+1) > toler^2,
   s=s(1:n)/x(n+1);
   x=x(1:n)/x(n+1);
   disp('Find an approximate complemetarity or maximal-feasible solution');
 else
   x=x(1:n)/s(n+1);
   s=s(1:n)/s(n+1);
   disp('The problem is infeasible');
 end;
 return
%
%  This program solves the monotone LCP.
%
%  See splcp.m
%
%  set eta = 1-gamma and gamma < 1.
%