diff --git a/matlab/lmmcp/lmmcp.m b/matlab/lmmcp/lmmcp.m
index 90a8f4195..f6c0864b0 100644
--- a/matlab/lmmcp/lmmcp.m
+++ b/matlab/lmmcp/lmmcp.m
@@ -22,9 +22,10 @@ function [x,FVAL,EXITFLAG,OUTPUT,JACOB] = lmmcp(FUN,x,lb,ub,options,varargin)
 %      preprocess : activate preprocessor for phase I (default = 1)
 %      presteps   : number of iterations in phase I (default = 20)
 %  Termination parameters
-%      epsilon2   : termination value of the merit function (default = 1E-16)
-%      MaxIter    : maximum number of iterations (default = 500)
+%      MaxIter    : Maximum number of iterations (default = 500)
 %      tmin       : safeguard stepsize (default = 1E-12)
+%      TolFun     : Termination tolerance on the function value, a positive 
+%                   scalar (default = sqrt(eps))
 %  Stepsize parameters
 %      m          : number of previous function values to use in the nonmonotone
 %                   line search rule (default = 10)
@@ -85,7 +86,7 @@ function [x,FVAL,EXITFLAG,OUTPUT,JACOB] = lmmcp(FUN,x,lb,ub,options,varargin)
 %            e-mail: kanzow@mathematik.uni-wuerzburg.de
 %                    petra@mathematik.uni-wuerzburg.de
 
