COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpapi08.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

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