optimization: Improved the karmarkar function 54/2854/5
Michael Baudin [Tue, 4 Jan 2011 15:37:10 +0000 (16:37 +0100)]
Change-Id: I7414ae7ac9c14242a4467ca8f25fad9008fec74b

16 files changed:
SEP/SEP_49_Linear_Programming.odt [new file with mode: 0644]
scilab/CHANGES_5.3.X
scilab/modules/optimization/help/en_US/karmarkar.xml
scilab/modules/optimization/macros/karmarkar.sci
scilab/modules/optimization/tests/nonreg_tests/bug_8719.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/nonreg_tests/bug_8719.tst [new file with mode: 0644]
scilab/modules/optimization/tests/nonreg_tests/bug_8720.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/nonreg_tests/bug_8720.tst [new file with mode: 0644]
scilab/modules/optimization/tests/nonreg_tests/bug_8726.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/nonreg_tests/bug_8726.tst [new file with mode: 0644]
scilab/modules/optimization/tests/nonreg_tests/bug_8727.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/nonreg_tests/bug_8727.tst [new file with mode: 0644]
scilab/modules/optimization/tests/nonreg_tests/bug_8775.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/nonreg_tests/bug_8775.tst [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/karmarkar.dia.ref [new file with mode: 0644]
scilab/modules/optimization/tests/unit_tests/karmarkar.tst [new file with mode: 0644]

diff --git a/SEP/SEP_49_Linear_Programming.odt b/SEP/SEP_49_Linear_Programming.odt
new file mode 100644 (file)
index 0000000..107c4ed
Binary files /dev/null and b/SEP/SEP_49_Linear_Programming.odt differ
index 5c7b22c..17a43ff 100644 (file)
@@ -196,6 +196,33 @@ Scilab:
 * Scilab startup script of the binary is now automatically synchronized with
   the one from the source tree.
 
+Optimization:
+=============
+
+* Improved the karmarkar linear optimization solver: 
+  * Let the user configure an output function.
+  * Put the number of iterations as an output argument.
+  * Put x0 as an optional argument.
+  * Manage inequality constraints.
+  * Manage bounds.
+  * Compute the exitflag.
+  * Compute the Lagrange multipliers.
+  * Provide documentation for the rtolf and gam parameters.
+  * Provide more examples.
+  * Provide documentation for the stopping rule.
+
+* bug 7193 fixed - The karmarkar help page did not document the 
+                   parameters eps, gamma, and crit.
+
+* bug 8719 fixed - The karmarkar function printed unwanted messages.
+
+* bug 8720 fixed - The karmarkar function stoped too early in the iterations.
+
+* bug 8726 fixed - The karmarkar function produced a division-by-zero error.
+
+* bug 8727 fixed - The karmarkar function requires the initial guess x0.
+
+* bug 8775 fixed - The karmarkar function did not detect unbounded problems.
 
 Bug fixes:
 ==========
index 9be4d4a..91ded5f 100644 (file)
@@ -2,6 +2,7 @@
 <!--
  * Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
  * Copyright (C) 2008 - INRIA
+ * Copyright (C) 2010 - DIGITEO - Michael Baudin
  * 
  * This file must be used under the terms of the CeCILL.
  * This source file is licensed as described in the file COPYING, which
   <refnamediv>
     <refname>karmarkar</refname>
 
-    <refpurpose>karmarkar algorithm</refpurpose>
+    <refpurpose>Solves a linear optimization problem.</refpurpose>
   </refnamediv>
 
   <refsynopsisdiv>
     <title>Calling Sequence</title>
 
-    <synopsis>[x1]=karmarkar(a,b,c,x0)</synopsis>
+    <synopsis>
+      xopt=karmarkar(Aeq,beq,c)
+      xopt=karmarkar(Aeq,beq,c,x0)
+      xopt=karmarkar(Aeq,beq,c,x0,rtolf)
+      xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam)
+      xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter)
+      xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun)
+      xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun,A,b)
+      xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun,A,b,lb)
+      xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun,A,b,lb,ub)
+      [xopt,fopt] = karmarkar(...)
+      [xopt,fopt,exitflag] = karmarkar(...)
+      [xopt,fopt,exitflag,iter] = karmarkar(...)
+      [xopt,fopt,exitflag,iter,yopt] = karmarkar(...)
+    </synopsis>
   </refsynopsisdiv>
 
   <refsection>
 
     <variablelist>
       <varlistentry>
-        <term>a</term>
+        <term>Aeq</term>
+        <listitem>
+          <para>
+            a n-by-p matrix of doubles, where <literal>n</literal> is the number
+            of linear equality constraints
+            and <literal>p</literal> is the number of unknowns,
+            the matrix in the linear equality constraints.
+          </para>
+        </listitem>
+      </varlistentry>
 
+      <varlistentry>
+        <term>beq</term>
         <listitem>
-          <para>matrix (n,p)</para>
+          <para>
+            a n-by-1 matrix of doubles,
+            the right hand side of the linear equality constraint.
+          </para>
         </listitem>
       </varlistentry>
 
       <varlistentry>
-        <term>b</term>
+        <term>c</term>
+        <listitem>
+          <para>
+            a p-by-1 matrix of doubles,
+            the linear part of the objective function.
+          </para>
+        </listitem>
+      </varlistentry>
 
+      <varlistentry>
+        <term>x0</term>
         <listitem>
-          <para>n - vector</para>
+          <para>
+            a p-by-1 matrix of doubles, the initial guess (default x0=[]).
+            If x0==[], the karmarkar function automatically computes a feasible initial
+            guess.
+          </para>
         </listitem>
       </varlistentry>
 
       <varlistentry>
-        <term>c</term>
+        <term>rtolf</term>
+        <listitem>
+          <para>
+            a 1-by-1 matrix of doubles,
+            a relative tolerance on <literal>f(x)=c'*x</literal> (default rtolf=1.d-5).
+          </para>
+        </listitem>
+      </varlistentry>
 
+      <varlistentry>
+        <term>gam</term>
         <listitem>
-          <para>p - vector</para>
+          <para>
+            a 1-by-1 matrix of doubles, the scaling factor (default gam=0.5).
+            The scaling factor must satisfy <literal>0&lt;gam&lt;1</literal>.
+          </para>
         </listitem>
       </varlistentry>
 
       <varlistentry>
-        <term>x0</term>
+        <term>maxiter</term>
+        <listitem>
+          <para>
+            a 1-by-1 matrix of floating point integers, the maximum number of iterations (default maxiter=200).
+            The maximum number of iterations must be greater than 1.
+          </para>
+        </listitem>
+      </varlistentry>
 
+      <varlistentry>
+        <term>outfun</term>
         <listitem>
-          <para>initial vector</para>
+          <para>
+            a function or a list, the output function. See below for details.
+          </para>
         </listitem>
       </varlistentry>
 
       <varlistentry>
-        <term>eps</term>
+        <term>A</term>
+        <listitem>
+          <para>
+            a ni-by-p matrix of doubles, the matrix of linear inequality constraints.
+          </para>
+        </listitem>
+      </varlistentry>
 
+      <varlistentry>
+        <term>b</term>
         <listitem>
-          <para>threshold (default value : 1.d-5)</para>
+          <para>
+            a ni-by-1 matrix of doubles, the right-hand side of linear inequality constraints.
+          </para>
         </listitem>
       </varlistentry>
 
       <varlistentry>
-        <term>gamma</term>
+        <term>lb</term>
+        <listitem>
+          <para>
+            a p-by-1 matrix of doubles, the lower bounds.
+          </para>
+        </listitem>
+      </varlistentry>
+
+      <varlistentry>
+        <term>ub</term>
+        <listitem>
+          <para>
+            a p-by-1 matrix of doubles, the upper bounds.
+          </para>
+        </listitem>
+      </varlistentry>
 
+      <varlistentry>
+        <term>xopt</term>
         <listitem>
-          <para>descent step <literal>0&lt;gamma&lt;1</literal> , default
-          value : 1/4</para>
+          <para>a p-by-1 matrix of doubles, the solution.</para>
         </listitem>
       </varlistentry>
 
       <varlistentry>
-        <term>x1</term>
+        <term>fopt</term>
+        <listitem>
+          <para>
+            a 1-by-1 matrix of doubles, the objective function
+            value at optimum, i.e. <literal>fopt=c'*xopt</literal>.
+          </para>
+        </listitem>
+      </varlistentry>
 
+      <varlistentry>
+        <term>exitflag</term>
         <listitem>
-          <para>solution</para>
+          <para>
+            a 1-by-1 matrix of floating point integers, the status of execution.
+            See below for details.
+          </para>
         </listitem>
       </varlistentry>
 
       <varlistentry>
-        <term>crit</term>
+        <term>iter</term>
+        <listitem>
+          <para>
+            a 1-by-1 matrix of floating point integers, the number of iterations.
+          </para>
+        </listitem>
+      </varlistentry>
 
+      <varlistentry>
+        <term>yopt</term>
         <listitem>
-          <para>value of c'*x1</para>
+          <para>
+            a <literal>struct</literal> containing the dual solution.
+            The structure yopt has four fields : ineqlin, eqlin, upper, lower.
+            The field yopt.ineqlin is the Lagrange multiplier for the inequality constraints.
+            The field yopt.eqlin is the Lagrange multiplier for the equality constraints.
+            The field yopt.upper is the Lagrange multiplier for the upper bounds.
+            The field yopt.lower is the Lagrange multiplier for the lower bounds.
+          </para>
         </listitem>
       </varlistentry>
+
     </variablelist>
   </refsection>
 
   <refsection>
     <title>Description</title>
 
