glpapi08.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

/* glpapi08.c (interior-point method routines) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glpapi.h"
#include "glpipm.h"

/***********************************************************************
*  NAME
*
*  glp_interior - solve LP problem with the interior-point method
*
*  SYNOPSIS
*
*  int glp_interior(glp_prob *lp, const void *parm);
*
*  DESCRIPTION
*
*  The routine glp_interior is a driver to the LP solver based on the
*  interior-point method.
*
*  The parameter parm is reserved for use in the future and must be
*  specified as NULL.
*
*  Currently this routine implements an easy variant of the primal-dual
*  interior-point method based on Mehrotra's technique.
*
*  This routine transforms the original LP problem to an equivalent LP
*  problem in the standard formulation (all constraints are equalities,
*  all variables are non-negative), calls the routine ipm_main to solve
*  the transformed problem, and then transforms an obtained solution to
*  the solution of the original problem.
*
*  RETURNS
*
*  0  The LP problem instance has been successfully solved. This code
*     does not necessarily mean that the solver has found optimal
*     solution. It only means that the solution process was successful.
*
*  GLP_EFAIL
*     The problem has no rows/columns.
*
*  GLP_ENOFEAS
*     The problem has no feasible (primal/dual) solution.
*
*  GLP_ENOCVG
*     Very slow convergence or divergence.
*
*  GLP_EITLIM
*     Iteration limit exceeded.
*
*  GLP_EINSTAB
*     Numerical instability on solving Newtonian system. */

struct dsa
{     /* working area used by glp_interior routine */
      LPX *lp;
      /* original LP instance to be solved */
      int orig_m;
      /* number of rows in original LP */
      int orig_n;
      /* number of columns in original LP */
      int *ref; /* int ref[1+orig_m+orig_n]; */
      /* this array contains references from components of original LP
         to components of transformed LP */
      /*--------------------------------------------------------------*/
      /* following components define transformed LP in standard form:
         minimize    c * x + c[0]
         subject to  A * x = b,   x >= 0 */
      int m;
      /* number of rows (constraints) in transformed LP */
      int n;
      /* number of columns (variables) in transformed LP */
      int size;
      /* size of arrays ia, ja, ar (enlarged automatically) */
      int ne;
      /* number of non-zeros in transformed constraint matrix */
      int *ia; /* int ia[1+size]; */
      int *ja; /* int ja[1+size]; */
      double *ar; /* double ar[1+size]; */
      /* transformed constraint matrix in storage-by-indices format
         which has m rows and n columns */
      double *b; /* double b[1+m]; */
      /* right-hand sides; b[0] is not used */
      double *c; /* double c[1+n]; */
      /* objective coefficients; c[0] is constant term */
      /*--------------------------------------------------------------*/
      /* following arrays define solution of transformed LP computed by
         interior-point solver */
      double *x; /* double x[1+n]; */
      /* values of primal variables */
      double *y; /* double y[1+m]; */
      /* values of dual variables for equality constraints */
      double *z; /* double z[1+n]; */
      /* values of dual variables for non-negativity constraints */
};

static void calc_mn(struct dsa *dsa)
{     /* calculate the number of constraints (m) and variables (n) for
         transformed LP */
      LPX *lp = dsa->lp;
      int orig_m = dsa->orig_m;
      int orig_n = dsa->orig_n;
      int i, j, m = 0, n = 0;
      for (i = 1; i <= orig_m; i++)
      {  switch (lpx_get_row_type(lp, i))
         {  case LPX_FR:                  break;
            case LPX_LO: m++; n++;        break;
            case LPX_UP: m++; n++;        break;
            case LPX_DB: m += 2; n += 2;  break;
            case LPX_FX: m++;             break;
            default: xassert(lp != lp);
         }
      }
      for (j = 1; j <= orig_n; j++)
      {  switch (lpx_get_col_type(lp, j))
         {  case LPX_FR: n += 2;          break;
            case LPX_LO: n++;             break;
            case LPX_UP: n++;             break;
            case LPX_DB: m++; n += 2;     break;
            case LPX_FX:                  break;
            default: xassert(lp != lp);
         }
      }
      dsa->m = m;
      dsa->n = n;
      return;
}

