glpapi08.c

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

               break;
            case LPX_UP:
               /* source: -inf < x <= ub */
               /* result: x = ub - x', x' >= 0 */
               col_pval[k] = ub - x[j];
               break;
            case LPX_DB:
               /* source: lb <= x <= ub */
               /* result: x = lb + x', x' + x'' = ub - lb */
               col_pval[k] = lb + x[j];
               break;
            case LPX_FX:
               /* source: x = lb */
               /* result: just substitute */
               col_pval[k] = lb;
               break;
            default:
               xassert(type != type);
         }
      }
      /* compute primal values of auxiliary variables */
      /* xR = A * xS */
      ind = xcalloc(1+orig_n, sizeof(int));
      val = xcalloc(1+orig_n, sizeof(double));
      for (k = 1; k <= orig_m; k++)
      {  rii = glp_get_rii(lp, k);
         temp = 0.0;
         len = lpx_get_mat_row(lp, k, ind, val);
         for (t = 1; t <= len; t++)
         {  sjj = glp_get_sjj(lp, ind[t]);
            temp += (rii * val[t] * sjj) * col_pval[ind[t]];
         }
         row_pval[k] = temp;
      }
      xfree(ind);
      xfree(val);
      /* compute dual values of auxiliary variables */
      for (k = 1; k <= orig_m; k++)
      {  type = lpx_get_row_type(lp, k);
         i = ref[k];
         switch (type)
         {  case LPX_FR:
               xassert(i == 0);
               row_dval[k] = 0.0;
               break;
            case LPX_LO:
            case LPX_UP:
            case LPX_DB:
            case LPX_FX:
               xassert(1 <= i && i <= m);
               row_dval[k] = (dir == LPX_MIN ? +1.0 : -1.0) * y[i];
               break;
            default:
               xassert(type != type);
         }
      }
      /* compute dual values of structural variables */
      /* dS = cS - A' * (dR - cR) */
      ind = xcalloc(1+orig_m, sizeof(int));
      val = xcalloc(1+orig_m, sizeof(double));
      for (k = 1; k <= orig_n; k++)
      {  sjj = glp_get_sjj(lp, k);
         temp = lpx_get_obj_coef(lp, k) / sjj;
         len = lpx_get_mat_col(lp, k, ind, val);
         for (t = 1; t <= len; t++)
         {  rii = glp_get_rii(lp, ind[t]);
            temp -= (rii * val[t] * sjj) * row_dval[ind[t]];
         }
         col_dval[k] = temp;
      }
      xfree(ind);
      xfree(val);
      /* unscale solution of original LP */
      for (i = 1; i <= orig_m; i++)
      {  rii = glp_get_rii(lp, i);
         row_pval[i] /= rii;
         row_dval[i] *= rii;
      }
      for (j = 1; j <= orig_n; j++)
      {  sjj = glp_get_sjj(lp, j);
         col_pval[j] *= sjj;
         col_dval[j] /= sjj;
      }
      return;
}

