1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpapi08.c	Sun Nov 06 20:59:10 2011 +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, 2011 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 */