lemon-project-template-glpk

annotate deps/glpk/src/glpapi07.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 /* glpapi07.c (exact simplex solver) */
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 "glpapi.h"
alpar@9 26 #include "glpssx.h"
alpar@9 27
alpar@9 28 /***********************************************************************
alpar@9 29 * NAME
alpar@9 30 *
alpar@9 31 * glp_exact - solve LP problem in exact arithmetic
alpar@9 32 *
alpar@9 33 * SYNOPSIS
alpar@9 34 *
alpar@9 35 * int glp_exact(glp_prob *lp, const glp_smcp *parm);
alpar@9 36 *
alpar@9 37 * DESCRIPTION
alpar@9 38 *
alpar@9 39 * The routine glp_exact is a tentative implementation of the primal
alpar@9 40 * two-phase simplex method based on exact (rational) arithmetic. It is
alpar@9 41 * similar to the routine glp_simplex, however, for all internal
alpar@9 42 * computations it uses arithmetic of rational numbers, which is exact
alpar@9 43 * in mathematical sense, i.e. free of round-off errors unlike floating
alpar@9 44 * point arithmetic.
alpar@9 45 *
alpar@9 46 * Note that the routine glp_exact uses inly two control parameters
alpar@9 47 * passed in the structure glp_smcp, namely, it_lim and tm_lim.
alpar@9 48 *
alpar@9 49 * RETURNS
alpar@9 50 *
alpar@9 51 * 0 The LP problem instance has been successfully solved. This code
alpar@9 52 * does not necessarily mean that the solver has found optimal
alpar@9 53 * solution. It only means that the solution process was successful.
alpar@9 54 *
alpar@9 55 * GLP_EBADB
alpar@9 56 * Unable to start the search, because the initial basis specified
alpar@9 57 * in the problem object is invalid--the number of basic (auxiliary
alpar@9 58 * and structural) variables is not the same as the number of rows in
alpar@9 59 * the problem object.
alpar@9 60 *
alpar@9 61 * GLP_ESING
alpar@9 62 * Unable to start the search, because the basis matrix correspodning
alpar@9 63 * to the initial basis is exactly singular.
alpar@9 64 *
alpar@9 65 * GLP_EBOUND
alpar@9 66 * Unable to start the search, because some double-bounded variables
alpar@9 67 * have incorrect bounds.
alpar@9 68 *
alpar@9 69 * GLP_EFAIL
alpar@9 70 * The problem has no rows/columns.
alpar@9 71 *
alpar@9 72 * GLP_EITLIM
alpar@9 73 * The search was prematurely terminated, because the simplex
alpar@9 74 * iteration limit has been exceeded.
alpar@9 75 *
alpar@9 76 * GLP_ETMLIM
alpar@9 77 * The search was prematurely terminated, because the time limit has
alpar@9 78 * been exceeded. */
alpar@9 79
alpar@9 80 static void set_d_eps(mpq_t x, double val)
alpar@9 81 { /* convert double val to rational x obtaining a more adequate
alpar@9 82 fraction than provided by mpq_set_d due to allowing a small
alpar@9 83 approximation error specified by a given relative tolerance;
alpar@9 84 for example, mpq_set_d would give the following
alpar@9 85 1/3 ~= 0.333333333333333314829616256247391... ->
alpar@9 86 -> 6004799503160661/18014398509481984
alpar@9 87 while this routine gives exactly 1/3 */
alpar@9 88 int s, n, j;
alpar@9 89 double f, p, q, eps = 1e-9;
alpar@9 90 mpq_t temp;
alpar@9 91 xassert(-DBL_MAX <= val && val <= +DBL_MAX);
alpar@9 92 #if 1 /* 30/VII-2008 */
alpar@9 93 if (val == floor(val))
alpar@9 94 { /* if val is integral, do not approximate */
alpar@9 95 mpq_set_d(x, val);
alpar@9 96 goto done;
alpar@9 97 }
alpar@9 98 #endif
alpar@9 99 if (val > 0.0)
alpar@9 100 s = +1;
alpar@9 101 else if (val < 0.0)
alpar@9 102 s = -1;
alpar@9 103 else
alpar@9 104 { mpq_set_si(x, 0, 1);
alpar@9 105 goto done;
alpar@9 106 }
alpar@9 107 f = frexp(fabs(val), &n);
alpar@9 108 /* |val| = f * 2^n, where 0.5 <= f < 1.0 */
alpar@9 109 fp2rat(f, 0.1 * eps, &p, &q);
alpar@9 110 /* f ~= p / q, where p and q are integers */
alpar@9 111 mpq_init(temp);
alpar@9 112 mpq_set_d(x, p);
alpar@9 113 mpq_set_d(temp, q);
alpar@9 114 mpq_div(x, x, temp);
alpar@9 115 mpq_set_si(temp, 1, 1);
alpar@9 116 for (j = 1; j <= abs(n); j++)
alpar@9 117 mpq_add(temp, temp, temp);
alpar@9 118 if (n > 0)
alpar@9 119 mpq_mul(x, x, temp);
alpar@9 120 else if (n < 0)
alpar@9 121 mpq_div(x, x, temp);
alpar@9 122 mpq_clear(temp);
alpar@9 123 if (s < 0) mpq_neg(x, x);
alpar@9 124 /* check that the desired tolerance has been attained */
alpar@9 125 xassert(fabs(val - mpq_get_d(x)) <= eps * (1.0 + fabs(val)));
alpar@9 126 done: return;
alpar@9 127 }
alpar@9 128
alpar@9 129 static void load_data(SSX *ssx, LPX *lp)
alpar@9 130 { /* load LP problem data into simplex solver workspace */
alpar@9 131 int m = ssx->m;
alpar@9 132 int n = ssx->n;
alpar@9 133 int nnz = ssx->A_ptr[n+1]-1;
alpar@9 134 int j, k, type, loc, len, *ind;
alpar@9 135 double lb, ub, coef, *val;
alpar@9 136 xassert(lpx_get_num_rows(lp) == m);
alpar@9 137 xassert(lpx_get_num_cols(lp) == n);
alpar@9 138 xassert(lpx_get_num_nz(lp) == nnz);
alpar@9 139 /* types and bounds of rows and columns */
alpar@9 140 for (k = 1; k <= m+n; k++)
alpar@9 141 { if (k <= m)
alpar@9 142 { type = lpx_get_row_type(lp, k);
alpar@9 143 lb = lpx_get_row_lb(lp, k);
alpar@9 144 ub = lpx_get_row_ub(lp, k);
alpar@9 145 }
alpar@9 146 else
alpar@9 147 { type = lpx_get_col_type(lp, k-m);
alpar@9 148 lb = lpx_get_col_lb(lp, k-m);
alpar@9 149 ub = lpx_get_col_ub(lp, k-m);
alpar@9 150 }
alpar@9 151 switch (type)
alpar@9 152 { case LPX_FR: type = SSX_FR; break;
alpar@9 153 case LPX_LO: type = SSX_LO; break;
alpar@9 154 case LPX_UP: type = SSX_UP; break;
alpar@9 155 case LPX_DB: type = SSX_DB; break;
alpar@9 156 case LPX_FX: type = SSX_FX; break;
alpar@9 157 default: xassert(type != type);
alpar@9 158 }
alpar@9 159 ssx->type[k] = type;
alpar@9 160 set_d_eps(ssx->lb[k], lb);
alpar@9 161 set_d_eps(ssx->ub[k], ub);
alpar@9 162 }
alpar@9 163 /* optimization direction */
alpar@9 164 switch (lpx_get_obj_dir(lp))
alpar@9 165 { case LPX_MIN: ssx->dir = SSX_MIN; break;
alpar@9 166 case LPX_MAX: ssx->dir = SSX_MAX; break;
alpar@9 167 default: xassert(lp != lp);
alpar@9 168 }
alpar@9 169 /* objective coefficients */
alpar@9 170 for (k = 0; k <= m+n; k++)
alpar@9 171 { if (k == 0)
alpar@9 172 coef = lpx_get_obj_coef(lp, 0);
alpar@9 173 else if (k <= m)
alpar@9 174 coef = 0.0;
alpar@9 175 else
alpar@9 176 coef = lpx_get_obj_coef(lp, k-m);
alpar@9 177 set_d_eps(ssx->coef[k], coef);
alpar@9 178 }
alpar@9 179 /* constraint coefficients */
alpar@9 180 ind = xcalloc(1+m, sizeof(int));
alpar@9 181 val = xcalloc(1+m, sizeof(double));
alpar@9 182 loc = 0;
alpar@9 183 for (j = 1; j <= n; j++)
alpar@9 184 { ssx->A_ptr[j] = loc+1;
alpar@9 185 len = lpx_get_mat_col(lp, j, ind, val);
alpar@9 186 for (k = 1; k <= len; k++)
alpar@9 187 { loc++;
alpar@9 188 ssx->A_ind[loc] = ind[k];
alpar@9 189 set_d_eps(ssx->A_val[loc], val[k]);
alpar@9 190 }
alpar@9 191 }
alpar@9 192 xassert(loc == nnz);
alpar@9 193 xfree(ind);
alpar@9 194 xfree(val);
alpar@9 195 return;
alpar@9 196 }
alpar@9 197
alpar@9 198 static int load_basis(SSX *ssx, LPX *lp)
alpar@9 199 { /* load current LP basis into simplex solver workspace */
alpar@9 200 int m = ssx->m;
alpar@9 201 int n = ssx->n;
alpar@9 202 int *type = ssx->type;
alpar@9 203 int *stat = ssx->stat;
alpar@9 204 int *Q_row = ssx->Q_row;
alpar@9 205 int *Q_col = ssx->Q_col;
alpar@9 206 int i, j, k;
alpar@9 207 xassert(lpx_get_num_rows(lp) == m);
alpar@9 208 xassert(lpx_get_num_cols(lp) == n);
alpar@9 209 /* statuses of rows and columns */
alpar@9 210 for (k = 1; k <= m+n; k++)
alpar@9 211 { if (k <= m)
alpar@9 212 stat[k] = lpx_get_row_stat(lp, k);
alpar@9 213 else
alpar@9 214 stat[k] = lpx_get_col_stat(lp, k-m);
alpar@9 215 switch (stat[k])
alpar@9 216 { case LPX_BS:
alpar@9 217 stat[k] = SSX_BS;
alpar@9 218 break;
alpar@9 219 case LPX_NL:
alpar@9 220 stat[k] = SSX_NL;
alpar@9 221 xassert(type[k] == SSX_LO || type[k] == SSX_DB);
alpar@9 222 break;
alpar@9 223 case LPX_NU:
alpar@9 224 stat[k] = SSX_NU;
alpar@9 225 xassert(type[k] == SSX_UP || type[k] == SSX_DB);
alpar@9 226 break;
alpar@9 227 case LPX_NF:
alpar@9 228 stat[k] = SSX_NF;
alpar@9 229 xassert(type[k] == SSX_FR);
alpar@9 230 break;
alpar@9 231 case LPX_NS:
alpar@9 232 stat[k] = SSX_NS;
alpar@9 233 xassert(type[k] == SSX_FX);
alpar@9 234 break;
alpar@9 235 default:
alpar@9 236 xassert(stat != stat);
alpar@9 237 }
alpar@9 238 }
alpar@9 239 /* build permutation matix Q */
alpar@9 240 i = j = 0;
alpar@9 241 for (k = 1; k <= m+n; k++)
alpar@9 242 { if (stat[k] == SSX_BS)
alpar@9 243 { i++;
alpar@9 244 if (i > m) return 1;
alpar@9 245 Q_row[k] = i, Q_col[i] = k;
alpar@9 246 }
alpar@9 247 else
alpar@9 248 { j++;
alpar@9 249 if (j > n) return 1;
alpar@9 250 Q_row[k] = m+j, Q_col[m+j] = k;
alpar@9 251 }
alpar@9 252 }
alpar@9 253 xassert(i == m && j == n);
alpar@9 254 return 0;
alpar@9 255 }
alpar@9 256
alpar@9 257 int glp_exact(glp_prob *lp, const glp_smcp *parm)
alpar@9 258 { glp_smcp _parm;
alpar@9 259 SSX *ssx;
alpar@9 260 int m = lpx_get_num_rows(lp);
alpar@9 261 int n = lpx_get_num_cols(lp);
alpar@9 262 int nnz = lpx_get_num_nz(lp);
alpar@9 263 int i, j, k, type, pst, dst, ret, *stat;
alpar@9 264 double lb, ub, *prim, *dual, sum;
alpar@9 265 if (parm == NULL)
alpar@9 266 parm = &_parm, glp_init_smcp((glp_smcp *)parm);
alpar@9 267 /* check control parameters */
alpar@9 268 if (parm->it_lim < 0)
alpar@9 269 xerror("glp_exact: it_lim = %d; invalid parameter\n",
alpar@9 270 parm->it_lim);
alpar@9 271 if (parm->tm_lim < 0)
alpar@9 272 xerror("glp_exact: tm_lim = %d; invalid parameter\n",
alpar@9 273 parm->tm_lim);
alpar@9 274 /* the problem must have at least one row and one column */
alpar@9 275 if (!(m > 0 && n > 0))
alpar@9 276 { xprintf("glp_exact: problem has no rows/columns\n");
alpar@9 277 return GLP_EFAIL;
alpar@9 278 }
alpar@9 279 #if 1
alpar@9 280 /* basic solution is currently undefined */
alpar@9 281 lp->pbs_stat = lp->dbs_stat = GLP_UNDEF;
alpar@9 282 lp->obj_val = 0.0;
alpar@9 283 lp->some = 0;
alpar@9 284 #endif
alpar@9 285 /* check that all double-bounded variables have correct bounds */
alpar@9 286 for (k = 1; k <= m+n; k++)
alpar@9 287 { if (k <= m)
alpar@9 288 { type = lpx_get_row_type(lp, k);
alpar@9 289 lb = lpx_get_row_lb(lp, k);
alpar@9 290 ub = lpx_get_row_ub(lp, k);
alpar@9 291 }
alpar@9 292 else
alpar@9 293 { type = lpx_get_col_type(lp, k-m);
alpar@9 294 lb = lpx_get_col_lb(lp, k-m);
alpar@9 295 ub = lpx_get_col_ub(lp, k-m);
alpar@9 296 }
alpar@9 297 if (type == LPX_DB && lb >= ub)
alpar@9 298 { xprintf("glp_exact: %s %d has invalid bounds\n",
alpar@9 299 k <= m ? "row" : "column", k <= m ? k : k-m);
alpar@9 300 return GLP_EBOUND;
alpar@9 301 }
alpar@9 302 }
alpar@9 303 /* create the simplex solver workspace */
alpar@9 304 xprintf("glp_exact: %d rows, %d columns, %d non-zeros\n",
alpar@9 305 m, n, nnz);
alpar@9 306 #ifdef HAVE_GMP
alpar@9 307 xprintf("GNU MP bignum library is being used\n");
alpar@9 308 #else
alpar@9 309 xprintf("GLPK bignum module is being used\n");
alpar@9 310 xprintf("(Consider installing GNU MP to attain a much better perf"
alpar@9 311 "ormance.)\n");
alpar@9 312 #endif
alpar@9 313 ssx = ssx_create(m, n, nnz);
alpar@9 314 /* load LP problem data into the workspace */
alpar@9 315 load_data(ssx, lp);
alpar@9 316 /* load current LP basis into the workspace */
alpar@9 317 if (load_basis(ssx, lp))
alpar@9 318 { xprintf("glp_exact: initial LP basis is invalid\n");
alpar@9 319 ret = GLP_EBADB;
alpar@9 320 goto done;
alpar@9 321 }
alpar@9 322 /* inherit some control parameters from the LP object */
alpar@9 323 #if 0
alpar@9 324 ssx->it_lim = lpx_get_int_parm(lp, LPX_K_ITLIM);
alpar@9 325 ssx->it_cnt = lpx_get_int_parm(lp, LPX_K_ITCNT);
alpar@9 326 ssx->tm_lim = lpx_get_real_parm(lp, LPX_K_TMLIM);
alpar@9 327 #else
alpar@9 328 ssx->it_lim = parm->it_lim;
alpar@9 329 ssx->it_cnt = lp->it_cnt;
alpar@9 330 ssx->tm_lim = (double)parm->tm_lim / 1000.0;
alpar@9 331 #endif
alpar@9 332 ssx->out_frq = 5.0;
alpar@9 333 ssx->tm_beg = xtime();
alpar@9 334 ssx->tm_lag = xlset(0);
alpar@9 335 /* solve LP */
alpar@9 336 ret = ssx_driver(ssx);
alpar@9 337 /* copy back some statistics to the LP object */
alpar@9 338 #if 0
alpar@9 339 lpx_set_int_parm(lp, LPX_K_ITLIM, ssx->it_lim);
alpar@9 340 lpx_set_int_parm(lp, LPX_K_ITCNT, ssx->it_cnt);
alpar@9 341 lpx_set_real_parm(lp, LPX_K_TMLIM, ssx->tm_lim);
alpar@9 342 #else
alpar@9 343 lp->it_cnt = ssx->it_cnt;
alpar@9 344 #endif
alpar@9 345 /* analyze the return code */
alpar@9 346 switch (ret)
alpar@9 347 { case 0:
alpar@9 348 /* optimal solution found */
alpar@9 349 ret = 0;
alpar@9 350 pst = LPX_P_FEAS, dst = LPX_D_FEAS;
alpar@9 351 break;
alpar@9 352 case 1:
alpar@9 353 /* problem has no feasible solution */
alpar@9 354 ret = 0;
alpar@9 355 pst = LPX_P_NOFEAS, dst = LPX_D_INFEAS;
alpar@9 356 break;
alpar@9 357 case 2:
alpar@9 358 /* problem has unbounded solution */
alpar@9 359 ret = 0;
alpar@9 360 pst = LPX_P_FEAS, dst = LPX_D_NOFEAS;
alpar@9 361 #if 1
alpar@9 362 xassert(1 <= ssx->q && ssx->q <= n);
alpar@9 363 lp->some = ssx->Q_col[m + ssx->q];
alpar@9 364 xassert(1 <= lp->some && lp->some <= m+n);
alpar@9 365 #endif
alpar@9 366 break;
alpar@9 367 case 3:
alpar@9 368 /* iteration limit exceeded (phase I) */
alpar@9 369 ret = GLP_EITLIM;
alpar@9 370 pst = LPX_P_INFEAS, dst = LPX_D_INFEAS;
alpar@9 371 break;
alpar@9 372 case 4:
alpar@9 373 /* iteration limit exceeded (phase II) */
alpar@9 374 ret = GLP_EITLIM;
alpar@9 375 pst = LPX_P_FEAS, dst = LPX_D_INFEAS;
alpar@9 376 break;
alpar@9 377 case 5:
alpar@9 378 /* time limit exceeded (phase I) */
alpar@9 379 ret = GLP_ETMLIM;
alpar@9 380 pst = LPX_P_INFEAS, dst = LPX_D_INFEAS;
alpar@9 381 break;
alpar@9 382 case 6:
alpar@9 383 /* time limit exceeded (phase II) */
alpar@9 384 ret = GLP_ETMLIM;
alpar@9 385 pst = LPX_P_FEAS, dst = LPX_D_INFEAS;
alpar@9 386 break;
alpar@9 387 case 7:
alpar@9 388 /* initial basis matrix is singular */
alpar@9 389 ret = GLP_ESING;
alpar@9 390 goto done;
alpar@9 391 default:
alpar@9 392 xassert(ret != ret);
alpar@9 393 }
alpar@9 394 /* obtain final basic solution components */
alpar@9 395 stat = xcalloc(1+m+n, sizeof(int));
alpar@9 396 prim = xcalloc(1+m+n, sizeof(double));
alpar@9 397 dual = xcalloc(1+m+n, sizeof(double));
alpar@9 398 for (k = 1; k <= m+n; k++)
alpar@9 399 { if (ssx->stat[k] == SSX_BS)
alpar@9 400 { i = ssx->Q_row[k]; /* x[k] = xB[i] */
alpar@9 401 xassert(1 <= i && i <= m);
alpar@9 402 stat[k] = LPX_BS;
alpar@9 403 prim[k] = mpq_get_d(ssx->bbar[i]);
alpar@9 404 dual[k] = 0.0;
alpar@9 405 }
alpar@9 406 else
alpar@9 407 { j = ssx->Q_row[k] - m; /* x[k] = xN[j] */
alpar@9 408 xassert(1 <= j && j <= n);
alpar@9 409 switch (ssx->stat[k])
alpar@9 410 { case SSX_NF:
alpar@9 411 stat[k] = LPX_NF;
alpar@9 412 prim[k] = 0.0;
alpar@9 413 break;
alpar@9 414 case SSX_NL:
alpar@9 415 stat[k] = LPX_NL;
alpar@9 416 prim[k] = mpq_get_d(ssx->lb[k]);
alpar@9 417 break;
alpar@9 418 case SSX_NU:
alpar@9 419 stat[k] = LPX_NU;
alpar@9 420 prim[k] = mpq_get_d(ssx->ub[k]);
alpar@9 421 break;
alpar@9 422 case SSX_NS:
alpar@9 423 stat[k] = LPX_NS;
alpar@9 424 prim[k] = mpq_get_d(ssx->lb[k]);
alpar@9 425 break;
alpar@9 426 default:
alpar@9 427 xassert(ssx != ssx);
alpar@9 428 }
alpar@9 429 dual[k] = mpq_get_d(ssx->cbar[j]);
alpar@9 430 }
alpar@9 431 }
alpar@9 432 /* and store them into the LP object */
alpar@9 433 pst = pst - LPX_P_UNDEF + GLP_UNDEF;
alpar@9 434 dst = dst - LPX_D_UNDEF + GLP_UNDEF;
alpar@9 435 for (k = 1; k <= m+n; k++)
alpar@9 436 stat[k] = stat[k] - LPX_BS + GLP_BS;
alpar@9 437 sum = lpx_get_obj_coef(lp, 0);
alpar@9 438 for (j = 1; j <= n; j++)
alpar@9 439 sum += lpx_get_obj_coef(lp, j) * prim[m+j];
alpar@9 440 lpx_put_solution(lp, 1, &pst, &dst, &sum,
alpar@9 441 &stat[0], &prim[0], &dual[0], &stat[m], &prim[m], &dual[m]);
alpar@9 442 xfree(stat);
alpar@9 443 xfree(prim);
alpar@9 444 xfree(dual);
alpar@9 445 done: /* delete the simplex solver workspace */
alpar@9 446 ssx_delete(ssx);
alpar@9 447 /* return to the application program */
alpar@9 448 return ret;
alpar@9 449 }
alpar@9 450
alpar@9 451 /* eof */