glpapi19.c

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

      switch (ret)
      {  case 0:
            /* optimal circulation found */
            ret = 0;
            break;
         case 1:
            /* no feasible circulation exists */
            ret = GLP_ENOPFS;
            break;
         case 2:
            /* integer overflow occured */
            ret = GLP_ERANGE;
            goto done;
         case 3:
            /* optimality test failed (logic error) */
            ret = GLP_EFAIL;
            goto done;
         default:
            xassert(ret != ret);
      }
      /* store solution components */
      /* (objective function = the total cost) */
      if (sol != NULL)
      {  temp = 0.0;
         for (k = 1; k <= na; k++)
            temp += (double)cost[k] * (double)x[k];
         *sol = temp;
      }
      /* (arc flows) */
      if (a_x >= 0)
      {  k = 0;
         for (i = 1; i <= G->nv; i++)
         {  v = G->v[i];
            for (a = v->out; a != NULL; a = a->t_next)
            {  temp = (double)x[++k];
               memcpy((char *)a->data + a_x, &temp, sizeof(double));
            }
         }
      }
      /* (node potentials = Lagrange multipliers) */
      if (v_pi >= 0)
      {  for (i = 1; i <= G->nv; i++)
         {  v = G->v[i];
            temp = - (double)pi[i];
            memcpy((char *)v->data + v_pi, &temp, sizeof(double));
         }
      }
done: /* free working arrays */
      xfree(tail);
      xfree(head);
      xfree(low);
      xfree(cap);
      xfree(cost);
      xfree(x);
      xfree(pi);
      return ret;
}

/***********************************************************************
*  NAME
*
*  glp_maxflow_lp - convert maximum flow problem to LP
*
*  SYNOPSIS
*
*  void glp_maxflow_lp(glp_prob *lp, glp_graph *G, int names, int s,
*     int t, int a_cap);
*
*  DESCRIPTION
*
*  The routine glp_maxflow_lp builds an LP problem, which corresponds
*  to the maximum flow problem on the specified network G. */

void glp_maxflow_lp(glp_prob *lp, glp_graph *G, int names, int s,
      int t, int a_cap)
{     glp_vertex *v;
      glp_arc *a;
      int i, j, type, ind[1+2];
      double cap, val[1+2];
      if (!(names == GLP_ON || names == GLP_OFF))
         xerror("glp_maxflow_lp: names = %d; invalid parameter\n",
            names);
      if (!(1 <= s && s <= G->nv))
         xerror("glp_maxflow_lp: s = %d; source node number out of rang"
            "e\n", s);
      if (!(1 <= t && t <= G->nv))
         xerror("glp_maxflow_lp: t = %d: sink node number out of range "
            "\n", t);
      if (s == t)
         xerror("glp_maxflow_lp: s = t = %d; source and sink nodes must"
            " be distinct\n", s);
      if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
         xerror("glp_maxflow_lp: a_cap = %d; invalid offset\n", a_cap);
      glp_erase_prob(lp);
      if (names) glp_set_prob_name(lp, G->name);
      glp_set_obj_dir(lp, GLP_MAX);
      glp_add_rows(lp, G->nv);
      for (i = 1; i <= G->nv; i++)
      {  v = G->v[i];
         if (names) glp_set_row_name(lp, i, v->name);
         if (i == s)
            type = GLP_LO;
         else if (i == t)
            type = GLP_UP;
         else
            type = GLP_FX;
         glp_set_row_bnds(lp, i, type, 0.0, 0.0);
      }
      if (G->na > 0) glp_add_cols(lp, G->na);
      for (i = 1, j = 0; i <= G->nv; i++)
      {  v = G->v[i];
         for (a = v->out; a != NULL; a = a->t_next)
         {  j++;
            if (names)
            {  char name[50+1];
               sprintf(name, "x[%d,%d]", a->tail->i, a->head->i);
               xassert(strlen(name) < sizeof(name));
               glp_set_col_name(lp, j, name);
            }
            if (a->tail->i != a->head->i)
            {  ind[1] = a->tail->i, val[1] = +1.0;
               ind[2] = a->head->i, val[2] = -1.0;
               glp_set_mat_col(lp, j, 2, ind, val);
            }
            if (a_cap >= 0)
               memcpy(&cap, (char *)a->data + a_cap, sizeof(double));
            else
               cap = 1.0;
            if (cap == DBL_MAX)
               type = GLP_LO;
            else if (cap != 0.0)
               type = GLP_DB;
            else
               type = GLP_FX;
            glp_set_col_bnds(lp, j, type, 0.0, cap);
            if (a->tail->i == s)
               glp_set_obj_coef(lp, j, +1.0);
            else if (a->head->i == s)
               glp_set_obj_coef(lp, j, -1.0);
         }
      }
      xassert(j == G->na);
      return;
}