int glp_interior(glp_prob *_lp, const void *parm)
{     /* easy-to-use driver to the interior-point method */
      struct dsa _dsa, *dsa = &_dsa;
      int ret, status;
      double *row_pval, *row_dval, *col_pval, *col_dval;
      if (parm != NULL)
         xerror("glp_interior: parm = %p; invalid parameter\n", parm);
      /* interior-point solution is currently undefined */
      _lp->ipt_stat = GLP_UNDEF;
      _lp->ipt_obj = 0.0;
      /* initialize working area */
      dsa->lp = _lp;
      dsa->orig_m = lpx_get_num_rows(dsa->lp);
      dsa->orig_n = lpx_get_num_cols(dsa->lp);
      dsa->ref = NULL;
      dsa->m = 0;
      dsa->n = 0;
      dsa->size = 0;
      dsa->ne = 0;
      dsa->ia = NULL;
      dsa->ja = NULL;
      dsa->ar = NULL;
      dsa->b = NULL;
      dsa->c = NULL;
      dsa->x = NULL;
      dsa->y = NULL;
      dsa->z = NULL;
      /* check if the problem is empty */
      if (!(dsa->orig_m > 0 && dsa->orig_n > 0))
      {  xprintf("glp_interior: problem has no rows and/or columns\n");
         ret = GLP_EFAIL;
         goto done;
      }
      /* issue warning about dense columns */
      if (dsa->orig_m >= 200)
      {  int j, len, ndc = 0;
         for (j = 1; j <= dsa->orig_n; j++)
         {  len = lpx_get_mat_col(dsa->lp, j, NULL, NULL);
            if ((double)len > 0.30 * (double)dsa->orig_m) ndc++;
         }
         if (ndc == 1)
            xprintf("glp_interior: WARNING: PROBLEM HAS ONE DENSE COLUM"
               "N\n");
         else if (ndc > 0)
            xprintf("glp_interior: WARNING: PROBLEM HAS %d DENSE COLUMN"
               "S\n", ndc);
      }
      /* determine dimension of transformed LP */
      xprintf("glp_interior: original LP problem has %d rows and %d col"
         "umns\n", dsa->orig_m, dsa->orig_n);
      calc_mn(dsa);
      xprintf("glp_interior: transformed LP problem has %d rows and %d "
         "columns\n", dsa->m, dsa->n);
      /* transform original LP to standard formulation */
      dsa->ref = xcalloc(1+dsa->orig_m+dsa->orig_n, sizeof(int));
      dsa->size = lpx_get_num_nz(dsa->lp);
      if (dsa->size < dsa->m) dsa->size = dsa->m;
      if (dsa->size < dsa->n) dsa->size = dsa->n;
      dsa->ne = 0;
      dsa->ia = xcalloc(1+dsa->size, sizeof(int));
      dsa->ja = xcalloc(1+dsa->size, sizeof(int));
      dsa->ar = xcalloc(1+dsa->size, sizeof(double));
      dsa->b = xcalloc(1+dsa->m, sizeof(double));
      dsa->c = xcalloc(1+dsa->n, sizeof(double));
      transform(dsa);
      /* solve the transformed LP problem */
      {  int *A_ptr = xcalloc(1+dsa->m+1, sizeof(int));
         int *A_ind = xcalloc(1+dsa->ne, sizeof(int));
         double *A_val = xcalloc(1+dsa->ne, sizeof(double));
         int i, k, pos, len;
         /* determine row lengths */
         for (i = 1; i <= dsa->m; i++)
            A_ptr[i] = 0;
         for (k = 1; k <= dsa->ne; k++)
            A_ptr[dsa->ia[k]]++;
         /* set up row pointers */
         pos = 1;
         for (i = 1; i <= dsa->m; i++)
            len = A_ptr[i], pos += len, A_ptr[i] = pos;
         A_ptr[dsa->m+1] = pos;
         xassert(pos - 1 == dsa->ne);
         /* build the matrix */
         for (k = 1; k <= dsa->ne; k++)
         {  pos = --A_ptr[dsa->ia[k]];
            A_ind[pos] = dsa->ja[k];
            A_val[pos] = dsa->ar[k];
         }
         xfree(dsa->ia), dsa->ia = NULL;
         xfree(dsa->ja), dsa->ja = NULL;
         xfree(dsa->ar), dsa->ar = NULL;
         dsa->x = xcalloc(1+dsa->n, sizeof(double));
         dsa->y = xcalloc(1+dsa->m, sizeof(double));
         dsa->z = xcalloc(1+dsa->n, sizeof(double));
         ret = ipm_main(dsa->m, dsa->n, A_ptr, A_ind, A_val, dsa->b,
            dsa->c, dsa->x, dsa->y, dsa->z);
         xfree(A_ptr);
         xfree(A_ind);
         xfree(A_val);
         xfree(dsa->b), dsa->b = NULL;
         xfree(dsa->c), dsa->c = NULL;
      }
      /* analyze return code reported by the solver */
      switch (ret)
      {  case 0:
            /* optimal solution found */
            ret = 0;
            break;
         case 1:
            /* problem has no feasible (primal or dual) solution */
            ret = GLP_ENOFEAS;
            break;
         case 2:
            /* no convergence */
            ret = GLP_ENOCVG;
            break;
         case 3:
            /* iteration limit exceeded */
            ret = GLP_EITLIM;
            break;
         case 4:
            /* numerical instability on solving Newtonian system */
            ret = GLP_EINSTAB;
            break;
         default:
            xassert(ret != ret);
      }
      /* recover solution of original LP */
      row_pval = xcalloc(1+dsa->orig_m, sizeof(double));
      row_dval = xcalloc(1+dsa->orig_m, sizeof(double));
      col_pval = xcalloc(1+dsa->orig_n, sizeof(double));
      col_dval = xcalloc(1+dsa->orig_n, sizeof(double));
      restore(dsa, row_pval, row_dval, col_pval, col_dval);
      xfree(dsa->ref), dsa->ref = NULL;
      xfree(dsa->x), dsa->x = NULL;
      xfree(dsa->y), dsa->y = NULL;
      xfree(dsa->z), dsa->z = NULL;
      /* store solution components into the problem object */
      status = (ret == 0 ? LPX_T_OPT : LPX_T_UNDEF);
      lpx_put_ipt_soln(dsa->lp, status, row_pval, row_dval, col_pval,
         col_dval);
      xfree(row_pval);
      xfree(row_dval);
      xfree(col_pval);
      xfree(col_dval);
done: /* return to the calling program */
      return ret;
}

