--- a/scripts/ChangeLog	Tue Apr 19 19:27:44 2005 +0000
+++ b/scripts/ChangeLog	Thu Apr 21 14:43:20 2005 +0000
@@ -1,3 +1,9 @@
+2005-04-21  John W. Eaton  <jwe@octave.org>
+
+	* optimization/glpk.m: Handle SENSE argument.
+
+	* optimization/qp.m, optimization/sqp.m: New files.
+
 2005-04-08  John W. Eaton  <jwe@octave.org>
 
 	* Makefile.in (clean, distclean, maintainer-clean):

--- a/scripts/optimization/glpk.m	Tue Apr 19 19:27:44 2005 +0000
+++ b/scripts/optimization/glpk.m	Thu Apr 21 14:43:20 2005 +0000
@@ -232,7 +232,7 @@
 ## percentage of @code{tolbnd} or @code{toldj}).
 ## 
 ## @item tolbnd (@code{LPX_K_TOLBND}, default: 10e-7)
-## Relative tolerance used to check ifthe current basic solution is primal
+## Relative tolerance used to check if the current basic solution is primal
 ## feasible.  It is not recommended that you change this parameter unless you
 ## have a detailed understanding of its purpose.
 ## 
@@ -246,7 +246,7 @@
 ## simplex table.  It is not recommended that you change this parameter unless you
 ## have a detailed understanding of its purpose.
 ## 
-## @item objll (@code{LPX_K_OBJLL)}, default: -DBL_MAX	
+## @item objll (@code{LPX_K_OBJLL}, default: -DBL_MAX)
 ## Lower limit of the objective function.  If on the phase II the objective
 ## function reaches this limit and continues decreasing, the solver stops
 ## the search.  This parameter is used in the dual simplex method only.
@@ -268,7 +268,7 @@
 ## output. Non-positive value means no delay.
 ## 
 ## @item tolint (@code{LPX_K_TOLINT}, default: 10e-5)
-## Relative tolerance used to check ifthe current basic solution is integer
+## Relative tolerance used to check if the current basic solution is integer
 ## feasible.  It is not recommended that you change this parameter unless
 ## you have a detailed understanding of its purpose.
 ## 
@@ -477,7 +477,7 @@
   ## 7) Vector with the type of variables
 
   if (nargin > 6)
-    if isempty (vartype)
+    if (isempty (vartype))
       vartype = repmat ("C", nx, 1);
     elseif (! ischar (vartype) || all (size (vartype) > 1)
 	    || length (vartype) != nx)
@@ -492,7 +492,23 @@
     vartype = repmat ("C", nx, 1);
   endif
 