static void new_coef(struct dsa *dsa, int i, int j, double aij)
{     /* add new element to the transformed constraint matrix */
      xassert(1 <= i && i <= dsa->m);
      xassert(1 <= j && j <= dsa->n);
      xassert(aij != 0.0);
      if (dsa->ne == dsa->size)
      {  int *ia = dsa->ia;
         int *ja = dsa->ja;
         double *ar = dsa->ar;
         dsa->size += dsa->size;
         dsa->ia = xcalloc(1+dsa->size, sizeof(int));
         memcpy(&dsa->ia[1], &ia[1], dsa->ne * sizeof(int));
         xfree(ia);
         dsa->ja = xcalloc(1+dsa->size, sizeof(int));
         memcpy(&dsa->ja[1], &ja[1], dsa->ne * sizeof(int));
         xfree(ja);
         dsa->ar = xcalloc(1+dsa->size, sizeof(double));
         memcpy(&dsa->ar[1], &ar[1], dsa->ne * sizeof(double));
         xfree(ar);
      }
      xassert(dsa->ne < dsa->size);
      dsa->ne++;
      dsa->ia[dsa->ne] = i;
      dsa->ja[dsa->ne] = j;
      dsa->ar[dsa->ne] = aij;
      return;
}

static void transform(struct dsa *dsa)
{     /* transform original LP to standard formulation */
      LPX *lp = dsa->lp;
      int orig_m = dsa->orig_m;
      int orig_n = dsa->orig_n;
      int *ref = dsa->ref;
      int m = dsa->m;
      int n = dsa->n;
      double *b = dsa->b;
      double *c = dsa->c;
      int i, j, k, type, t, ii, len, *ind;
      double lb, ub, coef, rii, sjj, *val;
      /* initialize components of transformed LP */
      dsa->ne = 0;
      for (i = 1; i <= m; i++) b[i] = 0.0;
      c[0] = lpx_get_obj_coef(lp, 0);
      for (j = 1; j <= n; j++) c[j] = 0.0;
      /* i and j are, respectively, ordinal number of current row and
         ordinal number of current column in transformed LP */
      i = j = 0;
      /* transform rows (auxiliary variables) */
      for (k = 1; k <= orig_m; k++)
      {  type = lpx_get_row_type(lp, k);
         rii = glp_get_rii(lp, k);
         lb = lpx_get_row_lb(lp, k) * rii;
         ub = lpx_get_row_ub(lp, k) * rii;
         switch (type)
         {  case LPX_FR:
               /* source: -inf < (L.F.) < +inf */
               /* result: ignore free row */
               ref[k] = 0;
               break;
            case LPX_LO:
               /* source: lb <= (L.F.) < +inf */
               /* result: (L.F.) - x' = lb, x' >= 0 */
               i++; j++;
               ref[k] = i;
               new_coef(dsa, i, j, -1.0);
               b[i] = lb;
               break;
            case LPX_UP:
               /* source: -inf < (L.F.) <= ub */
               /* result: (L.F.) + x' = ub, x' >= 0 */
               i++; j++;
               ref[k] = i;
               new_coef(dsa, i, j, +1.0);
               b[i] = ub;
               break;
            case LPX_DB:
               /* source: lb <= (L.F.) <= ub */
               /* result: (L.F.) - x' = lb, x' + x'' = ub - lb */
               i++; j++;
               ref[k] = i;
               new_coef(dsa, i, j, -1.0);
               b[i] = lb;
               i++;
               new_coef(dsa, i, j, +1.0);
               j++;
               new_coef(dsa, i, j, +1.0);
               b[i] = ub - lb;
               break;
            case LPX_FX:
               /* source: (L.F.) = lb */
               /* result: (L.F.) = lb */
               i++;
               ref[k] = i;
               b[i] = lb;
               break;
            default:
               xassert(type != type);
         }
      }
      /* transform columns (structural variables) */
      ind = xcalloc(1+orig_m, sizeof(int));
      val = xcalloc(1+orig_m, sizeof(double));
      for (k = 1; k <= orig_n; k++)
      {  type = lpx_get_col_type(lp, k);
         sjj = glp_get_sjj(lp, k);
         lb = lpx_get_col_lb(lp, k) / sjj;
         ub = lpx_get_col_ub(lp, k) / sjj;
         coef = lpx_get_obj_coef(lp, k) * sjj;
         len = lpx_get_mat_col(lp, k, ind, val);
         for (t = 1; t <= len; t++)
            val[t] *= (glp_get_rii(lp, ind[t]) * sjj);
         switch (type)
         {  case LPX_FR:
               /* source: -inf < x < +inf */
               /* result: x = x' - x'', x' >= 0, x'' >= 0 */
               j++;
               ref[orig_m+k] = j;
               for (t = 1; t <= len; t++)
               {  ii = ref[ind[t]];
                  if (ii != 0) new_coef(dsa, ii, j, +val[t]);
               }
               c[j] = +coef;
               j++;
               for (t = 1; t <= len; t++)
               {  ii = ref[ind[t]];
                  if (ii != 0) new_coef(dsa, ii, j, -val[t]);
               }
               c[j] = -coef;
               break;
            case LPX_LO:
               /* source: lb <= x < +inf */
               /* result: x = lb + x', x' >= 0 */
               j++;
               ref[orig_m+k] = j;
               for (t = 1; t <= len; t++)
               {  ii = ref[ind[t]];
                  if (ii != 0)
                  {  new_coef(dsa, ii, j, val[t]);
                     b[ii] -= val[t] * lb;
                  }
               }
               c[j] = +coef;
               c[0] += c[j] * lb;
               break;
            case LPX_UP:
               /* source: -inf < x <= ub */
               /* result: x = ub - x', x' >= 0 */
               j++;
               ref[orig_m+k] = j;
               for (t = 1; t <= len; t++)
               {  ii = ref[ind[t]];
                  if (ii != 0)
                  {  new_coef(dsa, ii, j, -val[t]);
                     b[ii] -= val[t] * ub;
                  }
               }
               c[j] = -coef;
               c[0] -= c[j] * ub;
               break;
            case LPX_DB:
               /* source: lb <= x <= ub */
               /* result: x = lb + x', x' + x'' = ub - lb */
               j++;
               ref[orig_m+k] = j;
               for (t = 1; t <= len; t++)
               {  ii = ref[ind[t]];
                  if (ii != 0)
                  {  new_coef(dsa, ii, j, val[t]);
                     b[ii] -= val[t] * lb;
                  }
               }
               c[j] = +coef;
               c[0] += c[j] * lb;
               i++;
               new_coef(dsa, i, j, +1.0);
               j++;
               new_coef(dsa, i, j, +1.0);
               b[i] = ub - lb;
               break;
            case LPX_FX:
               /* source: x = lb */
               /* result: just substitute */
               ref[orig_m+k] = 0;
               for (t = 1; t <= len; t++)
               {  ii = ref[ind[t]];
                  if (ii != 0) b[ii] -= val[t] * lb;
               }
               c[0] += coef * lb;
               break;
            default:
               xassert(type != type);
         }
      }
      xfree(ind);
      xfree(val);
      /* end of transformation */
      xassert(i == m && j == n);
      /* change the objective sign in case of maximization */
      if (lpx_get_obj_dir(lp) == LPX_MAX)
         for (j = 0; j <= n; j++) c[j] = -c[j];
      return;
}

