|
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 |
|
72 static 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 |
|
120 int 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); |
|
212 done: /* 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 |
|
236 void 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 |
|
261 int 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 |
|
280 double 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 |
|
302 double 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 |
|
327 double 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 |
|
352 double 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 |
|
378 double 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 */ |