lemon-project-template-glpk

annotate deps/glpk/src/glpapi08.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 /* glpapi08.c (interior-point 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 "glpapi.h"
alpar@9 26 #include "glpipm.h"
alpar@9 27 #include "glpnpp.h"
alpar@9 28
alpar@9 29 /***********************************************************************
alpar@9 30 * NAME
alpar@9 31 *
alpar@9 32 * glp_interior - solve LP problem with the interior-point method
alpar@9 33 *
alpar@9 34 * SYNOPSIS
alpar@9 35 *
alpar@9 36 * int glp_interior(glp_prob *P, const glp_iptcp *parm);
alpar@9 37 *
alpar@9 38 * The routine glp_interior is a driver to the LP solver based on the
alpar@9 39 * interior-point method.
alpar@9 40 *
alpar@9 41 * The interior-point solver has a set of control parameters. Values of
alpar@9 42 * the control parameters can be passed in a structure glp_iptcp, which
alpar@9 43 * the parameter parm points to.
alpar@9 44 *
alpar@9 45 * Currently this routine implements an easy variant of the primal-dual
alpar@9 46 * interior-point method based on Mehrotra's technique.
alpar@9 47 *
alpar@9 48 * This routine transforms the original LP problem to an equivalent LP
alpar@9 49 * problem in the standard formulation (all constraints are equalities,
alpar@9 50 * all variables are non-negative), calls the routine ipm_main to solve
alpar@9 51 * the transformed problem, and then transforms an obtained solution to
alpar@9 52 * the solution of the original problem.
alpar@9 53 *
alpar@9 54 * RETURNS
alpar@9 55 *
alpar@9 56 * 0 The LP problem instance has been successfully solved. This code
alpar@9 57 * does not necessarily mean that the solver has found optimal
alpar@9 58 * solution. It only means that the solution process was successful.
alpar@9 59 *
alpar@9 60 * GLP_EFAIL
alpar@9 61 * The problem has no rows/columns.
alpar@9 62 *
alpar@9 63 * GLP_ENOCVG
alpar@9 64 * Very slow convergence or divergence.
alpar@9 65 *
alpar@9 66 * GLP_EITLIM
alpar@9 67 * Iteration limit exceeded.
alpar@9 68 *
alpar@9 69 * GLP_EINSTAB
alpar@9 70 * Numerical instability on solving Newtonian system. */
alpar@9 71
alpar@9 72 static void transform(NPP *npp)
alpar@9 73 { /* transform LP to the standard formulation */
alpar@9 74 NPPROW *row, *prev_row;
alpar@9 75 NPPCOL *col, *prev_col;
alpar@9 76 for (row = npp->r_tail; row != NULL; row = prev_row)
alpar@9 77 { prev_row = row->prev;
alpar@9 78 if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
alpar@9 79 npp_free_row(npp, row);
alpar@9 80 else if (row->lb == -DBL_MAX)
alpar@9 81 npp_leq_row(npp, row);
alpar@9 82 else if (row->ub == +DBL_MAX)
alpar@9 83 npp_geq_row(npp, row);
alpar@9 84 else if (row->lb != row->ub)
alpar@9 85 { if (fabs(row->lb) < fabs(row->ub))
alpar@9 86 npp_geq_row(npp, row);
alpar@9 87 else
alpar@9 88 npp_leq_row(npp, row);
alpar@9 89 }
alpar@9 90 }
alpar@9 91 for (col = npp->c_tail; col != NULL; col = prev_col)
alpar@9 92 { prev_col = col->prev;
alpar@9 93 if (col->lb == -DBL_MAX && col->ub == +DBL_MAX)
alpar@9 94 npp_free_col(npp, col);
alpar@9 95 else if (col->lb == -DBL_MAX)
alpar@9 96 npp_ubnd_col(npp, col);
alpar@9 97 else if (col->ub == +DBL_MAX)
alpar@9 98 { if (col->lb != 0.0)
alpar@9 99 npp_lbnd_col(npp, col);
alpar@9 100 }
alpar@9 101 else if (col->lb != col->ub)
alpar@9 102 { if (fabs(col->lb) < fabs(col->ub))
alpar@9 103 { if (col->lb != 0.0)
alpar@9 104 npp_lbnd_col(npp, col);
alpar@9 105 }
alpar@9 106 else
alpar@9 107 npp_ubnd_col(npp, col);
alpar@9 108 npp_dbnd_col(npp, col);
alpar@9 109 }
alpar@9 110 else
alpar@9 111 npp_fixed_col(npp, col);
alpar@9 112 }
alpar@9 113 for (row = npp->r_head; row != NULL; row = row->next)
alpar@9 114 xassert(row->lb == row->ub);
alpar@9 115 for (col = npp->c_head; col != NULL; col = col->next)
alpar@9 116 xassert(col->lb == 0.0 && col->ub == +DBL_MAX);
alpar@9 117 return;
alpar@9 118 }
alpar@9 119
alpar@9 120 int glp_interior(glp_prob *P, const glp_iptcp *parm)
alpar@9 121 { glp_iptcp _parm;
alpar@9 122 GLPROW *row;
alpar@9 123 GLPCOL *col;
alpar@9 124 NPP *npp = NULL;
alpar@9 125 glp_prob *prob = NULL;
alpar@9 126 int i, j, ret;
alpar@9 127 /* check control parameters */
alpar@9 128 if (parm == NULL)
alpar@9 129 glp_init_iptcp(&_parm), parm = &_parm;
alpar@9 130 if (!(parm->msg_lev == GLP_MSG_OFF ||
alpar@9 131 parm->msg_lev == GLP_MSG_ERR ||
alpar@9 132 parm->msg_lev == GLP_MSG_ON ||
alpar@9 133 parm->msg_lev == GLP_MSG_ALL))
alpar@9 134 xerror("glp_interior: msg_lev = %d; invalid parameter\n",
alpar@9 135 parm->msg_lev);
alpar@9 136 if (!(parm->ord_alg == GLP_ORD_NONE ||
alpar@9 137 parm->ord_alg == GLP_ORD_QMD ||
alpar@9 138 parm->ord_alg == GLP_ORD_AMD ||
alpar@9 139 parm->ord_alg == GLP_ORD_SYMAMD))
alpar@9 140 xerror("glp_interior: ord_alg = %d; invalid parameter\n",
alpar@9 141 parm->ord_alg);
alpar@9 142 /* interior-point solution is currently undefined */
alpar@9 143 P->ipt_stat = GLP_UNDEF;
alpar@9 144 P->ipt_obj = 0.0;
alpar@9 145 /* check bounds of double-bounded variables */
alpar@9 146 for (i = 1; i <= P->m; i++)
alpar@9 147 { row = P->row[i];
alpar@9 148 if (row->type == GLP_DB && row->lb >= row->ub)
alpar@9 149 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9 150 xprintf("glp_interior: row %d: lb = %g, ub = %g; incorre"
alpar@9 151 "ct bounds\n", i, row->lb, row->ub);
alpar@9 152 ret = GLP_EBOUND;
alpar@9 153 goto done;
alpar@9 154 }
alpar@9 155 }
alpar@9 156 for (j = 1; j <= P->n; j++)
alpar@9 157 { col = P->col[j];
alpar@9 158 if (col->type == GLP_DB && col->lb >= col->ub)
alpar@9 159 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9 160 xprintf("glp_interior: column %d: lb = %g, ub = %g; inco"
alpar@9 161 "rrect bounds\n", j, col->lb, col->ub);
alpar@9 162 ret = GLP_EBOUND;
alpar@9 163 goto done;
alpar@9 164 }
alpar@9 165 }
alpar@9 166 /* transform LP to the standard formulation */
alpar@9 167 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9 168 xprintf("Original LP has %d row(s), %d column(s), and %d non-z"
alpar@9 169 "ero(s)\n", P->m, P->n, P->nnz);
alpar@9 170 npp = npp_create_wksp();
alpar@9 171 npp_load_prob(npp, P, GLP_OFF, GLP_IPT, GLP_ON);
alpar@9 172 transform(npp);
alpar@9 173 prob = glp_create_prob();
alpar@9 174 npp_build_prob(npp, prob);
alpar@9 175 if (parm->msg_lev >= GLP_MSG_ALL)
alpar@9 176 xprintf("Working LP has %d row(s), %d column(s), and %d non-ze"
alpar@9 177 "ro(s)\n", prob->m, prob->n, prob->nnz);
alpar@9 178 #if 1
alpar@9 179 /* currently empty problem cannot be solved */
alpar@9 180 if (!(prob->m > 0 && prob->n > 0))
alpar@9 181 { if (parm->msg_lev >= GLP_MSG_ERR)
alpar@9 182 xprintf("glp_interior: unable to solve empty problem\n");
alpar@9 183 ret = GLP_EFAIL;
alpar@9 184 goto done;
alpar@9 185 }
alpar@9 186 #endif
alpar@9 187 /* scale the resultant LP */
alpar@9 188 { ENV *env = get_env_ptr();
alpar@9 189 int term_out = env->term_out;
alpar@9 190 env->term_out = GLP_OFF;
alpar@9 191 glp_scale_prob(prob, GLP_SF_EQ);
alpar@9 192 env->term_out = term_out;
alpar@9 193 }
alpar@9 194 /* warn about dense columns */
alpar@9 195 if (parm->msg_lev >= GLP_MSG_ON && prob->m >= 200)
alpar@9 196 { int len, cnt = 0;
alpar@9 197 for (j = 1; j <= prob->n; j++)
alpar@9 198 { len = glp_get_mat_col(prob, j, NULL, NULL);
alpar@9 199 if ((double)len >= 0.20 * (double)prob->m) cnt++;
alpar@9 200 }
alpar@9 201 if (cnt == 1)
alpar@9 202 xprintf("WARNING: PROBLEM HAS ONE DENSE COLUMN\n");
alpar@9 203 else if (cnt > 0)
alpar@9 204 xprintf("WARNING: PROBLEM HAS %d DENSE COLUMNS\n", cnt);
alpar@9 205 }
alpar@9 206 /* solve the transformed LP */
alpar@9 207 ret = ipm_solve(prob, parm);
alpar@9 208 /* postprocess solution from the transformed LP */
alpar@9 209 npp_postprocess(npp, prob);
alpar@9 210 /* and store solution to the original LP */
alpar@9 211 npp_unload_sol(npp, P);
alpar@9 212 done: /* free working program objects */
alpar@9 213 if (npp != NULL) npp_delete_wksp(npp);
alpar@9 214 if (prob != NULL) glp_delete_prob(prob);
alpar@9 215 /* return to the application program */
alpar@9 216 return ret;
alpar@9 217 }
alpar@9 218
alpar@9 219 /***********************************************************************
alpar@9 220 * NAME
alpar@9 221 *
alpar@9 222 * glp_init_iptcp - initialize interior-point solver control parameters
alpar@9 223 *
alpar@9 224 * SYNOPSIS
alpar@9 225 *
alpar@9 226 * void glp_init_iptcp(glp_iptcp *parm);
alpar@9 227 *
alpar@9 228 * DESCRIPTION
alpar@9 229 *
alpar@9 230 * The routine glp_init_iptcp initializes control parameters, which are
alpar@9 231 * used by the interior-point solver, with default values.
alpar@9 232 *
alpar@9 233 * Default values of the control parameters are stored in the glp_iptcp
alpar@9 234 * structure, which the parameter parm points to. */
alpar@9 235
alpar@9 236 void glp_init_iptcp(glp_iptcp *parm)
alpar@9 237 { parm->msg_lev = GLP_MSG_ALL;
alpar@9 238 parm->ord_alg = GLP_ORD_AMD;
alpar@9 239 return;
alpar@9 240 }
alpar@9 241
alpar@9 242 /***********************************************************************
alpar@9 243 * NAME
alpar@9 244 *
alpar@9 245 * glp_ipt_status - retrieve status of interior-point solution
alpar@9 246 *
alpar@9 247 * SYNOPSIS
alpar@9 248 *
alpar@9 249 * int glp_ipt_status(glp_prob *lp);
alpar@9 250 *
alpar@9 251 * RETURNS
alpar@9 252 *
alpar@9 253 * The routine glp_ipt_status reports the status of solution found by
alpar@9 254 * the interior-point solver as follows:
alpar@9 255 *
alpar@9 256 * GLP_UNDEF - interior-point solution is undefined;
alpar@9 257 * GLP_OPT - interior-point solution is optimal;
alpar@9 258 * GLP_INFEAS - interior-point solution is infeasible;
alpar@9 259 * GLP_NOFEAS - no feasible solution exists. */
alpar@9 260
alpar@9 261 int glp_ipt_status(glp_prob *lp)
alpar@9 262 { int ipt_stat = lp->ipt_stat;
alpar@9 263 return ipt_stat;
alpar@9 264 }
alpar@9 265
alpar@9 266 /***********************************************************************
alpar@9 267 * NAME
alpar@9 268 *
alpar@9 269 * glp_ipt_obj_val - retrieve objective value (interior point)
alpar@9 270 *
alpar@9 271 * SYNOPSIS
alpar@9 272 *
alpar@9 273 * double glp_ipt_obj_val(glp_prob *lp);
alpar@9 274 *
alpar@9 275 * RETURNS
alpar@9 276 *
alpar@9 277 * The routine glp_ipt_obj_val returns value of the objective function
alpar@9 278 * for interior-point solution. */
alpar@9 279
alpar@9 280 double glp_ipt_obj_val(glp_prob *lp)
alpar@9 281 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 282 double z;
alpar@9 283 z = lp->ipt_obj;
alpar@9 284 /*if (cps->round && fabs(z) < 1e-9) z = 0.0;*/
alpar@9 285 return z;
alpar@9 286 }
alpar@9 287
alpar@9 288 /***********************************************************************
alpar@9 289 * NAME
alpar@9 290 *
alpar@9 291 * glp_ipt_row_prim - retrieve row primal value (interior point)
alpar@9 292 *
alpar@9 293 * SYNOPSIS
alpar@9 294 *
alpar@9 295 * double glp_ipt_row_prim(glp_prob *lp, int i);
alpar@9 296 *
alpar@9 297 * RETURNS
alpar@9 298 *
alpar@9 299 * The routine glp_ipt_row_prim returns primal value of the auxiliary
alpar@9 300 * variable associated with i-th row. */
alpar@9 301
alpar@9 302 double glp_ipt_row_prim(glp_prob *lp, int i)
alpar@9 303 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 304 double pval;
alpar@9 305 if (!(1 <= i && i <= lp->m))
alpar@9 306 xerror("glp_ipt_row_prim: i = %d; row number out of range\n",
alpar@9 307 i);
alpar@9 308 pval = lp->row[i]->pval;
alpar@9 309 /*if (cps->round && fabs(pval) < 1e-9) pval = 0.0;*/
alpar@9 310 return pval;
alpar@9 311 }
alpar@9 312
alpar@9 313 /***********************************************************************
alpar@9 314 * NAME
alpar@9 315 *
alpar@9 316 * glp_ipt_row_dual - retrieve row dual value (interior point)
alpar@9 317 *
alpar@9 318 * SYNOPSIS
alpar@9 319 *
alpar@9 320 * double glp_ipt_row_dual(glp_prob *lp, int i);
alpar@9 321 *
alpar@9 322 * RETURNS
alpar@9 323 *
alpar@9 324 * The routine glp_ipt_row_dual returns dual value (i.e. reduced cost)
alpar@9 325 * of the auxiliary variable associated with i-th row. */
alpar@9 326
alpar@9 327 double glp_ipt_row_dual(glp_prob *lp, int i)
alpar@9 328 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 329 double dval;
alpar@9 330 if (!(1 <= i && i <= lp->m))
alpar@9 331 xerror("glp_ipt_row_dual: i = %d; row number out of range\n",
alpar@9 332 i);
alpar@9 333 dval = lp->row[i]->dval;
alpar@9 334 /*if (cps->round && fabs(dval) < 1e-9) dval = 0.0;*/
alpar@9 335 return dval;
alpar@9 336 }
alpar@9 337
alpar@9 338 /***********************************************************************
alpar@9 339 * NAME
alpar@9 340 *
alpar@9 341 * glp_ipt_col_prim - retrieve column primal value (interior point)
alpar@9 342 *
alpar@9 343 * SYNOPSIS
alpar@9 344 *
alpar@9 345 * double glp_ipt_col_prim(glp_prob *lp, int j);
alpar@9 346 *
alpar@9 347 * RETURNS
alpar@9 348 *
alpar@9 349 * The routine glp_ipt_col_prim returns primal value of the structural
alpar@9 350 * variable associated with j-th column. */
alpar@9 351
alpar@9 352 double glp_ipt_col_prim(glp_prob *lp, int j)
alpar@9 353 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 354 double pval;
alpar@9 355 if (!(1 <= j && j <= lp->n))
alpar@9 356 xerror("glp_ipt_col_prim: j = %d; column number out of range\n"
alpar@9 357 , j);
alpar@9 358 pval = lp->col[j]->pval;
alpar@9 359 /*if (cps->round && fabs(pval) < 1e-9) pval = 0.0;*/
alpar@9 360 return pval;
alpar@9 361 }
alpar@9 362
alpar@9 363 /***********************************************************************
alpar@9 364 * NAME
alpar@9 365 *
alpar@9 366 * glp_ipt_col_dual - retrieve column dual value (interior point)
alpar@9 367 *
alpar@9 368 * SYNOPSIS
alpar@9 369 *
alpar@9 370 * #include "glplpx.h"
alpar@9 371 * double glp_ipt_col_dual(glp_prob *lp, int j);
alpar@9 372 *
alpar@9 373 * RETURNS
alpar@9 374 *
alpar@9 375 * The routine glp_ipt_col_dual returns dual value (i.e. reduced cost)
alpar@9 376 * of the structural variable associated with j-th column. */
alpar@9 377
alpar@9 378 double glp_ipt_col_dual(glp_prob *lp, int j)
alpar@9 379 { /*struct LPXCPS *cps = lp->cps;*/
alpar@9 380 double dval;
alpar@9 381 if (!(1 <= j && j <= lp->n))
alpar@9 382 xerror("glp_ipt_col_dual: j = %d; column number out of range\n"
alpar@9 383 , j);
alpar@9 384 dval = lp->col[j]->dval;
alpar@9 385 /*if (cps->round && fabs(dval) < 1e-9) dval = 0.0;*/
alpar@9 386 return dval;
alpar@9 387 }
alpar@9 388
alpar@9 389 /* eof */