glpssx02.c

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

         if (ssx->it_lim > 0) ssx->it_lim--;
         ssx->it_cnt++;
      }
      /* display final progress of the search */
      show_progress(ssx, 1);
      /* restore components of the original problem, which were changed
         by the routine */
      for (k = 1; k <= m+n; k++)
      {  type[k] = orig_type[k];
         mpq_set(lb[k], orig_lb[k]);
         mpq_clear(orig_lb[k]);
         mpq_set(ub[k], orig_ub[k]);
         mpq_clear(orig_ub[k]);
      }
      ssx->dir = orig_dir;
      for (k = 0; k <= m+n; k++)
      {  mpq_set(coef[k], orig_coef[k]);
         mpq_clear(orig_coef[k]);
      }
      xfree(orig_type);
      xfree(orig_lb);
      xfree(orig_ub);
      xfree(orig_coef);
      /* return to the calling program */
      return ret;
}

/*----------------------------------------------------------------------
// ssx_phase_II - find optimal solution.
//
// This routine implements phase II of the primal simplex method.
//
// On exit the routine returns one of the following codes:
//
// 0 - optimal solution found;
// 1 - problem has unbounded solution;
// 2 - iterations limit exceeded;
// 3 - time limit exceeded.
----------------------------------------------------------------------*/

int ssx_phase_II(SSX *ssx)
{     int ret;
      /* display initial progress of the search */
      show_progress(ssx, 2);
      /* main loop starts here */
      for (;;)
      {  /* display current progress of the search */
#if 0
         if (utime() - ssx->tm_lag >= ssx->out_frq - 0.001)
#else
         if (xdifftime(xtime(), ssx->tm_lag) >= ssx->out_frq - 0.001)
#endif
            show_progress(ssx, 2);
         /* check if the iterations limit has been exhausted */
         if (ssx->it_lim == 0)
         {  ret = 2;
            break;
         }
         /* check if the time limit has been exhausted */
#if 0
         if (ssx->tm_lim >= 0.0 && ssx->tm_lim <= utime() - ssx->tm_beg)
#else
         if (ssx->tm_lim >= 0.0 &&
             ssx->tm_lim <= xdifftime(xtime(), ssx->tm_beg))
#endif
         {  ret = 3;
            break;
         }
         /* choose non-basic variable xN[q] */
         ssx_chuzc(ssx);
         /* if xN[q] cannot be chosen, the current basic solution is
            dual feasible and therefore optimal */
         if (ssx->q == 0)
         {  ret = 0;
            break;
         }
         /* compute q-th column of the simplex table */
         ssx_eval_col(ssx);
         /* choose basic variable xB[p] */
         ssx_chuzr(ssx);
         /* if xB[p] cannot be chosen, the problem has no dual feasible
            solution (i.e. unbounded) */
         if (ssx->p == 0)
         {  ret = 1;
            break;
         }
         /* update values of basic variables */
         ssx_update_bbar(ssx);
         if (ssx->p > 0)
         {  /* compute p-th row of the inverse inv(B) */
            ssx_eval_rho(ssx);
            /* compute p-th row of the simplex table */
            ssx_eval_row(ssx);
            xassert(mpq_cmp(ssx->aq[ssx->p], ssx->ap[ssx->q]) == 0);
#if 0
            /* update simplex multipliers */
            ssx_update_pi(ssx);
#endif
            /* update reduced costs of non-basic variables */
            ssx_update_cbar(ssx);
         }
         /* jump to the adjacent vertex of the polyhedron */
         ssx_change_basis(ssx);
         /* one simplex iteration has been performed */
         if (ssx->it_lim > 0) ssx->it_lim--;
         ssx->it_cnt++;
      }
      /* display final progress of the search */
      show_progress(ssx, 2);
      /* return to the calling program */
      return ret;
}

/*----------------------------------------------------------------------
// ssx_driver - base driver to exact simplex method.
//
// This routine is a base driver to a version of the primal simplex
// method using exact (bignum) arithmetic.
//
// On exit the routine returns one of the following codes:
//
// 0 - optimal solution found;
// 1 - problem has no feasible solution;
// 2 - problem has unbounded solution;
// 3 - iterations limit exceeded (phase I);
// 4 - iterations limit exceeded (phase II);
// 5 - time limit exceeded (phase I);
// 6 - time limit exceeded (phase II);
// 7 - initial basis matrix is exactly singular.
----------------------------------------------------------------------*/

int ssx_driver(SSX *ssx)
{     int m = ssx->m;
      int *type = ssx->type;
      mpq_t *lb = ssx->lb;
      mpq_t *ub = ssx->ub;
      int *Q_col = ssx->Q_col;
      mpq_t *bbar = ssx->bbar;
      int i, k, ret;
      ssx->tm_beg = xtime();
      /* factorize the initial basis matrix */
      if (ssx_factorize(ssx))
      {  xprintf("Initial basis matrix is singular\n");
         ret = 7;
         goto done;
      }
      /* compute values of basic variables */
      ssx_eval_bbar(ssx);
      /* check if the initial basic solution is primal feasible */
      for (i = 1; i <= m; i++)
      {  int t;
         k = Q_col[i]; /* x[k] = xB[i] */
         t = type[k];
         if (t == SSX_LO || t == SSX_DB || t == SSX_FX)
         {  /* x[k] has lower bound */
            if (mpq_cmp(bbar[i], lb[k]) < 0)
            {  /* which is violated */
               break;
            }
         }
         if (t == SSX_UP || t == SSX_DB || t == SSX_FX)
         {  /* x[k] has upper bound */
            if (mpq_cmp(bbar[i], ub[k]) > 0)
            {  /* which is violated */
               break;
            }
         }
      }
      if (i > m)
      {  /* no basic variable violates its bounds */
         ret = 0;
         goto skip;
      }
      /* phase I: find primal feasible solution */
      ret = ssx_phase_I(ssx);
      switch (ret)
      {  case 0:
            ret = 0;
            break;
         case 1:
            xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
            ret = 1;
            break;
         case 2:
            xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
            ret = 3;
            break;
         case 3:
            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
            ret = 5;
            break;
         default:
            xassert(ret != ret);
      }
      /* compute values of basic variables (actually only the objective
         value needs to be computed) */
      ssx_eval_bbar(ssx);
skip: /* compute simplex multipliers */
      ssx_eval_pi(ssx);
      /* compute reduced costs of non-basic variables */
      ssx_eval_cbar(ssx);
      /* if phase I failed, do not start phase II */
      if (ret != 0) goto done;
      /* phase II: find optimal solution */
      ret = ssx_phase_II(ssx);
      switch (ret)
      {  case 0:
            xprintf("OPTIMAL SOLUTION FOUND\n");
            ret = 0;
            break;
         case 1:
            xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
            ret = 2;
            break;
         case 2:
            xprintf("ITERATIONS LIMIT EXCEEDED; SEARCH TERMINATED\n");
            ret = 4;
            break;
         case 3:
            xprintf("TIME LIMIT EXCEEDED; SEARCH TERMINATED\n");
            ret = 6;
            break;
         default:
            xassert(ret != ret);
      }
done: /* decrease the time limit by the spent amount of time */
      if (ssx->tm_lim >= 0.0)
#if 0
      {  ssx->tm_lim -= utime() - ssx->tm_beg;
#else
      {  ssx->tm_lim -= xdifftime(xtime(), ssx->tm_beg);
#endif
         if (ssx->tm_lim < 0.0) ssx->tm_lim = 0.0;
      }
      return ret;
}

/* eof */