Shell the following script in an empty directory to create the mfiles
associated with the paper  "On the resolution of monotone
complementarity problems" by Geiger and Kanzow
------------
sed > DJosephy.m <<'#####' 's/^-//'
-function DFx = DJosephy (n,x)
-%
-%   Computes the Jacobian matrix F'(x) for the Josephy example.
-%
-DFx=zeros(n,n);
-DFx=[6*x(1)+2*x(2) 2*x(1)+4*x(2) 1 3
-     4*x(1)+1 2*x(2) 3 2
-     6*x(1)+x(2) x(1)+4*x(2) 2 3
-     2*x(1) 6*x(2) 2 3];
-return
#####
sed > DKojima.m <<'#####' 's/^-//'
-function DFx = DKojima (n,x)
-%
-%   Computes the Jacobian matrix F'(x) for the Kojima example.
-%
-DFx=zeros(n,n);
-DFx=[6*x(1)+2*x(2) 2*x(1)+4*x(2) 1 3
-     4*x(1)+1 2*x(2) 10 2
-     6*x(1)+x(2) x(1)+4*x(2) 2 9
-     2*x(1) 6*x(2) 2 3];
-return
#####
sed > DMathies.m <<'#####' 's/^-//'
-function DFx = DMathies (n,x)
-%
-%   Computes the Jacobian matrix for the Mathiesen example
-%   (different parameter vectors).
-%
-alpha=0.75;
-b2=1;
-b3=0.5;
-alpha=0.75;
-b2=1;
-b3=2;
-alpha=0.9;
-b2=5;
-b3=3;
-DFx=zeros(n,n);
-DFx=[0 -1 1 1
-     1 alpha*(b2*x(3)+b3*x(4))/(x(2)^2) -alpha*b2/x(2) -alpha*b3/x(2)
-     -1 0 (1-alpha)*b3*x(4)/(x(3)^2) (alpha-1)*b3/x(3)
-     -1 0 0 0];
-return
#####
sed > DNCPHS34.m <<'#####' 's/^-//'
-function DFx = DNCPHS34 (n,x)
-%
-%   Computes the Jacobian of example HS 34, formulated as an NCP.
-%
-DFx=zeros(n,n);
-DFx=[x(4)*exp(x(1)) 0 0 exp(x(1)) 0 1 0 0
-     0 x(5)*exp(x(2)) 0 -1 exp(x(2)) 0 1 0
-     0 0 0 0 -1 0 0 1
-     -exp(x(1)) 1 0 0 0 0 0 0
-     0 -exp(x(2)) 1 0 0 0 0 0
-     -1 0 0 0 0 0 0 0
-     0 -1 0 0 0 0 0 0
-     0 0 -1 0 0 0 0 0];
-return
#####
sed > DNCPHS35.m <<'#####' 's/^-//'
-function DFx = DNCPHS35 (n,x)
-%
-%   Computes the Jacobian of example HS 35, formulated as an NCP.
-%
-DFx=zeros(n,n);
-DFx=[4 2 2 1
-     2 4 0 1
-     2 0 2 2
-     -1 -1 -2 0];
-return
#####
sed > DNCPHS66.m <<'#####' 's/^-//'
-function DFx = DNCPHS66 (n,x)
-%
-%   Computes the Jacobian for example HS 66, formulated as an NCP.
-%
-DFx=zeros(n,n);
-DFx=[x(4)*exp(x(1)) 0 0 exp(x(1)) 0 1 0 0
-     0 x(5)*exp(x(2)) 0 -1 exp(x(2)) 0 1 0
-     0 0 0 0 -1 0 0 1
-     -exp(x(1)) 1 0 0 0 0 0 0
-     0 -exp(x(2)) 1 0 0 0 0 0
-     -1 0 0 0 0 0 0 0
-     0 -1 0 0 0 0 0 0
-     0 0 -1 0 0 0 0 0];
-return
#####
sed > DNCPHS76.m <<'#####' 's/^-//'
-function DFx = DNCPHS76 (n,x)
-%
-%   Computes the Jacobian for example 76, formulated as an NCP.
-%
-DFx=zeros(n,n);
-DFx=[2 0 -1 0 1 3 0
-     0 1 0 0 2 1 -1
-     -1 0 2 1 1 2 -4
-     0 0 1 1 1 -1 0
-     -1 -2 -1 -1 0 0 0
-     -3 -1 -2 1 0 0 0
-     0 1 4 0 0 0 0];
-return
#####
sed > DaPhiF2.m <<'#####' 's/^-//'
-function DaPhix = DaPhiF2 (x,w)
-%
-%   Computes (mainly) the partial derivative of the squared Fischer
-%   merit function with respect to the first argument.
-%
-n=length(x);
-null=1e-8;
-DaPhix=(sqrt(x.^2+w.^2)-x-w).*(x./max(null*ones(n,1),sqrt(x.^2+w.^2))-ones(n,1));
-return
#####
sed > DaPhiMS.m <<'#####' 's/^-//'
-function DaPhix = DaPhiMS (x,w)
-%
-%   Computes the partial derivative of the implicit Lagrangian merit
-%   function by Mangasarian and Solodov with respect to the first
-%   argument.
-%
-n=length(x);
-alpha=1.1;
-DaPhix=w+(max(zeros(n,1),x-alpha*w)-x-alpha*max(zeros(n,1),w-alpha*x))/alpha;
-return
#####
sed > DbPhiF2.m <<'#####' 's/^-//'
-function DbPhix = DbPhiF2 (x,w)
-%
-%   Computes (mainly) the partial derivative of the squared Fischer
-%   merit function with respect to the second argument.
-%
-n=length(x);
-null=1e-8;
-DbPhix=(sqrt(x.^2+w.^2)-x-w).*(w./max(null*ones(n,1),sqrt(x.^2+w.^2))-ones(n,1));
-return
#####
sed > DbPhiMS.m <<'#####' 's/^-//'
-function DbPhix = DbPhiMS (x,w)
-%
-%   Computes the partial derivative of the implicit Lagrangian merit
-%   function by Mangasarian and Solodov with respect to the second
-%   argument.
-%
-n=length(x);
-alpha=1.1;
-DbPhix=x+(-alpha*max(zeros(n,1),x-alpha*w)+max(zeros(n,1),w-alpha*x)-w)/alpha;
-return
#####
sed > Josephy.m <<'#####' 's/^-//'
-function Fx = Josephy (n,x)
-%
-%   Computes the function value F(x) for the Josephy example.
-%
-Fx=zeros(4,1);
-Fx(1)=3*x(1)^2+2*x(1)*x(2)+2*x(2)^2+x(3)+3*x(4)-6;
-Fx(2)=2*x(1)^2+x(1)+x(2)^2+3*x(3)+2*x(4)-2;
-Fx(3)=3*x(1)^2+x(1)*x(2)+2*x(2)^2+2*x(3)+3*x(4)-1;
-Fx(4)=x(1)^2+3*x(2)^2+2*x(3)+3*x(4)-3;
-return
#####
sed > Kojima.m <<'#####' 's/^-//'
-function Fx = Kojima (n,x)
-%
-%   Computes the function value F(x) for the Kojima example.
-%
-Fx=zeros(4,1);
-Fx(1)=3*x(1)^2+2*x(1)*x(2)+2*x(2)^2+x(3)+3*x(4)-6;
-Fx(2)=2*x(1)^2+x(1)+x(2)^2+10*x(3)+2*x(4)-2;
-Fx(3)=3*x(1)^2+x(1)*x(2)+2*x(2)^2+2*x(3)+9*x(4)-9;
-Fx(4)=x(1)^2+3*x(2)^2+2*x(3)+3*x(4)-3;
-return
#####
sed > LMCOAP.m <<'#####' 's/^-//'
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-%
-%   Programm LMCOAP.M:
-%
-%   Contains the main program (Algorithm 5.1: Limited memory
-%   BFGS method) of the following paper:
-%
-%   Carl Geiger and Christian Kanzow:
-%   "On the Resolution of Monotone Complementarity Problems",
-%   to appear in "Computational Optimization and Applications".
-%
-%
-%   The user has to specify the following parameters:
-%
-%   n   : dimension of the test problem,
-%   m   : number if vector pairs to be stored by the limited memory BFGS
-%         method (usually 3 <= m <= 7),
-%   FName: specifies the test problem,
-%   Phi  : specifies the merit function; two possible choices:
-%         Phi='PhiF2' --> squared Fischer-function
-%         Phi='PhiMS' --> implicit Lagrangian by Mangasarian and Solodov
-%   x   : starting vector.
-%
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
-n=4;
-m=5;
-FName='Josephy';
-DFName=['D',FName];
-Phi='PhiF2';
-DaPhi=['Da',Phi];
-DbPhi=['Db',Phi];
-x=[1,1,1,1]';
-Fx=feval(FName,n,x);
-DFx=feval(DFName,n,x);
-Phix=feval(Phi,x,Fx);
-DaPhix=feval(DaPhi,x,Fx);
-DbPhix=feval(DbPhi,x,Fx);
-Psix=sum(Phix);
-DPsix=DaPhix+DFx'*DbPhix;
-noPsieval=1;
-noDPsieval=1;
-normD=norm(DPsix)
-eps=1e-5;
-it=0
-itmax=500;
-format long;
-format compact;
-d=-DPsix;
-alpha=zeros(1,m);
-rho=zeros(1,m);
-SS=zeros(n,m);
-YY=zeros(n,m);
-while (it<itmax) & (normD>eps)
-   [t,xneu,Psixneu,DPsixneu,noPsieval,noDPsieval]=...
-   StrWolfe(n,x,Psix,d,DPsix,FName,DFName,Phi,DaPhi,DbPhi,noPsieval,noDPsieval);
-   t
-   s=xneu-x;
-   y=DPsixneu-DPsix;
-   ys=y'*s;
-   mdach=min(it,m-1);
-   if mdach<(m-1)
-      rho(mdach+1)=1/ys;
-      SS(:,mdach+1)=s;
-      YY(:,mdach+1)=y;
-   else
-      for i=1:mdach
-         rho(i)=rho(i+1);
-         SS(:,i)=SS(:,i+1);
-         YY(:,i)=YY(:,i+1);
-      end
-      rho(mdach+1)=1/ys;
-      SS(:,mdach+1)=s;
-      YY(:,mdach+1)=y;
-   end
-   q=DPsixneu;
-   for i=(mdach+1):-1:1
-      alpha(i)=rho(i)*(SS(:,i)'*q);
-      q=q-alpha(i)*YY(:,i);
-   end
-   % gamma=1;
-   gamma=ys/(y'*y);
-   r=gamma*q;
-   for i=1:(mdach+1)
-      beta=rho(i)*(YY(:,i)'*r);
-      r=r+(alpha(i)-beta)*SS(:,i);
-   end
-   d=-r;
-   it=it+1
-   x=xneu
-   Psix=Psixneu
-   DPsix=DPsixneu;
-   normD=norm(DPsix)
-end
-noPsieval
-noDPsieval
-Fx=feval(FName,n,x);
-res=norm(min(x,Fx))
-return
-
#####
sed > Mathies.m <<'#####' 's/^-//'
-function Fx = Mathies (n,x)
-%
-%   Computes the function value F(x) for the Mathiesen example
-%   (different parameter vectors).
-%
-alpha=0.75;
-b2=1;
-b3=0.5;
-alpha=0.75;
-b2=1;
-b3=2;
-alpha=0.9;
-b2=5;
-b3=3;
-Fx=zeros(n,1);
-Fx(1)=-x(2)+x(3)+x(4);
-Fx(2)=x(1)-alpha*(b2*x(3)+b3*x(4))/x(2);
-Fx(3)=b2-x(1)-(1-alpha)*(b2*x(3)+b3*x(4))/x(3);
-Fx(4)=b3-x(1);
-return
-
-Germany)
#####
sed > NCPHS34.m <<'#####' 's/^-//'
-function Fx = NCPHS34 (n,x)
-%
-%   Computes the function value of the example HS 34, formulated as an NCP.
-%
-Fx=zeros(n,1);
-Fx(1)=-1+x(4)*exp(x(1))+x(6);
-Fx(2)=-x(4)+x(5)*exp(x(2))+x(7);
-Fx(3)=-x(5)+x(8);
-Fx(4)=x(2)-exp(x(1));
-Fx(5)=x(3)-exp(x(2));
-Fx(6)=100-x(1);
-Fx(7)=100-x(2);
-Fx(8)=10-x(3);
-return
#####
sed > NCPHS35.m <<'#####' 's/^-//'
-function Fx = NCPHS35 (n,x)
-%
-%   Computes the function value of example HS 35, formulated as an NCP.
-%
-Fx=zeros(n,1);
-Fx(1)=-8+4*x(1)+2*x(2)+2*x(3)+x(4);
-Fx(2)=-6+4*x(2)+2*x(1)+x(4);
-Fx(3)=-4+2*x(3)+2*x(1)+2*x(4);
-Fx(4)=3-x(1)-x(2)-2*x(3);
-return
#####
sed > NCPHS66.m <<'#####' 's/^-//'
-function Fx = NCPHS66 (n,x)
-%
-%   Computes the function value of example HS 66, formulated as an NCP.
-%
-Fx=zeros(n,1);
-Fx(1)=-0.8+x(4)*exp(x(1))+x(6);
-Fx(2)=-x(4)+x(5)*exp(x(2))+x(7);
-Fx(3)=0.2-x(5)+x(8);
-Fx(4)=x(2)-exp(x(1));
-Fx(5)=x(3)-exp(x(2));
-Fx(6)=100-x(1);
-Fx(7)=100-x(2);
-Fx(8)=10-x(3);
-return
#####
sed > NCPHS76.m <<'#####' 's/^-//'
-function Fx = NCPHS76 (n,x)
-%
-%   Computes the function value of example HS 76, formulated as an NCP.
-%
-Fx=zeros(n,1);
-Fx(1)=2*x(1)-x(3)-1+x(5)+3*x(6);
-Fx(2)=x(2)-3+2*x(5)+x(6)-x(7);
-Fx(3)=2*x(3)-x(1)+x(4)+1+x(5)+2*x(6)-4*x(7);
-Fx(4)=x(4)+x(3)-1+x(5)-x(6);
-Fx(5)=5-x(1)-2*x(2)-x(3)-x(4);
-Fx(6)=4-3*x(1)-x(2)-2*x(3)+x(4);
-Fx(7)=x(2)+4*x(3)-1.5;
-return
#####
sed > PhiF2.m <<'#####' 's/^-//'
-function Phix = PhiF2 (x,w)
-%
-%   Computes the function value of the squared Fischer merit-function.
-%
-Phix=0.5*(sqrt(x.^2+w.^2)-x-w).^2;
-return
#####
sed > PhiMS.m <<'#####' 's/^-//'
-function Phix = PhiMS (x,w)
-%
-%   Computes the function value of the implicit Lagrangian merit
-%   function by Mangasarian and Solodov.
-%
-n=length(x);
-alpha=1.1;
-Phix=x.*w+((max(zeros(n,1),x-alpha*w)).^2-x.^2+...
-     (max(zeros(n,1),w-alpha*x)).^2-w.^2)/(2*alpha);
-return
#####
sed > StrWolfe.m <<'#####' 's/^-//'
-function [t,xneu,Psixneu,DPsixneu,noPsieval,noDPsieval]=...
-StrWolfe(n,x,Psix,d,DPsix,FName,DFName,Phi,DaPhi,DbPhi,noPsieval,noDPsieval)
-%
-%   This M-file computes a steplength t satisfying the strong
-%   Wolfe conditions using an algorithm described in the book
-%   by Fletcher (without using any interpolation- or extrapolation-
-%   strategies).
-%
-sigma=1e-4;
-beta=0.9;
-tau1=9;
-tau2=0.1;
-tau3=0.5;
-initslope=DPsix'*d;
-my=-Psix/(sigma*initslope);
-talt=0;
-Psixalt=Psix;
-t=min(1,my);
-xneu=x+t*d;
-stop=2;
-while stop==2
-   Fxneu=feval(FName,n,xneu);
-   Phixneu=feval(Phi,xneu,Fxneu);
-   Psixneu=sum(Phixneu);
-   noPsieval=noPsieval+1;
-   if (Psixneu>(Psix+sigma*t*initslope)) | (Psixneu>=Psixalt)
-      aa=talt;
-      bb=t;
-      Psixa=Psixalt;
-      stop=1;
-   else
-      DFxneu=feval(DFName,n,xneu);
-      DaPhixneu=feval(DaPhi,xneu,Fxneu);
-      DbPhixneu=feval(DbPhi,xneu,Fxneu);
-      DPsixneu=DaPhixneu+DFxneu'*DbPhixneu;
-      noDPsieval=noDPsieval+1;
-      newslope=DPsixneu'*d;
-      if abs(newslope)<=(-beta*initslope)
-         stop=0;
-      else
-         if newslope>=0
-            aa=t;
-            bb=talt;
-            Psixa=Psix;
-            stop=1;
-         else
-            if my<=(2*t-talt)
-               tneu=my;
-            else
-               tneu=0.5*(2*t-talt+min(my,t+tau1*(t-talt)));
-            end
-            talt=t;
-            Psixalt=Psixneu;
-            t=tneu;
-            xneu=x+t*d;
-         end
-      end
-   end
-end
-while stop==1
-   t=0.5*(aa+tau2*(bb-aa)+bb-tau3*(bb-aa));
-   xneu=x+t*d;
-   Fxneu=feval(FName,n,xneu);
-   Phixneu=feval(Phi,xneu,Fxneu);
-   Psixneu=sum(Phixneu);
-   noPsieval=noPsieval+1;
-   if (Psixneu>(Psix+sigma*t*initslope)) | (Psixneu>=Psixa)
-      bb=t;
-   else
-      DFxneu=feval(DFName,n,xneu);
-      DaPhixneu=feval(DaPhi,xneu,Fxneu);
-      DbPhixneu=feval(DbPhi,xneu,Fxneu);
-      DPsixneu=DaPhixneu+DFxneu'*DbPhixneu;
-      noDPsieval=noDPsieval+1;
-      newslope=DPsixneu'*d;
-      if abs(newslope)<=(-beta*initslope)
-         stop=0;
-      else
-         if (bb-aa)*newslope>=0
-            bb=aa;
-         end
-         aa=t;
-         Psixa=Psixneu;
-      end
-   end
-end
-return
#####
