1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpapi06.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,806 @@
     1.4 +/* glpapi06.c (simplex method routines) */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpios.h"
    1.29 +#include "glpnpp.h"
    1.30 +#include "glpspx.h"
    1.31 +
    1.32 +/***********************************************************************
    1.33 +*  NAME
    1.34 +*
    1.35 +*  glp_simplex - solve LP problem with the simplex method
    1.36 +*
    1.37 +*  SYNOPSIS
    1.38 +*
    1.39 +*  int glp_simplex(glp_prob *P, const glp_smcp *parm);
    1.40 +*
    1.41 +*  DESCRIPTION
    1.42 +*
    1.43 +*  The routine glp_simplex is a driver to the LP solver based on the
    1.44 +*  simplex method. This routine retrieves problem data from the
    1.45 +*  specified problem object, calls the solver to solve the problem
    1.46 +*  instance, and stores results of computations back into the problem
    1.47 +*  object.
    1.48 +*
    1.49 +*  The simplex solver has a set of control parameters. Values of the
    1.50 +*  control parameters can be passed in a structure glp_smcp, which the
    1.51 +*  parameter parm points to.
    1.52 +*
    1.53 +*  The parameter parm can be specified as NULL, in which case the LP
    1.54 +*  solver uses default settings.
    1.55 +*
    1.56 +*  RETURNS
    1.57 +*
    1.58 +*  0  The LP problem instance has been successfully solved. This code
    1.59 +*     does not necessarily mean that the solver has found optimal
    1.60 +*     solution. It only means that the solution process was successful.
    1.61 +*
    1.62 +*  GLP_EBADB
    1.63 +*     Unable to start the search, because the initial basis specified
    1.64 +*     in the problem object is invalid--the number of basic (auxiliary
    1.65 +*     and structural) variables is not the same as the number of rows in
    1.66 +*     the problem object.
    1.67 +*
    1.68 +*  GLP_ESING
    1.69 +*     Unable to start the search, because the basis matrix correspodning
    1.70 +*     to the initial basis is singular within the working precision.
    1.71 +*
    1.72 +*  GLP_ECOND
    1.73 +*     Unable to start the search, because the basis matrix correspodning
    1.74 +*     to the initial basis is ill-conditioned, i.e. its condition number
    1.75 +*     is too large.
    1.76 +*
    1.77 +*  GLP_EBOUND
    1.78 +*     Unable to start the search, because some double-bounded variables
    1.79 +*     have incorrect bounds.
    1.80 +*
    1.81 +*  GLP_EFAIL
    1.82 +*     The search was prematurely terminated due to the solver failure.
    1.83 +*
    1.84 +*  GLP_EOBJLL
    1.85 +*     The search was prematurely terminated, because the objective
    1.86 +*     function being maximized has reached its lower limit and continues
    1.87 +*     decreasing (dual simplex only).
    1.88 +*
    1.89 +*  GLP_EOBJUL
    1.90 +*     The search was prematurely terminated, because the objective
    1.91 +*     function being minimized has reached its upper limit and continues
    1.92 +*     increasing (dual simplex only).
    1.93 +*
    1.94 +*  GLP_EITLIM
    1.95 +*     The search was prematurely terminated, because the simplex
    1.96 +*     iteration limit has been exceeded.
    1.97 +*
    1.98 +*  GLP_ETMLIM
    1.99 +*     The search was prematurely terminated, because the time limit has
   1.100 +*     been exceeded.
   1.101 +*
   1.102 +*  GLP_ENOPFS
   1.103 +*     The LP problem instance has no primal feasible solution (only if
   1.104 +*     the LP presolver is used).
   1.105 +*
   1.106 +*  GLP_ENODFS
   1.107 +*     The LP problem instance has no dual feasible solution (only if the
   1.108 +*     LP presolver is used). */
   1.109 +
   1.110 +static void trivial_lp(glp_prob *P, const glp_smcp *parm)
   1.111 +{     /* solve trivial LP which has empty constraint matrix */
   1.112 +      GLPROW *row;
   1.113 +      GLPCOL *col;
   1.114 +      int i, j;
   1.115 +      double p_infeas, d_infeas, zeta;
   1.116 +      P->valid = 0;
   1.117 +      P->pbs_stat = P->dbs_stat = GLP_FEAS;
   1.118 +      P->obj_val = P->c0;
   1.119 +      P->some = 0;
   1.120 +      p_infeas = d_infeas = 0.0;
   1.121 +      /* make all auxiliary variables basic */
   1.122 +      for (i = 1; i <= P->m; i++)
   1.123 +      {  row = P->row[i];
   1.124 +         row->stat = GLP_BS;
   1.125 +         row->prim = row->dual = 0.0;
   1.126 +         /* check primal feasibility */
   1.127 +         if (row->type == GLP_LO || row->type == GLP_DB ||
   1.128 +             row->type == GLP_FX)
   1.129 +         {  /* row has lower bound */
   1.130 +            if (row->lb > + parm->tol_bnd)
   1.131 +            {  P->pbs_stat = GLP_NOFEAS;
   1.132 +               if (P->some == 0 && parm->meth != GLP_PRIMAL)
   1.133 +                  P->some = i;
   1.134 +            }
   1.135 +            if (p_infeas < + row->lb)
   1.136 +               p_infeas = + row->lb;
   1.137 +         }
   1.138 +         if (row->type == GLP_UP || row->type == GLP_DB ||
   1.139 +             row->type == GLP_FX)
   1.140 +         {  /* row has upper bound */
   1.141 +            if (row->ub < - parm->tol_bnd)
   1.142 +            {  P->pbs_stat = GLP_NOFEAS;
   1.143 +               if (P->some == 0 && parm->meth != GLP_PRIMAL)
   1.144 +                  P->some = i;
   1.145 +            }
   1.146 +            if (p_infeas < - row->ub)
   1.147 +               p_infeas = - row->ub;
   1.148 +         }
   1.149 +      }
   1.150 +      /* determine scale factor for the objective row */
   1.151 +      zeta = 1.0;
   1.152 +      for (j = 1; j <= P->n; j++)
   1.153 +      {  col = P->col[j];
   1.154 +         if (zeta < fabs(col->coef)) zeta = fabs(col->coef);
   1.155 +      }
   1.156 +      zeta = (P->dir == GLP_MIN ? +1.0 : -1.0) / zeta;
   1.157 +      /* make all structural variables non-basic */
   1.158 +      for (j = 1; j <= P->n; j++)
   1.159 +      {  col = P->col[j];
   1.160 +         if (col->type == GLP_FR)
   1.161 +            col->stat = GLP_NF, col->prim = 0.0;
   1.162 +         else if (col->type == GLP_LO)
   1.163 +lo:         col->stat = GLP_NL, col->prim = col->lb;
   1.164 +         else if (col->type == GLP_UP)
   1.165 +up:         col->stat = GLP_NU, col->prim = col->ub;
   1.166 +         else if (col->type == GLP_DB)
   1.167 +         {  if (zeta * col->coef > 0.0)
   1.168 +               goto lo;
   1.169 +            else if (zeta * col->coef < 0.0)
   1.170 +               goto up;
   1.171 +            else if (fabs(col->lb) <= fabs(col->ub))
   1.172 +               goto lo;
   1.173 +            else
   1.174 +               goto up;
   1.175 +         }
   1.176 +         else if (col->type == GLP_FX)
   1.177 +            col->stat = GLP_NS, col->prim = col->lb;
   1.178 +         col->dual = col->coef;
   1.179 +         P->obj_val += col->coef * col->prim;
   1.180 +         /* check dual feasibility */
   1.181 +         if (col->type == GLP_FR || col->type == GLP_LO)
   1.182 +         {  /* column has no upper bound */
   1.183 +            if (zeta * col->dual < - parm->tol_dj)
   1.184 +            {  P->dbs_stat = GLP_NOFEAS;
   1.185 +               if (P->some == 0 && parm->meth == GLP_PRIMAL)
   1.186 +                  P->some = P->m + j;
   1.187 +            }
   1.188 +            if (d_infeas < - zeta * col->dual)
   1.189 +               d_infeas = - zeta * col->dual;
   1.190 +         }
   1.191 +         if (col->type == GLP_FR || col->type == GLP_UP)
   1.192 +         {  /* column has no lower bound */
   1.193 +            if (zeta * col->dual > + parm->tol_dj)
   1.194 +            {  P->dbs_stat = GLP_NOFEAS;
   1.195 +               if (P->some == 0 && parm->meth == GLP_PRIMAL)
   1.196 +                  P->some = P->m + j;
   1.197 +            }
   1.198 +            if (d_infeas < + zeta * col->dual)
   1.199 +               d_infeas = + zeta * col->dual;
   1.200 +         }
   1.201 +      }
   1.202 +      /* simulate the simplex solver output */
   1.203 +      if (parm->msg_lev >= GLP_MSG_ON && parm->out_dly == 0)
   1.204 +      {  xprintf("~%6d: obj = %17.9e  infeas = %10.3e\n", P->it_cnt,
   1.205 +            P->obj_val, parm->meth == GLP_PRIMAL ? p_infeas : d_infeas);
   1.206 +      }
   1.207 +      if (parm->msg_lev >= GLP_MSG_ALL && parm->out_dly == 0)
   1.208 +      {  if (P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS)
   1.209 +            xprintf("OPTIMAL SOLUTION FOUND\n");
   1.210 +         else if (P->pbs_stat == GLP_NOFEAS)
   1.211 +            xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
   1.212 +         else if (parm->meth == GLP_PRIMAL)
   1.213 +            xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
   1.214 +         else
   1.215 +            xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n");
   1.216 +      }
   1.217 +      return;
   1.218 +}
   1.219 +
   1.220 +static int solve_lp(glp_prob *P, const glp_smcp *parm)
   1.221 +{     /* solve LP directly without using the preprocessor */
   1.222 +      int ret;
   1.223 +      if (!glp_bf_exists(P))
   1.224 +      {  ret = glp_factorize(P);
   1.225 +         if (ret == 0)
   1.226 +            ;
   1.227 +         else if (ret == GLP_EBADB)
   1.228 +         {  if (parm->msg_lev >= GLP_MSG_ERR)
   1.229 +               xprintf("glp_simplex: initial basis is invalid\n");
   1.230 +         }
   1.231 +         else if (ret == GLP_ESING)
   1.232 +         {  if (parm->msg_lev >= GLP_MSG_ERR)
   1.233 +               xprintf("glp_simplex: initial basis is singular\n");
   1.234 +         }
   1.235 +         else if (ret == GLP_ECOND)
   1.236 +         {  if (parm->msg_lev >= GLP_MSG_ERR)
   1.237 +               xprintf(
   1.238 +                  "glp_simplex: initial basis is ill-conditioned\n");
   1.239 +         }
   1.240 +         else
   1.241 +            xassert(ret != ret);
   1.242 +         if (ret != 0) goto done;
   1.243 +      }
   1.244 +      if (parm->meth == GLP_PRIMAL)
   1.245 +         ret = spx_primal(P, parm);
   1.246 +      else if (parm->meth == GLP_DUALP)
   1.247 +      {  ret = spx_dual(P, parm);
   1.248 +         if (ret == GLP_EFAIL && P->valid)
   1.249 +            ret = spx_primal(P, parm);
   1.250 +      }
   1.251 +      else if (parm->meth == GLP_DUAL)
   1.252 +         ret = spx_dual(P, parm);
   1.253 +      else
   1.254 +         xassert(parm != parm);
   1.255 +done: return ret;
   1.256 +}
   1.257 +
   1.258 +static int preprocess_and_solve_lp(glp_prob *P, const glp_smcp *parm)
   1.259 +{     /* solve LP using the preprocessor */
   1.260 +      NPP *npp;
   1.261 +      glp_prob *lp = NULL;
   1.262 +      glp_bfcp bfcp;
   1.263 +      int ret;
   1.264 +      if (parm->msg_lev >= GLP_MSG_ALL)
   1.265 +         xprintf("Preprocessing...\n");
   1.266 +      /* create preprocessor workspace */
   1.267 +      npp = npp_create_wksp();
   1.268 +      /* load original problem into the preprocessor workspace */
   1.269 +      npp_load_prob(npp, P, GLP_OFF, GLP_SOL, GLP_OFF);
   1.270 +      /* process LP prior to applying primal/dual simplex method */
   1.271 +      ret = npp_simplex(npp, parm);
   1.272 +      if (ret == 0)
   1.273 +         ;
   1.274 +      else if (ret == GLP_ENOPFS)
   1.275 +      {  if (parm->msg_lev >= GLP_MSG_ALL)
   1.276 +            xprintf("PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION\n");
   1.277 +      }
   1.278 +      else if (ret == GLP_ENODFS)
   1.279 +      {  if (parm->msg_lev >= GLP_MSG_ALL)
   1.280 +            xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n");
   1.281 +      }
   1.282 +      else
   1.283 +         xassert(ret != ret);
   1.284 +      if (ret != 0) goto done;
   1.285 +      /* build transformed LP */
   1.286 +      lp = glp_create_prob();
   1.287 +      npp_build_prob(npp, lp);
   1.288 +      /* if the transformed LP is empty, it has empty solution, which
   1.289 +         is optimal */
   1.290 +      if (lp->m == 0 && lp->n == 0)
   1.291 +      {  lp->pbs_stat = lp->dbs_stat = GLP_FEAS;
   1.292 +         lp->obj_val = lp->c0;
   1.293 +         if (parm->msg_lev >= GLP_MSG_ON && parm->out_dly == 0)
   1.294 +         {  xprintf("~%6d: obj = %17.9e  infeas = %10.3e\n", P->it_cnt,
   1.295 +               lp->obj_val, 0.0);
   1.296 +         }
   1.297 +         if (parm->msg_lev >= GLP_MSG_ALL)
   1.298 +            xprintf("OPTIMAL SOLUTION FOUND BY LP PREPROCESSOR\n");
   1.299 +         goto post;
   1.300 +      }
   1.301 +      if (parm->msg_lev >= GLP_MSG_ALL)
   1.302 +      {  xprintf("%d row%s, %d column%s, %d non-zero%s\n",
   1.303 +            lp->m, lp->m == 1 ? "" : "s", lp->n, lp->n == 1 ? "" : "s",
   1.304 +            lp->nnz, lp->nnz == 1 ? "" : "s");
   1.305 +      }
   1.306 +      /* inherit basis factorization control parameters */
   1.307 +      glp_get_bfcp(P, &bfcp);
   1.308 +      glp_set_bfcp(lp, &bfcp);
   1.309 +      /* scale the transformed problem */
   1.310 +      {  ENV *env = get_env_ptr();
   1.311 +         int term_out = env->term_out;
   1.312 +         if (!term_out || parm->msg_lev < GLP_MSG_ALL)
   1.313 +            env->term_out = GLP_OFF;
   1.314 +         else
   1.315 +            env->term_out = GLP_ON;
   1.316 +         glp_scale_prob(lp, GLP_SF_AUTO);
   1.317 +         env->term_out = term_out;
   1.318 +      }
   1.319 +      /* build advanced initial basis */
   1.320 +      {  ENV *env = get_env_ptr();
   1.321 +         int term_out = env->term_out;
   1.322 +         if (!term_out || parm->msg_lev < GLP_MSG_ALL)
   1.323 +            env->term_out = GLP_OFF;
   1.324 +         else
   1.325 +            env->term_out = GLP_ON;
   1.326 +         glp_adv_basis(lp, 0);
   1.327 +         env->term_out = term_out;
   1.328 +      }
   1.329 +      /* solve the transformed LP */
   1.330 +      lp->it_cnt = P->it_cnt;
   1.331 +      ret = solve_lp(lp, parm);
   1.332 +      P->it_cnt = lp->it_cnt;
   1.333 +      /* only optimal solution can be postprocessed */
   1.334 +      if (!(ret == 0 && lp->pbs_stat == GLP_FEAS && lp->dbs_stat ==
   1.335 +            GLP_FEAS))
   1.336 +      {  if (parm->msg_lev >= GLP_MSG_ERR)
   1.337 +            xprintf("glp_simplex: unable to recover undefined or non-op"
   1.338 +               "timal solution\n");
   1.339 +         if (ret == 0)
   1.340 +         {  if (lp->pbs_stat == GLP_NOFEAS)
   1.341 +               ret = GLP_ENOPFS;
   1.342 +            else if (lp->dbs_stat == GLP_NOFEAS)
   1.343 +               ret = GLP_ENODFS;
   1.344 +            else
   1.345 +               xassert(lp != lp);
   1.346 +         }
   1.347 +         goto done;
   1.348 +      }
   1.349 +post: /* postprocess solution from the transformed LP */
   1.350 +      npp_postprocess(npp, lp);
   1.351 +      /* the transformed LP is no longer needed */
   1.352 +      glp_delete_prob(lp), lp = NULL;
   1.353 +      /* store solution to the original problem */
   1.354 +      npp_unload_sol(npp, P);
   1.355 +      /* the original LP has been successfully solved */
   1.356 +      ret = 0;
   1.357 +done: /* delete the transformed LP, if it exists */
   1.358 +      if (lp != NULL) glp_delete_prob(lp);
   1.359 +      /* delete preprocessor workspace */
   1.360 +      npp_delete_wksp(npp);
   1.361 +      return ret;
   1.362 +}
   1.363 +
   1.364 +int glp_simplex(glp_prob *P, const glp_smcp *parm)
   1.365 +{     /* solve LP problem with the simplex method */
   1.366 +      glp_smcp _parm;
   1.367 +      int i, j, ret;
   1.368 +      /* check problem object */
   1.369 +      if (P == NULL || P->magic != GLP_PROB_MAGIC)
   1.370 +         xerror("glp_simplex: P = %p; invalid problem object\n", P);
   1.371 +      if (P->tree != NULL && P->tree->reason != 0)
   1.372 +         xerror("glp_simplex: operation not allowed\n");
   1.373 +      /* check control parameters */
   1.374 +      if (parm == NULL)
   1.375 +         parm = &_parm, glp_init_smcp((glp_smcp *)parm);
   1.376 +      if (!(parm->msg_lev == GLP_MSG_OFF ||
   1.377 +            parm->msg_lev == GLP_MSG_ERR ||
   1.378 +            parm->msg_lev == GLP_MSG_ON  ||
   1.379 +            parm->msg_lev == GLP_MSG_ALL ||
   1.380 +            parm->msg_lev == GLP_MSG_DBG))
   1.381 +         xerror("glp_simplex: msg_lev = %d; invalid parameter\n",
   1.382 +            parm->msg_lev);
   1.383 +      if (!(parm->meth == GLP_PRIMAL ||
   1.384 +            parm->meth == GLP_DUALP  ||
   1.385 +            parm->meth == GLP_DUAL))
   1.386 +         xerror("glp_simplex: meth = %d; invalid parameter\n",
   1.387 +            parm->meth);
   1.388 +      if (!(parm->pricing == GLP_PT_STD ||
   1.389 +            parm->pricing == GLP_PT_PSE))
   1.390 +         xerror("glp_simplex: pricing = %d; invalid parameter\n",
   1.391 +            parm->pricing);
   1.392 +      if (!(parm->r_test == GLP_RT_STD ||
   1.393 +            parm->r_test == GLP_RT_HAR))
   1.394 +         xerror("glp_simplex: r_test = %d; invalid parameter\n",
   1.395 +            parm->r_test);
   1.396 +      if (!(0.0 < parm->tol_bnd && parm->tol_bnd < 1.0))
   1.397 +         xerror("glp_simplex: tol_bnd = %g; invalid parameter\n",
   1.398 +            parm->tol_bnd);
   1.399 +      if (!(0.0 < parm->tol_dj && parm->tol_dj < 1.0))
   1.400 +         xerror("glp_simplex: tol_dj = %g; invalid parameter\n",
   1.401 +            parm->tol_dj);
   1.402 +      if (!(0.0 < parm->tol_piv && parm->tol_piv < 1.0))
   1.403 +         xerror("glp_simplex: tol_piv = %g; invalid parameter\n",
   1.404 +            parm->tol_piv);
   1.405 +      if (parm->it_lim < 0)
   1.406 +         xerror("glp_simplex: it_lim = %d; invalid parameter\n",
   1.407 +            parm->it_lim);
   1.408 +      if (parm->tm_lim < 0)
   1.409 +         xerror("glp_simplex: tm_lim = %d; invalid parameter\n",
   1.410 +            parm->tm_lim);
   1.411 +      if (parm->out_frq < 1)
   1.412 +         xerror("glp_simplex: out_frq = %d; invalid parameter\n",
   1.413 +            parm->out_frq);
   1.414 +      if (parm->out_dly < 0)
   1.415 +         xerror("glp_simplex: out_dly = %d; invalid parameter\n",
   1.416 +            parm->out_dly);
   1.417 +      if (!(parm->presolve == GLP_ON || parm->presolve == GLP_OFF))
   1.418 +         xerror("glp_simplex: presolve = %d; invalid parameter\n",
   1.419 +            parm->presolve);
   1.420 +      /* basic solution is currently undefined */
   1.421 +      P->pbs_stat = P->dbs_stat = GLP_UNDEF;
   1.422 +      P->obj_val = 0.0;
   1.423 +      P->some = 0;
   1.424 +      /* check bounds of double-bounded variables */
   1.425 +      for (i = 1; i <= P->m; i++)
   1.426 +      {  GLPROW *row = P->row[i];
   1.427 +         if (row->type == GLP_DB && row->lb >= row->ub)
   1.428 +         {  if (parm->msg_lev >= GLP_MSG_ERR)
   1.429 +               xprintf("glp_simplex: row %d: lb = %g, ub = %g; incorrec"
   1.430 +                  "t bounds\n", i, row->lb, row->ub);
   1.431 +            ret = GLP_EBOUND;
   1.432 +            goto done;
   1.433 +         }
   1.434 +      }
   1.435 +      for (j = 1; j <= P->n; j++)
   1.436 +      {  GLPCOL *col = P->col[j];
   1.437 +         if (col->type == GLP_DB && col->lb >= col->ub)
   1.438 +         {  if (parm->msg_lev >= GLP_MSG_ERR)
   1.439 +               xprintf("glp_simplex: column %d: lb = %g, ub = %g; incor"
   1.440 +                  "rect bounds\n", j, col->lb, col->ub);
   1.441 +            ret = GLP_EBOUND;
   1.442 +            goto done;
   1.443 +         }
   1.444 +      }
   1.445 +      /* solve LP problem */
   1.446 +      if (parm->msg_lev >= GLP_MSG_ALL)
   1.447 +      {  xprintf("GLPK Simplex Optimizer, v%s\n", glp_version());
   1.448 +         xprintf("%d row%s, %d column%s, %d non-zero%s\n",
   1.449 +            P->m, P->m == 1 ? "" : "s", P->n, P->n == 1 ? "" : "s",
   1.450 +            P->nnz, P->nnz == 1 ? "" : "s");
   1.451 +      }
   1.452 +      if (P->nnz == 0)
   1.453 +         trivial_lp(P, parm), ret = 0;
   1.454 +      else if (!parm->presolve)
   1.455 +         ret = solve_lp(P, parm);
   1.456 +      else
   1.457 +         ret = preprocess_and_solve_lp(P, parm);
   1.458 +done: /* return to the application program */
   1.459 +      return ret;
   1.460 +}
   1.461 +
   1.462 +/***********************************************************************
   1.463 +*  NAME
   1.464 +*
   1.465 +*  glp_init_smcp - initialize simplex method control parameters
   1.466 +*
   1.467 +*  SYNOPSIS
   1.468 +*
   1.469 +*  void glp_init_smcp(glp_smcp *parm);
   1.470 +*
   1.471 +*  DESCRIPTION
   1.472 +*
   1.473 +*  The routine glp_init_smcp initializes control parameters, which are
   1.474 +*  used by the simplex solver, with default values.
   1.475 +*
   1.476 +*  Default values of the control parameters are stored in a glp_smcp
   1.477 +*  structure, which the parameter parm points to. */
   1.478 +
   1.479 +void glp_init_smcp(glp_smcp *parm)
   1.480 +{     parm->msg_lev = GLP_MSG_ALL;
   1.481 +      parm->meth = GLP_PRIMAL;
   1.482 +      parm->pricing = GLP_PT_PSE;
   1.483 +      parm->r_test = GLP_RT_HAR;
   1.484 +      parm->tol_bnd = 1e-7;
   1.485 +      parm->tol_dj = 1e-7;
   1.486 +      parm->tol_piv = 1e-10;
   1.487 +      parm->obj_ll = -DBL_MAX;
   1.488 +      parm->obj_ul = +DBL_MAX;
   1.489 +      parm->it_lim = INT_MAX;
   1.490 +      parm->tm_lim = INT_MAX;
   1.491 +      parm->out_frq = 500;
   1.492 +      parm->out_dly = 0;
   1.493 +      parm->presolve = GLP_OFF;
   1.494 +      return;
   1.495 +}
   1.496 +
   1.497 +/***********************************************************************
   1.498 +*  NAME
   1.499 +*
   1.500 +*  glp_get_status - retrieve generic status of basic solution
   1.501 +*
   1.502 +*  SYNOPSIS
   1.503 +*
   1.504 +*  int glp_get_status(glp_prob *lp);
   1.505 +*
   1.506 +*  RETURNS
   1.507 +*
   1.508 +*  The routine glp_get_status reports the generic status of the basic
   1.509 +*  solution for the specified problem object as follows:
   1.510 +*
   1.511 +*  GLP_OPT    - solution is optimal;
   1.512 +*  GLP_FEAS   - solution is feasible;
   1.513 +*  GLP_INFEAS - solution is infeasible;
   1.514 +*  GLP_NOFEAS - problem has no feasible solution;
   1.515 +*  GLP_UNBND  - problem has unbounded solution;
   1.516 +*  GLP_UNDEF  - solution is undefined. */
   1.517 +
   1.518 +int glp_get_status(glp_prob *lp)
   1.519 +{     int status;
   1.520 +      status = glp_get_prim_stat(lp);
   1.521 +      switch (status)
   1.522 +      {  case GLP_FEAS:
   1.523 +            switch (glp_get_dual_stat(lp))
   1.524 +            {  case GLP_FEAS:
   1.525 +                  status = GLP_OPT;
   1.526 +                  break;
   1.527 +               case GLP_NOFEAS:
   1.528 +                  status = GLP_UNBND;
   1.529 +                  break;
   1.530 +               case GLP_UNDEF:
   1.531 +               case GLP_INFEAS:
   1.532 +                  status = status;
   1.533 +                  break;
   1.534 +               default:
   1.535 +                  xassert(lp != lp);
   1.536 +            }
   1.537 +            break;
   1.538 +         case GLP_UNDEF:
   1.539 +         case GLP_INFEAS:
   1.540 +         case GLP_NOFEAS:
   1.541 +            status = status;
   1.542 +            break;
   1.543 +         default:
   1.544 +            xassert(lp != lp);
   1.545 +      }
   1.546 +      return status;
   1.547 +}
   1.548 +
   1.549 +/***********************************************************************
   1.550 +*  NAME
   1.551 +*
   1.552 +*  glp_get_prim_stat - retrieve status of primal basic solution
   1.553 +*
   1.554 +*  SYNOPSIS
   1.555 +*
   1.556 +*  int glp_get_prim_stat(glp_prob *lp);
   1.557 +*
   1.558 +*  RETURNS
   1.559 +*
   1.560 +*  The routine glp_get_prim_stat reports the status of the primal basic
   1.561 +*  solution for the specified problem object as follows:
   1.562 +*
   1.563 +*  GLP_UNDEF  - primal solution is undefined;
   1.564 +*  GLP_FEAS   - primal solution is feasible;
   1.565 +*  GLP_INFEAS - primal solution is infeasible;
   1.566 +*  GLP_NOFEAS - no primal feasible solution exists. */
   1.567 +
   1.568 +int glp_get_prim_stat(glp_prob *lp)
   1.569 +{     int pbs_stat = lp->pbs_stat;
   1.570 +      return pbs_stat;
   1.571 +}
   1.572 +
   1.573 +/***********************************************************************
   1.574 +*  NAME
   1.575 +*
   1.576 +*  glp_get_dual_stat - retrieve status of dual basic solution
   1.577 +*
   1.578 +*  SYNOPSIS
   1.579 +*
   1.580 +*  int glp_get_dual_stat(glp_prob *lp);
   1.581 +*
   1.582 +*  RETURNS
   1.583 +*
   1.584 +*  The routine glp_get_dual_stat reports the status of the dual basic
   1.585 +*  solution for the specified problem object as follows:
   1.586 +*
   1.587 +*  GLP_UNDEF  - dual solution is undefined;
   1.588 +*  GLP_FEAS   - dual solution is feasible;
   1.589 +*  GLP_INFEAS - dual solution is infeasible;
   1.590 +*  GLP_NOFEAS - no dual feasible solution exists. */
   1.591 +
   1.592 +int glp_get_dual_stat(glp_prob *lp)
   1.593 +{     int dbs_stat = lp->dbs_stat;
   1.594 +      return dbs_stat;
   1.595 +}
   1.596 +
   1.597 +/***********************************************************************
   1.598 +*  NAME
   1.599 +*
   1.600 +*  glp_get_obj_val - retrieve objective value (basic solution)
   1.601 +*
   1.602 +*  SYNOPSIS
   1.603 +*
   1.604 +*  double glp_get_obj_val(glp_prob *lp);
   1.605 +*
   1.606 +*  RETURNS
   1.607 +*
   1.608 +*  The routine glp_get_obj_val returns value of the objective function
   1.609 +*  for basic solution. */
   1.610 +
   1.611 +double glp_get_obj_val(glp_prob *lp)
   1.612 +{     /*struct LPXCPS *cps = lp->cps;*/
   1.613 +      double z;
   1.614 +      z = lp->obj_val;
   1.615 +      /*if (cps->round && fabs(z) < 1e-9) z = 0.0;*/
   1.616 +      return z;
   1.617 +}
   1.618 +
   1.619 +/***********************************************************************
   1.620 +*  NAME
   1.621 +*
   1.622 +*  glp_get_row_stat - retrieve row status
   1.623 +*
   1.624 +*  SYNOPSIS
   1.625 +*
   1.626 +*  int glp_get_row_stat(glp_prob *lp, int i);
   1.627 +*
   1.628 +*  RETURNS
   1.629 +*
   1.630 +*  The routine glp_get_row_stat returns current status assigned to the
   1.631 +*  auxiliary variable associated with i-th row as follows:
   1.632 +*
   1.633 +*  GLP_BS - basic variable;
   1.634 +*  GLP_NL - non-basic variable on its lower bound;
   1.635 +*  GLP_NU - non-basic variable on its upper bound;
   1.636 +*  GLP_NF - non-basic free (unbounded) variable;
   1.637 +*  GLP_NS - non-basic fixed variable. */
   1.638 +
   1.639 +int glp_get_row_stat(glp_prob *lp, int i)
   1.640 +{     if (!(1 <= i && i <= lp->m))
   1.641 +         xerror("glp_get_row_stat: i = %d; row number out of range\n",
   1.642 +            i);
   1.643 +      return lp->row[i]->stat;
   1.644 +}
   1.645 +
   1.646 +/***********************************************************************
   1.647 +*  NAME
   1.648 +*
   1.649 +*  glp_get_row_prim - retrieve row primal value (basic solution)
   1.650 +*
   1.651 +*  SYNOPSIS
   1.652 +*
   1.653 +*  double glp_get_row_prim(glp_prob *lp, int i);
   1.654 +*
   1.655 +*  RETURNS
   1.656 +*
   1.657 +*  The routine glp_get_row_prim returns primal value of the auxiliary
   1.658 +*  variable associated with i-th row. */
   1.659 +
   1.660 +double glp_get_row_prim(glp_prob *lp, int i)
   1.661 +{     /*struct LPXCPS *cps = lp->cps;*/
   1.662 +      double prim;
   1.663 +      if (!(1 <= i && i <= lp->m))
   1.664 +         xerror("glp_get_row_prim: i = %d; row number out of range\n",
   1.665 +            i);
   1.666 +      prim = lp->row[i]->prim;
   1.667 +      /*if (cps->round && fabs(prim) < 1e-9) prim = 0.0;*/
   1.668 +      return prim;
   1.669 +}
   1.670 +
   1.671 +/***********************************************************************
   1.672 +*  NAME
   1.673 +*
   1.674 +*  glp_get_row_dual - retrieve row dual value (basic solution)
   1.675 +*
   1.676 +*  SYNOPSIS
   1.677 +*
   1.678 +*  double glp_get_row_dual(glp_prob *lp, int i);
   1.679 +*
   1.680 +*  RETURNS
   1.681 +*
   1.682 +*  The routine glp_get_row_dual returns dual value (i.e. reduced cost)
   1.683 +*  of the auxiliary variable associated with i-th row. */
   1.684 +
   1.685 +double glp_get_row_dual(glp_prob *lp, int i)
   1.686 +{     /*struct LPXCPS *cps = lp->cps;*/
   1.687 +      double dual;
   1.688 +      if (!(1 <= i && i <= lp->m))
   1.689 +         xerror("glp_get_row_dual: i = %d; row number out of range\n",
   1.690 +            i);
   1.691 +      dual = lp->row[i]->dual;
   1.692 +      /*if (cps->round && fabs(dual) < 1e-9) dual = 0.0;*/
   1.693 +      return dual;
   1.694 +}
   1.695 +
   1.696 +/***********************************************************************
   1.697 +*  NAME
   1.698 +*
   1.699 +*  glp_get_col_stat - retrieve column status
   1.700 +*
   1.701 +*  SYNOPSIS
   1.702 +*
   1.703 +*  int glp_get_col_stat(glp_prob *lp, int j);
   1.704 +*
   1.705 +*  RETURNS
   1.706 +*
   1.707 +*  The routine glp_get_col_stat returns current status assigned to the
   1.708 +*  structural variable associated with j-th column as follows:
   1.709 +*
   1.710 +*  GLP_BS - basic variable;
   1.711 +*  GLP_NL - non-basic variable on its lower bound;
   1.712 +*  GLP_NU - non-basic variable on its upper bound;
   1.713 +*  GLP_NF - non-basic free (unbounded) variable;
   1.714 +*  GLP_NS - non-basic fixed variable. */
   1.715 +
   1.716 +int glp_get_col_stat(glp_prob *lp, int j)
   1.717 +{     if (!(1 <= j && j <= lp->n))
   1.718 +         xerror("glp_get_col_stat: j = %d; column number out of range\n"
   1.719 +            , j);
   1.720 +      return lp->col[j]->stat;
   1.721 +}
   1.722 +
   1.723 +/***********************************************************************
   1.724 +*  NAME
   1.725 +*
   1.726 +*  glp_get_col_prim - retrieve column primal value (basic solution)
   1.727 +*
   1.728 +*  SYNOPSIS
   1.729 +*
   1.730 +*  double glp_get_col_prim(glp_prob *lp, int j);
   1.731 +*
   1.732 +*  RETURNS
   1.733 +*
   1.734 +*  The routine glp_get_col_prim returns primal value of the structural
   1.735 +*  variable associated with j-th column. */
   1.736 +
   1.737 +double glp_get_col_prim(glp_prob *lp, int j)
   1.738 +{     /*struct LPXCPS *cps = lp->cps;*/
   1.739 +      double prim;
   1.740 +      if (!(1 <= j && j <= lp->n))
   1.741 +         xerror("glp_get_col_prim: j = %d; column number out of range\n"
   1.742 +            , j);
   1.743 +      prim = lp->col[j]->prim;
   1.744 +      /*if (cps->round && fabs(prim) < 1e-9) prim = 0.0;*/
   1.745 +      return prim;
   1.746 +}
   1.747 +
   1.748 +/***********************************************************************
   1.749 +*  NAME
   1.750 +*
   1.751 +*  glp_get_col_dual - retrieve column dual value (basic solution)
   1.752 +*
   1.753 +*  SYNOPSIS
   1.754 +*
   1.755 +*  double glp_get_col_dual(glp_prob *lp, int j);
   1.756 +*
   1.757 +*  RETURNS
   1.758 +*
   1.759 +*  The routine glp_get_col_dual returns dual value (i.e. reduced cost)
   1.760 +*  of the structural variable associated with j-th column. */
   1.761 +
   1.762 +double glp_get_col_dual(glp_prob *lp, int j)
   1.763 +{     /*struct LPXCPS *cps = lp->cps;*/
   1.764 +      double dual;
   1.765 +      if (!(1 <= j && j <= lp->n))
   1.766 +         xerror("glp_get_col_dual: j = %d; column number out of range\n"
   1.767 +            , j);
   1.768 +      dual = lp->col[j]->dual;
   1.769 +      /*if (cps->round && fabs(dual) < 1e-9) dual = 0.0;*/
   1.770 +      return dual;
   1.771 +}
   1.772 +
   1.773 +/***********************************************************************
   1.774 +*  NAME
   1.775 +*
   1.776 +*  glp_get_unbnd_ray - determine variable causing unboundedness
   1.777 +*
   1.778 +*  SYNOPSIS
   1.779 +*
   1.780 +*  int glp_get_unbnd_ray(glp_prob *lp);
   1.781 +*
   1.782 +*  RETURNS
   1.783 +*
   1.784 +*  The routine glp_get_unbnd_ray returns the number k of a variable,
   1.785 +*  which causes primal or dual unboundedness. If 1 <= k <= m, it is
   1.786 +*  k-th auxiliary variable, and if m+1 <= k <= m+n, it is (k-m)-th
   1.787 +*  structural variable, where m is the number of rows, n is the number
   1.788 +*  of columns in the problem object. If such variable is not defined,
   1.789 +*  the routine returns 0.
   1.790 +*
   1.791 +*  COMMENTS
   1.792 +*
   1.793 +*  If it is not exactly known which version of the simplex solver
   1.794 +*  detected unboundedness, i.e. whether the unboundedness is primal or
   1.795 +*  dual, it is sufficient to check the status of the variable reported
   1.796 +*  with the routine glp_get_row_stat or glp_get_col_stat. If the
   1.797 +*  variable is non-basic, the unboundedness is primal, otherwise, if
   1.798 +*  the variable is basic, the unboundedness is dual (the latter case
   1.799 +*  means that the problem has no primal feasible dolution). */
   1.800 +
   1.801 +int glp_get_unbnd_ray(glp_prob *lp)
   1.802 +{     int k;
   1.803 +      k = lp->some;
   1.804 +      xassert(k >= 0);
   1.805 +      if (k > lp->m + lp->n) k = 0;
   1.806 +      return k;
   1.807 +}
   1.808 +
   1.809 +/* eof */