-  ## 8) Parameters vector
+  ## 8) Sense of optimization
+
+  if (nargin > 7)
+    if (isempty (sense))
+      sense = 1;
+    elseif (ischar (sense) || all (size (sense) > 1) || ! isreal (sense))
+      error ("SENSE must be an integer value");
+    elseif (sense >= 0)
+      sense = 1;
+    else
+      sense = -1;
+    endif
+  else
+    sense = 1;
+  endif
+
+  ## 9) Parameters vector
 
   if (nargin > 8)
     if (! isstruct (param))

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/optimization/qp.m	Thu Apr 21 14:43:20 2005 +0000
@@ -0,0 +1,280 @@
+## Copyright (C) 2000, 2001, 2004, 2005 Gabriele Pannocchia.
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+## 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{lambda}] =} qp (@var{x0}, @var{H}, @var{q}, @var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb}, @var{Ain}, @var{A_ub})
+## Solve the quadratic program
+## @ifinfo
+##
+## @example
+##      min 0.5 x'*H*x + x'*q
+##       x
+## @end example
+##
+## @end ifinfo
+## @iftex
+## @tex
+## @end tex
+## @end iftex
+## subject to
+## @ifinfo
+##
+## @example
+##      A*x = b
+##      lb <= x <= ub
+##      A_lb <= Ain*x <= A_ub
+## @end example
+## @end ifinfo
+## @iftex
+## @tex
+## @end tex
+## @end iftex
+##
+## @noindent
+## using a null-space active-set method.
+##
+## Any bound (@var{A}, @var{b}, @var{lb}, @var{ub}, @var{A_lb},
+## @var{A_ub}) may be set to the empty matrix (@code{[]}) if not
+## present.  If the initial guess is feasible the algorithm is faster.
+##
+## The value @var{info} is a structure with the following fields:
+## @table @code
+## @item solveiter
+## The number of iterations required to find the solution.
+## @item info
+## An integer indicating the status of the solution, as follows:
+## @table @asis
+## @item 0
+## The problem is feasible and convex.  Global solution found.
+## @item 1
+## The problem is not convex.  Local solution found.
+## @item 2
+## The problem is not convex and unbounded.
+## @item 3
+## Maximum number of iterations reached.
+## @item 6
+## The problem is infeasible.
+## @end table
+## @end table
+## @end deftypefn
+
+function [x, obj, INFO, lambda] = qp (x0, H, q, A, b, lb, ub, A_lb, A_in, A_ub)
+
+  if (nargin == 5 || nargin == 7 || nargin == 10)
+
+    ## Checking the quadratic penalty
+    n = issquare (H);
+    if (n == 0)
+      error ("qp: quadratic penalty matrix not square");
+    endif
+
+    n1 = issymmetric (H);
+    if (n1 == 0)
+      ## warning ("qp: quadratic penalty matrix not symmetric");
+      H = (H + H')/2;
+    endif
+
+    ## Checking the initial guess (if empty it is resized to the
+    ## right dimension and filled with 0)
+    if (isempty (x0))
+      x0 = zeros (n, 1);
+    elseif (length (x0) != n)
+      error ("qp: the initial guess has incorrect length");
+    endif
+
+    ## Linear penalty.
+    if (length (q) != n)
+      error ("qp: the linear term has incorrect length");
+    endif
+
+    ## Equality constraint matrices
+    if (isempty (A) || isempty(b))
+      n_eq = 0;
+      A = zeros (n_eq, n);
+      b = zeros (n_eq, 1);
+    else
+      [n_eq, n1] = size (A);
+      if (n1 != n)
+	error ("qp: equality constraint matrix has incorrect column dimension");
+      endif
+      if (length (b) != n_eq)
+	error ("qp: equality constraint matrix and vector have inconsistent dimension");
+      endif
+    endif
+
+    ## Bound constraints
+    Ain = zeros (0, n);
+    bin = zeros (0, 1);
+    n_in = 0;
+    if (nargin > 5)
+      if (! isempty (lb))
+	if (length(lb) != n)
+	  error ("qp: lower bound has incorrect length");
+	else
+	  Ain = [Ain; eye(n)];
+	  bin = [bin; lb];
+	endif
+      endif
+
+      if (! isempty (ub))
+	if (length (ub) != n)
+	  error ("qp: upper bound has incorrect length");
+	else
+	  Ain = [Ain; -eye(n)];
+	  bin = [bin; -ub];
+	endif
+      endif
+    endif
+
+    ## Inequality constraints
+    if (nargin > 7)
+      [dimA_in, n1] = size (A_in);
+      if (n1 != n)
+	error ("qp: inequality constraint matrix has incorrect column dimension");
+      else
+	if (! isempty (A_lb))
+	  if (length (A_lb) != dimA_in)
+	    error ("qp: inequality constraint matrix and lower bound vector inconsistent");
+	  else
+	    Ain = [Ain; A_in];
+	    bin = [bin; A_lb];
+	  endif
+	endif
+	if (! isempty (A_ub))
+	  if (length (A_ub) != dimA_in)
+	    error ("qp: inequality constraint matrix and upper bound vector inconsistent");
+	  else
+	    Ain = [Ain; -A_in];
+	    bin = [bin; -A_ub];
+	  endif
+	endif
+      endif
+    endif
+    n_in = length (bin);
+
+    ## Now we should have the following QP:
+    ##
+    ##   min_x  0.5*x'*H*x + x'*q
+    ##   s.t.   A*x = b
+    ##          Ain*x >= bin
+
+    ## Discard inequality constraints that have -Inf bounds since those
+    ## will never be active.
+    idx = isinf (bin) & bin < 0;
+    bin(idx) = [];
+    Ain(idx,:) = [];
+
+    ## Check if the initial guess is feasible.
+    rtol = sqrt (eps);
+
+    eq_infeasible = (n_eq > 0 && norm (A*x0-b) > rtol*(1+norm (b)));
+    in_infeasible = (n_in > 0 && any (Ain*x0-bin < -rtol*(1+norm (bin))));
+
+    info = 0;
+    if (eq_infeasible || in_infeasible)
+      ## The initial guess is not feasible.
+      ## First define xbar that is feasible with respect to the equality
+      ## constraints.
+      if (eq_infeasible)
+	if (rank (A) < n_eq)
+	  error ("qp: equality constraint matrix must be full row rank")
+	endif
+	xbar = pinv (A) * b;
+      else
+	xbar = x0;
+      endif
+
+      ## Check if xbar is feasible with respect to the inequality
+      ## constraints also.
+      if (n_in > 0)
+	res = Ain * xbar - bin;
+	if (any (res < -rtol * (1 + norm (bin))))
+	  ## xbar is not feasible with respect to the inequality
+	  ## constraints.  Compute a step in the null space of the
+	  ## equality constraints, by solving a QP.  If the slack is
+	  ## small, we have a feasible initial guess.  Otherwise, the
+	  ## problem is infeasible.
+	  if (n_eq > 0)
+	    Z = null (A);
+	    if (isempty (Z))
+	      ## The problem is infeasible because A is square and full
+	      ## rank, but xbar is not feasible.
+	      info = 6;
+	    endif
+	  endif
+
+	  if (info != 6)
+            ## Solve an LP with additional slack variables to find
+	    ## a feasible starting point.
+	    gamma = eye (n_in);
+	    if (n_eq > 0)
+	      Atmp = [Ain*Z, gamma];
+	      btmp = -res;
+	    else
+	      Atmp = [Ain, gamma];
+	      btmp = bin;
+	    endif
+	    ctmp = [zeros(n-n_eq, 1); ones(n_in, 1)];
+	    lb = [-Inf*ones(n-n_eq,1); zeros(n_in,1)];
+	    ub = [];
+	    ctype = repmat ("L", n_in, 1);
+	    [P, dummy, status] = glpk (ctmp, Atmp, btmp, lb, ub, ctype);
+
+	    if ((status == 180 || status == 181 || status == 151)
+		&& all (abs (P(n-n_eq+1:end)) < rtol * (1 + norm (btmp))))
+	      ## We found a feasible starting point
+	      if (n_eq > 0)
+		x0 = xbar + Z*P(1:n-n_eq)
+	      else
+		x0 = P(1:n);
+              endif
+	    else
+	      ## The problem is infeasible
+	      info = 6;
+	    endif
+	  endif
+	else
+	  ## xbar is feasible.  We use it a starting point.
+	  x0 = xbar;
+	endif
+      else
+	## xbar is feasible.  We use it a starting point.
+	x0 = xbar;
+      endif
+    endif
+
+    if (info == 0)
+      ## The initial (or computed) guess is feasible.
+      ## We call the solver.
+      maxit = 100;
+      [x, lambda, info, iter] = __qp__ (x0, H, q, A, b, Ain, bin, maxit);
+    else
+      iter = 0;
+      x = x0;
+      lambda = [];
+    endif
+    obj = 0.5 * x' * H * x + q' * x;
+    INFO.solveiter = iter;
+    INFO.info = info;
+
+  else
+    usage ("[x, obj, info, lambda] = qp (x0, H, q, A, b, lb, ub, A_lb, Ain, A_ub)");
+  endif
+
+endfunction

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/scripts/optimization/sqp.m	Thu Apr 21 14:43:20 2005 +0000
@@ -0,0 +1,603 @@
+## Copyright (C) 2005 John W. Eaton
+##
+## This file is part of Octave.
+##
+## Octave is free software; you can redistribute it and/or modify it
+## under the terms of the GNU General Public License as published by
+## the Free Software Foundation; either version 2, or (at your option)
+## any later version.
+##
+## Octave is distributed in the hope that it will be useful, but
+## WITHOUT ANY WARRANTY; without even the implied warranty of
+## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+## General Public License for more details.
+##
+## You should have received a copy of the GNU General Public License
+## along with Octave; see the file COPYING.  If not, write to the Free
+## Software Foundation, 59 Temple Place - Suite 330, Boston, MA
+## 02111-1307, USA.
+
+## -*- texinfo -*-
+## @deftypefn {Function File} {[@var{x}, @var{obj}, @var{info}, @var{iter}, @var{nf}, @var{lambda}] =} sqp (@var{x}, @var{phi}, @var{g}, @var{h})
+## Solve the nonlinear program
+## @ifinfo
+##
+## @example
+##      min phi (x)
+##       x
+## @end example
+##
+## @end ifinfo
+## @iftex
+## @tex
+## @end tex
+## @end iftex
+## subject to
+## @ifinfo
+##
+## @example
+##      g(x)  = 0
+##      h(x) >= 0
+## @end example
+## @end ifinfo
+## @iftex
+## @tex
+## @end tex
+## @end iftex
+##
+## @noindent
+## using a successive quadratic programming method.
+##
+## The first argument is the initial guess for the vector @var{x}.
+##
+## The second argument is a function handle pointing to the ojective
+## function.  The objective function must be of the form
+##
+## @example
+##      y = phi (x)
+## @end example
+##
+## @noindent
+## in which @var{x} is a vector and @var{y} is a scalar.
+##
+## The second argument may also be a 2- or 3-element cell array of
+## function handles.  The first element should point to the objective
+## function, the second should point to a function that computes the
+## gradient of the objective function, and the third should point to a
+## function to compute the hessian of the objective function.  If the
+## gradient function is not supplied, the gradient is computed by finite
+## differences.  If the hessian function is not supplied, a BFGS update
+## formula is used to approximate the hessian.
+##
+## If supplied, the gradient function must be of the form
+##
+## @example
+##      g = gradient (x)
+## @end example
+##
+## @noindent
+## in which @var{x} is a vector and @var{g} is a vector.
+##
+## If supplied, the hessian function must be of the form
+##
+## @example
+##      h = hessian (x)
+## @end example
+##
+## @noindent
+## in which @var{x} is a vector and @var{h} is a matrix.
+##
+## The third and fourth arguments are function handles pointing to
+## functions that compute the equality constraints and the inequality
+## constraints, respectively.
+##
+## If your problem does not have equality (or inequality) constraints,
+## you may pass an empty matrix for @var{cef} (or @var{cif}).
+##
+## If supplied, the equality and inequality constraint functions must be
+## of the form
+##
+## @example
+##      r = f (x)
+## @end example
+##
+## @noindent
+## in which @var{x} is a vector and @var{r} is a vector.
+## 
+## The third and fourth arguments may also be 2-element cell arrays of
+## function handles.  The first element should point to the constraint
+## function and the second should point to a function that computes the
+## gradient of the constraint function:
+##
+## @example
+##                 [ d f(x)   d f(x)        d f(x) ]
+##     transpose ( [ ------   -----   ...   ------ ] )
+##                 [  dx_1     dx_2          dx_N  ]
+## @end example
+##
+## Here is an example of calling @code{sqp}:
+##
+## @example
+##      function r = g (x)
+##        r = [ sumsq(x)-10; x(2)*x(3)-5*x(4)*x(5); x(1)^3+x(2)^3+1];
+##      endfunction
+##
+##      function obj = phi (x)
+##        obj = exp(prod(x)) - 0.5*(x(1)^3+x(2)^3+1)^2;
+##      endfunction
+##
+##      x0 = [-1.8; 1.7; 1.9; -0.8; -0.8];
+##
+##      [x, obj, info, iter, nf, lambda] = sqp (x0, @@phi, @@g, [])
+##
+##      x =
+##      
+##        -1.71714
+##         1.59571
+##         1.82725
+##        -0.76364
+##        -0.76364
+##      
+##      obj = 0.053950
+##      info = 101
+##      iter = 8
+##      nf = 10
+##      lambda =
+##      
+##        -0.0401627
+##         0.0379578
+##        -0.0052227
+## @end example
+##
+## The value returned in @var{info} may be one of the following:
+## @table @asis
+## @item 101
+## The algorithm terminated because the norm of the last step was less
+## than @code{tol * norm (x))} (the value of tol is currently fixed at
+## @code{sqrt (eps)}---edit @file{sqp.m} to modify this value.
+## @item 102
+## The BFGS update failed.
+## @item 103
+## The maximum number of iterations was reached (the maximum number of
+## allowed iterations is currently fixed at 100---edit @file{sqp.m} to
+## increase this value).
+## @end table
+## @end deftypefn
+## @seealso{qp}
+
+function [x, obj, info, iter, nf, lambda] = sqp (x, objf, cef, cif)
+
+  global nfun;
+  global __sqp_obj_fun__;
+  global __sqp_ce_fun__;
+  global __sqp_ci_fun__;
+
+  if (nargin >= 2 && nargin <= 4)
+
+    ## Choose an initial NxN symmetric positive definite Hessan
+    ## approximation B.
+
+    n = length (x);
+
+    B = eye (n, n);
+
+    ## Evaluate objective function, constraints, and gradients at initial
+    ## value of x.
+    ##
+    ## obj_fun
+    ## obj_grad
+    ## ce_fun  -- equality constraint functions
+    ## ci_fun  -- inequality constraint functions
+    ## A == [grad_{x_1} cx_fun, grad_{x_2} cx_fun, ..., grad_{x_n} cx_fun]^T
+
+    obj_grd = @fd_obj_grd;
+    have_hess = 0;
+    if (iscell (objf))
+      if (length (objf) > 0)
+	__sqp_obj_fun__ = obj_fun = objf{1};
+	if (length (objf) > 1)
+	  obj_grd = objf{2};
+	  if (length (objf) > 2)
+	    obj_hess = objf{3};
+	    have_hess = 1;
+	    B = feval (obj_hess, x);
+	  endif
+	endif
+      else
+	error ("sqp: invalid objective function");
+      endif
+    else
+      __sqp_obj_fun__ = obj_fun = objf;
+    endif
+
+    ce_fun = @empty_cf;
+    ce_grd = @empty_jac;
+    if (nargin > 2)
+      ce_grd = @fd_ce_jac;
+      if (iscell (cef))
+	if (length (cef) > 0)
+	  __sqp_ce_fun__ = ce_fun = cef{1};
+	  if (length (cef) > 1)
+	    ce_grd = cef{2};
+	  endif
+	else
+	  error ("sqp: invalid equality constraint function");
+	endif
+      elseif (! isempty (cef))
+	ce_fun = cef;
+      endif
+    endif
+    __sqp_ce_fun__ = ce_fun;
+
+    ci_fun = @empty_cf;
+    ci_grd = @empty_jac;
+    if (nargin > 3)
+      ci_grd = @fd_ci_jac;
+      if (iscell (cif))
+	if (length (cif) > 0)
+	  __sqp_ci_fun__ = ci_fun = cif{1};
+	  if (length (cif) > 1)
+	    ci_grd = cif{2};
+	  endif
+	else
+	  error ("sqp: invalid equality constraint function");
+	endif
+      elseif (! isempty (cif))
+	ci_fun = cif;
+      endif
+    endif
+    __sqp_ci_fun__ = ci_fun;
+
+    iter_max = 100;
+
+    iter = 0;
+
+    obj = feval (obj_fun, x);
+    nfun = 1;
+
+    c = feval (obj_grd, x);
+
+    ce = feval (ce_fun, x);
+    F = feval (ce_grd, x);
+
+    ci = feval (ci_fun, x);
+    C = feval (ci_grd, x);
+
+    A = [F; C];
+
+    ## Choose an initial lambda (x is provided by the caller).
+
+    lambda = 100 * ones (rows (A), 1);
+
+    qp_iter = 1;
+    alpha = 1;
+
+    ## report ();
+
+    ## report (iter, qp_iter, alpha, nfun, obj);
+
+    while (++iter < iter_max)
+
+      ## Check convergence.  This is just a simple check on the first
+      ## order necessary conditions.
+
+      ## IDX is the indices of the active inequality constraints.
+
+      nr_f = rows (F);
+
+      lambda_e = lambda((1:nr_f)');
+      lambda_i = lambda((nr_f+1:end)');
+
+      con = [ce; ci];
+
+      t0 = norm (c - A' * lambda);
+      t1 = norm (ce);
+      t2 = all (ci >= 0);
+      t3 = all (lambda_i >= 0);
+      t4 = norm (lambda .* con);
+
+      tol = sqrt (eps);
+
+      if (t2 && t3 && max ([t0; t1; t4]) < tol)
+	break;
+      endif
+
+      ## Compute search direction p by solving QP.
+
+      g = -ce;
+      d = -ci;
+
+      ## Discard inequality constraints that have -Inf bounds since those
+      ## will never be active.
+      idx = isinf (d) & d < 0;
+      d(idx) = [];
+      C(idx,:) = [];
+
+      [p, obj_qp, INFO, lambda] = qp (x, B, c, F, g, [], [], d, C,
+				      Inf * ones (size (d)));
+
+      info = INFO.info;
+
+      ## Check QP solution and attempt to recover if it has failed.
+
+      ## Choose mu such that p is a descent direction for the chosen
+      ## merit function phi.
+
+      [x_new, alpha, obj_new] = linesearch_L1 (x, p, obj_fun, obj_grd,
+					       ce_fun, ci_fun, lambda, obj);
+
+      ## Evaluate objective function, constraints, and gradients at
+      ## x_new.
+
+      c_new = feval (obj_grd, x_new);
+
+      ce_new = feval (ce_fun, x_new);
+      F_new = feval (ce_grd, x_new);
+
+      ci_new = feval (ci_fun, x_new);
+      C_new = feval (ci_grd, x_new);
+
+      A_new = [F_new; C_new];
+
+      ## Set
+      ##
+      ## s = alpha * p
+      ## y = grad_x L (x_new, lambda) - grad_x L (x, lambda})
+
+      y = c_new - c;
+
+      if (! isempty (A))
+	t = ((A_new - A)'*lambda);
+	y -= t;
+      endif
+
+      delx = x_new - x;
+
+      if (norm (delx) < tol * norm (x))
+	info = 101;
+	break;
+      endif
+
+      if (have_hess)
+
+	B = feval (obj_hess, x);
+
+      else
+
+	## Update B using a quasi-Newton formula.
+
+	delxt = delx';
+
+	## Damped BFGS.  Or maybe we would actually want to use the Hessian
+	## of the Lagrangian, computed directly.
+
+	d1 = delxt*B*delx;
+
+	t1 = 0.2 * d1;
+	t2 = delxt*y;
+
+	if (t2 < t1)
+	  theta = 0.8*d1/(d1 - t2);
+	else
+	  theta = 1;
+	endif
+
+	r = theta*y + (1-theta)*B*delx;
+
+	d2 = delxt*r;
+
+	if (d1 == 0 || d2 == 0)
+	  info = 102;
+	  break;
+	endif
+
+	B = B - B*delx*delxt*B/d1 + r*r'/d2;
+
+      endif
+
+      x = x_new;
+
+      obj = obj_new;
+
+      c = c_new;
+
+      ce = ce_new;
+      F = F_new;
+
+      ci = ci_new;
+      C = C_new;
+
+      A = A_new;
+
+      ## report (iter, qp_iter, alpha, nfun, obj);
+
+    endwhile
+
+    if (iter >= iter_max)
+      info = 103;
+    endif
+
+    nf = nfun;
+
+  else
+
+    usage ("[x, obj, info, iter, nf, lambda] = sqp (x, objf, cef, cif)");
+
+  endif
+
+### endfunction
+
+
+function [merit, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu)
+
+  global nfun;
+
+  ce = feval (ce_fun, x);
+  ci = feval (ci_fun, x);
+
+  idx = ci < 0;
+
+  con = [ce; ci(idx)];
+
+  if (isempty (obj))
+    obj = feval (obj_fun, x);
+    nfun++;
+  endif
+
+  merit = obj;
+  t = norm (con, 1) / mu;
+
+  if (! isempty (t))
+    merit += t;
+  endif
+
+### endfunction
+
+
+function [x_new, alpha, obj] = linesearch_L1 (x, p, obj_fun, obj_grd,
+					      ce_fun, ci_fun, lambda, obj)
+
+  ## Choose parameters
+  ##
+  ## eta in the range (0, 0.5)
+  ## tau in the range (0, 1)
+
+  eta = 0.25;
+  tau = 0.5;
+
+  delta_bar = sqrt (eps);
+
+  if (isempty (lambda))
+    mu = 1 / delta_bar;
+  else
+    mu = 1 / (norm (lambda, Inf) + delta_bar);
+  endif
+
+  alpha = 1;
+
+  c = feval (obj_grd, x);
+  ce = feval (ce_fun, x);
+
+  [phi_x_mu, obj] = phi_L1 (obj, obj_fun, ce_fun, ci_fun, x, mu);
+
+  D_phi_x_mu = c' * p;
+  d = feval (ci_fun, x);
+  ## only those elements of d corresponding
+  ## to violated constraints should be included.
+  idx = d < 0;
+  t = - norm ([ce; d(idx)], 1) / mu;
+  if (! isempty (t))
+    D_phi_x_mu += t;
+  endif
+
+  while (1)
+    [p1, obj] = phi_L1 ([], obj_fun, ce_fun, ci_fun, x+alpha*p, mu);
+    p2 = phi_x_mu+eta*alpha*D_phi_x_mu;
+    if (p1 > p2)
+      ## Reset alpha = tau_alpha * alpha for some tau_alpha in the
+      ## range (0, tau).
+      tau_alpha = 0.9 * tau;  ## ??
+      alpha = tau_alpha * alpha;
+    else
+      break;
+    endif
+  endwhile
+
+  ## Set x_new = x + alpha * p;
+
+  x_new = x + alpha * p;
+
+### endfunction
+
+
+function report (iter, qp_iter, alpha, nfun, obj)
+
+  if (nargin == 0)
+    printf ("  Itn ItQP     Step  Nfun     Objective\n");
+  else
+    printf ("%5d %4d %8.1g %5d %13.6e\n", iter, qp_iter, alpha, nfun, obj);
+  endif
+
+### endfunction
+
+
+function grd = fdgrd (f, x)
+
+  if (! isempty (f))
+    y0 = feval (f, x);
+    nx = length (x);
+    grd = zeros (nx, 1);
+    deltax = sqrt (eps);
+    for i = 1:nx
+      t = x(i);
+      x(i) += deltax;
+      grd(i) = (feval (f, x) - y0) / deltax;
+      x(i) = t;
+    endfor
+  else
+    grd = zeros (0, 1);
+  endif
+
+### endfunction
+
+
+function jac = fdjac (f, x)
+
+  if (! isempty (f))
+    y0 = feval (f, x);
+    nf = length (y0);
+    nx = length (x);
+    jac = zeros (nf, nx);
+    deltax = sqrt (eps);
+    for i = 1:nx
+      t = x(i);
+      x(i) += deltax;
+      jac(:,i) = (feval (f, x) - y0) / deltax;
+      x(i) = t;
+    endfor
+  else
+    jac = zeros  (0, nx);
+  endif
+
+### endfunction
+
+
+function grd = fd_obj_grd (x)
+
+  global __sqp_obj_fun__;
+
+  grd = fdgrd (__sqp_obj_fun__, x);
+
+### endfunction
+
+
+function res = empty_cf (x)
+
+  res = zeros (0, 1);
+
+### endfunction
+
+
+function res = empty_jac (x)
+
+  res = zeros (0, length (x));
+
+### endfunction
+
+
+function jac = fd_ce_jac (x)
+
+  global __sqp_ce_fun__;
+
+  jac = fdjac (__sqp_ce_fun__, x);
+
+### endfunction
+
+
+function jac = fd_ci_jac (x)
+
+  global __sqp_ci_fun__;
+
+  jac = fdjac (__sqp_ci_fun__, x);
+
+### endfunction