-    <para>Computes <literal>x</literal> which minimizes</para>
+    <para>
+      Computes the solution of linear programming problems.
+    </para>
+
+    <para>
+      This function has the two following modes.
+    </para>
+
+    <itemizedlist>
+      <listitem>
+        <para>
+          If no inequality constraints and no bound is given
+          (i.e. if <literal>(A==[] &#38; lb==[] &#38; ub==[])</literal>),
+          the function ensures that the variable is nonnegative.
+        </para>
+      </listitem>
+      <listitem>
+        <para>
+          If any inequality constraints or any bound is given
+          (i.e. if <literal>(A&#60;&#62;[] or lb&#60;&#62;[] or ub&#60;&#62;[])</literal>),
+          the function takes into account for this inequality or bound
+          (and does not ensure that the variable is nonnegative).
+        </para>
+      </listitem>
+    </itemizedlist>
+
+    <para>
+      If no inequality constraints and no bound is given
+      (i.e. if <literal>(A==[] &#38; lb==[] &#38; ub==[])</literal>),
+      solves the linear optimization problem:
+    </para>
+
+    <para>
+      <latex>
+        \begin{eqnarray}
+        \begin{array}{l}
+        \textrm{minimize } c^T \cdot x\\
+        A_{eq} x = b_{eq}\\
+        x \geq 0
+        \end{array}
+        \end{eqnarray}
+      </latex>
+    </para>
+
+    <para>
+      If any inequality constraints or any bound is given
+      (i.e. if <literal>(A&#60;&#62;[] | lb&#60;&#62;[] | ub&#60;&#62;[])</literal>),
+      solves the linear optimization problem:
+    </para>
+
+    <para>
+      <latex>
+        \begin{eqnarray}
+        \begin{array}{l}
+        \textrm{minimize } c^T \cdot x\\
+        A_{eq} x = b_{eq}\\
+        A x \leq b\\
+        \ell_b \leq x \leq u_b
+        \end{array}
+        \end{eqnarray}
+      </latex>
+    </para>
+
+    <para>
+      Any optional parameter equal to the empty matrix <literal>[]</literal> is replaced by
+      its default value.
+    </para>
+
+    <para>
+      The <literal>exitflag</literal> parameter allows to know why the algorithm
+      terminated.
+    </para>
+
+    <itemizedlist>
+
+      <listitem>
+        <para>
+          <literal>exitflag = 1</literal> if algorithm converged.
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          <literal>exitflag = 0</literal> if maximum number of iterations was reached.
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          <literal>exitflag = -1</literal> if no feasible point was found
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          <literal>exitflag = -2</literal> if problem is unbounded.
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          <literal>exitflag = -3</literal> if search direction became zero.
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          <literal>exitflag = -4</literal> if algorithm stopped on user's request.
+        </para>
+      </listitem>
+
+    </itemizedlist>
+
+    <para>
+      The output function <literal>outfun</literal> must have header
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+stop = outfun ( xopt , optimValues , state )
+]]>
+    </programlisting>
+
+    <para>
+      where
+      <literal>xopt</literal> is a p-by-1 matrix of double representing
+      the current point,
+      <literal>optimValues</literal> is a <literal>struct</literal>,
+      <literal>state</literal> is a 1-by-1 matrix of strings.
+      Here, <literal>stop</literal> is a 1-by-1 matrix of booleans,
+      which is <literal>%t</literal> if the algorithm must stop.
+    </para>
+
+    <para>
+      <literal>optimValues</literal> has the following fields:
+    </para>
+
+    <itemizedlist>
+
+      <listitem>
+        <para>
+          <literal>optimValues.funccount</literal>: a 1-by-1 matrix of floating point integers,
+          the number of function evaluations
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          <literal>optimValues.fval</literal>: a 1-by-1 matrix of doubles, the best function value
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          <literal>optimValues.iteration</literal>: a 1-by-1 matrix of floating point integers,
+          the current iteration number
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          <literal>optimValues.procedure</literal>: a 1-by-1 matrix of strings, the type of step performed.
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          <literal>optimValues.dualgap</literal>: a 1-by-1 matrix of doubles, the duality gap, i.e.
+          <literal>abs(yopt'*beq - fopt)</literal>.
+          At optimum, the duality gap is zero.
+        </para>
+      </listitem>
+
+    </itemizedlist>
+
+    <para>
+      The <literal>optimValues.procedure</literal> field can have the following values.
+    </para>
+
+    <itemizedlist>
+
+      <listitem>
+        <para>
+          If <literal>optimValues.procedure="x0"</literal>, then the algorithm is computing the initial
+          feasible guess <literal>x0</literal> (phase 1).
+        </para>
+      </listitem>
+
+      <listitem>
+        <para>
+          If <literal>optimValues.procedure="x*"</literal>, then the algorithm is computing the solution
+          of the linear program (phase 2).
+        </para>
+      </listitem>
+
+    </itemizedlist>
+
+
+    <para>
+      It might happen that the output function requires additionnal arguments to be evaluated.
+      In this case, we can use the following feature.
+      The function <literal>outfun</literal> can also be the list <literal>(outf,a1,a2,...)</literal>.
+      In this case <literal>outf</literal>, the first element in the list, must have the header:
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+       stop = outf ( xopt , optimValues , state , a1 , a2 , ... )
+]]>
+    </programlisting>
+
+    <para>
+      where the input arguments <literal>a1, a2, ...</literal>
+      will be automatically be appended at the end of the calling sequence.
+    </para>
 
-    <informalequation>
-      <mediaobject>
-        <imageobject>
-          <imagedata align="center" fileref="../mml/karmarkar_equation_1.mml" />
-        </imageobject>
-      </mediaobject>
-    </informalequation>
   </refsection>
 
   <refsection>
-    <title>Examples</title>
+    <title>Stopping rule</title>
+
+    <para>
+      The stopping rule is based on the number of iterations,
+      the relative tolerance on the function value, the duality gap and 
+      the user's output function.
+    </para>
+    
+    <para>
+      In both the phase 1 and phase 2 of the algorithm, we check the
+      duality gap and the boolean:
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+dualgap > 1.e5 * dualgapmin
+]]>
+</programlisting>
+
+    <para>
+      where <literal>dualgap</literal> is the absolute value of the duality gap, and 
+      <literal>dualgapmin</literal> is the minimum absolute duality gap measured during 
+      this step of the algorithm.
+      The duality gap is computed from 
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+dualgap = abs(yopt'*beq - fopt)
+]]>
+</programlisting>
+
+    <para>
+      where <literal>yopt</literal> is the Lagrange multiplier, <literal>beq</literal>
+      is the right hand side of the inequality constraints and <literal>fopt</literal>
+      is the minimum function value attained in the current phase.
+    </para>
+
+    <para>
+      During the second phase of the algorithm (i.e. once <literal>x0</literal> is
+      determined, the termination condition for the function value is based on the boolean:
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+    abs(fprev-fopt)<=rtolf*abs(fprev)
+]]>
+    </programlisting>
+
+    <para>
+      where <literal>fprev</literal> is the previous function value,
+      <literal>fopt</literal> is the current function value and
+      <literal>rtolf</literal> is the relative tolerance on the function value.
+    </para>
+
+  </refsection>
 
-    <programlisting role="example"><![CDATA[ 
-n=10;p=20;
-a=rand(n,p);
+  <refsection>
+    <title>Implementation notes</title>
+
+    <para>
+      The implementation is based on the primal affine scaling algorithm, as
+      discovered by Dikin in 1967, and then re-discovered by Barnes and Vanderbei et al in 1986.
+    </para>
+
+    <para>
+      If the scaling factor <literal>gam</literal> is closer to 1 (say <literal>gam=0.99</literal>,
+      for example), then the number of iterations may be lower.
+      Tsuchiya and Muramatsu proved that if an optimal solution exists, then, for any
+      <literal>gam</literal> lower than 2/3, the sequence converges to a point in the interior point of the optimal face.
+      Dikin proved convergence with <literal>gam=1/2</literal>.
+      Mascarenhas found two examples where the parameter <literal>gam=0.999</literal> lets the
+      algorithm converge to a vertex which is not optimal, if the initial guess is chosen
+      appropriately.
+    </para>
+
+  </refsection>
+
+  <refsection>
+    <title>Example #1</title>
+
+    <para>
+      In the following example, we solve a linear optimization problem with 2 linear
+      equality constraints and 3 unknowns.
+      The linear optimization problem is
+    </para>
+
+    <para>
+      <latex>
+        \begin{eqnarray}
+        \begin{array}{l}
+        \textrm{minimize } -x_1 -x_2\\
+        x_1 - x_2 = 0\\
+        x_1 + x_2 + x_3 = 2\\
+        x \geq 0
+        \end{array}
+        \end{eqnarray}
+      </latex>
+    </para>
+
+    <para>
+      The following script solves the problem.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+    Aeq = [
+    1 -1 0
+    1  1 1
+    ];
+    beq = [0;2];
+    c = [-1;-1;0];
+    x0 = [0.1;0.1;1.8];
+    [xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c)
+    xstar=[1 1 0]'
+]]>
+    </programlisting>
+
+    <para>
+      The previous script produces the following output.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+-->[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c)
+ yopt  =
+   ineqlin: [0x0 constant]
+   eqlin: [2x1 constant]
+   lower: [3x1 constant]
+   upper: [3x1 constant]
+ iter  =
+    68.  
+ exitflag  =
+    1.  
+ fopt  =
+  - 1.9999898  
+ xopt  =
+    0.9999949  
+    0.9999949  
+    0.0000102  
+]]>
+    </programlisting>
+
+    <para>
+      We can explore the Lagrange multipliers by detailing
+      each field of the yopt structure.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+-->yopt.ineqlin
+ ans  =
+     []
+-->yopt.eqlin
+ ans  =
+  - 6.483D-17  
+    1.         
+-->yopt.lower
+ ans  =
+  - 2.070D-10  
+  - 2.070D-10  
+    1.         
+-->yopt.upper
+ ans  =
+    0.  
+    0.  
+    0.  
+]]>
+    </programlisting>
+
+    <para>
+      We can as well give the initial guess x0, as in the following script.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+    Aeq = [
+    1 -1 0
+    1  1 1
+    ];
+    beq = [0;2];
+    c = [-1;-1;0];
+    x0 = [0.1;0.1;1.8];
+    [xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0)
+]]>
+    </programlisting>
+
+    <para>
+      In the case where we need more precision, we can reduce the
+      relative tolerance on the function value.
+      In general, reducing the tolerance increases the number of iterations.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+-->[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,[],1.e-5);
+-->disp([fopt iter])
+  - 1.9999898    68.  
+-->[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,[],1.e-7);
+-->disp([fopt iter])
+  - 1.9999998    74.  
+-->[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,[],1.e-9);
+-->disp([fopt iter])
+  - 2.    78.  
+]]>
+    </programlisting>
+
+  </refsection>
+
+  <refsection>
+    <title>Example #2</title>
+
+    <para>
+      In the following example, we solve a linear optimization problem with 10 random linear
+      equality constraints and 20 unknowns.
+      The initial guess is chosen at random in the [0,1]^p range.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+n=10;
+p=20;
+Aeq=rand(n,p);
 c=rand(p,1);
-x0=abs(rand(p,1));
-b=a*x0;
-x1=karmarkar(a,b,c,x0);
- ]]></programlisting>
+x0=rand(p,1);
+beq=Aeq*x0;
+xopt=karmarkar(Aeq,beq,c,x0);
+// Check constraints
+norm(Aeq*xopt-beq)
+]]>
+    </programlisting>
+  </refsection>
+
+  <refsection>
+    <title>Inequality constraints</title>
+
+    <para>
+      Consider the following linear program with inequality constraints.
+    </para>
+
+    <para>
+      <latex>
+        \begin{eqnarray}
+        \begin{array}{l}
+        \textrm{minimize } - 20 x_1 - 24 x_2\\
+        3 x_1 + 6 x_2 \leq 60 \\
+        4 x_1 + 2 x_2 \leq 32 \\
+        \end{array}
+        \end{eqnarray}
+      </latex>
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+c = [-20 -24]';
+A = [
+3 6
+4 2
+];
+b = [60 32]';
+xopt=karmarkar([],[],c,[],[],[],[],[],A,b)
+]]>
+    </programlisting>
+
+    <para>
+      The previous script produces the following output.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+-->xopt=karmarkar([],[],c,[],[],[],[],[],A,b)
+ xopt  =
+    3.9999125  
+    7.9999912  
+]]>
+    </programlisting>
+
   </refsection>