int glp_maxflow_ffalg(glp_graph *G, int s, int t, int a_cap,
      double *sol, int a_x, int v_cut)
{     /* find maximal flow with Ford-Fulkerson algorithm */
      glp_vertex *v;
      glp_arc *a;
      int nv, na, i, k, flag, *tail, *head, *cap, *x, ret;
      char *cut;
      double temp;
      if (!(1 <= s && s <= G->nv))
         xerror("glp_maxflow_ffalg: s = %d; source node number out of r"
            "ange\n", s);
      if (!(1 <= t && t <= G->nv))
         xerror("glp_maxflow_ffalg: t = %d: sink node number out of ran"
            "ge\n", t);
      if (s == t)
         xerror("glp_maxflow_ffalg: s = t = %d; source and sink nodes m"
            "ust be distinct\n", s);
      if (a_cap >= 0 && a_cap > G->a_size - (int)sizeof(double))
         xerror("glp_maxflow_ffalg: a_cap = %d; invalid offset\n",
            a_cap);
      if (v_cut >= 0 && v_cut > G->v_size - (int)sizeof(int))
         xerror("glp_maxflow_ffalg: v_cut = %d; invalid offset\n",
            v_cut);
      /* allocate working arrays */
      nv = G->nv;
      na = G->na;
      tail = xcalloc(1+na, sizeof(int));
      head = xcalloc(1+na, sizeof(int));
      cap = xcalloc(1+na, sizeof(int));
      x = xcalloc(1+na, sizeof(int));
      if (v_cut < 0)
         cut = NULL;
      else
         cut = xcalloc(1+nv, sizeof(char));
      /* copy the flow network */
      k = 0;
      for (i = 1; i <= G->nv; i++)
      {  v = G->v[i];
         for (a = v->out; a != NULL; a = a->t_next)
         {  k++;
            tail[k] = a->tail->i;
            head[k] = a->head->i;
            if (tail[k] == head[k])
            {  ret = GLP_EDATA;
               goto done;
            }
            if (a_cap >= 0)
               memcpy(&temp, (char *)a->data + a_cap, sizeof(double));
            else
               temp = 1.0;
            if (!(0.0 <= temp && temp <= (double)INT_MAX &&
                  temp == floor(temp)))
            {  ret = GLP_EDATA;
               goto done;
            }
            cap[k] = (int)temp;
         }
      }
      xassert(k == na);
      /* find maximal flow in the flow network */
      ffalg(nv, na, tail, head, s, t, cap, x, cut);
      ret = 0;
      /* store solution components */
      /* (objective function = total flow through the network) */
      if (sol != NULL)
      {  temp = 0.0;
         for (k = 1; k <= na; k++)
         {  if (tail[k] == s)
               temp += (double)x[k];
            else if (head[k] == s)
               temp -= (double)x[k];
         }
         *sol = temp;
      }
      /* (arc flows) */
      if (a_x >= 0)
      {  k = 0;
         for (i = 1; i <= G->nv; i++)
         {  v = G->v[i];
            for (a = v->out; a != NULL; a = a->t_next)
            {  temp = (double)x[++k];
               memcpy((char *)a->data + a_x, &temp, sizeof(double));
            }
         }
      }
      /* (node flags) */
      if (v_cut >= 0)
      {  for (i = 1; i <= G->nv; i++)
         {  v = G->v[i];
            flag = cut[i];
            memcpy((char *)v->data + v_cut, &flag, sizeof(int));
         }
      }
done: /* free working arrays */
      xfree(tail);
      xfree(head);
      xfree(cap);
      xfree(x);
      if (cut != NULL) xfree(cut);
      return ret;
}

/* eof */