-% ------------Initialization----------------
+%% Initialization
 defaultopt = struct(...
     'beta',       0.55,...
     'Big',        1e10,...
@@ -93,7 +94,6 @@ defaultopt = struct(...
     'deltamin',   1,...
     'Display',    'none',...
     'epsilon1',   1e-6,...
-    'epsilon2',   1e-16,...
     'eta',        0.95,...
     'kwatch',     20,...
     'lambda1',    0.1,...
@@ -106,6 +106,7 @@ defaultopt = struct(...
     'sigma1',     0.5,...
     'sigma2',     2,...
     'tmin',       1e-12,...
+    'TolFun',     sqrt(eps),...
     'watchdog',   1);
 
 if nargin < 4
@@ -122,21 +123,22 @@ else
   options = catstruct(defaultopt,options);
 end
 
+warning('off','MATLAB:rankDeficientMatrix')
 
 switch options.Display
-    case {'off','none'}
-        verbosity = 0;
-    case {'iter','iter-detailed'}
-        verbosity = 2;
-    case {'final','final-detailed'}
-        verbosity = 1;
-    otherwise
-        verbosity = 0;
+  case {'off','none'}
+    verbosity = 0;
+  case {'iter','iter-detailed'}
+    verbosity = 2;
+  case {'final','final-detailed'}
+    verbosity = 1;
+  otherwise
+    verbosity = 0;
 end
 
 % parameter settings
 eps1 = options.epsilon1;
-eps2 = options.epsilon2;
+eps2 = 0.5*options.TolFun^2;
 null = options.null;
 Big  = options.Big;
 
@@ -212,14 +214,14 @@ aux(1) = Psix;
 MaxPsi = Psix;
 
 if watchdog==1
-   kbest        = k;
-   xbest        = x;
-   Phibest      = Phix;
-   Psibest      = Psix;
-   DPhibest     = DPhix;
-   DPsibest     = DPsix;
-   normDPsibest = normDPsix;
-end;
+  kbest        = k;
+  xbest        = x;
+  Phibest      = Phix;
+  Psibest      = Psix;
+  DPhibest     = DPhix;
+  DPsibest     = DPsix;
+  normDPsibest = normDPsix;
+end
 
 % initial output
 if verbosity > 1
@@ -229,76 +231,74 @@ if verbosity > 1
   fprintf('%4.0f %24.5e %24.5e\n',k,Psix,normDPsix);
 end
 
-%
-%   Preprocessor using local method
-%
+%% Preprocessor using local method
 
 if preprocess==1
 
   if verbosity > 1
     disp('************************** Preprocessor ****************************')
   end
+  
+  normpLM=1;
+  while (k < presteps) && (Psix > eps2) && (normpLM>null)
+    k = k+1;
+    
+    % choice of Levenberg-Marquardt parameter, note that we do not use
+    % the condition estimator for large-scale problems, although this
+    % may cause numerical problems in some examples
 
-   normpLM=1;
-   while (k < presteps) && (Psix > eps2) && (normpLM>null)
-      k=k+1;
-
-      % choice of Levenberg-Marquardt parameter, note that we do not use
-      % the condition estimator for large-scale problems, although this
-      % may cause numerical problems in some examples
-
-      i=0;
-      mu=0;
-      if n<100
-         i=1;
-         mu=1e-16;
-         if condest(DPhix'*DPhix)>1e25
-            mu=1e-6/(k+1);
-         end
+    i  = false;
+    mu = 0;
+    if n<100
+      i = true;
+      mu = 1e-16;
+      if condest(DPhix'*DPhix)>1e25
+        mu = 1e-6/(k+1);
       end
-      if i==1
-         pLM= [ DPhix ; sqrt(mu)*speye(n)]\[-Phix;sparse(n,1)];
-      else
-         pLM=-DPhix\Phix;
-      end
-      normpLM=norm(pLM);
+    end
+    if i
+      pLM =  [DPhix; sqrt(mu)*speye(n)]\[-Phix; zeros(n,1)];
+    else
+      pLM = -DPhix\Phix;
+    end
+    normpLM = norm(pLM);
+    
+    % compute the projected Levenberg-Marquard step onto box Xk
+    lbnew = max(min(lb-x,0),-delta);
+    ubnew = min(max(ub-x,0),delta);
+    d     = max(lbnew,min(pLM,ubnew));
+    xnew  = x+d;
 
-      % compute the projected Levenberg-Marquard step onto box Xk
-      lbnew=max(min(lb-x,0),-delta);
-      ubnew=min(max(ub-x,0),delta);
-      d=max(lbnew,min(pLM,ubnew));
-      xnew=x+d;
+    % function evaluations etc.
+    [Fxnew,DFxnew] = feval(FUN,xnew,varargin{:});
+    Phixnew        = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
+    Psixnew        = 0.5*(Phixnew'*Phixnew);
+    normPhixnew    = norm(Phixnew);
+    
+    % update of delta
+    if normPhixnew<=eta*normPhix
+      delta = max(deltamin,sigma2*delta);
+    elseif normPhixnew>5*eta*normPhix
+      delta = max(deltamin,sigma1*delta);
+    end
 
-      % function evaluations etc.
-      [Fxnew,DFxnew] = feval(FUN,xnew,varargin{:});
-      Phixnew=Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
-      Psixnew=0.5*(Phixnew'*Phixnew);
-      normPhixnew=norm(Phixnew);
-
-      % update of delta
-      if normPhixnew<=eta*normPhix
-         delta=max(deltamin,sigma2*delta);
-      elseif normPhixnew>5*eta*normPhix
-         delta=max(deltamin,sigma1*delta);
-      end
-
-      % update
-      x=xnew;
-      Fx=Fxnew;
-      DFx=DFxnew;
-      Phix=Phixnew;
-      Psix=Psixnew;
-      normPhix=normPhixnew;
-      DPhix=DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset);
-      DPsix=DPhix'*Phix;
-      normDPsix=norm(DPsix,inf);
-
-      % output at each iteration
-      t=1;
-      if verbosity > 1
-        fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t);
-      end
-   end
+    % update
+    x         = xnew;
+    Fx        = Fxnew;
+    DFx       = DFxnew;
+    Phix      = Phixnew;
+    Psix      = Psixnew;
+    normPhix  = normPhixnew;
+    DPhix     = DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset);
+    DPsix     = DPhix'*Phix;
+    normDPsix = norm(DPsix,inf);
+    
+    % output at each iteration
+    t=1;
+    if verbosity > 1
+      fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t);
+    end
+  end
 end
 
 % terminate program or redefine current iterate as original initial point
@@ -323,12 +323,10 @@ elseif preprocess==1 && Psix>=eps2
   if verbosity > 1
     disp('******************** Restart with initial point ********************')
     fprintf('%4.0f %24.5e %24.5e\n',k_main,Psix0,normDPsix0);
-   end
+  end
 end
 
-%
-%   Main algorithm
-%
+%%   Main algorithm
 
 if verbosity > 1
   disp('************************** Main program ****************************')
@@ -336,96 +334,96 @@ end
 
 while (k < kmax) && (Psix > eps2)
 
-   % choice of Levenberg-Marquardt parameter, note that we do not use
-   % the condition estimator for large-scale problems, although this
-   % may cause numerical problems in some examples
+  % choice of Levenberg-Marquardt parameter, note that we do not use
+  % the condition estimator for large-scale problems, although this
+  % may cause numerical problems in some examples
 
-   i=0;
-   if n<100
-      i=1;
-      mu=1e-16;
-      if condest(DPhix'*DPhix)>1e25
-         mu=1e-1/(k+1);
-      end
-   end
+  i = false;
+  if n<100
+    i  = true;
+    mu = 1e-16;
+    if condest(DPhix'*DPhix)>1e25
+      mu = 1e-1/(k+1);
+    end
+  end
+  
+  % compute a Levenberg-Marquard direction
 
-   % compute a Levenberg-Marquard direction
+  if i
+    d = [DPhix; sqrt(mu)*speye(n)]\[-Phix; zeros(n,1)];
+  else
+    d = -DPhix\Phix;
+  end
 
-   if i==1
-      d= [ DPhix ; sqrt(mu)*eye(n)]\[-Phix;zeros(n,1)];
-   else
-      d=-DPhix\Phix;
-   end
+  % computation of steplength t using the nonmonotone Armijo-rule
+  % starting with the 6-th iteration
 
-   % computation of steplength t using the nonmonotone Armijo-rule
-   % starting with the 6-th iteration
+  % computation of steplength t using the monotone Armijo-rule if
+  % d is a 'good' descent direction or k<=5
 
-   % computation of steplength t using the monotone Armijo-rule if
-   % d is a 'good' descent direction or k<=5
+  t       = 1;
+  xnew    = x+d;
+  Fxnew   = feval(FUN,xnew,varargin{:});
+  Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
+  Psixnew = 0.5*(Phixnew'*Phixnew);
+  const   = sigma*DPsix'*d;
+  
+  while (Psixnew > MaxPsi + const*t)  && (t > tmin)
+    t       = t*beta;
+    xnew    = x+t*d;
+    Fxnew   = feval(FUN,xnew,varargin{:});
+    Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
+    Psixnew = 0.5*(Phixnew'*Phixnew);
+  end
 
-   t       = 1;
-   xnew    = x+d;
-   Fxnew   = feval(FUN,xnew,varargin{:});
-   Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
-   Psixnew = 0.5*(Phixnew'*Phixnew);
-   const   = sigma*DPsix'*d;
-
-   while (Psixnew > MaxPsi + const*t)  && (t > tmin)
-      t       = t*beta;
-      xnew    = x+t*d;
-      Fxnew   = feval(FUN,xnew,varargin{:});
-      Phixnew = Phi(xnew,Fxnew,lb,ub,lambda1,lambda2,n,Indexset);
-      Psixnew = 0.5*(Phixnew'*Phixnew);
-   end
-
-   % updatings
-   x=xnew;
-   Fx=Fxnew;
-   Phix=Phixnew;
-   Psix=Psixnew;
-   [~,DFx]=feval(FUN,x,varargin{:});
-   DPhix=DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset);
-   DPsix=DPhix'*Phix;
-   normDPsix=norm(DPsix);
-   k=k+1;
-   k_main=k_main+1;
-
-   if k_main<=5
-      aux(mod(k_main,m)+1)=Psix;
+  % updatings
+  x         = xnew;
+  Fx        = Fxnew;
+  Phix      = Phixnew;
+  Psix      = Psixnew;
+  [~,DFx]   = feval(FUN,x,varargin{:});
+  DPhix     = DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset);
+  DPsix     = DPhix'*Phix;
+  normDPsix = norm(DPsix);
+  k         = k+1;
+  k_main    = k_main+1;
+  
+  if k_main<=5
+    aux(mod(k_main,m)+1) = Psix;
+    MaxPsi               = Psix;
+  else
+    aux(mod(k_main,m)+1) = Psix;
+    MaxPsi               = max(aux);
+  end
+  
+  % updatings for the watchdog strategy
+  if watchdog ==1
+    if Psix<Psibest
+      kbest        = k;
+      xbest        = x;
+      Phibest      = Phix;
+      Psibest      = Psix;
+      DPhibest     = DPhix;
+      DPsibest     = DPsix;
+      normDPsibest = normDPsix;
+    elseif k-kbest>kwatch
+      x=xbest;
+      Phix=Phibest;
+      Psix=Psibest;
+      DPhix=DPhibest;
+      DPsix=DPsibest;
+      normDPsix=normDPsibest;
       MaxPsi=Psix;
-   else
-      aux(mod(k_main,m)+1) = Psix;
-      MaxPsi               = max(aux);
-   end;
+    end
+  end
 
-   % updatings for the watchdog strategy
-   if watchdog ==1
-      if Psix<Psibest
-         kbest        = k;
-         xbest        = x;
-         Phibest      = Phix;
-         Psibest      = Psix;
-         DPhibest     = DPhix;
-         DPsibest     = DPsix;
-         normDPsibest = normDPsix;
-      elseif k-kbest>kwatch
-         x=xbest;
-         Phix=Phibest;
-         Psix=Psibest;
-         DPhix=DPhibest;
-         DPsix=DPsibest;
-         normDPsix=normDPsibest;
-         MaxPsi=Psix;
-      end;
-   end;
+  if verbosity > 1
+    % output at each iteration
+    fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t);
+  end
+end
 
-   if verbosity > 1
-     % output at each iteration
-     fprintf('%4.0f %24.5e %24.5e %11.7g\n',k,Psix,normDPsix,t);
-   end
-end;
-
-% final output
+%% Final output
 if Psix<=eps2
   EXITFLAG = 1;
   if verbosity > 0, disp('Approximate solution found.'); end
@@ -446,9 +444,10 @@ OUTPUT.Psix       = Psix;
 OUTPUT.normDPsix  = normDPsix;
 JACOB             = DFx;
 
-% Subfunctions
+%% Subfunctions
 
 function y = Phi(x,Fx,lb,ub,lambda1,lambda2,n,Indexset)
+%% PHI
 
 y           = zeros(2*n,1);
 phi_u       = sqrt((ub-x).^2+Fx.^2)-ub+x+Fx;
@@ -472,7 +471,7 @@ y([LZ; I3]) = lambda2*(max(0,x(I3)-lb(I3)).*max(0,Fx(I3))+max(0,ub(I3)-x(I3)).*m
 
 
 function H = DPhi(x,Fx,DFx,lb,ub,lambda1,lambda2,n,Indexset)
-% DPHI evaluates an element of the C-subdifferential of operator Phi
+%% DPHI evaluates an element of the C-subdifferential of operator Phi
 
 null       = 1e-8;
 beta_l     = zeros(n,1);
@@ -521,14 +520,14 @@ Da(I1b)    = (x(I1b)-lb(I1b))./denom1(I1b)-1;
 Db(I1b)    = Fx(I1b)./denom1(I1b)-1;
 I1b        = Indexset==1 & beta_l~=0;
 if any(I1b)
-  Da(I1b)    = z(I1b)./denom2(I1b)-1;
-  Db(I1b)    = (DFx(I1b,:)*z)./denom2(I1b)-1;
+  Da(I1b)  = z(I1b)./denom2(I1b)-1;
+  Db(I1b)  = (DFx(I1b,:)*z)./denom2(I1b)-1;
 end
 
-I1a        = I(Indexset==1 & alpha_l==1);
+I1a         = I(Indexset==1 & alpha_l==1);
 if any(I1a)
-  H2(I1a,:)  = repmat(x(I1a)-lb(I1a),1,n).*DFx(I1a,:)+sparse(1:length(I1a),I1a, ...
-                                                    Fx(I1a),length(I1a),n,length(I1a));
+  H2(I1a,:) = bsxfun(@times,x(I1a)-lb(I1a),DFx(I1a,:))+...
+              sparse(1:length(I1a),I1a,Fx(I1a),length(I1a),n,length(I1a));
 end
 
 I2         = Indexset==2;
@@ -544,14 +543,14 @@ Da(I2b)    = (ub(I2b)-x(I2b))./denom1(I2b)-1;
 Db(I2b)    = -Fx(I2b)./denom1(I2b)-1;
 I2b        = Indexset==2 & beta_u~=0;
 if any(I2b)
-  Da(I2b)    = -z(I2b)./denom2(I2b)-1;
-  Db(I2b)    = -(DFx(I2b,:)*z)./denom2(I2b)-1;
+  Da(I2b)  = -z(I2b)./denom2(I2b)-1;
+  Db(I2b)  = -(DFx(I2b,:)*z)./denom2(I2b)-1;
 end
 
-I2a        = I(Indexset==2 & alpha_u==1);
+I2a         = I(Indexset==2 & alpha_u==1);
 if any(I2a)
-  H2(I2a,:)  = repmat(x(I2a)-ub(I2a),1,n).*DFx(I2a,:)+sparse(1:length(I2a),I2a, ...
-                                                    Fx(I2a),length(I2a),n,length(I2a));
+  H2(I2a,:) = bsxfun(@times,x(I2a)-ub(I2a),DFx(I2a,:))+...
+              sparse(1:length(I2a),I2a,Fx(I2a),length(I2a),n,length(I2a));
 end
 
 I3         = Indexset==3;
@@ -593,23 +592,23 @@ end
 Da(I3)     = ai(I3)+bi(I3).*ci(I3);
 Db(I3)     = bi(I3).*di(I3);
 
-I3a        = I(Indexset==3 & alpha_l==1 & alpha_u==1);
+I3a         = I(Indexset==3 & alpha_l==1 & alpha_u==1);
 if any(I3a)
-  H2(I3a,:)  = repmat(-lb(I3a)-ub(I3a)+2*x(I3a),1,n).*DFx(I3a,:)+...
-      2*sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
+  H2(I3a,:) = bsxfun(@times,-lb(I3a)-ub(I3a)+2*x(I3a),DFx(I3a,:))+...
+              2*sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
 end
-I3a        = I(Indexset==3 & alpha_l==1 & alpha_u~=1);
+I3a         = I(Indexset==3 & alpha_l==1 & alpha_u~=1);
 if any(I3a)
-  H2(I3a,:)  = repmat(x(I3a)-lb(I3a),1,n).*DFx(I3a,:)+...
-      sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
+  H2(I3a,:) = bsxfun(@times,x(I3a)-lb(I3a),DFx(I3a,:))+...
+              sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
 end
-I3a        = I(Indexset==3 & alpha_l~=1 & alpha_u==1);
+I3a         = I(Indexset==3 & alpha_l~=1 & alpha_u==1);
 if any(I3a)
-  H2(I3a,:)  = repmat(x(I3a)-ub(I3a),1,n).*DFx(I3a,:)+...
-      sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
+  H2(I3a,:) = bsxfun(@times,x(I3a)-ub(I3a),DFx(I3a,:))+...
+              sparse(1:length(I3a),I3a,Fx(I3a),length(I3a),n,length(I3a));
 end
 
-%H1         = sparse(1:n,1:n,Da,n,n,n)+Db(:,ones(n,1)).*DFx;
-H1         = bsxfun(@times,Db,DFx);
-H1         = spdiags(diag(H1)+Da,0,H1);
-H          = [lambda1*H1; lambda2*H2];
+H1 = bsxfun(@times,Db,DFx);
+H1 = spdiags(diag(H1)+Da,0,H1);
+
+H  = [lambda1*H1; lambda2*H2];