/***********************************************************************
*  NAME
*
*  glp_ipt_status - retrieve status of interior-point solution
*
*  SYNOPSIS
*
*  int glp_ipt_status(glp_prob *lp);
*
*  RETURNS
*
*  The routine glp_ipt_status reports the status of solution found by
*  the interior-point solver as follows:
*
*  GLP_UNDEF  - interior-point solution is undefined;
*  GLP_OPT    - interior-point solution is optimal. */

int glp_ipt_status(glp_prob *lp)
{     int ipt_stat = lp->ipt_stat;
      return ipt_stat;
}

/***********************************************************************
*  NAME
*
*  glp_ipt_obj_val - retrieve objective value (interior point)
*
*  SYNOPSIS
*
*  double glp_ipt_obj_val(glp_prob *lp);
*
*  RETURNS
*
*  The routine glp_ipt_obj_val returns value of the objective function
*  for interior-point solution. */

double glp_ipt_obj_val(glp_prob *lp)
{     struct LPXCPS *cps = lp->cps;
      double z;
      z = lp->ipt_obj;
      if (cps->round && fabs(z) < 1e-9) z = 0.0;
      return z;
}

/***********************************************************************
*  NAME
*
*  glp_ipt_row_prim - retrieve row primal value (interior point)
*
*  SYNOPSIS
*
*  double glp_ipt_row_prim(glp_prob *lp, int i);
*
*  RETURNS
*
*  The routine glp_ipt_row_prim returns primal value of the auxiliary
*  variable associated with i-th row. */

double glp_ipt_row_prim(glp_prob *lp, int i)
{     struct LPXCPS *cps = lp->cps;
      double pval;
      if (!(1 <= i && i <= lp->m))
         xerror("glp_ipt_row_prim: i = %d; row number out of range\n",
            i);
      pval = lp->row[i]->pval;
      if (cps->round && fabs(pval) < 1e-9) pval = 0.0;
      return pval;
}

/***********************************************************************
*  NAME
*
*  glp_ipt_row_dual - retrieve row dual value (interior point)
*
*  SYNOPSIS
*
*  double glp_ipt_row_dual(glp_prob *lp, int i);
*
*  RETURNS
*
*  The routine glp_ipt_row_dual returns dual value (i.e. reduced cost)
*  of the auxiliary variable associated with i-th row. */

double glp_ipt_row_dual(glp_prob *lp, int i)
{     struct LPXCPS *cps = lp->cps;
      double dval;
      if (!(1 <= i && i <= lp->m))
         xerror("glp_ipt_row_dual: i = %d; row number out of range\n",
            i);
      dval = lp->row[i]->dval;
      if (cps->round && fabs(dval) < 1e-9) dval = 0.0;
      return dval;
}

/***********************************************************************
*  NAME
*
*  glp_ipt_col_prim - retrieve column primal value (interior point)
*
*  SYNOPSIS
*
*  double glp_ipt_col_prim(glp_prob *lp, int j);
*
*  RETURNS
*
*  The routine glp_ipt_col_prim returns primal value of the structural
*  variable associated with j-th column. */

double glp_ipt_col_prim(glp_prob *lp, int j)
{     struct LPXCPS *cps = lp->cps;
      double pval;
      if (!(1 <= j && j <= lp->n))
         xerror("glp_ipt_col_prim: j = %d; column number out of range\n"
            , j);
      pval = lp->col[j]->pval;
      if (cps->round && fabs(pval) < 1e-9) pval = 0.0;
      return pval;
}

/***********************************************************************
*  NAME
*
*  glp_ipt_col_dual - retrieve column dual value (interior point)
*
*  SYNOPSIS
*
*  #include "glplpx.h"
*  double glp_ipt_col_dual(glp_prob *lp, int j);
*
*  RETURNS
*
*  The routine glp_ipt_col_dual returns dual value (i.e. reduced cost)
*  of the structural variable associated with j-th column. */

double glp_ipt_col_dual(glp_prob *lp, int j)
{     struct LPXCPS *cps = lp->cps;
      double dval;
      if (!(1 <= j && j <= lp->n))
         xerror("glp_ipt_col_dual: j = %d; column number out of range\n"
            , j);
      dval = lp->col[j]->dval;
      if (cps->round && fabs(dval) < 1e-9) dval = 0.0;
      return dval;
}

/* eof */