+
+  <refsection>
+    <title>With bounds</title>
+
+    <para>
+      Consider the following linear optimization problem.
+    </para>
+
+    <para>
+      <latex>
+        \begin{eqnarray}
+        \begin{array}{l}
+        \textrm{minimize } 2 x_1 + 5 x_2 - 2.5 x_3\\
+        x_1 +   S4 x_3 \leq 5 \\
+        E2 x1 - x2 - x3 \leq 0 \\
+        -2 \leq x_1 \leq 2 \\
+        1 \leq x_2 \\
+        0 \leq x_3 \leq 3 \\
+        \end{array}
+        \end{eqnarray}
+      </latex>
+    </para>
+    <para>
+      where <literal>S4 = sin(pi/4)/4</literal> and <literal>E2 = exp(2)</literal>.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+c = [ 2; 5; -2.5];
+S4 = sin(%pi/4)/4;
+E2 = exp(2);
+A = [
+1  0 S4
+E2 -1 -1
+];
+b = [ 5; 0];
+lb = [ -2; 1   ; 0 ];
+ub = [  2; %inf; 3 ];
+xstar = [-2;1;3];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb,ub)
+]]>
+    </programlisting>
+
+    <para>
+      The previous script produces the following output.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+-->[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb,ub)
+    yopt  =
+    ineqlin: [2x1 constant]
+    eqlin: [0x0 constant]
+    lower: [3x1 constant]
+    upper: [3x1 constant]
+    iter  =
+    76.
+    exitflag  =
+    1.
+    fopt  =
+    - 6.4999482
+    xopt  =
+    - 1.9999914
+    1.0000035
+    2.9999931
+]]>
+    </programlisting>
+
+  </refsection>
+
+  <refsection>
+    <title>Configuring an output function</title>
+
+    <para>
+      It may be useful to configure a callback, so that we can customize the
+      printed messages or create a plot.
+      Consider the following linear optimization problem.
+    </para>
+
+    <para>
+      <latex>
+        \begin{eqnarray}
+        \begin{array}{l}
+        \textrm{minimize } - x_1 - x_2\\
+        2 w x_1 + x_2 \leq 1+w^2, \quad w=0.0, 0.1, 0.2, ..., 1.0\\
+        x_1,x_2 \geq 0
+        \end{array}
+        \end{eqnarray}
+      </latex>
+    </para>
+
+    <para>
+      The following output function plots the current point
+      and prints the iteration number, the value of the objective
+      function.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+function stop = myoutputfunction ( xopt , optimValues , state )
+    localmsg = gettext("Iteration #%3.0f, state=%s, procedure=%s, fopt=%10.3e, x=[%s], dualgap=%10.3e\n")
+    xstr = strcat(msprintf("%10.3e\n",xopt)'," ")
+    mprintf(localmsg,optimValues.iteration,state,optimValues.procedure,optimValues.fval,xstr,optimValues.dualgap)
+    plot(xopt(1),xopt(2),"bo")
+    stop = %f
+endfunction
+]]>
+    </programlisting>
+
+    <para>
+      The following script defines the optimization problem and
+      runs the optimization.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+  n=11;
+  A = [2*linspace(0,1,n)',ones(n,1)];
+  b = 1 + linspace(0,1,n)'.^2;
+  c=[-1;-1];
+  // Plot the constraints
+  scf();
+  for i = 1 : n
+    plot(linspace(0,1,100),b(i)-A(i,1)*linspace(0,1,100),"b-")
+  end
+  // Run the optimization
+  xopt=karmarkar([],[],c,[],[],[],[],myoutputfunction,A,b);
+  // Plot the starting and ending points
+  plot(xopt(1),xopt(2),"k*")
+]]>
+    </programlisting>
+
+    <para>
+      The previous script produces the following output and creates a graphics.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+-->xopt=karmarkar([],[],c,[],[],[],[],myoutputfunction,A,b)
+Iteration #  0, state=init, procedure=x0, fopt=1.000e+000, x=[0.000e+000 0.000e+000], dualgap=Inf
+Iteration #  0, state=init, procedure=x0, fopt=1.000e+000, x=[0.000e+000 0.000e+000], dualgap=Inf
+Iteration #  1, state=init, procedure=x0, fopt=5.000e-001, x=[2.201e-001 -4.313e-002], dualgap=3.676e-001
+Iteration #  2, state=init, procedure=x0, fopt=2.500e-001, x=[3.283e-001 -6.512e-002], dualgap=2.140e-001
+Iteration #  3, state=init, procedure=x0, fopt=1.250e-001, x=[3.822e-001 -7.634e-002], dualgap=1.161e-001
+Iteration #  4, state=init, procedure=x0, fopt=6.250e-002, x=[4.090e-001 -8.202e-002], dualgap=6.033e-002
+Iteration #  5, state=init, procedure=x0, fopt=3.125e-002, x=[4.225e-001 -8.488e-002], dualgap=3.072e-002
+Iteration #  6, state=init, procedure=x0, fopt=1.563e-002, x=[4.292e-001 -8.632e-002], dualgap=1.549e-002
+[...]
+Iteration # 50, state=init, procedure=x0, fopt=8.882e-016, x=[4.359e-001 -8.775e-002], dualgap=8.882e-016
+Iteration # 51, state=init, procedure=x0, fopt=4.441e-016, x=[4.359e-001 -8.775e-002], dualgap=4.441e-016
+Iteration # 52, state=init, procedure=x0, fopt=2.220e-016, x=[4.359e-001 -8.775e-002], dualgap=2.220e-016
+Iteration # 52, state=init, procedure=x*, fopt=-3.481e-001, x=[4.359e-001 -8.775e-002], dualgap=Inf
+Iteration # 52, state=iter, procedure=x*, fopt=-3.481e-001, x=[4.359e-001 -8.775e-002], dualgap=Inf
+Iteration # 53, state=iter, procedure=x*, fopt=-7.927e-001, x=[5.249e-001 2.678e-001], dualgap=5.098e-001
+[...]
+Iteration # 65, state=iter, procedure=x*, fopt=-1.250e+000, x=[5.005e-001 7.494e-001], dualgap=1.258e-004
+Iteration # 66, state=iter, procedure=x*, fopt=-1.250e+000, x=[5.005e-001 7.494e-001], dualgap=5.941e-005
+Iteration # 67, state=iter, procedure=x*, fopt=-1.250e+000, x=[5.005e-001 7.495e-001], dualgap=2.882e-005
+Iteration # 68, state=iter, procedure=x*, fopt=-1.250e+000, x=[5.005e-001 7.495e-001], dualgap=1.418e-005
+Iteration # 69, state=done, procedure=x*, fopt=-1.250e+000, x=[5.005e-001 7.495e-001], dualgap=7.035e-006
+ xopt  =
+    0.5005127  
+    0.7494803  
+]]>
+    </programlisting>
+
+  </refsection>
+
+
+  <refsection>
+    <title>Infeasible problem</title>
+
+    <para>
+      Consider the following linear optimization problem.
+      It is extracted from "Linear Programming in Matlab"
+      Ferris, Mangasarian, Wright, 2008, Chapter 3, "The Simplex Method", Exercise 3-4-2 1.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+    // An infeasible problem.
+    // Minimize -3 x1 + x2
+    //  - x1 -   x2 >= -2
+    //  2 x1 + 2 x2 >= 10
+    // x >= 0
+    c = [-3;1];
+    A=[
+    -1 -1
+    2  2
+    ];
+    A=-A;
+    b=[-2;10];
+    b=-b;
+    lb=[0;0];
+    [xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb)
+]]>
+    </programlisting>
+
+    <para>
+      The previous script produces the following output.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+-->[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb)
+    yopt  =
+
+    ineqlin: [0x0 constant]
+    eqlin: [0x0 constant]
+    lower: [0x0 constant]
+    upper: [0x0 constant]
+    iter  =
+    40.
+    exitflag  =
+    - 1.
+    fopt  =
+    []
+    xopt  =
+    []
+]]>
+    </programlisting>
+
+  </refsection>
+
+
+  <refsection>
+    <title>Unbounded problem</title>
+
+    <para>
+      Consider the following linear optimization problem.
+      It is extracted from "Linear and Nonlinear Optimization"
+      Griva, Nash, Sofer, 2009, Chapter 5, "The Simplex Method", Example 5.3.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+  // An unbounded problem.
+  // Minimize -x1 - 2 x2
+  //  -1 x1 +   x2 <= 2
+  //  -2 x1 +   x2 <= 1
+  // x >= 0
+  c = [-1;-2];
+  A=[
+  -1  1
+  -2  1
+  ];
+  b=[2;1];
+  lb=[0;0];
+  [xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb)
+]]>
+    </programlisting>
+
+    <para>
+      The previous script produces the following output.
+      Notice that the algorithm detects that the duality gap has increased much
+      more than expected, indicating a failure of the algorithm to find an
+      optimal point.
+    </para>
+
+    <programlisting role="example">
+      <![CDATA[ 
+-->[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb)
+    yopt  =
+    ineqlin: [2x1 constant]
+    eqlin: [0x0 constant]
+    lower: [2x1 constant]
+    upper: [2x1 constant]
+    iter  =
+    59.
+    exitflag  =
+    - 2.
+    fopt  =
+    - 45374652.
+    xopt  =
+    15124883.
+    15124885.
+]]>
+    </programlisting>
+
+  </refsection>
+
+
+  <refsection>
+    <title>References</title>
+    <para>
+      "Iterative solution of problems of linear and quadratic programming",
+      Dikin,
+      Doklady Akademii Nausk SSSR, Vol. 174, pp. 747-748, 1967
+    </para>
+
+    <para>
+      "A New Polynomial Time Algorithm for Linear Programming",
+      Narendra Karmarkar, Combinatorica, Vol 4, nr. 4, p. 373–395, 1984.
+    </para>
+
+    <para>
+      "A variation on Karmarkar’s algorithm for solving linear programming problems,
+      Earl R. Barnes, Mathematical Programming, Volume 36, Number 2, 174-182, 1986.
+    </para>
+
+    <para>
+      "A modification of karmarkar's linear programming algorithm",
+      Robert J. Vanderbei, Marc S. Meketon and Barry A. Freedman,
+      Algorithmica, Volume 1, Numbers 1-4, 395-407, 1986.
+    </para>
+
+    <para>
+      "Practical Optimization: Algorithms and Engineering Applications",
+      Andreas Antoniou, Wu-Sheng Lu, Springer, 2007,
+      Chapter 12, "Linear Programming Part II: Interior Point Methods".
+    </para>
+
+    <para>
+      "Global Convergence of a Long-Step Affine Scaling Algorithm for Degenerate Linear Programming Problems",
+      Takashi Tsuchiya and Masakazu Muramatsu,
+      SIAM J. Optim. Volume 5, Issue 3, pp. 525-551 (August 1995)
+    </para>
+
+    <para>
+      "The convergence of dual variables",
+      Dikin,
+      Tech. report, Siberian Energy Institute, Russia, 1991
+    </para>
+
+    <para>
+      "The Affine Scaling Algorithm Fails for Stepsize 0.999",
+      Walter F. Mascarenhas, SIAM J. Optim. Volume 7, Issue 1, pp. 34-46 (1997)
+    </para>
+
+    <para>
+      "A Primal-Dual Exterior Point Algorithm For Linear Programming Problems"
+      Nikolaos Samaras, Angelo Sifaleras, Charalampos Triantafyllidis
+      Yugoslav Journal of Operations Research
+      Vol 19 (2009), Number 1, 123-132
+    </para>
+
+  </refsection>
+
 </refentry>
+
index 6ea8278..e8ad0d0 100644 (file)
@@ -1,5 +1,6 @@
 // Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
 // Copyright (C) INRIA
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
 // 
 // This file must be used under the terms of the CeCILL.
 // This source file is licensed as described in the file COPYING, which
 // are also available at    
 // http://www.cecill.info/licences/Licence_CeCILL_V2-en.txt
 //
-function [x1,crit]=karmarkar(a,b,c,x0,eps,Gamma)
-maxiter=200;
-epsdef=1.d-5
+
+// xopt=karmarkar(Aeq,beq,c)
+// xopt=karmarkar(Aeq,beq,c,x0)
+// xopt=karmarkar(Aeq,beq,c,x0,rtolf)
+// xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam)
+// xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter)
+// xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun)
+// xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun,A,b)
+// xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun,A,b,lb)
+// xopt=karmarkar(Aeq,beq,c,x0,rtolf,gam,maxiter,outfun,A,b,lb,ub)
+// [xopt,fopt] = karmarkar(...)
+// [xopt,fopt,exitflag] = karmarkar(...)
+// [xopt,fopt,exitflag,iter] = karmarkar(...)
+// [xopt,fopt,exitflag,iter,yopt] = karmarkar(...)
 //
