lemon-project-template-glpk
comparison deps/glpk/src/glpapi08.c @ 9:33de93886c88
Import GLPK 4.47
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 20:59:10 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:4da6898eee40 |
---|---|
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, 2011 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 */ |