glplpp01.c

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

      /* save some information about the original problem */
      lpp->orig_m = lpx_get_num_rows(orig);
      lpp->orig_n = lpx_get_num_cols(orig);
      lpp->orig_nnz = lpx_get_num_nz(orig);
      lpp->orig_dir = lpx_get_obj_dir(orig);
      /* allocate working arrays */
      c = xcalloc(1+lpp->orig_n, sizeof(double));
      ndx = xcalloc(1+lpp->orig_n, sizeof(int));
      val = xcalloc(1+lpp->orig_n, sizeof(double));
      /* auxiliary variables (i.e. rows) in the original problem may
         have non-zero objective coefficients; so, we substitute these
         auxiliary variables into the objective function in order that
         it depends only on structural variables (i.e. columns); the
         resultant vector of objective coefficients is accumulated in
         the working array c */
      for (j = 1; j <= lpp->orig_n; j++)
         c[j] = lpx_get_obj_coef(orig, j);
      for (i = 1; i <= lpp->orig_m; i++)
      {  /* obtain an objective coefficient at i-th row */
#if 0
         temp = lpx_get_row_coef(orig, i);
#else
         temp = 0.0;
#endif
         /* substitute i-th row into the objective function */
         if (temp != 0.0)
         {  len = lpx_get_mat_row(orig, i, ndx, val);
            for (t = 1; t <= len; t++) c[ndx[t]] += val[t] * temp;
         }
      }
      /* copy rows of the original problem into the workspace; each
         row created in the workspace is assigned a reference number,
         which is its ordinal number in the original problem */
      for (i = 1; i <= lpp->orig_m; i++)
      {  lpx_get_row_bnds(orig, i, &typx, &lb, &ub);
         if (typx == LPX_FR || typx == LPX_UP) lb = -DBL_MAX;
         if (typx == LPX_FR || typx == LPX_LO) ub = +DBL_MAX;
         if (typx == LPX_FX) ub = lb;
         lpp_add_row(lpp, lb, ub);
      }
      /* copy columns of the original problem into the workspace; each
         column created in the workspace is assigned a reference number,
         which its ordinal number in the original problem */
      for (j = 1; j <= lpp->orig_n; j++)
      {  lpx_get_col_bnds(orig, j, &typx, &lb, &ub);
         if (typx == LPX_FR || typx == LPX_UP) lb = -DBL_MAX;
         if (typx == LPX_FR || typx == LPX_LO) ub = +DBL_MAX;
         if (typx == LPX_FX) ub = lb;
         lpp_add_col(lpp, lb, ub, c[j]);
      }
      /* copy the constant term of the original objective function */
      lpp->c0 = lpx_get_obj_coef(orig, 0);
      /* if the original problem is maximization, change the sign of
         the objective function, because the transformed problem to be
         processed by the presolver must be minimization */
      if (lpp->orig_dir == LPX_MAX)
      {  for (col = lpp->col_ptr; col != NULL; col = col->next)
            col->c = - col->c;
         lpp->c0 = - lpp->c0;
      }
      /* build an auxiliary array to map column ordinal numbers to the
         corresponding pointers */
      xassert(sizeof(LPPCOL *) <= sizeof(double));
      map = (LPPCOL **)c;
      for (col = lpp->col_ptr; col != NULL; col = col->next)
         map[col->j] = col;
      /* copy the original constraint matrix into the workspace */
      for (row = lpp->row_ptr; row != NULL; row = row->next)
#if 1
      {  len = lpx_get_mat_row(orig, row->i, ndx, val);
         for (t = 1; t <= len; t++)
            lpp_add_aij(lpp, row, map[ndx[t]], val[t]);
      }
#else /* 27/XI-2003 (the problem persists) */
      {  double big, eps;
         len = lpx_get_mat_row(orig, row->i, ndx, val);
         big = 0.0;
         for (t = 1; t <= len; t++)
            if (big < fabs(val[t])) big = fabs(val[t]);
         eps = 1e-10 * big;
         for (t = 1; t <= len; t++)
         {  if (fabs(val[t]) < eps) continue;
            lpp_add_aij(lpp, row, map[ndx[t]], val[t]);
         }
      }
#endif
      /* free working arrays */
      xfree(c);
      xfree(ndx);
      xfree(val);
      return;
}