static void restore(struct dsa *dsa, double row_pval[],
      double row_dval[], double col_pval[], double col_dval[])
{     /* restore solution of original LP */
      LPX *lp = dsa->lp;
      int orig_m = dsa->orig_m;
      int orig_n = dsa->orig_n;
      int *ref = dsa->ref;
      int m = dsa->m;
      double *x = dsa->x;
      double *y = dsa->y;
      int dir = lpx_get_obj_dir(lp);
      int i, j, k, type, t, len, *ind;
      double lb, ub, rii, sjj, temp, *val;
      /* compute primal values of structural variables */
      for (k = 1; k <= orig_n; k++)
      {  j = ref[orig_m+k];
         type = lpx_get_col_type(lp, k);
         sjj = glp_get_sjj(lp, k);
         lb = lpx_get_col_lb(lp, k) / sjj;
         ub = lpx_get_col_ub(lp, k) / sjj;
         switch (type)
         {  case LPX_FR:
               /* source: -inf < x < +inf */
               /* result: x = x' - x'', x' >= 0, x'' >= 0 */
               col_pval[k] = x[j] - x[j+1];
               break;
            case LPX_LO:
               /* source: lb <= x < +inf */
               /* result: x = lb + x', x' >= 0 */
               col_pval[k] = lb + x[j];