alpar@1: /* glpapi06.c (simplex method routines) */
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  This code is part of GLPK (GNU Linear Programming Kit).
alpar@1: *
alpar@1: *  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@1: *  2009, 2010 Andrew Makhorin, Department for Applied Informatics,
alpar@1: *  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@1: *  E-mail: <mao@gnu.org>.
alpar@1: *
alpar@1: *  GLPK is free software: you can redistribute it and/or modify it
alpar@1: *  under the terms of the GNU General Public License as published by
alpar@1: *  the Free Software Foundation, either version 3 of the License, or
alpar@1: *  (at your option) any later version.
alpar@1: *
alpar@1: *  GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@1: *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@1: *  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@1: *  License for more details.
alpar@1: *
alpar@1: *  You should have received a copy of the GNU General Public License
alpar@1: *  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@1: ***********************************************************************/
alpar@1: 
alpar@1: #include "glpios.h"
alpar@1: #include "glpnpp.h"
alpar@1: #include "glpspx.h"
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_simplex - solve LP problem with the simplex method
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_simplex(glp_prob *P, const glp_smcp *parm);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine glp_simplex is a driver to the LP solver based on the
alpar@1: *  simplex method. This routine retrieves problem data from the
alpar@1: *  specified problem object, calls the solver to solve the problem
alpar@1: *  instance, and stores results of computations back into the problem
alpar@1: *  object.
alpar@1: *
alpar@1: *  The simplex solver has a set of control parameters. Values of the
alpar@1: *  control parameters can be passed in a structure glp_smcp, which the
alpar@1: *  parameter parm points to.
alpar@1: *
alpar@1: *  The parameter parm can be specified as NULL, in which case the LP
alpar@1: *  solver uses default settings.
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  0  The LP problem instance has been successfully solved. This code
alpar@1: *     does not necessarily mean that the solver has found optimal
alpar@1: *     solution. It only means that the solution process was successful.
alpar@1: *
alpar@1: *  GLP_EBADB
alpar@1: *     Unable to start the search, because the initial basis specified
alpar@1: *     in the problem object is invalid--the number of basic (auxiliary
alpar@1: *     and structural) variables is not the same as the number of rows in
alpar@1: *     the problem object.
alpar@1: *
alpar@1: *  GLP_ESING
alpar@1: *     Unable to start the search, because the basis matrix correspodning
alpar@1: *     to the initial basis is singular within the working precision.
alpar@1: *
alpar@1: *  GLP_ECOND
alpar@1: *     Unable to start the search, because the basis matrix correspodning
alpar@1: *     to the initial basis is ill-conditioned, i.e. its condition number
alpar@1: *     is too large.
alpar@1: *
alpar@1: *  GLP_EBOUND
alpar@1: *     Unable to start the search, because some double-bounded variables
alpar@1: *     have incorrect bounds.
alpar@1: *
alpar@1: *  GLP_EFAIL
alpar@1: *     The search was prematurely terminated due to the solver failure.
alpar@1: *
alpar@1: *  GLP_EOBJLL
alpar@1: *     The search was prematurely terminated, because the objective
alpar@1: *     function being maximized has reached its lower limit and continues
alpar@1: *     decreasing (dual simplex only).
alpar@1: *
alpar@1: *  GLP_EOBJUL
alpar@1: *     The search was prematurely terminated, because the objective
alpar@1: *     function being minimized has reached its upper limit and continues
alpar@1: *     increasing (dual simplex only).
alpar@1: *
alpar@1: *  GLP_EITLIM
alpar@1: *     The search was prematurely terminated, because the simplex
alpar@1: *     iteration limit has been exceeded.
alpar@1: *
alpar@1: *  GLP_ETMLIM
alpar@1: *     The search was prematurely terminated, because the time limit has
alpar@1: *     been exceeded.
alpar@1: *
alpar@1: *  GLP_ENOPFS
alpar@1: *     The LP problem instance has no primal feasible solution (only if
alpar@1: *     the LP presolver is used).
alpar@1: *
alpar@1: *  GLP_ENODFS
alpar@1: *     The LP problem instance has no dual feasible solution (only if the
alpar@1: *     LP presolver is used). */
alpar@1: 
alpar@1: static void trivial_lp(glp_prob *P, const glp_smcp *parm)
alpar@1: {     /* solve trivial LP which has empty constraint matrix */
alpar@1:       GLPROW *row;
alpar@1:       GLPCOL *col;
alpar@1:       int i, j;
alpar@1:       double p_infeas, d_infeas, zeta;
alpar@1:       P->valid = 0;
alpar@1:       P->pbs_stat = P->dbs_stat = GLP_FEAS;
alpar@1:       P->obj_val = P->c0;
alpar@1:       P->some = 0;
alpar@1:       p_infeas = d_infeas = 0.0;
alpar@1:       /* make all auxiliary variables basic */
alpar@1:       for (i = 1; i <= P->m; i++)
alpar@1:       {  row = P->row[i];
alpar@1:          row->stat = GLP_BS;
alpar@1:          row->prim = row->dual = 0.0;
alpar@1:          /* check primal feasibility */
alpar@1:          if (row->type == GLP_LO || row->type == GLP_DB ||
alpar@1:              row->type == GLP_FX)
alpar@1:          {  /* row has lower bound */
alpar@1:             if (row->lb > + parm->tol_bnd)
alpar@1:             {  P->pbs_stat = GLP_NOFEAS;
alpar@1:                if (P->some == 0 && parm->meth != GLP_PRIMAL)
alpar@1:                   P->some = i;
alpar@1:             }
alpar@1:             if (p_infeas < + row->lb)
alpar@1:                p_infeas = + row->lb;
alpar@1:          }
alpar@1:          if (row->type == GLP_UP || row->type == GLP_DB ||
alpar@1:              row->type == GLP_FX)
alpar@1:          {  /* row has upper bound */
alpar@1:             if (row->ub < - parm->tol_bnd)
alpar@1:             {  P->pbs_stat = GLP_NOFEAS;
alpar@1:                if (P->some == 0 && parm->meth != GLP_PRIMAL)
alpar@1:                   P->some = i;
alpar@1:             }
alpar@1:             if (p_infeas < - row->ub)
alpar@1:                p_infeas = - row->ub;
alpar@1:          }
alpar@1:       }
alpar@1:       /* determine scale factor for the objective row */
alpar@1:       zeta = 1.0;
alpar@1:       for (j = 1; j <= P->n; j++)
alpar@1:       {  col = P->col[j];
alpar@1:          if (zeta < fabs(col->coef)) zeta = fabs(col->coef);
alpar@1:       }
alpar@1:       zeta = (P->dir == GLP_MIN ? +1.0 : -1.0) / zeta;
alpar@1:       /* make all structural variables non-basic */
alpar@1:       for (j = 1; j <= P->n; j++)
alpar@1:       {  col = P->col[j];
alpar@1:          if (col->type == GLP_FR)
alpar@1:             col->stat = GLP_NF, col->prim = 0.0;
alpar@1:          else if (col->type == GLP_LO)
alpar@1: lo:         col->stat = GLP_NL, col->prim = col->lb;
alpar@1:          else if (col->type == GLP_UP)
alpar@1: up:         col->stat = GLP_NU, col->prim = col->ub;
alpar@1:          else if (col->type == GLP_DB)
alpar@1:          {  if (zeta * col->coef > 0.0)
alpar@1:                goto lo;
alpar@1:             else if (zeta * col->coef < 0.0)
alpar@1:                goto up;
alpar@1:             else if (fabs(col->lb) <= fabs(col->ub))
alpar@1:                goto lo;
alpar@1:             else
alpar@1:                goto up;
alpar@1:          }
alpar@1:          else if (col->type == GLP_FX)
alpar@1:             col->stat = GLP_NS, col->prim = col->lb;
alpar@1:          col->dual = col->coef;
alpar@1:          P->obj_val += col->coef * col->prim;
alpar@1:          /* check dual feasibility */
alpar@1:          if (col->type == GLP_FR || col->type == GLP_LO)
alpar@1:          {  /* column has no upper bound */
alpar@1:             if (zeta * col->dual < - parm->tol_dj)
alpar@1:             {  P->dbs_stat = GLP_NOFEAS;
alpar@1:                if (P->some == 0 && parm->meth == GLP_PRIMAL)
alpar@1:                   P->some = P->m + j;
alpar@1:             }
alpar@1:             if (d_infeas < - zeta * col->dual)
alpar@1:                d_infeas = - zeta * col->dual;
alpar@1:          }
alpar@1:          if (col->type == GLP_FR || col->type == GLP_UP)
alpar@1:          {  /* column has no lower bound */
alpar@1:             if (zeta * col->dual > + parm->tol_dj)
alpar@1:             {  P->dbs_stat = GLP_NOFEAS;
alpar@1:                if (P->some == 0 && parm->meth == GLP_PRIMAL)
alpar@1:                   P->some = P->m + j;
alpar@1:             }
alpar@1:             if (d_infeas < + zeta * col->dual)
alpar@1:                d_infeas = + zeta * col->dual;
alpar@1:          }
alpar@1:       }
alpar@1:       /* simulate the simplex solver output */
alpar@1:       if (parm->msg_lev >= GLP_MSG_ON && parm->out_dly == 0)
alpar@1:       {  xprintf("~%6d: obj = %17.9e  infeas = %10.3e\n", P->it_cnt,
alpar@1:             P->obj_val, parm->meth == GLP_PRIMAL ? p_infeas : d_infeas);
alpar@1:       }
alpar@1:       if (parm->msg_lev >= GLP_MSG_ALL && parm->out_dly == 0)
alpar@1:       {  if (P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS)
alpar@1:             xprintf("OPTIMAL SOLUTION FOUND\n");
alpar@1:          else if (P->pbs_stat == GLP_NOFEAS)
alpar@1:             xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
alpar@1:          else if (parm->meth == GLP_PRIMAL)
alpar@1:             xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
alpar@1:          else
alpar@1:             xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n");
alpar@1:       }
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: static int solve_lp(glp_prob *P, const glp_smcp *parm)
alpar@1: {     /* solve LP directly without using the preprocessor */
alpar@1:       int ret;
alpar@1:       if (!glp_bf_exists(P))
alpar@1:       {  ret = glp_factorize(P);
alpar@1:          if (ret == 0)
alpar@1:             ;
alpar@1:          else if (ret == GLP_EBADB)
alpar@1:          {  if (parm->msg_lev >= GLP_MSG_ERR)
alpar@1:                xprintf("glp_simplex: initial basis is invalid\n");
alpar@1:          }
alpar@1:          else if (ret == GLP_ESING)
alpar@1:          {  if (parm->msg_lev >= GLP_MSG_ERR)
alpar@1:                xprintf("glp_simplex: initial basis is singular\n");
alpar@1:          }
alpar@1:          else if (ret == GLP_ECOND)
alpar@1:          {  if (parm->msg_lev >= GLP_MSG_ERR)
alpar@1:                xprintf(
alpar@1:                   "glp_simplex: initial basis is ill-conditioned\n");
alpar@1:          }
alpar@1:          else
alpar@1:             xassert(ret != ret);
alpar@1:          if (ret != 0) goto done;
alpar@1:       }
alpar@1:       if (parm->meth == GLP_PRIMAL)
alpar@1:          ret = spx_primal(P, parm);
alpar@1:       else if (parm->meth == GLP_DUALP)
alpar@1:       {  ret = spx_dual(P, parm);
alpar@1:          if (ret == GLP_EFAIL && P->valid)
alpar@1:             ret = spx_primal(P, parm);
alpar@1:       }
alpar@1:       else if (parm->meth == GLP_DUAL)
alpar@1:          ret = spx_dual(P, parm);
alpar@1:       else
alpar@1:          xassert(parm != parm);
alpar@1: done: return ret;
alpar@1: }
alpar@1: 
alpar@1: static int preprocess_and_solve_lp(glp_prob *P, const glp_smcp *parm)
alpar@1: {     /* solve LP using the preprocessor */
alpar@1:       NPP *npp;
alpar@1:       glp_prob *lp = NULL;
alpar@1:       glp_bfcp bfcp;
alpar@1:       int ret;
alpar@1:       if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:          xprintf("Preprocessing...\n");
alpar@1:       /* create preprocessor workspace */
alpar@1:       npp = npp_create_wksp();
alpar@1:       /* load original problem into the preprocessor workspace */
alpar@1:       npp_load_prob(npp, P, GLP_OFF, GLP_SOL, GLP_OFF);
alpar@1:       /* process LP prior to applying primal/dual simplex method */
alpar@1:       ret = npp_simplex(npp, parm);
alpar@1:       if (ret == 0)
alpar@1:          ;
alpar@1:       else if (ret == GLP_ENOPFS)
alpar@1:       {  if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:             xprintf("PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION\n");
alpar@1:       }
alpar@1:       else if (ret == GLP_ENODFS)
alpar@1:       {  if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:             xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n");
alpar@1:       }
alpar@1:       else
alpar@1:          xassert(ret != ret);
alpar@1:       if (ret != 0) goto done;
alpar@1:       /* build transformed LP */
alpar@1:       lp = glp_create_prob();
alpar@1:       npp_build_prob(npp, lp);
alpar@1:       /* if the transformed LP is empty, it has empty solution, which
alpar@1:          is optimal */
alpar@1:       if (lp->m == 0 && lp->n == 0)
alpar@1:       {  lp->pbs_stat = lp->dbs_stat = GLP_FEAS;
alpar@1:          lp->obj_val = lp->c0;
alpar@1:          if (parm->msg_lev >= GLP_MSG_ON && parm->out_dly == 0)
alpar@1:          {  xprintf("~%6d: obj = %17.9e  infeas = %10.3e\n", P->it_cnt,
alpar@1:                lp->obj_val, 0.0);
alpar@1:          }
alpar@1:          if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:             xprintf("OPTIMAL SOLUTION FOUND BY LP PREPROCESSOR\n");
alpar@1:          goto post;
alpar@1:       }
alpar@1:       if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:       {  xprintf("%d row%s, %d column%s, %d non-zero%s\n",
alpar@1:             lp->m, lp->m == 1 ? "" : "s", lp->n, lp->n == 1 ? "" : "s",
alpar@1:             lp->nnz, lp->nnz == 1 ? "" : "s");
alpar@1:       }
alpar@1:       /* inherit basis factorization control parameters */
alpar@1:       glp_get_bfcp(P, &bfcp);
alpar@1:       glp_set_bfcp(lp, &bfcp);
alpar@1:       /* scale the transformed problem */
alpar@1:       {  ENV *env = get_env_ptr();
alpar@1:          int term_out = env->term_out;
alpar@1:          if (!term_out || parm->msg_lev < GLP_MSG_ALL)
alpar@1:             env->term_out = GLP_OFF;
alpar@1:          else
alpar@1:             env->term_out = GLP_ON;
alpar@1:          glp_scale_prob(lp, GLP_SF_AUTO);
alpar@1:          env->term_out = term_out;
alpar@1:       }
alpar@1:       /* build advanced initial basis */
alpar@1:       {  ENV *env = get_env_ptr();
alpar@1:          int term_out = env->term_out;
alpar@1:          if (!term_out || parm->msg_lev < GLP_MSG_ALL)
alpar@1:             env->term_out = GLP_OFF;
alpar@1:          else
alpar@1:             env->term_out = GLP_ON;
alpar@1:          glp_adv_basis(lp, 0);
alpar@1:          env->term_out = term_out;
alpar@1:       }
alpar@1:       /* solve the transformed LP */
alpar@1:       lp->it_cnt = P->it_cnt;
alpar@1:       ret = solve_lp(lp, parm);
alpar@1:       P->it_cnt = lp->it_cnt;
alpar@1:       /* only optimal solution can be postprocessed */
alpar@1:       if (!(ret == 0 && lp->pbs_stat == GLP_FEAS && lp->dbs_stat ==
alpar@1:             GLP_FEAS))
alpar@1:       {  if (parm->msg_lev >= GLP_MSG_ERR)
alpar@1:             xprintf("glp_simplex: unable to recover undefined or non-op"
alpar@1:                "timal solution\n");
alpar@1:          if (ret == 0)
alpar@1:          {  if (lp->pbs_stat == GLP_NOFEAS)
alpar@1:                ret = GLP_ENOPFS;
alpar@1:             else if (lp->dbs_stat == GLP_NOFEAS)
alpar@1:                ret = GLP_ENODFS;
alpar@1:             else
alpar@1:                xassert(lp != lp);
alpar@1:          }
alpar@1:          goto done;
alpar@1:       }
alpar@1: post: /* postprocess solution from the transformed LP */
alpar@1:       npp_postprocess(npp, lp);
alpar@1:       /* the transformed LP is no longer needed */
alpar@1:       glp_delete_prob(lp), lp = NULL;
alpar@1:       /* store solution to the original problem */
alpar@1:       npp_unload_sol(npp, P);
alpar@1:       /* the original LP has been successfully solved */
alpar@1:       ret = 0;
alpar@1: done: /* delete the transformed LP, if it exists */
alpar@1:       if (lp != NULL) glp_delete_prob(lp);
alpar@1:       /* delete preprocessor workspace */
alpar@1:       npp_delete_wksp(npp);
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: int glp_simplex(glp_prob *P, const glp_smcp *parm)
alpar@1: {     /* solve LP problem with the simplex method */
alpar@1:       glp_smcp _parm;
alpar@1:       int i, j, ret;
alpar@1:       /* check problem object */
alpar@1:       if (P == NULL || P->magic != GLP_PROB_MAGIC)
alpar@1:          xerror("glp_simplex: P = %p; invalid problem object\n", P);
alpar@1:       if (P->tree != NULL && P->tree->reason != 0)
alpar@1:          xerror("glp_simplex: operation not allowed\n");
alpar@1:       /* check control parameters */
alpar@1:       if (parm == NULL)
alpar@1:          parm = &_parm, glp_init_smcp((glp_smcp *)parm);
alpar@1:       if (!(parm->msg_lev == GLP_MSG_OFF ||
alpar@1:             parm->msg_lev == GLP_MSG_ERR ||
alpar@1:             parm->msg_lev == GLP_MSG_ON  ||
alpar@1:             parm->msg_lev == GLP_MSG_ALL ||
alpar@1:             parm->msg_lev == GLP_MSG_DBG))
alpar@1:          xerror("glp_simplex: msg_lev = %d; invalid parameter\n",
alpar@1:             parm->msg_lev);
alpar@1:       if (!(parm->meth == GLP_PRIMAL ||
alpar@1:             parm->meth == GLP_DUALP  ||
alpar@1:             parm->meth == GLP_DUAL))
alpar@1:          xerror("glp_simplex: meth = %d; invalid parameter\n",
alpar@1:             parm->meth);
alpar@1:       if (!(parm->pricing == GLP_PT_STD ||
alpar@1:             parm->pricing == GLP_PT_PSE))
alpar@1:          xerror("glp_simplex: pricing = %d; invalid parameter\n",
alpar@1:             parm->pricing);
alpar@1:       if (!(parm->r_test == GLP_RT_STD ||
alpar@1:             parm->r_test == GLP_RT_HAR))
alpar@1:          xerror("glp_simplex: r_test = %d; invalid parameter\n",
alpar@1:             parm->r_test);
alpar@1:       if (!(0.0 < parm->tol_bnd && parm->tol_bnd < 1.0))
alpar@1:          xerror("glp_simplex: tol_bnd = %g; invalid parameter\n",
alpar@1:             parm->tol_bnd);
alpar@1:       if (!(0.0 < parm->tol_dj && parm->tol_dj < 1.0))
alpar@1:          xerror("glp_simplex: tol_dj = %g; invalid parameter\n",
alpar@1:             parm->tol_dj);
alpar@1:       if (!(0.0 < parm->tol_piv && parm->tol_piv < 1.0))
alpar@1:          xerror("glp_simplex: tol_piv = %g; invalid parameter\n",
alpar@1:             parm->tol_piv);
alpar@1:       if (parm->it_lim < 0)
alpar@1:          xerror("glp_simplex: it_lim = %d; invalid parameter\n",
alpar@1:             parm->it_lim);
alpar@1:       if (parm->tm_lim < 0)
alpar@1:          xerror("glp_simplex: tm_lim = %d; invalid parameter\n",
alpar@1:             parm->tm_lim);
alpar@1:       if (parm->out_frq < 1)
alpar@1:          xerror("glp_simplex: out_frq = %d; invalid parameter\n",
alpar@1:             parm->out_frq);
alpar@1:       if (parm->out_dly < 0)
alpar@1:          xerror("glp_simplex: out_dly = %d; invalid parameter\n",
alpar@1:             parm->out_dly);
alpar@1:       if (!(parm->presolve == GLP_ON || parm->presolve == GLP_OFF))
alpar@1:          xerror("glp_simplex: presolve = %d; invalid parameter\n",
alpar@1:             parm->presolve);
alpar@1:       /* basic solution is currently undefined */
alpar@1:       P->pbs_stat = P->dbs_stat = GLP_UNDEF;
alpar@1:       P->obj_val = 0.0;
alpar@1:       P->some = 0;
alpar@1:       /* check bounds of double-bounded variables */
alpar@1:       for (i = 1; i <= P->m; i++)
alpar@1:       {  GLPROW *row = P->row[i];
alpar@1:          if (row->type == GLP_DB && row->lb >= row->ub)
alpar@1:          {  if (parm->msg_lev >= GLP_MSG_ERR)
alpar@1:                xprintf("glp_simplex: row %d: lb = %g, ub = %g; incorrec"
alpar@1:                   "t bounds\n", i, row->lb, row->ub);
alpar@1:             ret = GLP_EBOUND;
alpar@1:             goto done;
alpar@1:          }
alpar@1:       }
alpar@1:       for (j = 1; j <= P->n; j++)
alpar@1:       {  GLPCOL *col = P->col[j];
alpar@1:          if (col->type == GLP_DB && col->lb >= col->ub)
alpar@1:          {  if (parm->msg_lev >= GLP_MSG_ERR)
alpar@1:                xprintf("glp_simplex: column %d: lb = %g, ub = %g; incor"
alpar@1:                   "rect bounds\n", j, col->lb, col->ub);
alpar@1:             ret = GLP_EBOUND;
alpar@1:             goto done;
alpar@1:          }
alpar@1:       }
alpar@1:       /* solve LP problem */
alpar@1:       if (parm->msg_lev >= GLP_MSG_ALL)
alpar@1:       {  xprintf("GLPK Simplex Optimizer, v%s\n", glp_version());
alpar@1:          xprintf("%d row%s, %d column%s, %d non-zero%s\n",
alpar@1:             P->m, P->m == 1 ? "" : "s", P->n, P->n == 1 ? "" : "s",
alpar@1:             P->nnz, P->nnz == 1 ? "" : "s");
alpar@1:       }
alpar@1:       if (P->nnz == 0)
alpar@1:          trivial_lp(P, parm), ret = 0;
alpar@1:       else if (!parm->presolve)
alpar@1:          ret = solve_lp(P, parm);
alpar@1:       else
alpar@1:          ret = preprocess_and_solve_lp(P, parm);
alpar@1: done: /* return to the application program */
alpar@1:       return ret;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_init_smcp - initialize simplex method control parameters
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  void glp_init_smcp(glp_smcp *parm);
alpar@1: *
alpar@1: *  DESCRIPTION
alpar@1: *
alpar@1: *  The routine glp_init_smcp initializes control parameters, which are
alpar@1: *  used by the simplex solver, with default values.
alpar@1: *
alpar@1: *  Default values of the control parameters are stored in a glp_smcp
alpar@1: *  structure, which the parameter parm points to. */
alpar@1: 
alpar@1: void glp_init_smcp(glp_smcp *parm)
alpar@1: {     parm->msg_lev = GLP_MSG_ALL;
alpar@1:       parm->meth = GLP_PRIMAL;
alpar@1:       parm->pricing = GLP_PT_PSE;
alpar@1:       parm->r_test = GLP_RT_HAR;
alpar@1:       parm->tol_bnd = 1e-7;
alpar@1:       parm->tol_dj = 1e-7;
alpar@1:       parm->tol_piv = 1e-10;
alpar@1:       parm->obj_ll = -DBL_MAX;
alpar@1:       parm->obj_ul = +DBL_MAX;
alpar@1:       parm->it_lim = INT_MAX;
alpar@1:       parm->tm_lim = INT_MAX;
alpar@1:       parm->out_frq = 500;
alpar@1:       parm->out_dly = 0;
alpar@1:       parm->presolve = GLP_OFF;
alpar@1:       return;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_status - retrieve generic status of basic solution
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_get_status(glp_prob *lp);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_status reports the generic status of the basic
alpar@1: *  solution for the specified problem object as follows:
alpar@1: *
alpar@1: *  GLP_OPT    - solution is optimal;
alpar@1: *  GLP_FEAS   - solution is feasible;
alpar@1: *  GLP_INFEAS - solution is infeasible;
alpar@1: *  GLP_NOFEAS - problem has no feasible solution;
alpar@1: *  GLP_UNBND  - problem has unbounded solution;
alpar@1: *  GLP_UNDEF  - solution is undefined. */
alpar@1: 
alpar@1: int glp_get_status(glp_prob *lp)
alpar@1: {     int status;
alpar@1:       status = glp_get_prim_stat(lp);
alpar@1:       switch (status)
alpar@1:       {  case GLP_FEAS:
alpar@1:             switch (glp_get_dual_stat(lp))
alpar@1:             {  case GLP_FEAS:
alpar@1:                   status = GLP_OPT;
alpar@1:                   break;
alpar@1:                case GLP_NOFEAS:
alpar@1:                   status = GLP_UNBND;
alpar@1:                   break;
alpar@1:                case GLP_UNDEF:
alpar@1:                case GLP_INFEAS:
alpar@1:                   status = status;
alpar@1:                   break;
alpar@1:                default:
alpar@1:                   xassert(lp != lp);
alpar@1:             }
alpar@1:             break;
alpar@1:          case GLP_UNDEF:
alpar@1:          case GLP_INFEAS:
alpar@1:          case GLP_NOFEAS:
alpar@1:             status = status;
alpar@1:             break;
alpar@1:          default:
alpar@1:             xassert(lp != lp);
alpar@1:       }
alpar@1:       return status;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_prim_stat - retrieve status of primal basic solution
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_get_prim_stat(glp_prob *lp);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_prim_stat reports the status of the primal basic
alpar@1: *  solution for the specified problem object as follows:
alpar@1: *
alpar@1: *  GLP_UNDEF  - primal solution is undefined;
alpar@1: *  GLP_FEAS   - primal solution is feasible;
alpar@1: *  GLP_INFEAS - primal solution is infeasible;
alpar@1: *  GLP_NOFEAS - no primal feasible solution exists. */
alpar@1: 
alpar@1: int glp_get_prim_stat(glp_prob *lp)
alpar@1: {     int pbs_stat = lp->pbs_stat;
alpar@1:       return pbs_stat;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_dual_stat - retrieve status of dual basic solution
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_get_dual_stat(glp_prob *lp);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_dual_stat reports the status of the dual basic
alpar@1: *  solution for the specified problem object as follows:
alpar@1: *
alpar@1: *  GLP_UNDEF  - dual solution is undefined;
alpar@1: *  GLP_FEAS   - dual solution is feasible;
alpar@1: *  GLP_INFEAS - dual solution is infeasible;
alpar@1: *  GLP_NOFEAS - no dual feasible solution exists. */
alpar@1: 
alpar@1: int glp_get_dual_stat(glp_prob *lp)
alpar@1: {     int dbs_stat = lp->dbs_stat;
alpar@1:       return dbs_stat;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_obj_val - retrieve objective value (basic solution)
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  double glp_get_obj_val(glp_prob *lp);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_obj_val returns value of the objective function
alpar@1: *  for basic solution. */
alpar@1: 
alpar@1: double glp_get_obj_val(glp_prob *lp)
alpar@1: {     /*struct LPXCPS *cps = lp->cps;*/
alpar@1:       double z;
alpar@1:       z = lp->obj_val;
alpar@1:       /*if (cps->round && fabs(z) < 1e-9) z = 0.0;*/
alpar@1:       return z;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_row_stat - retrieve row status
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_get_row_stat(glp_prob *lp, int i);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_row_stat returns current status assigned to the
alpar@1: *  auxiliary variable associated with i-th row as follows:
alpar@1: *
alpar@1: *  GLP_BS - basic variable;
alpar@1: *  GLP_NL - non-basic variable on its lower bound;
alpar@1: *  GLP_NU - non-basic variable on its upper bound;
alpar@1: *  GLP_NF - non-basic free (unbounded) variable;
alpar@1: *  GLP_NS - non-basic fixed variable. */
alpar@1: 
alpar@1: int glp_get_row_stat(glp_prob *lp, int i)
alpar@1: {     if (!(1 <= i && i <= lp->m))
alpar@1:          xerror("glp_get_row_stat: i = %d; row number out of range\n",
alpar@1:             i);
alpar@1:       return lp->row[i]->stat;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_row_prim - retrieve row primal value (basic solution)
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  double glp_get_row_prim(glp_prob *lp, int i);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_row_prim returns primal value of the auxiliary
alpar@1: *  variable associated with i-th row. */
alpar@1: 
alpar@1: double glp_get_row_prim(glp_prob *lp, int i)
alpar@1: {     /*struct LPXCPS *cps = lp->cps;*/
alpar@1:       double prim;
alpar@1:       if (!(1 <= i && i <= lp->m))
alpar@1:          xerror("glp_get_row_prim: i = %d; row number out of range\n",
alpar@1:             i);
alpar@1:       prim = lp->row[i]->prim;
alpar@1:       /*if (cps->round && fabs(prim) < 1e-9) prim = 0.0;*/
alpar@1:       return prim;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_row_dual - retrieve row dual value (basic solution)
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  double glp_get_row_dual(glp_prob *lp, int i);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_row_dual returns dual value (i.e. reduced cost)
alpar@1: *  of the auxiliary variable associated with i-th row. */
alpar@1: 
alpar@1: double glp_get_row_dual(glp_prob *lp, int i)
alpar@1: {     /*struct LPXCPS *cps = lp->cps;*/
alpar@1:       double dual;
alpar@1:       if (!(1 <= i && i <= lp->m))
alpar@1:          xerror("glp_get_row_dual: i = %d; row number out of range\n",
alpar@1:             i);
alpar@1:       dual = lp->row[i]->dual;
alpar@1:       /*if (cps->round && fabs(dual) < 1e-9) dual = 0.0;*/
alpar@1:       return dual;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_col_stat - retrieve column status
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_get_col_stat(glp_prob *lp, int j);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_col_stat returns current status assigned to the
alpar@1: *  structural variable associated with j-th column as follows:
alpar@1: *
alpar@1: *  GLP_BS - basic variable;
alpar@1: *  GLP_NL - non-basic variable on its lower bound;
alpar@1: *  GLP_NU - non-basic variable on its upper bound;
alpar@1: *  GLP_NF - non-basic free (unbounded) variable;
alpar@1: *  GLP_NS - non-basic fixed variable. */
alpar@1: 
alpar@1: int glp_get_col_stat(glp_prob *lp, int j)
alpar@1: {     if (!(1 <= j && j <= lp->n))
alpar@1:          xerror("glp_get_col_stat: j = %d; column number out of range\n"
alpar@1:             , j);
alpar@1:       return lp->col[j]->stat;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_col_prim - retrieve column primal value (basic solution)
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  double glp_get_col_prim(glp_prob *lp, int j);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_col_prim returns primal value of the structural
alpar@1: *  variable associated with j-th column. */
alpar@1: 
alpar@1: double glp_get_col_prim(glp_prob *lp, int j)
alpar@1: {     /*struct LPXCPS *cps = lp->cps;*/
alpar@1:       double prim;
alpar@1:       if (!(1 <= j && j <= lp->n))
alpar@1:          xerror("glp_get_col_prim: j = %d; column number out of range\n"
alpar@1:             , j);
alpar@1:       prim = lp->col[j]->prim;
alpar@1:       /*if (cps->round && fabs(prim) < 1e-9) prim = 0.0;*/
alpar@1:       return prim;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_col_dual - retrieve column dual value (basic solution)
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  double glp_get_col_dual(glp_prob *lp, int j);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_col_dual returns dual value (i.e. reduced cost)
alpar@1: *  of the structural variable associated with j-th column. */
alpar@1: 
alpar@1: double glp_get_col_dual(glp_prob *lp, int j)
alpar@1: {     /*struct LPXCPS *cps = lp->cps;*/
alpar@1:       double dual;
alpar@1:       if (!(1 <= j && j <= lp->n))
alpar@1:          xerror("glp_get_col_dual: j = %d; column number out of range\n"
alpar@1:             , j);
alpar@1:       dual = lp->col[j]->dual;
alpar@1:       /*if (cps->round && fabs(dual) < 1e-9) dual = 0.0;*/
alpar@1:       return dual;
alpar@1: }
alpar@1: 
alpar@1: /***********************************************************************
alpar@1: *  NAME
alpar@1: *
alpar@1: *  glp_get_unbnd_ray - determine variable causing unboundedness
alpar@1: *
alpar@1: *  SYNOPSIS
alpar@1: *
alpar@1: *  int glp_get_unbnd_ray(glp_prob *lp);
alpar@1: *
alpar@1: *  RETURNS
alpar@1: *
alpar@1: *  The routine glp_get_unbnd_ray returns the number k of a variable,
alpar@1: *  which causes primal or dual unboundedness. If 1 <= k <= m, it is
alpar@1: *  k-th auxiliary variable, and if m+1 <= k <= m+n, it is (k-m)-th
alpar@1: *  structural variable, where m is the number of rows, n is the number
alpar@1: *  of columns in the problem object. If such variable is not defined,
alpar@1: *  the routine returns 0.
alpar@1: *
alpar@1: *  COMMENTS
alpar@1: *
alpar@1: *  If it is not exactly known which version of the simplex solver
alpar@1: *  detected unboundedness, i.e. whether the unboundedness is primal or
alpar@1: *  dual, it is sufficient to check the status of the variable reported
alpar@1: *  with the routine glp_get_row_stat or glp_get_col_stat. If the
alpar@1: *  variable is non-basic, the unboundedness is primal, otherwise, if
alpar@1: *  the variable is basic, the unboundedness is dual (the latter case
alpar@1: *  means that the problem has no primal feasible dolution). */
alpar@1: 
alpar@1: int glp_get_unbnd_ray(glp_prob *lp)
alpar@1: {     int k;
alpar@1:       k = lp->some;
alpar@1:       xassert(k >= 0);
alpar@1:       if (k > lp->m + lp->n) k = 0;
alpar@1:       return k;
alpar@1: }
alpar@1: 
alpar@1: /* eof */