/*----------------------------------------------------------------------
-- lpp_append_tqe - append new transformation queue entry.
--
-- *Synopsis*
--
-- #include "glplpp.h"
-- void *lpp_append_tqe(LPP *lpp, int type, int size);
--
-- *Description*
--
-- The routine lpp_append_tqe appends a new transformation queue entry
-- to the list.
--
-- The parameter type is the entry type.
--
-- The parameter size is the size of a specific part, in bytes.
--
-- *Returns*
--
-- The routine returns a pointer to a specific part, which is allocated
-- and attached to the entry. */

void *lpp_append_tqe(LPP *lpp, int type, int size)
{     LPPTQE *tqe;
      tqe = dmp_get_atomv(lpp->tqe_pool, sizeof(LPPTQE));
      tqe->type = type;
      tqe->info = dmp_get_atomv(lpp->tqe_pool, size);
      tqe->next = lpp->tqe_list;
      lpp->tqe_list = tqe;
      return tqe->info;
}

/*----------------------------------------------------------------------
-- lpp_build_prob - build resultant problem.
--
-- *Synopsis*
--
-- #include "glplpp.h"
-- LPX *lpp_build_prob(LPP *lpp);
--
-- *Description*
--
-- The routine lpp_build_prob converts the resultant LP problem from an
-- internal format, in which the problem is stored in the workspace, to
-- the LPX format.
--
-- *Returns*
--
-- The routine returns a pointer to the LP problem object. */

LPX *lpp_build_prob(LPP *lpp)
{     LPX *prob;
      LPPROW *row;
      LPPCOL *col;
      int i, j, typx;
      /* count number of rows and columns in the resultant problem */
      lpp->m = lpp->n = 0;
      for (row = lpp->row_ptr; row != NULL; row = row->next) lpp->m++;
      for (col = lpp->col_ptr; col != NULL; col = col->next) lpp->n++;
      /* allocate two arrays to save reference numbers assigned to rows
         and columns of the resultant problem */
      lpp->row_ref = xcalloc(1+lpp->m, sizeof(int));
      lpp->col_ref = xcalloc(1+lpp->n, sizeof(int));
      /* create LP problem object */
      prob = lpx_create_prob();
      /* the resultant problem should have the same optimization sense
         as the original problem */
      lpx_set_obj_dir(prob, lpp->orig_dir);
      /* set the constant term of the objective function */
      lpx_set_obj_coef(prob, 0,
         lpp->orig_dir == LPX_MIN ? + lpp->c0 : - lpp->c0);
      /* create rows of the resultant problem */
#if 0 /* 03/VII-2008 */
      xassert(lpp->m > 0);
#endif
      if (lpp->m > 0)
         lpx_add_rows(prob, lpp->m);
      for (i = 1, row = lpp->row_ptr; i <= lpp->m; i++, row = row->next)
      {  xassert(row != NULL);
         lpp->row_ref[i] = row->i;
         row->i = i;
         if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
            typx = LPX_FR;
         else if (row->ub == +DBL_MAX)
            typx = LPX_LO;
         else if (row->lb == -DBL_MAX)
            typx = LPX_UP;
         else if (row->lb != row->ub)
            typx = LPX_DB;
         else
            typx = LPX_FX;
         lpx_set_row_bnds(prob, i, typx, row->lb, row->ub);
      }
      xassert(row == NULL);
      /* create columns of the resultant problem */
#if 0 /* 03/VII-2008 */
      xassert(lpp->n > 0);
#endif
      if (lpp->n > 0)
         lpx_add_cols(prob, lpp->n);
      for (j = 1, col = lpp->col_ptr; j <= lpp->n; j++, col = col->next)
      {  xassert(col != NULL);
         lpp->col_ref[j] = col->j;
         col->j = j;
         if (col->lb == -DBL_MAX && col->ub == +DBL_MAX)
            typx = LPX_FR;
         else if (col->ub == +DBL_MAX)
            typx = LPX_LO;
         else if (col->lb == -DBL_MAX)
            typx = LPX_UP;
         else if (col->lb != col->ub)
            typx = LPX_DB;
         else
            typx = LPX_FX;
         lpx_set_col_bnds(prob, j, typx, col->lb, col->ub);
         lpx_set_obj_coef(prob, j,
            lpp->orig_dir == LPX_MIN ? + col->c : - col->c);
      }
      xassert(col == NULL);
      /* create the constraint matrix of the resultant problem */
#if 0
      info.lpp = lpp;
      info.row = NULL;
      info.aij = NULL;
      lpx_load_mat(prob, &info, next_aij);
#else
      {  LPPAIJ *aij;
         int len, *ind;
         double *val;
         ind = xcalloc(1+lpp->n, sizeof(int));
         val = xcalloc(1+lpp->n, sizeof(double));
         for (row = lpp->row_ptr; row != NULL; row = row->next)
         {  len = 0;
            for (aij = row->ptr; aij != NULL; aij = aij->r_next)
               len++, ind[len] = aij->col->j, val[len] = aij->val;
            lpx_set_mat_row(prob, row->i, len, ind, val);
         }
         xfree(ind);
         xfree(val);
      }
#endif
      /* count number of non-zeros in the resultant problem */
      lpp->nnz = lpx_get_num_nz(prob);
      /* internal data structures that represnts the resultant problem
         are no longer needed, so free them */
      dmp_delete_pool(lpp->row_pool), lpp->row_pool = NULL;
      dmp_delete_pool(lpp->col_pool), lpp->col_pool = NULL;
      dmp_delete_pool(lpp->aij_pool), lpp->aij_pool = NULL;
      lpp->row_ptr = NULL, lpp->col_ptr = NULL;
      lpp->row_que = NULL, lpp->col_que = NULL;
      /* return a pointer to the built LP problem object */
      return prob;
}