-[lhs,rhs]=argn(0)
-select rhs
-case 4 then
-  Gamma=1/4,eps=epsdef
-case 5 then
-  Gamma=1/4
-case 6,
-else
- error('[x1,crit]=karmarkar(a,b,c,x0 [,eps [,Gamma]])')
-end
-//verification des donnees
-[n,p]=size(a)
-w=size(b)
-if w(1)<>n then 
-  error(msprintf(gettext("%s: Wrong size for input argument #%d."),'karmarkar',2));
-end
-w=size(c)
-if w(1)<>p then
-  error(msprintf(gettext("%s: Wrong size for input argument #%d."),'karmarkar',1));
-end
-w=size(x0)
-if w(1)<>p then
-  error(msprintf(gettext("%s: Wrong size for input argument #%d."),'karmarkar',4));
-end
-if min(x0)<0|norm(a*x0-b)>eps then 
-  error(msprintf(gettext("%s: x0 is not feasible."),'karmarkar'));
-end
-//
-x1=x0;tc=c';
-crit=tc*x1;
-test=eps+1
-count=0
-while test>eps&count<=maxiter
-    count=count+1;
-//    ax=a*diag(x1);
-    ax=a.*(ones(size(a,1),1)*x1');
-    xc=x1.*c;
-    y=ax'\xc;
-//    y=(ax*ax')\(ax*xc)
-    d=-xc+ax'*y;
-    dk=x1.*d;
-    if min(dk)>0 then 
-      error(msprintf(gettext("%s: Unbounded problem!"),'karmarkar'));
-    end
-    alpha=-Gamma/min(d);
-    test=alpha*(norm(d)**2)/max(1,abs(crit));
-    x1=x1+alpha*dk;
-    crit=tc*x1;
-    write(%io(2),[count,crit,test],'(f3.0,3x,e10.3,3x,e10.3)')
-end
-  
+// exitflag = 1 if algorithm converged.
+// exitflag = 0 if maximum number of iterations was reached.
+// exitflag = -1 if no feasible point was found
+// exitflag = -2 if problem is unbounded.
+// exitflag = -3 if search direction became zero.
+// exitflag = -4 if algorithm stopped on user's request.
+// exitflag = -5 if duality gap became too large
+// exitflag = -%inf on internal error.
+
+function [xopt,fopt,exitflag,iter,yopt]=karmarkar(varargin)
+    function argin = argindefault ( rhs , vararglist , ivar , default )
+        // Returns the value of the input argument #ivar.
+        // If this argument was not provided, or was equal to the 
+        // empty matrix, returns the default value.
+        if ( rhs < ivar ) then
+            argin = default
+        else
+            if ( vararglist(ivar) <> [] ) then
+                argin = vararglist(ivar)
+            else
+                argin = default
+            end
+        end
+    endfunction
+
+    //
+    [lhs,rhs]=argn(0)
+    if ( rhs<3 | rhs>12 | rhs==9 ) then
+        error(msprintf(gettext("%s: Wrong number of input arguments: %d to %d expected.\n"),"karmarkar",3,12));
+    end
+    //
+    Aeq = varargin(1)
+    beq = varargin(2)
+    c = varargin(3)
+    x0 = argindefault ( rhs , varargin , 4 , [] )
+    rtolf = argindefault ( rhs , varargin , 5 , 1.d-5 )
+    gam = argindefault ( rhs , varargin , 6 , 1/2 )
+    maxiter = argindefault ( rhs , varargin , 7 , 200 )
+    __karmarkar_outfun__ = argindefault ( rhs , varargin , 8 , [] )
+    A = argindefault ( rhs , varargin , 9 , [] )
+    b = argindefault ( rhs , varargin , 10 , [] )
+    lb = argindefault ( rhs , varargin , 11 , [] )
+    ub = argindefault ( rhs , varargin , 12 , [] )
+    //
+    // Check input arguments
+    //
+    //
+    // Check type
+    //
+    if ( typeof(Aeq) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",1));
+    end
+    if ( typeof(beq) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",2));
+    end
+    if ( typeof(c) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",3));
+    end
+    if ( typeof(x0) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",4));
+    end
+    if ( typeof(rtolf) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",5));
+    end
+    if ( typeof(gam) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",6));
+    end
+    if ( typeof(maxiter) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",7));
+    end
+    if ( __karmarkar_outfun__ <> [] ) then
+        if ( and(typeof(__karmarkar_outfun__) <> ["function" "list" "fptr"] ) ) then
+            error(msprintf(gettext("%s: Wrong type for input argument #%d: %s expected.\n"),"karmarkar",8,"function or list"));
+        end
+        if ( typeof(__karmarkar_outfun__)=="list" ) then
+            if ( and(typeof(__karmarkar_outfun__(1))<>["function" "fptr"] ) ) then
+                error(msprintf(gettext("%s: Wrong type for input argument #%d: %s expected.\n"),"karmarkar",8,"function"));
+            end
+        end
+    end
+    if ( typeof(A) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",9));
+    end
+    if ( typeof(b) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",10));
+    end
+    if ( typeof(lb) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",11));
+    end
+    if ( typeof(ub) <> "constant" ) then
+        error(msprintf(gettext("%s: Wrong type for input argument #%d: Real matrix expected.\n"),"karmarkar",12));
+    end
+    //
+    // Check size
+    //
+    [ne,pe]=size(Aeq)
+    [ni,pi]=size(A)
+    plb=size(lb,"*")
+    pub=size(ub,"*")
+    if ( Aeq <> [] ) then
+        p = pe
+    elseif ( A <> [] ) then
+        p = pi
+    elseif ( lb <> [] ) then
+        p = plb
+    elseif ( ub <> [] ) then
+        p = pub
+    else
+        error(msprintf(gettext("%s: Either one of Aeq, A, lb or ub must be non-empty."),"karmarkar"));
+    end
+    // The case where Aeq==[] and A==[] is treated below.
+    if ( ( Aeq <> [] & beq == [] ) | ( Aeq == [] & beq <> [] ) ) then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",2));
+    end
+    if ( Aeq <> [] ) then
+    if ( or ( size(beq) <> [ne 1] ) ) then 
+        error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",2));
+    end
+    end
+    if ( or ( size(c) <> [p 1] ) ) then 
+        error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",3));
+    end
+    if ( x0 <> [] ) then
+        if ( or ( size(x0) <> [p 1] ) ) then 
+            error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",4));
+        end
+    end
+    if ( or ( size(rtolf) <> [1 1] ) ) then 
+        error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",5));
+    end
+    if ( or ( size(gam) <> [1 1] ) ) then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",6));
+    end
+    if ( or ( size(maxiter) <> [1 1] ) ) then 
+        error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",7));
+    end
+    if ( A <> [] ) then
+        if ( Aeq <> [] ) then
+            if ( pi <> pe ) then 
+                error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",9));
+            end
+        end
+    end
+    if ( ( A <> [] & b == [] ) | ( A == [] & b <> [] ) ) then
+        error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",10));
+    end
+    if ( b <> [] ) then
+        if ( size(b) <> [ni 1] ) then 
+            error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",10));
+        end
+    end
+    if ( lb <> [] ) then
+        if ( or ( size(lb) <> [p 1] ) ) then 
+            error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",11));
+        end
+    end
+    if ( ub <> [] ) then
+        if ( or ( size(ub) <> [p 1] ) ) then 
+            error(msprintf(gettext("%s: Wrong size for input argument #%d."),"karmarkar",12));
+        end
+    end
+    //
+    // Check content
+    //
+    if ( rtolf < 0 | rtolf > 1 ) then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d. rtolf must be in [0,1]."),"karmarkar",5));
+    end
+    if ( gam < 0 | gam > 1 ) then
+        error(msprintf(gettext("%s: Wrong value for input argument #%d. gam must be in [0,1]."),"karmarkar",6));
+    end
+    if ( maxiter < 1 ) then 
+        error(msprintf(gettext("%s: Wrong value for input argument #%d. maxiter must be greater than 1."),"karmarkar",7));
+    end
+    if ( floor(maxiter) <> maxiter ) then 
+        error(msprintf(gettext("%s: Wrong value for input argument #%d. maxiter must be a floating point integer."),"karmarkar",7));
+    end
+    if ( lb == [] & ub == [] ) then
+    if ( x0 <> [] ) then
+        if ( min(x0)<0 ) then
+            error(msprintf(gettext("%s: Wrong value for input argument #%d. x0 is not positive."),"karmarkar",4));
+        end
+    end
+    end
+    if ( lb <> [] & ub <> [] ) then 
+    if ( or ( ub < lb ) ) then 
+        error(msprintf(gettext("%s: Wrong value for input argument #%d. One entry of the upper bound ub is greater than the lower bound lb."),"karmarkar",12));
+    end
+    end
+    if ( lb <> [] & x0 <> [] ) then
+    if ( or ( lb > x0 ) ) then 
+        error(msprintf(gettext("%s: Wrong value for input argument #%d. x0 lower than lower bound lb."),"karmarkar",12));
+    end
+    end
+    if ( ub <> [] & x0 <> [] ) then
+    if ( or ( ub < x0 ) ) then 
+        error(msprintf(gettext("%s: Wrong value for input argument #%d. x0 greater than upper bound ub."),"karmarkar",12));
+    end
+    end
+    //
+    // Proceed
+    //
+    // Transform the general LP into standard form
+    [AAeq,bbeq,cc,xx0,pinit,newposvars] = karmarkar_preprocess ( Aeq , beq , c , A , b , lb , ub , x0 )
+    //
+    // Solve the standard form LP
+    outfun = list(__karmarkar_outfun__ , pinit , newposvars )
+    [xxopt,fopt,exitflag,iter,yyopt] = karmarkar_findStandardLP ( AAeq , bbeq , cc , xx0 , rtolf , gam , maxiter , outfun )
+    //
+    // Extract the solution from the initial problem
+    [xopt,fopt,yopt] = karmarkar_postprocess ( Aeq , beq , c , A , b , lb , ub , pinit , xxopt , yyopt , exitflag )
+endfunction
+
+function [AAeq,bbeq,cc,xx0,pinit,newposvars] = karmarkar_preprocess ( Aeq , beq , c , A , b , lb , ub , x0 )
+    // Transform the general LP into a standard LP.
+    //
+    // Parameters
+    // Aeq : a ne-by-p matrix of doubles, the matrix of linear equality constraints
+    // beq : a ne-by-1 matrix of doubles, the right hand side of linear equality constraints
+    // A : a ni-by-p matrix of doubles, the matrix of linear inequality constraints
+    // b : a ni-by-1 matrix of doubles, the right hand side of linear inequality constraints
+    // c : a p-by-1 matrix of doubles, the linear objective
+    // x0 : a p-by-1 matrix of doubles, the initial guess. If x0=[], this means that no initial guess is provided.
+    // AAeq : a m-by-q matrix of doubles, the matrix of linear equality constraints. We have m = ne + ni and q = p + ni.
+    // bbeq : a m-by-1 matrix of doubles, the right hand side of linear equality constraints
+    // cc : a p-by-1 matrix of doubles, the linear objective
+    // xx0 : a q-by-1 matrix of doubles, the initial guess. If xx0=[], this means that no initial guess is provided.
+    // pinit : the initial value of p. The solution of the original general LP is xs(1:pinit), where xs is the solution of the modified standard LP.
+    // newposvars : a 1-by-1 matrix of booleans, %f if no positive variables have been introduced, %t if positive variables have been introduced.
+    //
+    // Description
+    // Transform the general linear program into a standard linear program.
+    //
+    // More precisely, if A <> [] transform the general linear program (L.P.):
+    //
+    // min c'*x
+    // A * x <= b
+    // Aeq * x = beq
+    //
+    // into the standard LP:
+    //
+    // min cc'*z
+    // AAeq * z = bbeq
+    // z >= 0
+    //
+    // where z is the new unknown, positive, which may contain slack variables.
+    //
+    // An unrestricted variable xi is transformed into x = xp - xn, where xp,xn >= 0.
+    // This turns the inequality constraint Ax <= b into [A -A][xp;xn] <= b.
+    // The inequality constraint is turned into an equality constraint by introducing 
+    // slack variables si.
+    // This transforms the inequality into: [A -A][xp;xn] + s = b which can be written :
+    //
+    // [A -A I][xp;xn;s] = b
+    //
+    // Therefore, the number of variables increases from pinit to 2*pinit + ni.
+    // The initial variable is x(1:pinit)-x(pinit+1:2*pinit).
+    // The same happens for the equality constraints which is turned from 
+    // Aeq*x=b to [Aeq -Aeq 0][xp;xn;s]=b.
+    // The cost function is turned from c'*x into [c;-c;0]'[xp;xn].
+    //
+    newposvars = %f
+    //
+    [ne,pe]=size(Aeq)
+    [ni,pi]=size(A)
+    plb=size(lb,"*")
+    pub=size(ub,"*")
+    if ( Aeq <> [] ) then
+        pinit = pe
+    elseif ( A <> [] ) then
+        pinit = pi
+    elseif ( lb <> [] ) then
+        pinit = plb
+    elseif ( ub <> [] ) then
+        pinit = pub
+    end
+    //
+    if ( A == [] & lb == [] & ub == [] ) then
+        AAeq = Aeq
+        bbeq = beq
+        cc = c
+        xx0 = x0
+        return
+    end
+    //
+    // Process the lower bound : x >= lb.
+    // Move it as an inequality constraint: -I * x <= -lb.
+    // TODO : process the case lb <> [] by shifting x (instead of introducing positive variables).
+    if ( lb <> [] ) then
+        if ( A == [] ) then
+            A = -eye(plb,plb)
+            b = -lb
+        else
+            A = [A;-eye(plb,plb)]
+            b = [b;-lb]
+        end
+    end
+    //
+    // Process the upper bound : x <= ub.
+    // Move it as an inequality constraint : I * x <= ub.
+    if ( ub <> [] ) then
+        if ( A == [] ) then
+            A = eye(pub,pub)
+            b = ub
+        else
+            A = [A;eye(pub,pub)]
+            b = [b;ub]
+        end
+    end
+    //
+    // Remove constraints where b(i) = %inf.
+    // Such a A(i,:)*x <= %inf = b(i) will be satisfied anyway, but 
+    // may cause failures in the algorithm.
+    iinf = find(b == %inf)
+    b(iinf) = []
+    A(iinf,:) = []
+    //
+    // Create the map from the initial constraints to the final constraints.
+    //
+    // Update the number of inequalities,
+    // given that the bounds have been updated.
+    [ni,pi]=size(A)
+    //
+    // Initialize AAeq, bbeq, cc and xx0.
+    //
+    // If ni inequality constraints are given, transform the problem by 
+    // adding pinit positive variables and ni slack variables.
+    // The inequality is Ax <= b.
+    if ( A <> [] ) then
+        //
+        // Create the matrix 
+        // AAeq = [
+        //     Aeq -Aeq 0
+        //     A   -A   I
+        // ]
+        AAeq (1:ne+ni,1:2*pinit+ni) = zeros(ne+ni,2*pinit+ni)
+        if ( Aeq <> [] ) then
+            AAeq (1:ne,1:pinit) = Aeq
+            AAeq (1:ne,pinit+1:2*pinit) = -Aeq
+        end
+        AAeq (ne+1:ne+ni,1:pinit) = A
+        AAeq (ne+1:ne+ni,pinit+1:2*pinit) = -A
+        AAeq (ne+1:ne+ni,2*pinit+1:2*pinit+ni) = eye(ni,ni)
+        bbeq = [beq;b]
+        cc = [c;-c;zeros(ni,1)]
+        newposvars = %t
+        if ( x0 == [] ) then
+            xx0 = []
+        else
+            xx0 = zeros(2*pinit+ni,1)
+            xx0(1:pinit) = max(x0,0)
+            xx0(pinit+1:2*pinit) = -min(x0,0)
+            s = b - A*x0
+            if ( min(s)<0 ) then
+              error(msprintf(gettext("%s: Wrong value for input argument #%d. x0 does not satisfy the inequality constraints."),"karmarkar",4));
+            end
+            xx0(2*pinit+1:2*pinit+ni) = s
+        end
+    end
+endfunction
+
+function [xopt,fopt,yopt] = karmarkar_postprocess ( Aeq , beq , c , A , b , lb , ub , pinit , xxopt , yyopt , exitflag )
+    // Transform the solution of the standard LP into the solution of the original general LP.
+    //
+    // Extract the solution from the initial problem
+    if ( pinit < size(xxopt,"*") ) then
+        xopt = xxopt(1:pinit) - xxopt(pinit+1:2*pinit)
+    else
+        xopt = xxopt
+    end
+    //
+    // Extract the dual solution.
+    [ne,pe]=size(Aeq)
+    [ni,pi]=size(A)
+    plb=size(lb,"*")
+    pub=size(ub,"*")
+    // 
+    // Initialize
+    if ( yyopt == [] ) then
+        yopt.ineqlin = []
+        yopt.eqlin = []
+        yopt.lower = []
+        yopt.upper = []
+    else
+        yopt.ineqlin = []
+        yopt.eqlin = []
+        yopt.lower = zeros(pinit,1)
+        yopt.upper = zeros(pinit,1)
+    end
+    //
+    // Update depending on the presence of the options.
+    kstart = 1
+    if ( ne > 0 ) then
+      kstop = kstart + ne - 1
+      yopt.eqlin = -yyopt(kstart:kstop)
+      kstart = kstop + 1
+    end
+    if ( ni > 0 ) then
+      noninf = find(b<>%inf)
+      kinf = find(b==%inf)
+      kstop = kstart + ni - 1 - size(kinf,"*")
+      yopt.ineqlin(noninf) = -yyopt(kstart:kstop)
+      kstart = kstop + 1
+    end
+    if ( ni == 0 & plb == 0 & pub == 0 ) then
+      yopt.lower = c - Aeq'*yyopt
+    elseif ( plb > 0 ) then
+      noninf = find(lb<>%inf)
+      kinf = find(lb==%inf)
+      kstop = kstart + plb - 1 - size(kinf,"*")
+      yopt.lower(noninf) = -yyopt(kstart:kstop)
+      kstart = kstop + 1
+    end
+    if ( pub > 0 ) then
+      noninf = find(ub<>%inf)
+      kinf = find(ub==%inf)
+      kstop = kstart + pub - 1 - size(kinf,"*")
+      yopt.upper(noninf) = -yyopt(kstart:kstop)
+      kstart = kstop + 1
+    end
+endfunction
+
+function [xopt,fopt,exitflag,iter,yopt] = karmarkar_findStandardLP ( Aeq , beq , c , x0 , rtolf , gam , maxiter , outfun )
+    // Solves a linear problem in standard form (searches for x0 if necessary).
+    //
+    // Parameters
+    // pinit : a 1-by-1 matrix of floating point integers, the number of parameters before the introduction of slack variables.
+    //
+    // Given x0, uses a primal affine scaling (P.A.S.) algorithm to iteratively
+    // find the solution of a linear program in standard form:
+    //
+    // min c'*x
+    // Aeq*x = beq
+    // x >= 0
+    //
+    // If x0 is the empty matrix, compute a strictly feasible point.
+    //
+    [ne,p]=size(Aeq)
+    if ( x0 == [] ) then
+        //
+        // Compute a strictly feasible point.
+        //
+        // Add an extra variable x(p+1) and solve :
+        //
+        // min z(p+1)
+        // AAeq*z = beq
+        // z >= 0
+        //
+        // where z=[x(1) x(2) ... x(p) x(p+1)].
+        // AAeq has the same number of rows, but one column more than Aeq.
+        // The initial guess for the modified problem is z0=[1 1 ... 1]=e(n+1).
+        // The last column of AAeq is beq-Aeq*e(n) so that the initial guess z0 satisfies AAeq*z0=beq.
+        //
+        // References
+        // "A Primal-Dual Exterior Point Algorithm For Linear Programming Problems"
+        // Nikolaos Samaras, Angelo Sifaleras, Charalampos Triantafyllidis
+        // Yugoslav Journal of Operations Research
+        // Vol 19 (2009), Number 1, 123-132
+        //
+        // "A modification of karmarkar's linear programming algorithm",
+        // Robert J. Vanderbei, Marc S. Meketon and Barry A. Freedman,
+        // Algorithmica, Volume 1, Numbers 1-4, 395-407, 1986.
+        //
+        s = beq - Aeq*ones(p,1)
+        AAeq = [Aeq,s]
+        cc = [zeros(p,1);1]
+        xx0 = ones(p+1,1)
+        step = 1
+        xfeasmax = %eps
+        iterstart = 0
+        [xopt1,fopt1,exitflag1,iter1,yopt1] = karmarkar_findxopt ( AAeq , beq , cc , xx0 , rtolf , gam , maxiter , outfun , step , xfeasmax , iterstart )
+        if ( exitflag1 <> 1 ) then
+            // The algorithm did not converge.
+            xopt = []
+            fopt = []
+            exitflag = -1
+            iter = iter1
+            yopt = []
+            return
+        end
+        x0 = xopt1(1:p)
+    else
+        iter1 = 0
+    end
+    //
+    // Check that the initial guess is feasible.
+    if ( x0 <> [] ) then
+        if ( norm(Aeq*x0-beq)>sqrt(%eps)*norm(beq) ) then
+            error(msprintf(gettext("%s: Wrong value for input argument #%d. x0 does not satisfy the equality constraints."),"karmarkar",4));
+        end
+    end
+    //
+    // Find the optimum x
+    step = 2
+    xfeasmax = %nan
+    iterstart = iter1
+    [xopt,fopt,exitflag,iter,yopt] = karmarkar_findxopt ( Aeq , beq , c , x0 , rtolf , gam , maxiter , outfun , step , xfeasmax , iterstart )
+endfunction
+
+function [xopt,fopt,exitflag,iter,yopt] = karmarkar_findxopt ( Aeq , beq , c , x0 , rtolf , gam , maxiter , outfun , step , xfeasmax , iterstart )
+    // Solves a linear problem in standard form, given x0.
+    //
+    // Parameters
+    // step : a 1-by-1 matrix of floating point integers, the kind of algorithm performed. We have step=1 during the search for a feasible point x0, and step=2 during the search for x*.
+    // xfeasmax : a 1-by-1 matrix of double, the maximum value of the feasibility variable, during the search for a feasible point x0.
+    // pinit : a 1-by-1 matrix of floating point integers, the number of parameters before the introduction of slack variables.
+    // iterstart : the initial number of iterations, including the iterations of the previous steps
+    //
+    // Description
+    // Given x0, uses a primal affine scaling (P.A.S.) algorithm to iteratively
+    // find the solution of a linear program in standard form:
+    //
+    // min c'*x
+    // Aeq*x = beq
+    // x >= 0
+    //
+    // We assume that x0 is strictly feasible, i.e. || Aeq*x0-beq || is small and x0 > 0.
+    //
+    // If step = 1, we stop when xopt($) is below the feasibility threshold xfeasmax.
+    // If step = 2, we stop when the objective function does not vary anymore, or the maximum number 
+    // of iterations exceeds the maximum, or the users asks to.
+    //
+    // exitflag = 1 if algorithm converged.
+    // exitflag = 0 if maximum number of iterations was reached.
+    // exitflag = -2 if problem is unbounded.
+    // exitflag = -3 if search direction became zero.
+    // exitflag = -4 if step=2 and algorithm stopped on user's request.
+    // exitflag = -%inf on internal error.
+    //
+    // References
+    //      "A variation on Karmarkar’s algorithm for solving linear programming problems,
+    //      Earl R. Barnes, Mathematical Programming, Volume 36, Number 2, 174-182, 1986.
+    //      
+    //      "A modification of karmarkar's linear programming algorithm",
+    //      Robert J. Vanderbei, Marc S. Meketon and Barry A. Freedman,
+    //      Algorithmica, Volume 1, Numbers 1-4, 395-407, 1986.
+    //      
+    //      "Practical Optimization: Algorithms and Engineering Applications",
+    //      Andreas Antoniou, Wu-Sheng Lu, Springer, 2007,
+    //      Chapter 12, "Linear Programming Part II: Interior Point Methods".
+    //
+    [ne,p]=size(Aeq)
+    xopt=x0
+    yopt = []
+    tc=c'
+    fopt=tc*xopt
+    dualgapmin = %inf
+    dualgap = %inf
+    funccount = 1
+    fprev = fopt+1
+    iter=iterstart
+    s = zeros(p,1)
+    stop = %f
+    exitflag = -%inf
+    firstloop = %t
+    //
+    state = "init"
+    if ( step == 1 ) then
+        procedure = "x0"
+    else
+        procedure = "x*"
+    end
+    optimValues = struct(...
+    "funccount" , funccount , ...
+    "fval" , fopt , ...
+    "iteration" , iter , ...
+    "procedure" , procedure, ...
+    "dualgap" , dualgap ...
+    );
+    stop = karmarkar_outfunDriver ( xopt , optimValues , state , outfun )
+    //
+    while ( %t )
+        if ( iter >= maxiter ) then
+            exitflag = 0
+            break
+        end
+        if ( step == 1 ) then
+            if ( xopt($) <= xfeasmax ) then
+                exitflag = 1
+                break
+            end
+        end
+        if ( step == 2 ) then
+            if ( abs(fprev-fopt)<=rtolf*abs(fprev) ) then
+                exitflag = 1
+                break
+            end
+        end
+        //
+        // Compute the duality gap
+        if ( dualgap > 1.e5 * dualgapmin ) then
+            // Unbounded problem.
+            exitflag = -2
+            break
+        end
+        //
+        // Calls back the output function
+        if ( step == 1 ) then
+            state = "init"
+        else
+            state = "iter"
+        end
+        if ( step == 1 ) then
+            procedure = "x0"
+        else
+            procedure = "x*"
+        end
+        optimValues = struct(...
+        "funccount" , funccount , ...
+        "fval" , fopt , ...
+        "iteration" , iter , ...
+        "procedure" , procedure, ...
+        "dualgap" , dualgap ...
+        );
+        stop = karmarkar_outfunDriver ( xopt , optimValues , state , outfun )
+        if ( stop ) then
+            exitflag = -4
+            break
+        end
+        iter=iter+1
+        // Compute B as B = Aeq*X where X = diag(xopt).
+        // The following method is equivalent, but faster.
+        xt = xopt'
+        B = Aeq.*xt(ones(ne,1),:)
+        v = xopt.*c
+        // y = inv(B*B') * (B*v) i.e. y is the solution of (B*B') y = B*v. 
+        // This implies that y is the solution of (B')*y = v, i.e. 
+        yopt = B'\v
+        p = -v+B'*yopt
+        if ( min(p)==0 ) then
+            exitflag = -3
+            break
+        end
+        d = xopt.*p
+        if ( min(d)>0 ) then 
+            // Unbounded problem.
+            exitflag = -2
+            break
+        end
+        alpha = -gam / min(p)
+        s = alpha*d
+        xopt=xopt+s
+        fprev = fopt
+        fopt=tc*xopt
+        funccount = funccount + 1
+        //
+        // Compute the duality gap
+        dualgap = abs(yopt'*beq - fopt)
+        if ( firstloop ) then
+            dualgapmin = dualgap
+        else
+            if ( dualgapmin > dualgap ) then
+                dualgapmin = dualgap
+            end
+        end
+        firstloop = %f
+    end
+    //
+    if ( step == 1 ) then
+        state = "init"
+    else
+        state = "done"
+    end
+    if ( step == 1 ) then
+        procedure = "x0"
+    else
+        procedure = "x*"
+    end
+    optimValues = struct(...
+    "funccount" , funccount , ...
+    "fval" , fopt , ...
+    "iteration" , iter , ...
+    "procedure" , procedure, ...
+    "dualgap" , dualgap ...
+    );
+    stop = karmarkar_outfunDriver ( xopt , optimValues , state , outfun )
+endfunction
+
+function stop = karmarkar_outfunDriver ( xopt , optimValues , state , outfun )
+    //
+    // The driver for the output function.
+    // outfun : a list where the first item is the output function.
+    __karmarkar_outfun__ = outfun (1) 
+    pinit = outfun (2)
+    newposvars = outfun (3)
+    //
+    // Calls back the user's output function, if required.
+    // Reduce the size of the vectors, to take into account for potential slack variables.
+    if ( __karmarkar_outfun__ <> [] ) then
+        if ( newposvars ) then
+            xopt = xopt(1:pinit) - xopt(pinit+1:2*pinit)
+        else
+            xopt = xopt(1:pinit)
+        end
+        cbktype = typeof( __karmarkar_outfun__ )
+        if ( cbktype == "list" ) then
+            __karmarkar_outfun__f_ = __karmarkar_outfun__ (1)
+            stop = __karmarkar_outfun__f_ ( xopt , optimValues , state , __karmarkar_outfun__ (2:$))
+        elseif ( or(cbktype == ["function" "fptr"] ) ) then
+            stop = __karmarkar_outfun__ ( xopt , optimValues , state )
+        end
+    else
+        stop = %f
+    end
 endfunction
+
+
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8719.dia.ref b/scilab/modules/optimization/tests/nonreg_tests/bug_8719.dia.ref
new file mode 100644 (file)
index 0000000..ed80d6d
--- /dev/null
@@ -0,0 +1,44 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8719 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8719
+//
+// <-- Short Description -->
+// The karmarkar function prints unwanted messages.
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+c = [-20 -24 0 0]';
+a = [
+3 6 1 0
+4 2 0 1
+];
+b = [60 32]';
+expected = [4 8 0 0]';
+x0 = [
+4.1128205
+7.7333333
+1.2615385
+0.0820513
+];
+xopt=karmarkar(a,b,c,x0);
+assert_close ( xopt , expected , 1.e-3 );
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8719.tst b/scilab/modules/optimization/tests/nonreg_tests/bug_8719.tst
new file mode 100644 (file)
index 0000000..1cd14eb
--- /dev/null
@@ -0,0 +1,47 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8719 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8719
+//
+// <-- Short Description -->
+// The karmarkar function prints unwanted messages.
+//
+
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+
+c = [-20 -24 0 0]';
+a = [
+3 6 1 0
+4 2 0 1
+];
+b = [60 32]';
+expected = [4 8 0 0]';
+x0 = [
+4.1128205
+7.7333333
+1.2615385
+0.0820513
+];
+xopt=karmarkar(a,b,c,x0);
+assert_close ( xopt , expected , 1.e-3 );
+
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8720.dia.ref b/scilab/modules/optimization/tests/nonreg_tests/bug_8720.dia.ref
new file mode 100644 (file)
index 0000000..510d98e
--- /dev/null
@@ -0,0 +1,39 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8720 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8720
+//
+// <-- Short Description -->
+// The karmarkar function may stop too early in the iterations.
+//
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+c = 1.e-20 * [-20 -24 0 0]';
+a = [
+3 6 1 0
+4 2 0 1
+];
+b = [60 32]';
+x0 = [4.1128205  7.7333333  1.2615385  0.0820513]';
+expected = [4 8 0 0]';
+xopt = karmarkar(a,b,c,x0);
+assert_close ( xopt , expected , 1.e-4 );
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8720.tst b/scilab/modules/optimization/tests/nonreg_tests/bug_8720.tst
new file mode 100644 (file)
index 0000000..02276bf
--- /dev/null
@@ -0,0 +1,42 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8720 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8720
+//
+// <-- Short Description -->
+// The karmarkar function may stop too early in the iterations.
+//
+
+function flag = assert_close ( computed, expected, epsilon )
+  if expected==0.0 then
+    shift = norm(computed-expected);
+  else
+    shift = norm(computed-expected)/norm(expected);
+  end
+  if shift < epsilon then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+
+c = 1.e-20 * [-20 -24 0 0]';
+a = [
+3 6 1 0
+4 2 0 1
+];
+b = [60 32]';
+x0 = [4.1128205  7.7333333  1.2615385  0.0820513]';
+expected = [4 8 0 0]';
+xopt = karmarkar(a,b,c,x0);
+assert_close ( xopt , expected , 1.e-4 );
+
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8726.dia.ref b/scilab/modules/optimization/tests/nonreg_tests/bug_8726.dia.ref
new file mode 100644 (file)
index 0000000..2fc6bac
--- /dev/null
@@ -0,0 +1,34 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8726 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8726
+//
+// <-- Short Description -->
+// The karmarkar function may produce a division-by-zero error.
+//
+function flag = MY_assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+Aeq = [
+1 -1 0
+1  1 1
+];
+beq = [0;2];
+c = [-1;-1;0];
+x0 = [1;1;0];
+xopt=karmarkar(Aeq,beq,c,x0);
+xexpected = [1;1;0];
+MY_assert_equal ( xopt , xexpected );
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8726.tst b/scilab/modules/optimization/tests/nonreg_tests/bug_8726.tst
new file mode 100644 (file)
index 0000000..e1e15ab
--- /dev/null
@@ -0,0 +1,39 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8726 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8726
+//
+// <-- Short Description -->
+// The karmarkar function may produce a division-by-zero error.
+//
+
+function flag = MY_assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+
+
+Aeq = [
+1 -1 0
+1  1 1
+];
+beq = [0;2];
+c = [-1;-1;0];
+x0 = [1;1;0];
+xopt=karmarkar(Aeq,beq,c,x0);
+xexpected = [1;1;0];
+MY_assert_equal ( xopt , xexpected );
+
+
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8727.dia.ref b/scilab/modules/optimization/tests/nonreg_tests/bug_8727.dia.ref
new file mode 100644 (file)
index 0000000..a087b86
--- /dev/null
@@ -0,0 +1,65 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8727 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8727
+//
+// <-- Short Description -->
+// The karmarkar function requires the initial guess x0.
+//
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+    if ( expected == computed ) then
+      flag = 1
+      return
+    end
+    if ( expected==0.0 ) then
+        shift = norm(computed-expected);
+    else
+        shift = norm(computed-expected)/norm(expected);
+    end
+    if shift < epsilon then
+        flag = 1;
+    else
+        flag = 0;
+    end
+    if flag <> 1 then bugmes();quit;end
+endfunction
+n=11;
+Aeq = zeros(n,n+2);
+Aeq(:,1) = 2*linspace(0,1,n)';
+Aeq(:,2) = ones(11,1);
+Aeq(1:n,3:n+2) = eye(n,n);
+beq = 1 + linspace(0,1,n)'.^2;
+c=[-1;-1;zeros(n,1)];
+xopt=karmarkar(Aeq,beq,c);
+xstar  = [ 
+    0.5041940  
+    0.7457937  
+    0.2542063  
+    0.1633675  
+    0.0925287  
+    0.0416899  
+    0.0108511  
+    0.0000123  
+    0.0091735  
+    0.0383347  
+    0.0874959  
+    0.1566571  
+    0.2458183  
+];
+assert_close ( xopt , xstar , 1.e-6 );
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8727.tst b/scilab/modules/optimization/tests/nonreg_tests/bug_8727.tst
new file mode 100644 (file)
index 0000000..cececd5
--- /dev/null
@@ -0,0 +1,70 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8727 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8727
+//
+// <-- Short Description -->
+// The karmarkar function requires the initial guess x0.
+//
+
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+    if ( expected == computed ) then
+      flag = 1
+      return
+    end
+    if ( expected==0.0 ) then
+        shift = norm(computed-expected);
+    else
+        shift = norm(computed-expected)/norm(expected);
+    end
+    if shift < epsilon then
+        flag = 1;
+    else
+        flag = 0;
+    end
+    if flag <> 1 then pause,end
+endfunction
+
+
+n=11;
+Aeq = zeros(n,n+2);
+Aeq(:,1) = 2*linspace(0,1,n)';
+Aeq(:,2) = ones(11,1);
+Aeq(1:n,3:n+2) = eye(n,n);
+beq = 1 + linspace(0,1,n)'.^2;
+c=[-1;-1;zeros(n,1)];
+xopt=karmarkar(Aeq,beq,c);
+
+xstar  = [ 
+    0.5041940  
+    0.7457937  
+    0.2542063  
+    0.1633675  
+    0.0925287  
+    0.0416899  
+    0.0108511  
+    0.0000123  
+    0.0091735  
+    0.0383347  
+    0.0874959  
+    0.1566571  
+    0.2458183  
+];
+assert_close ( xopt , xstar , 1.e-6 );
+
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8775.dia.ref b/scilab/modules/optimization/tests/nonreg_tests/bug_8775.dia.ref
new file mode 100644 (file)
index 0000000..77299fb
--- /dev/null
@@ -0,0 +1,34 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8775 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8775
+//
+// <-- Short Description -->
+// The karmarkar function might diverge toward a non-optimal point.
+//
+function flag = MY_assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then bugmes();quit;end
+endfunction
+// An unbounded problem.
+Aeq = [
+ 2 -2 -1 1 0
+-1 -4  1 0 1
+];
+beq = [-1;-1];
+c = [2;9;3;0;0];
+x0 = [0.2;0.7;1;1;1];
+[xopt,fopt,exitflag]=karmarkar(Aeq,beq,c,x0,0,0.999);
+MY_assert_equal ( exitflag , -2 );
diff --git a/scilab/modules/optimization/tests/nonreg_tests/bug_8775.tst b/scilab/modules/optimization/tests/nonreg_tests/bug_8775.tst
new file mode 100644 (file)
index 0000000..81b86c8
--- /dev/null
@@ -0,0 +1,37 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// <-- Non-regression test for bug 8775 -->
+//
+// <-- Bugzilla URL -->
+// http://bugzilla.scilab.org/show_bug.cgi?id=8775
+//
+// <-- Short Description -->
+// The karmarkar function might diverge toward a non-optimal point.
+//
+
+function flag = MY_assert_equal ( computed , expected )
+  if computed==expected then
+    flag = 1;
+  else
+    flag = 0;
+  end
+  if flag <> 1 then pause,end
+endfunction
+
+// An unbounded problem.
+Aeq = [
+ 2 -2 -1 1 0
+-1 -4  1 0 1
+];
+beq = [-1;-1];
+c = [2;9;3;0;0];
+x0 = [0.2;0.7;1;1;1];
+[xopt,fopt,exitflag]=karmarkar(Aeq,beq,c,x0,0,0.999);
+MY_assert_equal ( exitflag , -2 );
+
diff --git a/scilab/modules/optimization/tests/unit_tests/karmarkar.dia.ref b/scilab/modules/optimization/tests/unit_tests/karmarkar.dia.ref
new file mode 100644 (file)
index 0000000..f5f42dd
--- /dev/null
@@ -0,0 +1,648 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+// <-- JVM NOT MANDATORY -->
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+    if ( expected == computed ) then
+      flag = 1
+      return
+    end
+    if ( expected==0.0 ) then
+        shift = norm(computed-expected);
+    else
+        shift = norm(computed-expected)/norm(expected);
+    end
+    if shift < epsilon then
+        flag = 1;
+    else
+        flag = 0;
+    end
+    if flag <> 1 then bugmes();quit;end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+    if computed==expected then
+        flag = 1;
+    else
+        flag = 0;
+    end
+    if flag <> 1 then bugmes();quit;end
+endfunction
+// With slack variables:
+//
+// Minimize -20.x1 - 24.x2 such as:
+// 3.x1 + 6.x2 + x3 = 60
+// 4.x1 + 2.x2 + x4 = 32
+// x >= 0
+c = [-20 -24 0 0]';
+Aeq = [
+3 6 1 0
+4 2 0 1
+];
+beq = [60 32]';
+xexpected = [4 8 0 0]';
+fexpected = -272;
+x0 = [
+4.1128205
+7.7333333
+1.2615385
+0.0820513
+];
+[xopt,fopt]=karmarkar(Aeq,beq,c,x0);
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+//
+// Configure the relative tolerance
+rtolf=1.e-6;
+[xopt,fopt]=karmarkar(Aeq,beq,c,x0,rtolf);
+assert_close ( xopt , xexpected , 1.e-4 );
+assert_close ( fopt , fexpected , 1.e-5 );
+//
+// Configure the gamma
+gam = 0.1;
+[xopt,fopt]=karmarkar(Aeq,beq,c,x0,[],gam);
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+////////////////////////////////////////////////////////////
+//
+// Check new API (from Scilab v5.3.x).
+//
+// Check exit flag
+[xopt,fopt,exitflag]=karmarkar(Aeq,beq,c,x0);
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+assert_equal ( exitflag , 1 );
+//
+// Check number of iterations
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,x0);
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter>10 , %t );
+//
+// Check dual solution
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0);
+lambda.ineqlin = [];
+lambda.eqlin = [28/9;8/3];
+lambda.upper = [0;0;0;0];
+lambda.lower = [0;0;28/9;8/3];
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter>10 , %t );
+assert_equal ( yopt.ineqlin , lambda.ineqlin );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-8 );
+assert_close ( yopt.lower , lambda.lower , 1.e-8 );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Check number of iterations, with default options
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,x0,[],[],10);
+assert_equal ( exitflag , 0 );
+assert_equal ( iter , 10 );
+//
+// Check output function
+function stop = myoutputfunction ( xopt , optimValues , state )
+    localmsg = gettext("Iteration #%3.0f, state=%s, procedure=%s, fopt=%10.3e, x=[%s], dualgap=%10.3e\n")
+    xstr = strcat(msprintf("%10.3e\n",xopt)'," ")
+    teststring = sprintf(localmsg,optimValues.iteration,state,optimValues.procedure,optimValues.fval,xstr,optimValues.dualgap)
+    stop = %f
+endfunction
+xopt=karmarkar(Aeq,beq,c,x0,[],[],[],myoutputfunction);
+assert_close ( xopt , xexpected , 1.e-3 );
+//
+// Check output function, without initial guess
+xopt=karmarkar(Aeq,beq,c,[],[],[],[],myoutputfunction);
+assert_close ( xopt , xexpected , 1.e-3 );
+//
+// Check that the output function can stop the algorithm
+function stop = myoutputfunctionStop ( xopt , optimValues , state )
+    stop = (iter >= 7)
+endfunction
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,x0,[],[],[],myoutputfunctionStop);
+assert_close ( xopt , xexpected , 1.e-2 );
+assert_close ( fopt , fexpected , 1.e-3 );
+assert_equal ( exitflag , -4 );
+assert_equal ( iter , 7 );
+//
+// Check output function with additionnal arguments
+function stop = myoutputfunction2 ( xopt , optimValues , state , myAeq , mybeq , myc )
+    localmsg = gettext("Iteration #%3.0f, fopt=%10.3e, state=%s, ||Ax-beq||=%.3e\n")
+    teststring = sprintf(localmsg,optimValues.iteration,optimValues.fval,state,norm(myAeq*xopt-mybeq))
+    stop = %f
+endfunction
+xopt=karmarkar(Aeq,beq,c,x0,[],[],[],list(myoutputfunction2,Aeq,beq,c));
+assert_close ( xopt , xexpected , 1.e-3 );
+//
+// References
+// "Practical Optimization", Antoniou, Lu, 2007
+// Chapter 11, "Linear Programming Part I: The simplex method",
+// Example 11.9, p. 361
+// Chapter 12, "Linear Programming Part II: Interior point methods",
+// Example 12.2, p.382
+// 
+// Minimize 2.x1 + 9.x2 + 3.x3
+// -2.x1 + 2.x2 + x3 - x4 = 1
+//    x1 + 4.x2 - x3 - x5 = 1
+// x >= 0
+Aeq = [
+-2 2 1 -1 0
+1 4 -1 0 -1
+];
+beq = [1;1];
+c = [2;9;3;0;0];
+x0 = [0.2;0.7;1;1;1];
+gam = 0.9999;
+rtolf = 1.e-4;
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0,rtolf,gam,[],myoutputfunction);
+xstar = [0 1/3 1/3 0 0]';
+fstar = 4;
+lambda.ineqlin = [];
+lambda.eqlin = [-3.5;-0.5];
+lambda.upper = [0;0;0;0;0];
+lambda.lower = [8.5;0;0;3.5;0.5];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter , 4 );
+assert_equal ( yopt.ineqlin , lambda.ineqlin );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_close ( yopt.lower , lambda.lower , 1.e-6 );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Minimize -x1 -x2
+// x1 - x2      = 0
+// x1 + x2 + x3 = 2
+// x >= 0
+//
+// Let karmarkar find a feasible x0.
+Aeq = [
+1 -1 0
+1  1 1
+];
+beq = [0;2];
+c = [-1;-1;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c);
+xstar = [1;1;0];
+fstar = -2;
+lambda.ineqlin = [];
+lambda.eqlin = [0;1];
+lambda.upper = [0;0;0];
+lambda.lower = [0;0;1];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_equal ( yopt.ineqlin , lambda.ineqlin );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_close ( yopt.lower , lambda.lower , 1.e-6 );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Give a linear inequality A*x <= b.
+//
+// Minimize -x1 -x2
+// x1 - x2  = 0
+// x1 + x2 <= 2
+//
+// Give x0.
+Aeq = [
+1 -1
+];
+beq = 0;
+c = [-1;-1];
+A = [1 1];
+b = 2;
+x0 = [0.1;0.1];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0,[],[],[],[],A,b);
+xstar=[1 1]';
+fstar = c'*xstar;
+lambda.ineqlin = 1;
+lambda.eqlin = 0;
+lambda.upper = [0;0];
+lambda.lower = [0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-6 );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Give a linear inequality A*x <= b
+//
+// Minimize -x1 -x2
+// x1 - x2  = 0
+// x1 + x2 <= 2
+//
+// Do not give x0.
+Aeq = [
+1 -1
+];
+beq = 0;
+c = [-1;-1];
+A = [1 1];
+b = 2;
+function stop = myoutputfunction2 ( xopt , optimValues , state )
+    assert_equal ( size(xopt) , [2 1] );
+    assert_equal ( or(state==["init","iter","done"]) , %t );
+    assert_equal ( or(optimValues.procedure==["x0","x*"]) , %t );
+    stop = %f
+endfunction
+Attention: redéfinition de la fonction: myoutputfunction2       . Utilisez funcprot(0) pour éviter ce message
+
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],myoutputfunction2,A,b);
+xstar=[1 1]';
+fstar = c'*xstar;
+lambda.ineqlin = 1;
+lambda.eqlin = 0;
+lambda.upper = [0;0];
+lambda.lower = [0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-6 );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Minimize -20.x1 - 24.x2 such as:
+// 3.x1 + 6.x2 >= 60
+// 4.x1 + 2.x2 >= 32
+c = [-20 -24]';
+Aeq=[];
+beq=[];
+A = [
+3 6
+4 2
+];
+b = [60 32]';
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b);
+xstar = [4 8]';
+fstar = c'*xstar;
+lambda.ineqlin = [28/9;8/3];
+lambda.eqlin = [];
+lambda.upper = [0;0];
+lambda.lower = [0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-6 );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// References
+// "Practical Optimization", Antoniou, Lu, 2007
+// Chapter 11, "Linear Programming Part I: The simplex method",
+// Example 11.9, p. 361
+// Chapter 12, "Linear Programming Part II: Interior point methods",
+// Example 12.2, p.382
+// 
+// Minimize 2.x1 + 9.x2 + 3.x3
+// +2.x1 - 2.x2 - x3 <= -1
+//   -x1 - 4.x2 + x3 <= -1
+// x >= 0
+//
+// Give x0
+Aeq=[];
+beq=[];
+A = [
+ 2 -2 -1
+-1 -4  1
+-1  0  0
+ 0 -1  0
+ 0  0 -1
+];
+b = [-1;-1;0;0;0];
+c = [2;9;3];
+x0 = [0.2;0.7;1];
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,x0,[],[],[],[],A,b);
+xstar = [0 1/3 1/3]';
+fstar = c'*xstar;
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+//
+// Let x0 be found by the algorithm.
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b);
+xstar = [0 1/3 1/3]';
+fstar = c'*xstar;
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+/////////////////////////////////////////////////////
+// References
+// "Lipsol toolbox", Yin Zhang, Scilab port by Rubio Scola, example0.sce.
+//
+// Minimize 2.x1 + 5.x2 - 2.5.x3
+//    x1        S4 x3 <= 5
+// E2 x1 - x2    - x3 <= 0
+//    x1              <= 2
+//         x2         <= %inf
+//                 x3 <= 3
+//  - x1              <= 2
+//       - x2         <= -1
+//               - x3 <= 0
+// 
+// where
+// S4 = sin(pi/4)/4
+// E2 = exp(2)
+//   
+S4 = sin(%pi/4)/4;
+E2 = exp(2);
+c = [ 2; 5; -2.5];
+A = [ 
+       1  0 S4
+      E2 -1 -1
+       1  0  0
+       0  0  1
+      -1  0  0
+       0 -1  0
+       0  0 -1
+];
+b = [ 5; 0;2;3;2;-1;0];
+Aeq = [];
+beq = [];
+xstar = [-2;1;3];
+fstar = c'*xstar;
+lambda.ineqlin = [0;0;0;2.5;2;5;0];
+lambda.eqlin = [];
+lambda.upper = [0;0;0];
+lambda.lower = [0;0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b);
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_equal ( yopt.eqlin , lambda.eqlin );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Minimize 2.x1 + 9.x2 + 3.x3
+//  2 x1 - 2 x2 - x3 <= -1
+// -1 x1 - 4 x2 + x3 <= -1
+// - x1              <= -1
+// - x2              <= 0
+// - x3              <= 0
+//
+Aeq=[];
+beq=[];
+A = [
+ 2 -2 -1
+-1 -4  1
+-1  0  0
+ 0 -1  0
+ 0  0 -1
+];
+b = [-1;-1;-1;0;0];
+c = [2;9;3];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b);
+xstar = [1;0.5;2];
+fstar = c'*xstar;
+lambda.ineqlin = [3.5;0.5;8.5;0;0];
+lambda.eqlin = [];
+lambda.upper = [0;0;0];
+lambda.lower = [0;0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_equal ( yopt.eqlin , lambda.eqlin );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+//////////////////////////////////////////////////////////////////////////
+//
+// Set lower bound and do not give x0.
+Aeq=[];
+beq=[];
+A = [
+ 2 -2 -1
+-1 -4  1
+];
+b = [-1;-1];
+c = [2;9;3];
+lb = [1;0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b,lb);
+xstar = [1;0.5;2];
+fstar = c'*xstar;
+lambda.ineqlin = [3.5;0.5];
+lambda.eqlin = [];
+lambda.upper = [0;0;0];
+lambda.lower = [8.5;0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_equal ( yopt.eqlin , lambda.eqlin );
+assert_close ( yopt.lower , lambda.lower , 1.e-9 );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Set lower bound and give x0.
+Aeq=[];
+beq=[];
+A = [
+ 2 -2 -1
+-1 -4  1
+];
+b = [-1;-1];
+c = [2;9;3];
+x0 = [1.337848;0.885225;2.535279];
+lb = [1;0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0,[],[],[],[],A,b,lb);
+xstar = [1;0.5;2];
+fstar = c'*xstar;
+lambda.ineqlin = [3.5;0.5];
+lambda.eqlin = [];
+lambda.upper = [0;0;0];
+lambda.lower = [8.5;0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_close ( yopt.eqlin , lambda.eqlin , %eps );
+assert_close ( yopt.lower , lambda.lower , 1.e-9 );
+assert_close ( yopt.upper , lambda.upper , 1.e-9 );
+// References
+// LIPSOL is a set of Linear-programming Interior-Point SOLvers written 
+// by Yin Zhang. 
+// The original Matlab-based code has been adapted to Scilab 
+// by H. Rubio Scola. 
+//
+// Minimize 2 x1 + 5 x2 - 2.5 x3
+//    x1 +   S4 x3 <= 5
+// E2 x1 - x2 - x3 <= 0
+//  -2 <= x1 <= 2
+//   1 <= x2
+//   0 <= x3 <= 3
+// where:
+// S4 = sin(pi/4)/4
+// E2 = exp(2)
+//
+// Do not give x0.
+c = [ 2; 5; -2.5];
+S4 = sin(%pi/4)/4;
+E2 = exp(2);
+A = [ 
+   1  0 S4
+  E2 -1 -1
+];
+b = [ 5; 0];
+lb = [ -2; 1   ; 0 ];
+ub = [  2; %inf; 3 ];
+Aeq = [];
+beq = [];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b,lb,ub);
+xstar = [-2;1;3];
+fstar = c'*xstar;
+lambda.ineqlin = [0;0];
+lambda.eqlin = [];
+lambda.upper = [0;0;2.5];
+lambda.lower = [2;5;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_close ( yopt.eqlin , lambda.eqlin , %eps );
+assert_close ( yopt.lower , lambda.lower , 1.e-9 );
+assert_close ( yopt.upper , lambda.upper , 1.e-9 );
+//
+// An unbounded problem.
+c = [-20 -24]';
+A = [
+-3 -6
+-4 -2
+];
+b = [-60 -32]';
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b);
+assert_equal ( exitflag , -2 );
+//
+// "Linear Programming in Matlab"
+// Ferris, Mangasarian, Wright
+// 2008
+// Chapter 3, "The Simplex Method", Exercise 3-4-2 1.
+//
+// An infeasible problem.
+// Minimize -3 x1 + x2
+//  - x1 -   x2 >= -2
+//  2 x1 + 2 x2 >= 10
+// x >= 0
+c = [-3;1];
+A=[
+ -1 -1
+  2  2
+];
+A=-A;
+b=[-2;10];
+b=-b;
+lb=[0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb);
+assert_equal ( xopt , [] );
+assert_equal ( fopt , [] );
+assert_equal ( exitflag , -1 );
+assert_equal ( iter > 0 , %t );
+assert_equal ( yopt.ineqlin , [] );
+assert_equal ( yopt.eqlin , [] );
+assert_equal ( yopt.lower , [] );
+assert_equal ( yopt.upper , [] );
+//
+// "Linear Programming in Matlab"
+// Ferris, Mangasarian, Wright
+// 2008
+// Chapter 3, "The Simplex Method", Exercise 3-4-2 2.
+//
+// An unbounded problem.
+// Minimize -x1 + x2
+//  2 x1 -   x2 >= 1
+//    x1 + 2 x2 >= 2
+// x >= 0
+c = [-1;1];
+A=[
+ 2 -1
+ 1  2
+];
+A=-A;
+b=[1;2];
+b=-b;
+lb=[0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb);
+assert_equal ( exitflag , -2 );
+//
+// "Linear and Nonlinear Optimization"
+// Griva, Nash, Sofer
+// 2009
+// Chapter 5, "The Simplex Method", Example 5.3
+//
+// An unbounded problem.
+// Minimize -x1 - 2 x2
+//  -1 x1 +   x2 <= 2
+//  -2 x1 +   x2 <= 1
+// x >= 0
+c = [-1;-2];
+A=[
+ -1  1
+ -2  1
+];
+b=[2;1];
+lb=[0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb);
+assert_equal ( exitflag , -2 );
+//
+// "Linear and Nonlinear Optimization"
+// Griva, Nash, Sofer
+// 2009
+// Chapter 5, "The Simplex Method", Example 5.6
+//
+// An unfeasible problem.
+// Minimize -x1
+//  -  x1 -   x2 <= -6
+//   2 x1 + 3 x2 <=  4
+// x >= 0
+c = [-1;0];
+A=[
+ -1  -1
+  2   3
+];
+b=[-6;4];
+lb=[0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb);
+assert_equal ( exitflag , -1 );
+//
+// Example from the help page
+// Check that the output points remain feasible
+function stop = myoutputfunction3 ( xopt , optimValues , state , A , b )
+    assert_equal ( and(A*xopt<=b) , %t );
+    stop = %f
+endfunction
+n=11;
+A = [2*linspace(0,1,n)',ones(n,1)];
+b = 1 + linspace(0,1,n)'.^2;
+c=[-1;-1];
+xopt=karmarkar([],[],c,[],[],[],[],list(myoutputfunction3,A,b),A,b);
+xstar = [0.5005127;0.7494803];
+assert_close ( xopt , xstar , 1.e-4 );
diff --git a/scilab/modules/optimization/tests/unit_tests/karmarkar.tst b/scilab/modules/optimization/tests/unit_tests/karmarkar.tst
new file mode 100644 (file)
index 0000000..327caa5
--- /dev/null
@@ -0,0 +1,658 @@
+// =============================================================================
+// Scilab ( http://www.scilab.org/ ) - This file is part of Scilab
+// Copyright (C) 2010 - DIGITEO - Michael Baudin
+//
+//  This file is distributed under the same license as the Scilab package.
+// =============================================================================
+
+// <-- JVM NOT MANDATORY -->
+
+//
+// assert_close --
+//   Returns 1 if the two real matrices computed and expected are close,
+//   i.e. if the relative distance between computed and expected is lesser than epsilon.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_close ( computed, expected, epsilon )
+    if ( expected == computed ) then
+      flag = 1
+      return
+    end
+    if ( expected==0.0 ) then
+        shift = norm(computed-expected);
+    else
+        shift = norm(computed-expected)/norm(expected);
+    end
+    if shift < epsilon then
+        flag = 1;
+    else
+        flag = 0;
+    end
+    if flag <> 1 then pause,end
+endfunction
+//
+// assert_equal --
+//   Returns 1 if the two real matrices computed and expected are equal.
+// Arguments
+//   computed, expected : the two matrices to compare
+//   epsilon : a small number
+//
+function flag = assert_equal ( computed , expected )
+    if computed==expected then
+        flag = 1;
+    else
+        flag = 0;
+    end
+    if flag <> 1 then pause,end
+endfunction
+
+
+// With slack variables:
+//
+// Minimize -20.x1 - 24.x2 such as:
+// 3.x1 + 6.x2 + x3 = 60
+// 4.x1 + 2.x2 + x4 = 32
+// x >= 0
+c = [-20 -24 0 0]';
+Aeq = [
+3 6 1 0
+4 2 0 1
+];
+beq = [60 32]';
+xexpected = [4 8 0 0]';
+fexpected = -272;
+x0 = [
+4.1128205
+7.7333333
+1.2615385
+0.0820513
+];
+[xopt,fopt]=karmarkar(Aeq,beq,c,x0);
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+//
+// Configure the relative tolerance
+rtolf=1.e-6;
+[xopt,fopt]=karmarkar(Aeq,beq,c,x0,rtolf);
+assert_close ( xopt , xexpected , 1.e-4 );
+assert_close ( fopt , fexpected , 1.e-5 );
+//
+// Configure the gamma
+gam = 0.1;
+[xopt,fopt]=karmarkar(Aeq,beq,c,x0,[],gam);
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+////////////////////////////////////////////////////////////
+//
+// Check new API (from Scilab v5.3.x).
+//
+// Check exit flag
+[xopt,fopt,exitflag]=karmarkar(Aeq,beq,c,x0);
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+assert_equal ( exitflag , 1 );
+//
+// Check number of iterations
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,x0);
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter>10 , %t );
+//
+// Check dual solution
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0);
+lambda.ineqlin = [];
+lambda.eqlin = [28/9;8/3];
+lambda.upper = [0;0;0;0];
+lambda.lower = [0;0;28/9;8/3];
+assert_close ( xopt , xexpected , 1.e-3 );
+assert_close ( fopt , fexpected , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter>10 , %t );
+assert_equal ( yopt.ineqlin , lambda.ineqlin );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-8 );
+assert_close ( yopt.lower , lambda.lower , 1.e-8 );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Check number of iterations, with default options
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,x0,[],[],10);
+assert_equal ( exitflag , 0 );
+assert_equal ( iter , 10 );
+//
+// Check output function
+function stop = myoutputfunction ( xopt , optimValues , state )
+    localmsg = gettext("Iteration #%3.0f, state=%s, procedure=%s, fopt=%10.3e, x=[%s], dualgap=%10.3e\n")
+    xstr = strcat(msprintf("%10.3e\n",xopt)'," ")
+    teststring = sprintf(localmsg,optimValues.iteration,state,optimValues.procedure,optimValues.fval,xstr,optimValues.dualgap)
+    stop = %f
+endfunction
+xopt=karmarkar(Aeq,beq,c,x0,[],[],[],myoutputfunction);
+assert_close ( xopt , xexpected , 1.e-3 );
+//
+// Check output function, without initial guess
+xopt=karmarkar(Aeq,beq,c,[],[],[],[],myoutputfunction);
+assert_close ( xopt , xexpected , 1.e-3 );
+//
+// Check that the output function can stop the algorithm
+function stop = myoutputfunctionStop ( xopt , optimValues , state )
+    stop = (iter >= 7)
+endfunction
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,x0,[],[],[],myoutputfunctionStop);
+assert_close ( xopt , xexpected , 1.e-2 );
+assert_close ( fopt , fexpected , 1.e-3 );
+assert_equal ( exitflag , -4 );
+assert_equal ( iter , 7 );
+//
+// Check output function with additionnal arguments
+function stop = myoutputfunction2 ( xopt , optimValues , state , myAeq , mybeq , myc )
+    localmsg = gettext("Iteration #%3.0f, fopt=%10.3e, state=%s, ||Ax-beq||=%.3e\n")
+    teststring = sprintf(localmsg,optimValues.iteration,optimValues.fval,state,norm(myAeq*xopt-mybeq))
+    stop = %f
+endfunction
+xopt=karmarkar(Aeq,beq,c,x0,[],[],[],list(myoutputfunction2,Aeq,beq,c));
+assert_close ( xopt , xexpected , 1.e-3 );
+//
+// References
+// "Practical Optimization", Antoniou, Lu, 2007
+// Chapter 11, "Linear Programming Part I: The simplex method",
+// Example 11.9, p. 361
+// Chapter 12, "Linear Programming Part II: Interior point methods",
+// Example 12.2, p.382
+// 
+// Minimize 2.x1 + 9.x2 + 3.x3
+// -2.x1 + 2.x2 + x3 - x4 = 1
+//    x1 + 4.x2 - x3 - x5 = 1
+// x >= 0
+Aeq = [
+-2 2 1 -1 0
+1 4 -1 0 -1
+];
+beq = [1;1];
+c = [2;9;3;0;0];
+x0 = [0.2;0.7;1;1;1];
+gam = 0.9999;
+rtolf = 1.e-4;
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0,rtolf,gam,[],myoutputfunction);
+xstar = [0 1/3 1/3 0 0]';
+fstar = 4;
+lambda.ineqlin = [];
+lambda.eqlin = [-3.5;-0.5];
+lambda.upper = [0;0;0;0;0];
+lambda.lower = [8.5;0;0;3.5;0.5];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter , 4 );
+assert_equal ( yopt.ineqlin , lambda.ineqlin );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_close ( yopt.lower , lambda.lower , 1.e-6 );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Minimize -x1 -x2
+// x1 - x2      = 0
+// x1 + x2 + x3 = 2
+// x >= 0
+//
+// Let karmarkar find a feasible x0.
+Aeq = [
+1 -1 0
+1  1 1
+];
+beq = [0;2];
+c = [-1;-1;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c);
+xstar = [1;1;0];
+fstar = -2;
+lambda.ineqlin = [];
+lambda.eqlin = [0;1];
+lambda.upper = [0;0;0];
+lambda.lower = [0;0;1];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_equal ( yopt.ineqlin , lambda.ineqlin );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_close ( yopt.lower , lambda.lower , 1.e-6 );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Give a linear inequality A*x <= b.
+//
+// Minimize -x1 -x2
+// x1 - x2  = 0
+// x1 + x2 <= 2
+//
+// Give x0.
+Aeq = [
+1 -1
+];
+beq = 0;
+c = [-1;-1];
+A = [1 1];
+b = 2;
+x0 = [0.1;0.1];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0,[],[],[],[],A,b);
+xstar=[1 1]';
+fstar = c'*xstar;
+lambda.ineqlin = 1;
+lambda.eqlin = 0;
+lambda.upper = [0;0];
+lambda.lower = [0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-6 );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Give a linear inequality A*x <= b
+//
+// Minimize -x1 -x2
+// x1 - x2  = 0
+// x1 + x2 <= 2
+//
+// Do not give x0.
+Aeq = [
+1 -1
+];
+beq = 0;
+c = [-1;-1];
+A = [1 1];
+b = 2;
+function stop = myoutputfunction2 ( xopt , optimValues , state )
+    assert_equal ( size(xopt) , [2 1] );
+    assert_equal ( or(state==["init","iter","done"]) , %t );
+    assert_equal ( or(optimValues.procedure==["x0","x*"]) , %t );
+    stop = %f
+endfunction
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],myoutputfunction2,A,b);
+xstar=[1 1]';
+fstar = c'*xstar;
+lambda.ineqlin = 1;
+lambda.eqlin = 0;
+lambda.upper = [0;0];
+lambda.lower = [0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-6 );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Minimize -20.x1 - 24.x2 such as:
+// 3.x1 + 6.x2 >= 60
+// 4.x1 + 2.x2 >= 32
+c = [-20 -24]';
+Aeq=[];
+beq=[];
+A = [
+3 6
+4 2
+];
+b = [60 32]';
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b);
+xstar = [4 8]';
+fstar = c'*xstar;
+lambda.ineqlin = [28/9;8/3];
+lambda.eqlin = [];
+lambda.upper = [0;0];
+lambda.lower = [0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-6 );
+assert_close ( yopt.eqlin , lambda.eqlin , 1.e-6 );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// References
+// "Practical Optimization", Antoniou, Lu, 2007
+// Chapter 11, "Linear Programming Part I: The simplex method",
+// Example 11.9, p. 361
+// Chapter 12, "Linear Programming Part II: Interior point methods",
+// Example 12.2, p.382
+// 
+// Minimize 2.x1 + 9.x2 + 3.x3
+// +2.x1 - 2.x2 - x3 <= -1
+//   -x1 - 4.x2 + x3 <= -1
+// x >= 0
+//
+// Give x0
+Aeq=[];
+beq=[];
+A = [
+ 2 -2 -1
+-1 -4  1
+-1  0  0
+ 0 -1  0
+ 0  0 -1
+];
+b = [-1;-1;0;0;0];
+c = [2;9;3];
+x0 = [0.2;0.7;1];
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,x0,[],[],[],[],A,b);
+xstar = [0 1/3 1/3]';
+fstar = c'*xstar;
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+//
+// Let x0 be found by the algorithm.
+[xopt,fopt,exitflag,iter]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b);
+xstar = [0 1/3 1/3]';
+fstar = c'*xstar;
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+/////////////////////////////////////////////////////
+// References
+// "Lipsol toolbox", Yin Zhang, Scilab port by Rubio Scola, example0.sce.
+//
+// Minimize 2.x1 + 5.x2 - 2.5.x3
+//    x1        S4 x3 <= 5
+// E2 x1 - x2    - x3 <= 0
+//    x1              <= 2
+//         x2         <= %inf
+//                 x3 <= 3
+//  - x1              <= 2
+//       - x2         <= -1
+//               - x3 <= 0
+// 
+// where
+// S4 = sin(pi/4)/4
+// E2 = exp(2)
+//   
+
+S4 = sin(%pi/4)/4;
+E2 = exp(2);
+c = [ 2; 5; -2.5];
+A = [ 
+       1  0 S4
+      E2 -1 -1
+       1  0  0
+       0  0  1
+      -1  0  0
+       0 -1  0
+       0  0 -1
+];
+b = [ 5; 0;2;3;2;-1;0];
+Aeq = [];
+beq = [];
+xstar = [-2;1;3];
+fstar = c'*xstar;
+lambda.ineqlin = [0;0;0;2.5;2;5;0];
+lambda.eqlin = [];
+lambda.upper = [0;0;0];
+lambda.lower = [0;0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b);
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_equal ( yopt.eqlin , lambda.eqlin );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+
+//
+// Minimize 2.x1 + 9.x2 + 3.x3
+//  2 x1 - 2 x2 - x3 <= -1
+// -1 x1 - 4 x2 + x3 <= -1
+// - x1              <= -1
+// - x2              <= 0
+// - x3              <= 0
+//
+Aeq=[];
+beq=[];
+A = [
+ 2 -2 -1
+-1 -4  1
+-1  0  0
+ 0 -1  0
+ 0  0 -1
+];
+b = [-1;-1;-1;0;0];
+c = [2;9;3];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b);
+xstar = [1;0.5;2];
+fstar = c'*xstar;
+lambda.ineqlin = [3.5;0.5;8.5;0;0];
+lambda.eqlin = [];
+lambda.upper = [0;0;0];
+lambda.lower = [0;0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_equal ( yopt.eqlin , lambda.eqlin );
+assert_equal ( yopt.lower , lambda.lower );
+assert_equal ( yopt.upper , lambda.upper );
+
+//////////////////////////////////////////////////////////////////////////
+//
+// Set lower bound and do not give x0.
+Aeq=[];
+beq=[];
+A = [
+ 2 -2 -1
+-1 -4  1
+];
+b = [-1;-1];
+c = [2;9;3];
+lb = [1;0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b,lb);
+xstar = [1;0.5;2];
+fstar = c'*xstar;
+lambda.ineqlin = [3.5;0.5];
+lambda.eqlin = [];
+lambda.upper = [0;0;0];
+lambda.lower = [8.5;0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_equal ( yopt.eqlin , lambda.eqlin );
+assert_close ( yopt.lower , lambda.lower , 1.e-9 );
+assert_equal ( yopt.upper , lambda.upper );
+//
+// Set lower bound and give x0.
+Aeq=[];
+beq=[];
+A = [
+ 2 -2 -1
+-1 -4  1
+];
+b = [-1;-1];
+c = [2;9;3];
+x0 = [1.337848;0.885225;2.535279];
+lb = [1;0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,x0,[],[],[],[],A,b,lb);
+xstar = [1;0.5;2];
+fstar = c'*xstar;
+lambda.ineqlin = [3.5;0.5];
+lambda.eqlin = [];
+lambda.upper = [0;0;0];
+lambda.lower = [8.5;0;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_close ( yopt.eqlin , lambda.eqlin , %eps );
+assert_close ( yopt.lower , lambda.lower , 1.e-9 );
+assert_close ( yopt.upper , lambda.upper , 1.e-9 );
+// References
+// LIPSOL is a set of Linear-programming Interior-Point SOLvers written 
+// by Yin Zhang. 
+// The original Matlab-based code has been adapted to Scilab 
+// by H. Rubio Scola. 
+//
+// Minimize 2 x1 + 5 x2 - 2.5 x3
+//    x1 +   S4 x3 <= 5
+// E2 x1 - x2 - x3 <= 0
+//  -2 <= x1 <= 2
+//   1 <= x2
+//   0 <= x3 <= 3
+// where:
+// S4 = sin(pi/4)/4
+// E2 = exp(2)
+//
+// Do not give x0.
+c = [ 2; 5; -2.5];
+S4 = sin(%pi/4)/4;
+E2 = exp(2);
+A = [ 
+   1  0 S4
+  E2 -1 -1
+];
+b = [ 5; 0];
+lb = [ -2; 1   ; 0 ];
+ub = [  2; %inf; 3 ];
+Aeq = [];
+beq = [];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar(Aeq,beq,c,[],[],[],[],[],A,b,lb,ub);
+xstar = [-2;1;3];
+fstar = c'*xstar;
+lambda.ineqlin = [0;0];
+lambda.eqlin = [];
+lambda.upper = [0;0;2.5];
+lambda.lower = [2;5;0];
+assert_close ( xopt , xstar , 1.e-4 );
+assert_close ( fopt , fstar , 1.e-4 );
+assert_equal ( exitflag , 1 );
+assert_equal ( iter > 0 , %t );
+assert_close ( yopt.ineqlin , lambda.ineqlin , 1.e-9 );
+assert_close ( yopt.eqlin , lambda.eqlin , %eps );
+assert_close ( yopt.lower , lambda.lower , 1.e-9 );
+assert_close ( yopt.upper , lambda.upper , 1.e-9 );
+//
+// An unbounded problem.
+c = [-20 -24]';
+A = [
+-3 -6
+-4 -2
+];
+b = [-60 -32]';
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b);
+assert_equal ( exitflag , -2 );
+//
+// "Linear Programming in Matlab"
+// Ferris, Mangasarian, Wright
+// 2008
+// Chapter 3, "The Simplex Method", Exercise 3-4-2 1.
+//
+// An infeasible problem.
+// Minimize -3 x1 + x2
+//  - x1 -   x2 >= -2
+//  2 x1 + 2 x2 >= 10
+// x >= 0
+c = [-3;1];
+A=[
+ -1 -1
+  2  2
+];
+A=-A;
+b=[-2;10];
+b=-b;
+lb=[0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb);
+assert_equal ( xopt , [] );
+assert_equal ( fopt , [] );
+assert_equal ( exitflag , -1 );
+assert_equal ( iter > 0 , %t );
+assert_equal ( yopt.ineqlin , [] );
+assert_equal ( yopt.eqlin , [] );
+assert_equal ( yopt.lower , [] );
+assert_equal ( yopt.upper , [] );
+
+
+//
+// "Linear Programming in Matlab"
+// Ferris, Mangasarian, Wright
+// 2008
+// Chapter 3, "The Simplex Method", Exercise 3-4-2 2.
+//
+// An unbounded problem.
+// Minimize -x1 + x2
+//  2 x1 -   x2 >= 1
+//    x1 + 2 x2 >= 2
+// x >= 0
+c = [-1;1];
+A=[
+ 2 -1
+ 1  2
+];
+A=-A;
+b=[1;2];
+b=-b;
+lb=[0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb);
+assert_equal ( exitflag , -2 );
+
+//
+// "Linear and Nonlinear Optimization"
+// Griva, Nash, Sofer
+// 2009
+// Chapter 5, "The Simplex Method", Example 5.3
+//
+// An unbounded problem.
+// Minimize -x1 - 2 x2
+//  -1 x1 +   x2 <= 2
+//  -2 x1 +   x2 <= 1
+// x >= 0
+c = [-1;-2];
+A=[
+ -1  1
+ -2  1
+];
+b=[2;1];
+lb=[0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb);
+assert_equal ( exitflag , -2 );
+//
+// "Linear and Nonlinear Optimization"
+// Griva, Nash, Sofer
+// 2009
+// Chapter 5, "The Simplex Method", Example 5.6
+//
+// An unfeasible problem.
+// Minimize -x1
+//  -  x1 -   x2 <= -6
+//   2 x1 + 3 x2 <=  4
+// x >= 0
+c = [-1;0];
+A=[
+ -1  -1
+  2   3
+];
+b=[-6;4];
+lb=[0;0];
+[xopt,fopt,exitflag,iter,yopt]=karmarkar([],[],c,[],[],[],[],[],A,b,lb);
+assert_equal ( exitflag , -1 );
+//
+// Example from the help page
+// Check that the output points remain feasible
+function stop = myoutputfunction3 ( xopt , optimValues , state , A , b )
+    assert_equal ( and(A*xopt<=b) , %t );
+    stop = %f
+endfunction
+n=11;
+A = [2*linspace(0,1,n)',ones(n,1)];
+b = 1 + linspace(0,1,n)'.^2;
+c=[-1;-1];
+xopt=karmarkar([],[],c,[],[],[],[],list(myoutputfunction3,A,b),A,b);
+xstar = [0.5005127;0.7494803];
+assert_close ( xopt , xstar , 1.e-4 );
+
+