updating lmmcp.m from RECS

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,6 +123,7 @@ else
   options = catstruct(defaultopt,options);
 end
 
+warning('off','MATLAB:rankDeficientMatrix')
 
 switch options.Display
   case {'off','none'}
@@ -136,7 +138,7 @@ end
 
 % parameter settings
 eps1 = options.epsilon1;
-eps2 = options.epsilon2;
+eps2 = 0.5*options.TolFun^2;
 null = options.null;
 Big  = options.Big;
 
@@ -219,7 +221,7 @@ if watchdog==1
   DPhibest     = DPhix;
   DPsibest     = DPsix;
   normDPsibest = normDPsix;
-end;
+end
 
 % initial output
 if verbosity > 1
@@ -229,9 +231,7 @@ 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
 
@@ -247,17 +247,17 @@ if preprocess==1
     % the condition estimator for large-scale problems, although this
     % may cause numerical problems in some examples
 
-      i=0;
+    i  = false;
     mu = 0;
     if n<100
-         i=1;
+      i = true;
       mu = 1e-16;
       if condest(DPhix'*DPhix)>1e25
         mu = 1e-6/(k+1);
       end
     end
-      if i==1
-         pLM= [ DPhix ; sqrt(mu)*speye(n)]\[-Phix;sparse(n,1)];
+    if i
+      pLM =  [DPhix; sqrt(mu)*speye(n)]\[-Phix; zeros(n,1)];
     else
       pLM = -DPhix\Phix;
     end
@@ -326,9 +326,7 @@ elseif preprocess==1 && Psix>=eps2
   end
 end
 
-%
-%   Main algorithm
-%
+%%   Main algorithm
 
 if verbosity > 1
   disp('************************** Main program ****************************')
@@ -340,9 +338,9 @@ while (k < kmax) && (Psix > eps2)
   % the condition estimator for large-scale problems, although this
   % may cause numerical problems in some examples
 
-   i=0;
+  i = false;
   if n<100
-      i=1;
+    i  = true;
     mu = 1e-16;
     if condest(DPhix'*DPhix)>1e25
       mu = 1e-1/(k+1);
@@ -351,8 +349,8 @@ while (k < kmax) && (Psix > eps2)
   
   % compute a Levenberg-Marquard direction
 
-   if i==1
-      d= [ DPhix ; sqrt(mu)*eye(n)]\[-Phix;zeros(n,1)];
+  if i
+    d = [DPhix; sqrt(mu)*speye(n)]\[-Phix; zeros(n,1)];
   else
     d = -DPhix\Phix;
   end
@@ -396,7 +394,7 @@ while (k < kmax) && (Psix > eps2)
   else
     aux(mod(k_main,m)+1) = Psix;
     MaxPsi               = max(aux);
-   end;
+  end
   
   % updatings for the watchdog strategy
   if watchdog ==1
@@ -416,16 +414,16 @@ while (k < kmax) && (Psix > eps2)
       DPsix=DPsibest;
       normDPsix=normDPsibest;
       MaxPsi=Psix;
-      end;
-   end;
+    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;
+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);
@@ -527,8 +526,8 @@ end
 
 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;
@@ -550,8 +549,8 @@ end
 
 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;
@@ -595,21 +594,21 @@ Db(I3)     = bi(I3).*di(I3);
 
 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,:)+...
+  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);
 if any(I3a)
-  H2(I3a,:)  = repmat(x(I3a)-lb(I3a),1,n).*DFx(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);
 if any(I3a)
-  H2(I3a,:)  = repmat(x(I3a)-ub(I3a),1,n).*DFx(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];