/*----------------------------------------------------------------------
-- lpp_alloc_sol - allocate recovered solution segment.
--
-- *Synopsis*
--
-- #include "glplpp.h"
-- void lpp_alloc_sol(LPP *lpp);
--
-- *Description*
--
-- The routine lpp_alloc_sol allocates and initializes the recovered
-- solution segment. */

void lpp_alloc_sol(LPP *lpp)
{     int i, j;
      lpp->row_stat = xcalloc(1+lpp->nrows, sizeof(int));
      lpp->row_prim = xcalloc(1+lpp->nrows, sizeof(double));
      lpp->row_dual = xcalloc(1+lpp->nrows, sizeof(double));
      lpp->col_stat = xcalloc(1+lpp->ncols, sizeof(int));
      lpp->col_prim = xcalloc(1+lpp->ncols, sizeof(double));
      lpp->col_dual = xcalloc(1+lpp->ncols, sizeof(double));
      for (i = 1; i <= lpp->nrows; i++) lpp->row_stat[i] = 0;
      for (j = 1; j <= lpp->ncols; j++) lpp->col_stat[j] = 0;
      return;
}

/*----------------------------------------------------------------------
-- lpp_load_sol - load basic solution into LP presolver workspace.
--
-- *Synopsis*
--
-- #include "glplpp.h"
-- void lpp_load_sol(LPP *lpp, LPX *prob);
--
-- *Description*
--
-- The routine lpp_load_sol loads a basic solution of the resultant LP
-- problem into the LP presolver workspace. */

void lpp_load_sol(LPP *lpp, LPX *prob)
{     int i, j, ref, stat;
      double prim, dual;
      xassert(lpp->m == lpx_get_num_rows(prob));
      xassert(lpp->n == lpx_get_num_cols(prob));
      xassert(lpp->orig_dir == lpx_get_obj_dir(prob));
      xassert(lpx_get_status(prob) != LPX_UNDEF);
      for (i = 1; i <= lpp->m; i++)
      {  lpx_get_row_info(prob, i, &stat, &prim, &dual);
         ref = lpp->row_ref[i];
         xassert(1 <= ref && ref <= lpp->nrows);
         xassert(lpp->row_stat[ref] == 0);
         lpp->row_stat[ref] = stat;
         lpp->row_prim[ref] = prim;
         lpp->row_dual[ref] =
            (lpp->orig_dir == LPX_MIN ? + dual : - dual);
      }
      for (j = 1; j <= lpp->n; j++)
      {  lpx_get_col_info(prob, j, &stat, &prim, &dual);
         ref = lpp->col_ref[j];
         xassert(1 <= ref && ref <= lpp->ncols);
         xassert(lpp->col_stat[ref] == 0);
         lpp->col_stat[ref] = stat;
         lpp->col_prim[ref] = prim;
         lpp->col_dual[ref] =
            (lpp->orig_dir == LPX_MIN ? + dual : - dual);
      }
      xfree(lpp->row_ref), lpp->row_ref = NULL;
      xfree(lpp->col_ref), lpp->col_ref = NULL;
      return;
}

/*----------------------------------------------------------------------
-- lpp_unload_sol - unload basic solution from LP presolver workspace.
--
-- *Synopsis*
--
-- #include "glplpp.h"
-- void lpp_unload_sol(LPP *lpp, LPX *orig);
--
-- *Description*
--
-- The routine lpp_unload_sol unloads a recovered basic solution from
-- the LP presolver workspace into the original problem object, which
-- the parameter orig points to. */

