lemon-project-template-glpk

annotate deps/glpk/src/glpapi06.c @ 9:33de93886c88

find changesets by author, revision, files, or words in the commit message
Import GLPK 4.47
author Alpar Juttner <alpar@cs.elte.hu>
date Sun, 06 Nov 2011 20:59:10 +0100
parents
children
rev   line source
alpar@9 1 /* glpapi06.c (simplex method routines) */
alpar@9 2
alpar@9 3 /***********************************************************************
alpar@9 4 * This code is part of GLPK (GNU Linear Programming Kit).
alpar@9 5 *
alpar@9 6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
alpar@9 7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
alpar@9 8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
alpar@9 9 * E-mail: <mao@gnu.org>.
alpar@9 10 *
alpar@9 11 * GLPK is free software: you can redistribute it and/or modify it
alpar@9 12 * under the terms of the GNU General Public License as published by
alpar@9 13 * the Free Software Foundation, either version 3 of the License, or
alpar@9 14 * (at your option) any later version.
alpar@9 15 *
alpar@9 16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
alpar@9 17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
alpar@9 18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
alpar@9 19 * License for more details.
alpar@9 20 *
alpar@9 21 * You should have received a copy of the GNU General Public License
alpar@9 22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
alpar@9 23 ***********************************************************************/
alpar@9 24
alpar@9 25 #include "glpios.h"
alpar@9 26 #include "glpnpp.h"
alpar@9 27 #include "glpspx.h"
alpar@9 28
alpar@9 29 /***********************************************************************
alpar@9 30 * NAME
alpar@9 31 *
alpar@9 32 * glp_simplex - solve LP problem with the simplex method
alpar@9 33 *
alpar@9 34 * SYNOPSIS
alpar@9 35 *
alpar@9 36 * int glp_simplex(glp_prob *P, const glp_smcp *parm);
alpar@9 37 *
alpar@9 38 * DESCRIPTION
alpar@9 39 *
alpar@9 40 * The routine glp_simplex is a driver to the LP solver based on the
alpar@9 41 * simplex method. This routine retrieves problem data from the
alpar@9 42 * specified problem object, calls the solver to solve the problem
alpar@9 43 * instance, and stores results of computations back into the problem
alpar@9 44 * object.
alpar@9 45 *
alpar@9 46 * The simplex solver has a set of control parameters. Values of the
alpar@9 47 * control parameters can be passed in a structure glp_smcp, which the
alpar@9 48 * parameter parm points to.
alpar@9 49 *
alpar@9 50 * The parameter parm can be specified as NULL, in which case the LP
alpar@9 51 * solver uses default settings.
alpar@9 52 *
alpar@9 53 * RETURNS
alpar@9 54 *
alpar@9 55 * 0 The LP problem instance has been successfully solved. This code
alpar@9 56 * does not necessarily mean that the solver has found optimal
alpar@9 57 * solution. It only means that the solution process was successful.
alpar@9 58 *
alpar@9 59 * GLP_EBADB
alpar@9 60 * Unable to start the search, because the initial basis specified
alpar@9 61 * in the problem object is invalid--the number of basic (auxiliary
alpar@9 62 * and structural) variables is not the same as the number of rows in
alpar@9 63 * the problem object.
alpar@9 64 *
alpar@9 65 * GLP_ESING
alpar@9 66 * Unable to start the search, because the basis matrix correspodning
alpar@9 67 * to the initial basis is singular within the working precision.
alpar@9 68 *
alpar@9 69 * GLP_ECOND
alpar@9 70 * Unable to start the search, because the basis matrix correspodning
alpar@9 71 * to the initial basis is ill-conditioned, i.e. its condition number
alpar@9 72 * is too large.
alpar@9 73 *
alpar@9 74 * GLP_EBOUND
alpar@9 75 * Unable to start the search, because some double-bounded variables
alpar@9 76 * have incorrect bounds.
alpar@9 77 *
alpar@9 78 * GLP_EFAIL
alpar@9 79 * The search was prematurely terminated due to the solver failure.
alpar@9 80 *
alpar@9 81 * GLP_EOBJLL
alpar@9 82 * The search was prematurely terminated, because the objective
alpar@9 83 * function being maximized has reached its lower limit and continues
alpar@9 84 * decreasing (dual simplex only).
alpar@9 85 *
alpar@9 86 * GLP_EOBJUL
alpar@9 87 * The search was prematurely terminated, because the objective
alpar@9 88 * function being minimized has reached its upper limit and continues
alpar@9 89 * increasing (dual simplex only).
alpar@9 90 *
alpar@9 91 * GLP_EITLIM
alpar@9 92 * The search was prematurely terminated, because the simplex
alpar@9 93 * iteration limit has been exceeded.
alpar@9 94 *
alpar@9 95 * GLP_ETMLIM
alpar@9 96 * The search was prematurely terminated, because the time limit has
alpar@9 97 * been exceeded.
alpar@9 98 *
alpar@9 99 * GLP_ENOPFS
alpar@9 100 * The LP problem instance has no primal feasible solution (only if
alpar@9 101 * the LP presolver is used).
alpar@9 102 *
alpar@9 103 * GLP_ENODFS
alpar@9 104 * The LP problem instance has no dual feasible solution (only if the
alpar@9 105 * LP presolver is used). */
alpar@9 106
alpar@9 107 static void trivial_lp(glp_prob *P, const glp_smcp *parm)
alpar@9 108 { /* solve trivial LP which has empty constraint matrix */
alpar@9 109 GLPROW *row;
alpar@9 110 GLPCOL *col;
alpar@9 111 int i, j;
alpar@9 112 double p_infeas, d_infeas, zeta;
alpar@9 113 P->valid = 0;
alpar@9 114 P->pbs_stat = P->dbs_stat = GLP_FEAS;
alpar@9 115 P->obj_val = P->c0;
alpar@9 116 P->some = 0;
alpar@9 117 p_infeas = d_infeas = 0.0;
alpar@9 118 /* make all auxiliary variables basic */
alpar@9 119 for (i = 1; i <= P->m; i++)
alpar@9 120 { row = P->row[i];
alpar@9 121 row->stat = GLP_BS;
alpar@9 122 row->prim = row->dual = 0.0;
alpar@9 123 /* check primal feasibility */
alpar@9 124 if (row->type == GLP_LO || row->type == GLP_DB ||
alpar@9 125 row->type == GLP_FX)
alpar@9 126 { /* row has lower bound */
alpar@9 127 if (row->lb > + parm->tol_bnd)
alpar@9 128 { P->pbs_stat = GLP_NOFEAS;
alpar@9 129 if (P->some == 0 && parm->meth != GLP_PRIMAL)
alpar@9 130 P->some = i;
alpar@9 131 }
alpar@9 132 if (p_infeas < + row->lb)
alpar@9 133 p_infeas = + row->lb;
alpar@9 134 }
alpar@9 135 if (row->type == GLP_UP || row->type == GLP_DB ||
alpar@9 136 row->type == GLP_FX)
alpar@9 137 { /* row has upper bound */
alpar@9 138 if (row->ub < - parm->tol_bnd)
alpar@9 139 { P->pbs_stat = GLP_NOFEAS;
alpar@9 140 if (P->some == 0 && parm->meth != GLP_PRIMAL)
alpar@9 141 P->some = i;
alpar@9 142 }
alpar@9 143 if (p_infeas < - row->ub)
alpar@9 144 p_infeas = - row->ub;
alpar@9 145 }
alpar@9 146 }
alpar@9 147 /* determine scale factor for the objective row */
alpar@9 148 zeta = 1.0;
alpar@9 149 for (j = 1; j <= P->n; j++)
alpar@9 150 { col = P->col[j];
alpar@9 151 if (zeta < fabs(col->coef)) zeta = fabs(col->coef);
alpar@9 152 }
alpar@9 153 zeta = (P->dir == GLP_MIN ? +1.0 : -1.0) / zeta;
alpar@9 154 /* make all structural variables non-basic */
alpar@9 155 for (j = 1; j <= P->n; j++)
alpar@9 156 { col = P->col[j];
alpar@9 157 if (col->type == GLP_FR)
alpar@9 158 col->stat = GLP_NF, col->prim = 0.0;
alpar@9 159 else if (col->type == GLP_LO)
alpar@9 160 lo: col->stat = GLP_NL, col->prim = col->lb;
alpar@9 161 else if (col->type == GLP_UP)
alpar@9 162 up: col->stat = GLP_NU, col->prim = col->ub;
alpar@9 163 else if (col->type == GLP_DB)
alpar@9 164 { if (zeta * col->coef > 0.0)
alpar@9 165 goto lo;
alpar@9 166 else if (zeta * col->coef < 0.0)
alpar@9 167 goto up;
alpar@9 168 else if (fabs(col->lb) <= fabs(col->ub))
alpar@9 169 goto lo;
alpar@9 170 else
alpar@9 171 goto up;
alpar@9 172 }
alpar@9 173 else if (col->type == GLP_FX)
alpar@9 174 col->stat = GLP_NS, col->prim = col->lb;
alpar@9 175 col->dual = col->coef;
alpar@9 176 P->obj_val += col->coef * col->prim;
alpar@9 177 /* check dual feasibility */
alpar@9 178 if (col->type == GLP_FR || col->type == GLP_LO)
alpar@9 179 { /* column has no upper bound */
alpar@9 180 if (zeta * col->dual < - parm->tol_dj)
alpar@9 181 { P->dbs_stat = GLP_NOFEAS;
alpar@9 182 if (P->some == 0 && parm->meth == GLP_PRIMAL)
alpar@9 183 P->some = P->m + j;
alpar@9 184 }
alpar@9 185 if (d_infeas < - zeta * col->dual)
alpar@9 186 d_infeas = - zeta * col->dual;
alpar@9 187 }
alpar@9 188 if (col->type == GLP_FR || col->type == GLP_UP)
alpar@9 189 { /* column has no lower bound */
alpar@9 190 if (zeta * col->dual > + parm->tol_dj)
alpar@9 191 { P->dbs_stat = GLP_NOFEAS;
alpar@9 192 if (P->some == 0 && parm->meth == GLP_PRIMAL)
alpar@9 193 P->some = P->m + j;
alpar@9 194 }
alpar@9 195 if (d_infeas < + zeta * col->dual)
alpar@9 196 d_infeas = + zeta * col->dual;
alpar@9 197 }
alpar@9 198 }
alpar@9 199 /* simulate the simplex solver output */
alpar@9 200 if (parm->msg_lev >= GLP_MSG_ON && parm->out_dly == 0)
alpar@9 201 { xprintf("~%6d: obj = %17.9e infeas = %10.3e\n", P->it_cnt,
alpar@9 202 P->obj_val, parm->meth == GLP_PRIMAL ? p_infeas : d_infeas);
alpar@9 203 }
alpar@9 204 if (parm->msg_lev >= GLP_MSG_ALL && parm->out_dly == 0)
alpar@9 205 { if (P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS)
alpar@9 206 xprintf("OPTIMAL SOLUTION FOUND\n");
alpar@9 207 else if (P->pbs_stat == GLP_NOFEAS)
alpar@9 208 xprintf("PROBLEM HAS NO FEASIBLE SOLUTION\n");
alpar@9 209 else if (parm->meth == GLP_PRIMAL)
alpar@9 210 xprintf("PROBLEM HAS UNBOUNDED SOLUTION\n");
alpar@9 211 else
alpar@9 212 xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n");
alpar@9 213 }
alpar@9 214 return;
alpar@9 215 }
alpar@9 216
alpar@9 217 static int solve_lp(glp_prob *P, const glp_smcp *parm)
alpar@9 218 { /* solve LP directly without using the preprocessor */
alpar@9 219 int ret;
alpar@9 220 if (!glp_bf_exists(P))
alpar@9 221 { ret = glp_factorize(P);
alpar@9 222 if (ret == 0)
alpar@9 223 ;
alpar@9 224 else if (ret == GLP_EBADB)
alpar@9 225 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9 226 xprintf("glp_simplex: initial basis is invalid\n");
alpar@9 227 }
alpar@9 228 else if (ret == GLP_ESING)
alpar@9 229 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9 230 xprintf("glp_simplex: initial basis is singular\n");
alpar@9 231 }
alpar@9 232 else if (ret == GLP_ECOND)
alpar@9 233 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9 234 xprintf(
alpar@9 235 "glp_simplex: initial basis is ill-conditioned\n");
alpar@9 236 }
alpar@9 237 else
alpar@9 238 xassert(ret != ret);
alpar@9 239 if (ret != 0) goto done;
alpar@9 240 }
alpar@9 241 if (parm->meth == GLP_PRIMAL)
alpar@9 242 ret = spx_primal(P, parm);
alpar@9 243 else if (parm->meth == GLP_DUALP)
alpar@9 244 { ret = spx_dual(P, parm);
alpar@9 245 if (ret == GLP_EFAIL && P->valid)
alpar@9 246 ret = spx_primal(P, parm);
alpar@9 247 }
alpar@9 248 else if (parm->meth == GLP_DUAL)
alpar@9 249 ret = spx_dual(P, parm);
alpar@9 250 else
alpar@9 251 xassert(parm != parm);
alpar@9 252 done: return ret;
alpar@9 253 }
alpar@9 254
alpar@9 255 static int preprocess_and_solve_lp(glp_prob *P, const glp_smcp *parm)
alpar@9 256 { /* solve LP using the preprocessor */
alpar@9 257 NPP *npp;
alpar@9 258 glp_prob *lp = NULL;
alpar@9 259 glp_bfcp bfcp;
alpar@9 260 int ret;
alpar@9 261 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9 262 xprintf("Preprocessing...\n");
alpar@9 263 /* create preprocessor workspace */
alpar@9 264 npp = npp_create_wksp();
alpar@9 265 /* load original problem into the preprocessor workspace */
alpar@9 266 npp_load_prob(npp, P, GLP_OFF, GLP_SOL, GLP_OFF);
alpar@9 267 /* process LP prior to applying primal/dual simplex method */
alpar@9 268 ret = npp_simplex(npp, parm);
alpar@9 269 if (ret == 0)
alpar@9 270 ;
alpar@9 271 else if (ret == GLP_ENOPFS)
alpar@9 272 { if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9 273 xprintf("PROBLEM HAS NO PRIMAL FEASIBLE SOLUTION\n");
alpar@9 274 }
alpar@9 275 else if (ret == GLP_ENODFS)
alpar@9 276 { if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9 277 xprintf("PROBLEM HAS NO DUAL FEASIBLE SOLUTION\n");
alpar@9 278 }
alpar@9 279 else
alpar@9 280 xassert(ret != ret);
alpar@9 281 if (ret != 0) goto done;
alpar@9 282 /* build transformed LP */
alpar@9 283 lp = glp_create_prob();
alpar@9 284 npp_build_prob(npp, lp);
alpar@9 285 /* if the transformed LP is empty, it has empty solution, which
alpar@9 286 is optimal */
alpar@9 287 if (lp->m == 0 && lp->n == 0)
alpar@9 288 { lp->pbs_stat = lp->dbs_stat = GLP_FEAS;
alpar@9 289 lp->obj_val = lp->c0;
alpar@9 290 if (parm->msg_lev >= GLP_MSG_ON && parm->out_dly == 0)
alpar@9 291 { xprintf("~%6d: obj = %17.9e infeas = %10.3e\n", P->it_cnt,
alpar@9 292 lp->obj_val, 0.0);
alpar@9 293 }
alpar@9 294 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9 295 xprintf("OPTIMAL SOLUTION FOUND BY LP PREPROCESSOR\n");
alpar@9 296 goto post;
alpar@9 297 }
alpar@9 298 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9 299 { xprintf("%d row%s, %d column%s, %d non-zero%s\n",
alpar@9 300 lp->m, lp->m == 1 ? "" : "s", lp->n, lp->n == 1 ? "" : "s",
alpar@9 301 lp->nnz, lp->nnz == 1 ? "" : "s");
alpar@9 302 }
alpar@9 303 /* inherit basis factorization control parameters */
alpar@9 304 glp_get_bfcp(P, &bfcp);
alpar@9 305 glp_set_bfcp(lp, &bfcp);
alpar@9 306 /* scale the transformed problem */
alpar@9 307 { ENV *env = get_env_ptr();
alpar@9 308 int term_out = env->term_out;
alpar@9 309 if (!term_out || parm->msg_lev < GLP_MSG_ALL)
alpar@9 310 env->term_out = GLP_OFF;
alpar@9 311 else
alpar@9 312 env->term_out = GLP_ON;
alpar@9 313 glp_scale_prob(lp, GLP_SF_AUTO);
alpar@9 314 env->term_out = term_out;
alpar@9 315 }
alpar@9 316 /* build advanced initial basis */
alpar@9 317 { ENV *env = get_env_ptr();
alpar@9 318 int term_out = env->term_out;
alpar@9 319 if (!term_out || parm->msg_lev < GLP_MSG_ALL)
alpar@9 320 env->term_out = GLP_OFF;
alpar@9 321 else
alpar@9 322 env->term_out = GLP_ON;
alpar@9 323 glp_adv_basis(lp, 0);
alpar@9 324 env->term_out = term_out;
alpar@9 325 }
alpar@9 326 /* solve the transformed LP */
alpar@9 327 lp->it_cnt = P->it_cnt;
alpar@9 328 ret = solve_lp(lp, parm);
alpar@9 329 P->it_cnt = lp->it_cnt;
alpar@9 330 /* only optimal solution can be postprocessed */
alpar@9 331 if (!(ret == 0 && lp->pbs_stat == GLP_FEAS && lp->dbs_stat ==
alpar@9 332 GLP_FEAS))
alpar@9 333 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9 334 xprintf("glp_simplex: unable to recover undefined or non-op"
alpar@9 335 "timal solution\n");
alpar@9 336 if (ret == 0)
alpar@9 337 { if (lp->pbs_stat == GLP_NOFEAS)
alpar@9 338 ret = GLP_ENOPFS;
alpar@9 339 else if (lp->dbs_stat == GLP_NOFEAS)
alpar@9 340 ret = GLP_ENODFS;
alpar@9 341 else
alpar@9 342 xassert(lp != lp);
alpar@9 343 }
alpar@9 344 goto done;
alpar@9 345 }
alpar@9 346 post: /* postprocess solution from the transformed LP */
alpar@9 347 npp_postprocess(npp, lp);
alpar@9 348 /* the transformed LP is no longer needed */
alpar@9 349 glp_delete_prob(lp), lp = NULL;
alpar@9 350 /* store solution to the original problem */
alpar@9 351 npp_unload_sol(npp, P);
alpar@9 352 /* the original LP has been successfully solved */
alpar@9 353 ret = 0;
alpar@9 354 done: /* delete the transformed LP, if it exists */
alpar@9 355 if (lp != NULL) glp_delete_prob(lp);
alpar@9 356 /* delete preprocessor workspace */
alpar@9 357 npp_delete_wksp(npp);
alpar@9 358 return ret;
alpar@9 359 }
alpar@9 360
alpar@9 361 int glp_simplex(glp_prob *P, const glp_smcp *parm)
alpar@9 362 { /* solve LP problem with the simplex method */
alpar@9 363 glp_smcp _parm;
alpar@9 364 int i, j, ret;
alpar@9 365 /* check problem object */
alpar@9 366 if (P == NULL || P->magic != GLP_PROB_MAGIC)
alpar@9 367 xerror("glp_simplex: P = %p; invalid problem object\n", P);
alpar@9 368 if (P->tree != NULL && P->tree->reason != 0)
alpar@9 369 xerror("glp_simplex: operation not allowed\n");
alpar@9 370 /* check control parameters */
alpar@9 371 if (parm == NULL)
alpar@9 372 parm = &_parm, glp_init_smcp((glp_smcp *)parm);
alpar@9 373 if (!(parm->msg_lev == GLP_MSG_OFF ||
alpar@9 374 parm->msg_lev == GLP_MSG_ERR ||
alpar@9 375 parm->msg_lev == GLP_MSG_ON ||
alpar@9 376 parm->msg_lev == GLP_MSG_ALL ||
alpar@9 377 parm->msg_lev == GLP_MSG_DBG))
alpar@9 378 xerror("glp_simplex: msg_lev = %d; invalid parameter\n",
alpar@9 379 parm->msg_lev);
alpar@9 380 if (!(parm->meth == GLP_PRIMAL ||
alpar@9 381 parm->meth == GLP_DUALP ||
alpar@9 382 parm->meth == GLP_DUAL))
alpar@9 383 xerror("glp_simplex: meth = %d; invalid parameter\n",
alpar@9 384 parm->meth);
alpar@9 385 if (!(parm->pricing == GLP_PT_STD ||
alpar@9 386 parm->pricing == GLP_PT_PSE))
alpar@9 387 xerror("glp_simplex: pricing = %d; invalid parameter\n",
alpar@9 388 parm->pricing);
alpar@9 389 if (!(parm->r_test == GLP_RT_STD ||
alpar@9 390 parm->r_test == GLP_RT_HAR))
alpar@9 391 xerror("glp_simplex: r_test = %d; invalid parameter\n",
alpar@9 392 parm->r_test);
alpar@9 393 if (!(0.0 < parm->tol_bnd && parm->tol_bnd < 1.0))
alpar@9 394 xerror("glp_simplex: tol_bnd = %g; invalid parameter\n",
alpar@9 395 parm->tol_bnd);
alpar@9 396 if (!(0.0 < parm->tol_dj && parm->tol_dj < 1.0))
alpar@9 397 xerror("glp_simplex: tol_dj = %g; invalid parameter\n",
alpar@9 398 parm->tol_dj);
alpar@9 399 if (!(0.0 < parm->tol_piv && parm->tol_piv < 1.0))
alpar@9 400 xerror("glp_simplex: tol_piv = %g; invalid parameter\n",
alpar@9 401 parm->tol_piv);
alpar@9 402 if (parm->it_lim < 0)
alpar@9 403 xerror("glp_simplex: it_lim = %d; invalid parameter\n",
alpar@9 404 parm->it_lim);
alpar@9 405 if (parm->tm_lim < 0)
alpar@9 406 xerror("glp_simplex: tm_lim = %d; invalid parameter\n",
alpar@9 407 parm->tm_lim);
alpar@9 408 if (parm->out_frq < 1)
alpar@9 409 xerror("glp_simplex: out_frq = %d; invalid parameter\n",
alpar@9 410 parm->out_frq);
alpar@9 411 if (parm->out_dly < 0)
alpar@9 412 xerror("glp_simplex: out_dly = %d; invalid parameter\n",
alpar@9 413 parm->out_dly);
alpar@9 414 if (!(parm->presolve == GLP_ON || parm->presolve == GLP_OFF))
alpar@9 415 xerror("glp_simplex: presolve = %d; invalid parameter\n",
alpar@9 416 parm->presolve);
alpar@9 417 /* basic solution is currently undefined */
alpar@9 418 P->pbs_stat = P->dbs_stat = GLP_UNDEF;
alpar@9 419 P->obj_val = 0.0;
alpar@9 420 P->some = 0;
alpar@9 421 /* check bounds of double-bounded variables */
alpar@9 422 for (i = 1; i <= P->m; i++)
alpar@9 423 { GLPROW *row = P->row[i];
alpar@9 424 if (row->type == GLP_DB && row->lb >= row->ub)
alpar@9 425 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9 426 xprintf("glp_simplex: row %d: lb = %g, ub = %g; incorrec"
alpar@9 427 "t bounds\n", i, row->lb, row->ub);
alpar@9 428 ret = GLP_EBOUND;
alpar@9 429 goto done;
alpar@9 430 }
alpar@9 431 }
alpar@9 432 for (j = 1; j <= P->n; j++)
alpar@9 433 { GLPCOL *col = P->col[j];
alpar@9 434 if (col->type == GLP_DB && col->lb >= col->ub)
alpar@9 435 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9 436 xprintf("glp_simplex: column %d: lb = %g, ub = %g; incor"
alpar@9 437 "rect bounds\n", j, col->lb, col->ub);
alpar@9 438 ret = GLP_EBOUND;
alpar@9 439 goto done;
alpar@9 440 }
alpar@9 441 }
alpar@9 442 /* solve LP problem */
alpar@9 443 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9 444 { xprintf("GLPK Simplex Optimizer, v%s\n", glp_version());
alpar@9 445 xprintf("%d row%s, %d column%s, %d non-zero%s\n",
alpar@9 446 P->m, P->m == 1 ? "" : "s", P->n, P->n == 1 ? "" : "s",
alpar@9 447 P->nnz, P->nnz == 1 ? "" : "s");
alpar@9 448 }
alpar@9 449 if (P->nnz == 0)
alpar@9 450 trivial_lp(P, parm), ret = 0;
alpar@9 451 else if (!parm->presolve)
alpar@9 452 ret = solve_lp(P, parm);
alpar@9 453 else
alpar@9 454 ret = preprocess_and_solve_lp(P, parm);
alpar@9 455 done: /* return to the application program */
alpar@9 456 return ret;
alpar@9 457 }
alpar@9 458
alpar@9 459 /***********************************************************************
alpar@9 460 * NAME
alpar@9 461 *
alpar@9 462 * glp_init_smcp - initialize simplex method control parameters
alpar@9 463 *
alpar@9 464 * SYNOPSIS
alpar@9 465 *
alpar@9 466 * void glp_init_smcp(glp_smcp *parm);
alpar@9 467 *
alpar@9 468 * DESCRIPTION
alpar@9 469 *
alpar@9 470 * The routine glp_init_smcp initializes control parameters, which are
alpar@9 471 * used by the simplex solver, with default values.
alpar@9 472 *
alpar@9 473 * Default values of the control parameters are stored in a glp_smcp
alpar@9 474 * structure, which the parameter parm points to. */
alpar@9 475
alpar@9 476 void glp_init_smcp(glp_smcp *parm)
alpar@9 477 { parm->msg_lev = GLP_MSG_ALL;
alpar@9 478 parm->meth = GLP_PRIMAL;
alpar@9 479 parm->pricing = GLP_PT_PSE;
alpar@9 480 parm->r_test = GLP_RT_HAR;
alpar@9 481 parm->tol_bnd = 1e-7;
alpar@9 482 parm->tol_dj = 1e-7;
alpar@9 483 parm->tol_piv = 1e-10;
alpar@9 484 parm->obj_ll = -DBL_MAX;
alpar@9 485 parm->obj_ul = +DBL_MAX;
alpar@9 486 parm->it_lim = INT_MAX;
alpar@9 487 parm->tm_lim = INT_MAX;
alpar@9 488 parm->out_frq = 500;
alpar@9 489 parm->out_dly = 0;
alpar@9 490 parm->presolve = GLP_OFF;
alpar@9 491 return;
alpar@9 492 }
alpar@9 493
alpar@9 494 /***********************************************************************
alpar@9 495 * NAME
alpar@9 496 *
alpar@9 497 * glp_get_status - retrieve generic status of basic solution
alpar@9 498 *
alpar@9 499 * SYNOPSIS
alpar@9 500 *
alpar@9 501 * int glp_get_status(glp_prob *lp);
alpar@9 502 *
alpar@9 503 * RETURNS
alpar@9 504 *
alpar@9 505 * The routine glp_get_status reports the generic status of the basic
alpar@9 506 * solution for the specified problem object as follows:
alpar@9 507 *
alpar@9 508 * GLP_OPT - solution is optimal;
alpar@9 509 * GLP_FEAS - solution is feasible;
alpar@9 510 * GLP_INFEAS - solution is infeasible;
alpar@9 511 * GLP_NOFEAS - problem has no feasible solution;
alpar@9 512 * GLP_UNBND - problem has unbounded solution;
alpar@9 513 * GLP_UNDEF - solution is undefined. */
alpar@9 514
alpar@9 515 int glp_get_status(glp_prob *lp)
alpar@9 516 { int status;
alpar@9 517 status = glp_get_prim_stat(lp);
alpar@9 518 switch (status)
alpar@9 519 { case GLP_FEAS:
alpar@9 520 switch (glp_get_dual_stat(lp))
alpar@9 521 { case GLP_FEAS:
alpar@9 522 status = GLP_OPT;
alpar@9 523 break;
alpar@9 524 case GLP_NOFEAS:
alpar@9 525 status = GLP_UNBND;
alpar@9 526 break;
alpar@9 527 case GLP_UNDEF:
alpar@9 528 case GLP_INFEAS:
alpar@9 529 status = status;
alpar@9 530 break;
alpar@9 531 default:
alpar@9 532 xassert(lp != lp);
alpar@9 533 }
alpar@9 534 break;
alpar@9 535 case GLP_UNDEF:
alpar@9 536 case GLP_INFEAS:
alpar@9 537 case GLP_NOFEAS:
alpar@9 538 status = status;
alpar@9 539 break;
alpar@9 540 default:
alpar@9 541 xassert(lp != lp);
alpar@9 542 }
alpar@9 543 return status;
alpar@9 544 }
alpar@9 545
alpar@9 546 /***********************************************************************
alpar@9 547 * NAME
alpar@9 548 *
alpar@9 549 * glp_get_prim_stat - retrieve status of primal basic solution
alpar@9 550 *
alpar@9 551 * SYNOPSIS
alpar@9 552 *
alpar@9 553 * int glp_get_prim_stat(glp_prob *lp);
alpar@9 554 *
alpar@9 555 * RETURNS
alpar@9 556 *
alpar@9 557 * The routine glp_get_prim_stat reports the status of the primal basic
alpar@9 558 * solution for the specified problem object as follows:
alpar@9 559 *
alpar@9 560 * GLP_UNDEF - primal solution is undefined;
alpar@9 561 * GLP_FEAS - primal solution is feasible;
alpar@9 562 * GLP_INFEAS - primal solution is infeasible;
alpar@9 563 * GLP_NOFEAS - no primal feasible solution exists. */
alpar@9 564
alpar@9 565 int glp_get_prim_stat(glp_prob *lp)
alpar@9 566 { int pbs_stat = lp->pbs_stat;
alpar@9 567 return pbs_stat;
alpar@9 568 }
alpar@9 569
alpar@9 570 /***********************************************************************
alpar@9 571 * NAME
alpar@9 572 *
alpar@9 573 * glp_get_dual_stat - retrieve status of dual basic solution
alpar@9 574 *
alpar@9 575 * SYNOPSIS
alpar@9 576 *
alpar@9 577 * int glp_get_dual_stat(glp_prob *lp);
alpar@9 578 *
alpar@9 579 * RETURNS
alpar@9 580 *
alpar@9 581 * The routine glp_get_dual_stat reports the status of the dual basic
alpar@9 582 * solution for the specified problem object as follows:
alpar@9 583 *
alpar@9 584 * GLP_UNDEF - dual solution is undefined;
alpar@9 585 * GLP_FEAS - dual solution is feasible;
alpar@9 586 * GLP_INFEAS - dual solution is infeasible;
alpar@9 587 * GLP_NOFEAS - no dual feasible solution exists. */
alpar@9 588
alpar@9 589 int glp_get_dual_stat(glp_prob *lp)
alpar@9 590 { int dbs_stat = lp->dbs_stat;
alpar@9 591 return dbs_stat;
alpar@9 592 }
alpar@9 593
alpar@9 594 /***********************************************************************
alpar@9 595 * NAME
alpar@9 596 *
alpar@9 597 * glp_get_obj_val - retrieve objective value (basic solution)
alpar@9 598 *
alpar@9 599 * SYNOPSIS
alpar@9 600 *
alpar@9 601 * double glp_get_obj_val(glp_prob *lp);
alpar@9 602 *
alpar@9 603 * RETURNS
alpar@9 604 *
alpar@9 605 * The routine glp_get_obj_val returns value of the objective function
alpar@9 606 * for basic solution. */
alpar@9 607
alpar@9 608 double glp_get_obj_val(glp_prob *lp)
alpar@9 609 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 610 double z;
alpar@9 611 z = lp->obj_val;
alpar@9 612 /*if (cps->round && fabs(z) < 1e-9) z = 0.0;*/
alpar@9 613 return z;
alpar@9 614 }
alpar@9 615
alpar@9 616 /***********************************************************************
alpar@9 617 * NAME
alpar@9 618 *
alpar@9 619 * glp_get_row_stat - retrieve row status
alpar@9 620 *
alpar@9 621 * SYNOPSIS
alpar@9 622 *
alpar@9 623 * int glp_get_row_stat(glp_prob *lp, int i);
alpar@9 624 *
alpar@9 625 * RETURNS
alpar@9 626 *
alpar@9 627 * The routine glp_get_row_stat returns current status assigned to the
alpar@9 628 * auxiliary variable associated with i-th row as follows:
alpar@9 629 *
alpar@9 630 * GLP_BS - basic variable;
alpar@9 631 * GLP_NL - non-basic variable on its lower bound;
alpar@9 632 * GLP_NU - non-basic variable on its upper bound;
alpar@9 633 * GLP_NF - non-basic free (unbounded) variable;
alpar@9 634 * GLP_NS - non-basic fixed variable. */
alpar@9 635
alpar@9 636 int glp_get_row_stat(glp_prob *lp, int i)
alpar@9 637 { if (!(1 <= i && i <= lp->m))
alpar@9 638 xerror("glp_get_row_stat: i = %d; row number out of range\n",
alpar@9 639 i);
alpar@9 640 return lp->row[i]->stat;
alpar@9 641 }
alpar@9 642
alpar@9 643 /***********************************************************************
alpar@9 644 * NAME
alpar@9 645 *
alpar@9 646 * glp_get_row_prim - retrieve row primal value (basic solution)
alpar@9 647 *
alpar@9 648 * SYNOPSIS
alpar@9 649 *
alpar@9 650 * double glp_get_row_prim(glp_prob *lp, int i);
alpar@9 651 *
alpar@9 652 * RETURNS
alpar@9 653 *
alpar@9 654 * The routine glp_get_row_prim returns primal value of the auxiliary
alpar@9 655 * variable associated with i-th row. */
alpar@9 656
alpar@9 657 double glp_get_row_prim(glp_prob *lp, int i)
alpar@9 658 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 659 double prim;
alpar@9 660 if (!(1 <= i && i <= lp->m))
alpar@9 661 xerror("glp_get_row_prim: i = %d; row number out of range\n",
alpar@9 662 i);
alpar@9 663 prim = lp->row[i]->prim;
alpar@9 664 /*if (cps->round && fabs(prim) < 1e-9) prim = 0.0;*/
alpar@9 665 return prim;
alpar@9 666 }
alpar@9 667
alpar@9 668 /***********************************************************************
alpar@9 669 * NAME
alpar@9 670 *
alpar@9 671 * glp_get_row_dual - retrieve row dual value (basic solution)
alpar@9 672 *
alpar@9 673 * SYNOPSIS
alpar@9 674 *
alpar@9 675 * double glp_get_row_dual(glp_prob *lp, int i);
alpar@9 676 *
alpar@9 677 * RETURNS
alpar@9 678 *
alpar@9 679 * The routine glp_get_row_dual returns dual value (i.e. reduced cost)
alpar@9 680 * of the auxiliary variable associated with i-th row. */
alpar@9 681
alpar@9 682 double glp_get_row_dual(glp_prob *lp, int i)
alpar@9 683 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 684 double dual;
alpar@9 685 if (!(1 <= i && i <= lp->m))
alpar@9 686 xerror("glp_get_row_dual: i = %d; row number out of range\n",
alpar@9 687 i);
alpar@9 688 dual = lp->row[i]->dual;
alpar@9 689 /*if (cps->round && fabs(dual) < 1e-9) dual = 0.0;*/
alpar@9 690 return dual;
alpar@9 691 }
alpar@9 692
alpar@9 693 /***********************************************************************
alpar@9 694 * NAME
alpar@9 695 *
alpar@9 696 * glp_get_col_stat - retrieve column status
alpar@9 697 *
alpar@9 698 * SYNOPSIS
alpar@9 699 *
alpar@9 700 * int glp_get_col_stat(glp_prob *lp, int j);
alpar@9 701 *
alpar@9 702 * RETURNS
alpar@9 703 *
alpar@9 704 * The routine glp_get_col_stat returns current status assigned to the
alpar@9 705 * structural variable associated with j-th column as follows:
alpar@9 706 *
alpar@9 707 * GLP_BS - basic variable;
alpar@9 708 * GLP_NL - non-basic variable on its lower bound;
alpar@9 709 * GLP_NU - non-basic variable on its upper bound;
alpar@9 710 * GLP_NF - non-basic free (unbounded) variable;
alpar@9 711 * GLP_NS - non-basic fixed variable. */
alpar@9 712
alpar@9 713 int glp_get_col_stat(glp_prob *lp, int j)
alpar@9 714 { if (!(1 <= j && j <= lp->n))
alpar@9 715 xerror("glp_get_col_stat: j = %d; column number out of range\n"
alpar@9 716 , j);
alpar@9 717 return lp->col[j]->stat;
alpar@9 718 }
alpar@9 719
alpar@9 720 /***********************************************************************
alpar@9 721 * NAME
alpar@9 722 *
alpar@9 723 * glp_get_col_prim - retrieve column primal value (basic solution)
alpar@9 724 *
alpar@9 725 * SYNOPSIS
alpar@9 726 *
alpar@9 727 * double glp_get_col_prim(glp_prob *lp, int j);
alpar@9 728 *
alpar@9 729 * RETURNS
alpar@9 730 *
alpar@9 731 * The routine glp_get_col_prim returns primal value of the structural
alpar@9 732 * variable associated with j-th column. */
alpar@9 733
alpar@9 734 double glp_get_col_prim(glp_prob *lp, int j)
alpar@9 735 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 736 double prim;
alpar@9 737 if (!(1 <= j && j <= lp->n))
alpar@9 738 xerror("glp_get_col_prim: j = %d; column number out of range\n"
alpar@9 739 , j);
alpar@9 740 prim = lp->col[j]->prim;
alpar@9 741 /*if (cps->round && fabs(prim) < 1e-9) prim = 0.0;*/
alpar@9 742 return prim;
alpar@9 743 }
alpar@9 744
alpar@9 745 /***********************************************************************
alpar@9 746 * NAME
alpar@9 747 *
alpar@9 748 * glp_get_col_dual - retrieve column dual value (basic solution)
alpar@9 749 *
alpar@9 750 * SYNOPSIS
alpar@9 751 *
alpar@9 752 * double glp_get_col_dual(glp_prob *lp, int j);
alpar@9 753 *
alpar@9 754 * RETURNS
alpar@9 755 *
alpar@9 756 * The routine glp_get_col_dual returns dual value (i.e. reduced cost)
alpar@9 757 * of the structural variable associated with j-th column. */
alpar@9 758
alpar@9 759 double glp_get_col_dual(glp_prob *lp, int j)
alpar@9 760 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 761 double dual;
alpar@9 762 if (!(1 <= j && j <= lp->n))
alpar@9 763 xerror("glp_get_col_dual: j = %d; column number out of range\n"
alpar@9 764 , j);
alpar@9 765 dual = lp->col[j]->dual;
alpar@9 766 /*if (cps->round && fabs(dual) < 1e-9) dual = 0.0;*/
alpar@9 767 return dual;
alpar@9 768 }
alpar@9 769
alpar@9 770 /***********************************************************************
alpar@9 771 * NAME
alpar@9 772 *
alpar@9 773 * glp_get_unbnd_ray - determine variable causing unboundedness
alpar@9 774 *
alpar@9 775 * SYNOPSIS
alpar@9 776 *
alpar@9 777 * int glp_get_unbnd_ray(glp_prob *lp);
alpar@9 778 *
alpar@9 779 * RETURNS
alpar@9 780 *
alpar@9 781 * The routine glp_get_unbnd_ray returns the number k of a variable,
alpar@9 782 * which causes primal or dual unboundedness. If 1 <= k <= m, it is
alpar@9 783 * k-th auxiliary variable, and if m+1 <= k <= m+n, it is (k-m)-th
alpar@9 784 * structural variable, where m is the number of rows, n is the number
alpar@9 785 * of columns in the problem object. If such variable is not defined,
alpar@9 786 * the routine returns 0.
alpar@9 787 *
alpar@9 788 * COMMENTS
alpar@9 789 *
alpar@9 790 * If it is not exactly known which version of the simplex solver
alpar@9 791 * detected unboundedness, i.e. whether the unboundedness is primal or
alpar@9 792 * dual, it is sufficient to check the status of the variable reported
alpar@9 793 * with the routine glp_get_row_stat or glp_get_col_stat. If the
alpar@9 794 * variable is non-basic, the unboundedness is primal, otherwise, if
alpar@9 795 * the variable is basic, the unboundedness is dual (the latter case
alpar@9 796 * means that the problem has no primal feasible dolution). */
alpar@9 797
alpar@9 798 int glp_get_unbnd_ray(glp_prob *lp)
alpar@9 799 { int k;
alpar@9 800 k = lp->some;
alpar@9 801 xassert(k >= 0);
alpar@9 802 if (k > lp->m + lp->n) k = 0;
alpar@9 803 return k;
alpar@9 804 }
alpar@9 805
alpar@9 806 /* eof */