1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpapi08.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,389 @@
1.4 +/* glpapi08.c (interior-point method routines) */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpapi.h"
1.29 +#include "glpipm.h"
1.30 +#include "glpnpp.h"
1.31 +
1.32 +/***********************************************************************
1.33 +* NAME
1.34 +*
1.35 +* glp_interior - solve LP problem with the interior-point method
1.36 +*
1.37 +* SYNOPSIS
1.38 +*
1.39 +* int glp_interior(glp_prob *P, const glp_iptcp *parm);
1.40 +*
1.41 +* The routine glp_interior is a driver to the LP solver based on the
1.42 +* interior-point method.
1.43 +*
1.44 +* The interior-point solver has a set of control parameters. Values of
1.45 +* the control parameters can be passed in a structure glp_iptcp, which
1.46 +* the parameter parm points to.
1.47 +*
1.48 +* Currently this routine implements an easy variant of the primal-dual
1.49 +* interior-point method based on Mehrotra's technique.
1.50 +*
1.51 +* This routine transforms the original LP problem to an equivalent LP
1.52 +* problem in the standard formulation (all constraints are equalities,
1.53 +* all variables are non-negative), calls the routine ipm_main to solve
1.54 +* the transformed problem, and then transforms an obtained solution to
1.55 +* the solution of the original problem.
1.56 +*
1.57 +* RETURNS
1.58 +*
1.59 +* 0 The LP problem instance has been successfully solved. This code
1.60 +* does not necessarily mean that the solver has found optimal
1.61 +* solution. It only means that the solution process was successful.
1.62 +*
1.63 +* GLP_EFAIL
1.64 +* The problem has no rows/columns.
1.65 +*
1.66 +* GLP_ENOCVG
1.67 +* Very slow convergence or divergence.
1.68 +*
1.69 +* GLP_EITLIM
1.70 +* Iteration limit exceeded.
1.71 +*
1.72 +* GLP_EINSTAB
1.73 +* Numerical instability on solving Newtonian system. */
1.74 +
1.75 +static void transform(NPP *npp)
1.76 +{ /* transform LP to the standard formulation */
1.77 + NPPROW *row, *prev_row;
1.78 + NPPCOL *col, *prev_col;
1.79 + for (row = npp->r_tail; row != NULL; row = prev_row)
1.80 + { prev_row = row->prev;
1.81 + if (row->lb == -DBL_MAX && row->ub == +DBL_MAX)
1.82 + npp_free_row(npp, row);
1.83 + else if (row->lb == -DBL_MAX)
1.84 + npp_leq_row(npp, row);
1.85 + else if (row->ub == +DBL_MAX)
1.86 + npp_geq_row(npp, row);
1.87 + else if (row->lb != row->ub)
1.88 + { if (fabs(row->lb) < fabs(row->ub))
1.89 + npp_geq_row(npp, row);
1.90 + else
1.91 + npp_leq_row(npp, row);
1.92 + }
1.93 + }
1.94 + for (col = npp->c_tail; col != NULL; col = prev_col)
1.95 + { prev_col = col->prev;
1.96 + if (col->lb == -DBL_MAX && col->ub == +DBL_MAX)
1.97 + npp_free_col(npp, col);
1.98 + else if (col->lb == -DBL_MAX)
1.99 + npp_ubnd_col(npp, col);
1.100 + else if (col->ub == +DBL_MAX)
1.101 + { if (col->lb != 0.0)
1.102 + npp_lbnd_col(npp, col);
1.103 + }
1.104 + else if (col->lb != col->ub)
1.105 + { if (fabs(col->lb) < fabs(col->ub))
1.106 + { if (col->lb != 0.0)
1.107 + npp_lbnd_col(npp, col);
1.108 + }
1.109 + else
1.110 + npp_ubnd_col(npp, col);
1.111 + npp_dbnd_col(npp, col);
1.112 + }
1.113 + else
1.114 + npp_fixed_col(npp, col);
1.115 + }
1.116 + for (row = npp->r_head; row != NULL; row = row->next)
1.117 + xassert(row->lb == row->ub);
1.118 + for (col = npp->c_head; col != NULL; col = col->next)
1.119 + xassert(col->lb == 0.0 && col->ub == +DBL_MAX);
1.120 + return;
1.121 +}
1.122 +
1.123 +int glp_interior(glp_prob *P, const glp_iptcp *parm)
1.124 +{ glp_iptcp _parm;
1.125 + GLPROW *row;
1.126 + GLPCOL *col;
1.127 + NPP *npp = NULL;
1.128 + glp_prob *prob = NULL;
1.129 + int i, j, ret;
1.130 + /* check control parameters */
1.131 + if (parm == NULL)
1.132 + glp_init_iptcp(&_parm), parm = &_parm;
1.133 + if (!(parm->msg_lev == GLP_MSG_OFF ||
1.134 + parm->msg_lev == GLP_MSG_ERR ||
1.135 + parm->msg_lev == GLP_MSG_ON ||
1.136 + parm->msg_lev == GLP_MSG_ALL))
1.137 + xerror("glp_interior: msg_lev = %d; invalid parameter\n",
1.138 + parm->msg_lev);
1.139 + if (!(parm->ord_alg == GLP_ORD_NONE ||
1.140 + parm->ord_alg == GLP_ORD_QMD ||
1.141 + parm->ord_alg == GLP_ORD_AMD ||
1.142 + parm->ord_alg == GLP_ORD_SYMAMD))
1.143 + xerror("glp_interior: ord_alg = %d; invalid parameter\n",
1.144 + parm->ord_alg);
1.145 + /* interior-point solution is currently undefined */
1.146 + P->ipt_stat = GLP_UNDEF;
1.147 + P->ipt_obj = 0.0;
1.148 + /* check bounds of double-bounded variables */
1.149 + for (i = 1; i <= P->m; i++)
1.150 + { row = P->row[i];
1.151 + if (row->type == GLP_DB && row->lb >= row->ub)
1.152 + { if (parm->msg_lev >= GLP_MSG_ERR)
1.153 + xprintf("glp_interior: row %d: lb = %g, ub = %g; incorre"
1.154 + "ct bounds\n", i, row->lb, row->ub);
1.155 + ret = GLP_EBOUND;
1.156 + goto done;
1.157 + }
1.158 + }
1.159 + for (j = 1; j <= P->n; j++)
1.160 + { col = P->col[j];
1.161 + if (col->type == GLP_DB && col->lb >= col->ub)
1.162 + { if (parm->msg_lev >= GLP_MSG_ERR)
1.163 + xprintf("glp_interior: column %d: lb = %g, ub = %g; inco"
1.164 + "rrect bounds\n", j, col->lb, col->ub);
1.165 + ret = GLP_EBOUND;
1.166 + goto done;
1.167 + }
1.168 + }
1.169 + /* transform LP to the standard formulation */
1.170 + if (parm->msg_lev >= GLP_MSG_ALL)
1.171 + xprintf("Original LP has %d row(s), %d column(s), and %d non-z"
1.172 + "ero(s)\n", P->m, P->n, P->nnz);
1.173 + npp = npp_create_wksp();
1.174 + npp_load_prob(npp, P, GLP_OFF, GLP_IPT, GLP_ON);
1.175 + transform(npp);
1.176 + prob = glp_create_prob();
1.177 + npp_build_prob(npp, prob);
1.178 + if (parm->msg_lev >= GLP_MSG_ALL)
1.179 + xprintf("Working LP has %d row(s), %d column(s), and %d non-ze"
1.180 + "ro(s)\n", prob->m, prob->n, prob->nnz);
1.181 +#if 1
1.182 + /* currently empty problem cannot be solved */
1.183 + if (!(prob->m > 0 && prob->n > 0))
1.184 + { if (parm->msg_lev >= GLP_MSG_ERR)
1.185 + xprintf("glp_interior: unable to solve empty problem\n");
1.186 + ret = GLP_EFAIL;
1.187 + goto done;
1.188 + }
1.189 +#endif
1.190 + /* scale the resultant LP */
1.191 + { ENV *env = get_env_ptr();
1.192 + int term_out = env->term_out;
1.193 + env->term_out = GLP_OFF;
1.194 + glp_scale_prob(prob, GLP_SF_EQ);
1.195 + env->term_out = term_out;
1.196 + }
1.197 + /* warn about dense columns */
1.198 + if (parm->msg_lev >= GLP_MSG_ON && prob->m >= 200)
1.199 + { int len, cnt = 0;
1.200 + for (j = 1; j <= prob->n; j++)
1.201 + { len = glp_get_mat_col(prob, j, NULL, NULL);
1.202 + if ((double)len >= 0.20 * (double)prob->m) cnt++;
1.203 + }
1.204 + if (cnt == 1)
1.205 + xprintf("WARNING: PROBLEM HAS ONE DENSE COLUMN\n");
1.206 + else if (cnt > 0)
1.207 + xprintf("WARNING: PROBLEM HAS %d DENSE COLUMNS\n", cnt);
1.208 + }
1.209 + /* solve the transformed LP */
1.210 + ret = ipm_solve(prob, parm);
1.211 + /* postprocess solution from the transformed LP */
1.212 + npp_postprocess(npp, prob);
1.213 + /* and store solution to the original LP */
1.214 + npp_unload_sol(npp, P);
1.215 +done: /* free working program objects */
1.216 + if (npp != NULL) npp_delete_wksp(npp);
1.217 + if (prob != NULL) glp_delete_prob(prob);
1.218 + /* return to the application program */
1.219 + return ret;
1.220 +}
1.221 +
1.222 +/***********************************************************************
1.223 +* NAME
1.224 +*
1.225 +* glp_init_iptcp - initialize interior-point solver control parameters
1.226 +*
1.227 +* SYNOPSIS
1.228 +*
1.229 +* void glp_init_iptcp(glp_iptcp *parm);
1.230 +*
1.231 +* DESCRIPTION
1.232 +*
1.233 +* The routine glp_init_iptcp initializes control parameters, which are
1.234 +* used by the interior-point solver, with default values.
1.235 +*
1.236 +* Default values of the control parameters are stored in the glp_iptcp
1.237 +* structure, which the parameter parm points to. */
1.238 +
1.239 +void glp_init_iptcp(glp_iptcp *parm)
1.240 +{ parm->msg_lev = GLP_MSG_ALL;
1.241 + parm->ord_alg = GLP_ORD_AMD;
1.242 + return;
1.243 +}
1.244 +
1.245 +/***********************************************************************
1.246 +* NAME
1.247 +*
1.248 +* glp_ipt_status - retrieve status of interior-point solution
1.249 +*
1.250 +* SYNOPSIS
1.251 +*
1.252 +* int glp_ipt_status(glp_prob *lp);
1.253 +*
1.254 +* RETURNS
1.255 +*
1.256 +* The routine glp_ipt_status reports the status of solution found by
1.257 +* the interior-point solver as follows:
1.258 +*
1.259 +* GLP_UNDEF - interior-point solution is undefined;
1.260 +* GLP_OPT - interior-point solution is optimal;
1.261 +* GLP_INFEAS - interior-point solution is infeasible;
1.262 +* GLP_NOFEAS - no feasible solution exists. */
1.263 +
1.264 +int glp_ipt_status(glp_prob *lp)
1.265 +{ int ipt_stat = lp->ipt_stat;
1.266 + return ipt_stat;
1.267 +}
1.268 +
1.269 +/***********************************************************************
1.270 +* NAME
1.271 +*
1.272 +* glp_ipt_obj_val - retrieve objective value (interior point)
1.273 +*
1.274 +* SYNOPSIS
1.275 +*
1.276 +* double glp_ipt_obj_val(glp_prob *lp);
1.277 +*
1.278 +* RETURNS
1.279 +*
1.280 +* The routine glp_ipt_obj_val returns value of the objective function
1.281 +* for interior-point solution. */
1.282 +
1.283 +double glp_ipt_obj_val(glp_prob *lp)
1.284 +{ /*struct LPXCPS *cps = lp->cps;*/
1.285 + double z;
1.286 + z = lp->ipt_obj;
1.287 + /*if (cps->round && fabs(z) < 1e-9) z = 0.0;*/
1.288 + return z;
1.289 +}
1.290 +
1.291 +/***********************************************************************
1.292 +* NAME
1.293 +*
1.294 +* glp_ipt_row_prim - retrieve row primal value (interior point)
1.295 +*
1.296 +* SYNOPSIS
1.297 +*
1.298 +* double glp_ipt_row_prim(glp_prob *lp, int i);
1.299 +*
1.300 +* RETURNS
1.301 +*
1.302 +* The routine glp_ipt_row_prim returns primal value of the auxiliary
1.303 +* variable associated with i-th row. */
1.304 +
1.305 +double glp_ipt_row_prim(glp_prob *lp, int i)
1.306 +{ /*struct LPXCPS *cps = lp->cps;*/
1.307 + double pval;
1.308 + if (!(1 <= i && i <= lp->m))
1.309 + xerror("glp_ipt_row_prim: i = %d; row number out of range\n",
1.310 + i);
1.311 + pval = lp->row[i]->pval;
1.312 + /*if (cps->round && fabs(pval) < 1e-9) pval = 0.0;*/
1.313 + return pval;
1.314 +}
1.315 +
1.316 +/***********************************************************************
1.317 +* NAME
1.318 +*
1.319 +* glp_ipt_row_dual - retrieve row dual value (interior point)
1.320 +*
1.321 +* SYNOPSIS
1.322 +*
1.323 +* double glp_ipt_row_dual(glp_prob *lp, int i);
1.324 +*
1.325 +* RETURNS
1.326 +*
1.327 +* The routine glp_ipt_row_dual returns dual value (i.e. reduced cost)
1.328 +* of the auxiliary variable associated with i-th row. */
1.329 +
1.330 +double glp_ipt_row_dual(glp_prob *lp, int i)
1.331 +{ /*struct LPXCPS *cps = lp->cps;*/
1.332 + double dval;
1.333 + if (!(1 <= i && i <= lp->m))
1.334 + xerror("glp_ipt_row_dual: i = %d; row number out of range\n",
1.335 + i);
1.336 + dval = lp->row[i]->dval;
1.337 + /*if (cps->round && fabs(dval) < 1e-9) dval = 0.0;*/
1.338 + return dval;
1.339 +}
1.340 +
1.341 +/***********************************************************************
1.342 +* NAME
1.343 +*
1.344 +* glp_ipt_col_prim - retrieve column primal value (interior point)
1.345 +*
1.346 +* SYNOPSIS
1.347 +*
1.348 +* double glp_ipt_col_prim(glp_prob *lp, int j);
1.349 +*
1.350 +* RETURNS
1.351 +*
1.352 +* The routine glp_ipt_col_prim returns primal value of the structural
1.353 +* variable associated with j-th column. */
1.354 +
1.355 +double glp_ipt_col_prim(glp_prob *lp, int j)
1.356 +{ /*struct LPXCPS *cps = lp->cps;*/
1.357 + double pval;
1.358 + if (!(1 <= j && j <= lp->n))
1.359 + xerror("glp_ipt_col_prim: j = %d; column number out of range\n"
1.360 + , j);
1.361 + pval = lp->col[j]->pval;
1.362 + /*if (cps->round && fabs(pval) < 1e-9) pval = 0.0;*/
1.363 + return pval;
1.364 +}
1.365 +
1.366 +/***********************************************************************
1.367 +* NAME
1.368 +*
1.369 +* glp_ipt_col_dual - retrieve column dual value (interior point)
1.370 +*
1.371 +* SYNOPSIS
1.372 +*
1.373 +* #include "glplpx.h"
1.374 +* double glp_ipt_col_dual(glp_prob *lp, int j);
1.375 +*
1.376 +* RETURNS
1.377 +*
1.378 +* The routine glp_ipt_col_dual returns dual value (i.e. reduced cost)
1.379 +* of the structural variable associated with j-th column. */
1.380 +
1.381 +double glp_ipt_col_dual(glp_prob *lp, int j)
1.382 +{ /*struct LPXCPS *cps = lp->cps;*/
1.383 + double dval;
1.384 + if (!(1 <= j && j <= lp->n))
1.385 + xerror("glp_ipt_col_dual: j = %d; column number out of range\n"
1.386 + , j);
1.387 + dval = lp->col[j]->dval;
1.388 + /*if (cps->round && fabs(dval) < 1e-9) dval = 0.0;*/
1.389 + return dval;
1.390 +}
1.391 +
1.392 +/* eof */