void lpp_unload_sol(LPP *lpp, LPX *orig)
{     int i, j, k, m, n, typx, tagx, p_stat, d_stat;
      double sum;
      m = lpp->orig_m;
      n = lpp->orig_n;
      xassert(m == lpx_get_num_rows(orig));
      xassert(n == lpx_get_num_cols(orig));
      xassert(lpp->orig_dir == lpx_get_obj_dir(orig));
      /* check row and column statuses */
      xassert(m <= lpp->nrows);
      xassert(n <= lpp->ncols);
      for (k = 1; k <= m+n; k++)
      {  tagx = (k <= m ? lpp->row_stat[k] : lpp->col_stat[k-m]);
         if (tagx != LPX_BS)
         {  if (k <= m)
               lpx_get_row_bnds(orig, k, &typx, NULL, NULL);
            else
               lpx_get_col_bnds(orig, k-m, &typx, NULL, NULL);
            switch (typx)
            {  case LPX_FR:
                  xassert(tagx == LPX_NF);
                  break;
               case LPX_LO:
                  xassert(tagx == LPX_NL);
                  break;
               case LPX_UP:
                  xassert(tagx == LPX_NU);
                  break;
               case LPX_DB:
                  xassert(tagx == LPX_NL || tagx == LPX_NU);
                  break;
               case LPX_FX:
                  xassert(tagx == LPX_NS);
                  break;
               default:
                  xassert(orig != orig);
            }
         }
      }
      /* if the original problem is maximization, change signs of dual
         values */
      if (lpp->orig_dir == LPX_MAX)
      {  for (i = 1; i <= m; i++) lpp->row_dual[i] = -lpp->row_dual[i];
         for (j = 1; j <= n; j++) lpp->col_dual[j] = -lpp->col_dual[j];
      }
      /* store solution components into the original problem object (it
         is assumed that the recovered solution is optimal) */
      p_stat = d_stat = GLP_FEAS;
      for (i = 1; i <= m; i++)
         lpp->row_stat[i] = lpp->row_stat[i] - LPX_BS + GLP_BS;
      for (j = 1; j <= n; j++)
         lpp->col_stat[j] = lpp->col_stat[j] - LPX_BS + GLP_BS;
      sum = lpx_get_obj_coef(orig, 0);
      for (j = 1; j <= n; j++)
         sum += lpx_get_obj_coef(orig, j) * lpp->col_prim[j];
      lpx_put_solution(orig, 1, &p_stat, &d_stat, &sum,
         lpp->row_stat, lpp->row_prim, lpp->row_dual,
         lpp->col_stat, lpp->col_prim, lpp->col_dual);
      for (i = 1; i <= m; i++)
         lpp->row_stat[i] = lpp->row_stat[i] - GLP_BS + LPX_BS;
      for (j = 1; j <= n; j++)
         lpp->col_stat[j] = lpp->col_stat[j] - GLP_BS + LPX_BS;
      return;
}

/*----------------------------------------------------------------------
-- lpp_delete_wksp - delete LP presolver workspace.
--
-- *Synopsis*
--
-- #include "glplpp.h"
-- void lpp_delete_wksp(LPP *lpp);
--
-- *Description*
--
-- The routine lpp_delete_wksp deletes an LP presolver workspace, which
-- the parameter lpp points to, freeing all the memory allocated to this
-- object. */

void lpp_delete_wksp(LPP *lpp)
{     if (lpp->row_pool != NULL) dmp_delete_pool(lpp->row_pool);
      if (lpp->col_pool != NULL) dmp_delete_pool(lpp->col_pool);
      if (lpp->aij_pool != NULL) dmp_delete_pool(lpp->aij_pool);
      if (lpp->tqe_pool != NULL) dmp_delete_pool(lpp->tqe_pool);
      if (lpp->row_ref != NULL) xfree(lpp->row_ref);
      if (lpp->col_ref != NULL) xfree(lpp->col_ref);
      if (lpp->row_stat != NULL) xfree(lpp->row_stat);
      if (lpp->row_prim != NULL) xfree(lpp->row_prim);
      if (lpp->row_dual != NULL) xfree(lpp->row_dual);
      if (lpp->col_stat != NULL) xfree(lpp->col_stat);
      if (lpp->col_prim != NULL) xfree(lpp->col_prim);
      if (lpp->col_dual != NULL) xfree(lpp->col_dual);
      xfree(lpp);
      return;
}

/* eof */