1 | /* glpnpp02.c */ |
---|
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 "glpnpp.h" |
---|
26 | |
---|
27 | /*********************************************************************** |
---|
28 | * NAME |
---|
29 | * |
---|
30 | * npp_free_row - process free (unbounded) row |
---|
31 | * |
---|
32 | * SYNOPSIS |
---|
33 | * |
---|
34 | * #include "glpnpp.h" |
---|
35 | * void npp_free_row(NPP *npp, NPPROW *p); |
---|
36 | * |
---|
37 | * DESCRIPTION |
---|
38 | * |
---|
39 | * The routine npp_free_row processes row p, which is free (i.e. has |
---|
40 | * no finite bounds): |
---|
41 | * |
---|
42 | * -inf < sum a[p,j] x[j] < +inf. (1) |
---|
43 | * j |
---|
44 | * |
---|
45 | * PROBLEM TRANSFORMATION |
---|
46 | * |
---|
47 | * Constraint (1) cannot be active, so it is redundant and can be |
---|
48 | * removed from the original problem. |
---|
49 | * |
---|
50 | * Removing row p leads to removing a column of multiplier pi[p] for |
---|
51 | * this row in the dual system. Since row p has no bounds, pi[p] = 0, |
---|
52 | * so removing the column does not affect the dual solution. |
---|
53 | * |
---|
54 | * RECOVERING BASIC SOLUTION |
---|
55 | * |
---|
56 | * In solution to the original problem row p is inactive constraint, |
---|
57 | * so it is assigned status GLP_BS, and multiplier pi[p] is assigned |
---|
58 | * zero value. |
---|
59 | * |
---|
60 | * RECOVERING INTERIOR-POINT SOLUTION |
---|
61 | * |
---|
62 | * In solution to the original problem row p is inactive constraint, |
---|
63 | * so its multiplier pi[p] is assigned zero value. |
---|
64 | * |
---|
65 | * RECOVERING MIP SOLUTION |
---|
66 | * |
---|
67 | * None needed. */ |
---|
68 | |
---|
69 | struct free_row |
---|
70 | { /* free (unbounded) row */ |
---|
71 | int p; |
---|
72 | /* row reference number */ |
---|
73 | }; |
---|
74 | |
---|
75 | static int rcv_free_row(NPP *npp, void *info); |
---|
76 | |
---|
77 | void npp_free_row(NPP *npp, NPPROW *p) |
---|
78 | { /* process free (unbounded) row */ |
---|
79 | struct free_row *info; |
---|
80 | /* the row must be free */ |
---|
81 | xassert(p->lb == -DBL_MAX && p->ub == +DBL_MAX); |
---|
82 | /* create transformation stack entry */ |
---|
83 | info = npp_push_tse(npp, |
---|
84 | rcv_free_row, sizeof(struct free_row)); |
---|
85 | info->p = p->i; |
---|
86 | /* remove the row from the problem */ |
---|
87 | npp_del_row(npp, p); |
---|
88 | return; |
---|
89 | } |
---|
90 | |
---|
91 | static int rcv_free_row(NPP *npp, void *_info) |
---|
92 | { /* recover free (unbounded) row */ |
---|
93 | struct free_row *info = _info; |
---|
94 | if (npp->sol == GLP_SOL) |
---|
95 | npp->r_stat[info->p] = GLP_BS; |
---|
96 | if (npp->sol != GLP_MIP) |
---|
97 | npp->r_pi[info->p] = 0.0; |
---|
98 | return 0; |
---|
99 | } |
---|
100 | |
---|
101 | /*********************************************************************** |
---|
102 | * NAME |
---|
103 | * |
---|
104 | * npp_geq_row - process row of 'not less than' type |
---|
105 | * |
---|
106 | * SYNOPSIS |
---|
107 | * |
---|
108 | * #include "glpnpp.h" |
---|
109 | * void npp_geq_row(NPP *npp, NPPROW *p); |
---|
110 | * |
---|
111 | * DESCRIPTION |
---|
112 | * |
---|
113 | * The routine npp_geq_row processes row p, which is 'not less than' |
---|
114 | * inequality constraint: |
---|
115 | * |
---|
116 | * L[p] <= sum a[p,j] x[j] (<= U[p]), (1) |
---|
117 | * j |
---|
118 | * |
---|
119 | * where L[p] < U[p], and upper bound may not exist (U[p] = +oo). |
---|
120 | * |
---|
121 | * PROBLEM TRANSFORMATION |
---|
122 | * |
---|
123 | * Constraint (1) can be replaced by equality constraint: |
---|
124 | * |
---|
125 | * sum a[p,j] x[j] - s = L[p], (2) |
---|
126 | * j |
---|
127 | * |
---|
128 | * where |
---|
129 | * |
---|
130 | * 0 <= s (<= U[p] - L[p]) (3) |
---|
131 | * |
---|
132 | * is a non-negative surplus variable. |
---|
133 | * |
---|
134 | * Since in the primal system there appears column s having the only |
---|
135 | * non-zero coefficient in row p, in the dual system there appears a |
---|
136 | * new row: |
---|
137 | * |
---|
138 | * (-1) pi[p] + lambda = 0, (4) |
---|
139 | * |
---|
140 | * where (-1) is coefficient of column s in row p, pi[p] is multiplier |
---|
141 | * of row p, lambda is multiplier of column q, 0 is coefficient of |
---|
142 | * column s in the objective row. |
---|
143 | * |
---|
144 | * RECOVERING BASIC SOLUTION |
---|
145 | * |
---|
146 | * Status of row p in solution to the original problem is determined |
---|
147 | * by its status and status of column q in solution to the transformed |
---|
148 | * problem as follows: |
---|
149 | * |
---|
150 | * +--------------------------------------+------------------+ |
---|
151 | * | Transformed problem | Original problem | |
---|
152 | * +-----------------+--------------------+------------------+ |
---|
153 | * | Status of row p | Status of column s | Status of row p | |
---|
154 | * +-----------------+--------------------+------------------+ |
---|
155 | * | GLP_BS | GLP_BS | N/A | |
---|
156 | * | GLP_BS | GLP_NL | GLP_BS | |
---|
157 | * | GLP_BS | GLP_NU | GLP_BS | |
---|
158 | * | GLP_NS | GLP_BS | GLP_BS | |
---|
159 | * | GLP_NS | GLP_NL | GLP_NL | |
---|
160 | * | GLP_NS | GLP_NU | GLP_NU | |
---|
161 | * +-----------------+--------------------+------------------+ |
---|
162 | * |
---|
163 | * Value of row multiplier pi[p] in solution to the original problem |
---|
164 | * is the same as in solution to the transformed problem. |
---|
165 | * |
---|
166 | * 1. In solution to the transformed problem row p and column q cannot |
---|
167 | * be basic at the same time; otherwise the basis matrix would have |
---|
168 | * two linear dependent columns: unity column of auxiliary variable |
---|
169 | * of row p and unity column of variable s. |
---|
170 | * |
---|
171 | * 2. Though in the transformed problem row p is equality constraint, |
---|
172 | * it may be basic due to primal degenerate solution. |
---|
173 | * |
---|
174 | * RECOVERING INTERIOR-POINT SOLUTION |
---|
175 | * |
---|
176 | * Value of row multiplier pi[p] in solution to the original problem |
---|
177 | * is the same as in solution to the transformed problem. |
---|
178 | * |
---|
179 | * RECOVERING MIP SOLUTION |
---|
180 | * |
---|
181 | * None needed. */ |
---|
182 | |
---|
183 | struct ineq_row |
---|
184 | { /* inequality constraint row */ |
---|
185 | int p; |
---|
186 | /* row reference number */ |
---|
187 | int s; |
---|
188 | /* column reference number for slack/surplus variable */ |
---|
189 | }; |
---|
190 | |
---|
191 | static int rcv_geq_row(NPP *npp, void *info); |
---|
192 | |
---|
193 | void npp_geq_row(NPP *npp, NPPROW *p) |
---|
194 | { /* process row of 'not less than' type */ |
---|
195 | struct ineq_row *info; |
---|
196 | NPPCOL *s; |
---|
197 | /* the row must have lower bound */ |
---|
198 | xassert(p->lb != -DBL_MAX); |
---|
199 | xassert(p->lb < p->ub); |
---|
200 | /* create column for surplus variable */ |
---|
201 | s = npp_add_col(npp); |
---|
202 | s->lb = 0.0; |
---|
203 | s->ub = (p->ub == +DBL_MAX ? +DBL_MAX : p->ub - p->lb); |
---|
204 | /* and add it to the transformed problem */ |
---|
205 | npp_add_aij(npp, p, s, -1.0); |
---|
206 | /* create transformation stack entry */ |
---|
207 | info = npp_push_tse(npp, |
---|
208 | rcv_geq_row, sizeof(struct ineq_row)); |
---|
209 | info->p = p->i; |
---|
210 | info->s = s->j; |
---|
211 | /* replace the row by equality constraint */ |
---|
212 | p->ub = p->lb; |
---|
213 | return; |
---|
214 | } |
---|
215 | |
---|
216 | static int rcv_geq_row(NPP *npp, void *_info) |
---|
217 | { /* recover row of 'not less than' type */ |
---|
218 | struct ineq_row *info = _info; |
---|
219 | if (npp->sol == GLP_SOL) |
---|
220 | { if (npp->r_stat[info->p] == GLP_BS) |
---|
221 | { if (npp->c_stat[info->s] == GLP_BS) |
---|
222 | { npp_error(); |
---|
223 | return 1; |
---|
224 | } |
---|
225 | else if (npp->c_stat[info->s] == GLP_NL || |
---|
226 | npp->c_stat[info->s] == GLP_NU) |
---|
227 | npp->r_stat[info->p] = GLP_BS; |
---|
228 | else |
---|
229 | { npp_error(); |
---|
230 | return 1; |
---|
231 | } |
---|
232 | } |
---|
233 | else if (npp->r_stat[info->p] == GLP_NS) |
---|
234 | { if (npp->c_stat[info->s] == GLP_BS) |
---|
235 | npp->r_stat[info->p] = GLP_BS; |
---|
236 | else if (npp->c_stat[info->s] == GLP_NL) |
---|
237 | npp->r_stat[info->p] = GLP_NL; |
---|
238 | else if (npp->c_stat[info->s] == GLP_NU) |
---|
239 | npp->r_stat[info->p] = GLP_NU; |
---|
240 | else |
---|
241 | { npp_error(); |
---|
242 | return 1; |
---|
243 | } |
---|
244 | } |
---|
245 | else |
---|
246 | { npp_error(); |
---|
247 | return 1; |
---|
248 | } |
---|
249 | } |
---|
250 | return 0; |
---|
251 | } |
---|
252 | |
---|
253 | /*********************************************************************** |
---|
254 | * NAME |
---|
255 | * |
---|
256 | * npp_leq_row - process row of 'not greater than' type |
---|
257 | * |
---|
258 | * SYNOPSIS |
---|
259 | * |
---|
260 | * #include "glpnpp.h" |
---|
261 | * void npp_leq_row(NPP *npp, NPPROW *p); |
---|
262 | * |
---|
263 | * DESCRIPTION |
---|
264 | * |
---|
265 | * The routine npp_leq_row processes row p, which is 'not greater than' |
---|
266 | * inequality constraint: |
---|
267 | * |
---|
268 | * (L[p] <=) sum a[p,j] x[j] <= U[p], (1) |
---|
269 | * j |
---|
270 | * |
---|
271 | * where L[p] < U[p], and lower bound may not exist (L[p] = +oo). |
---|
272 | * |
---|
273 | * PROBLEM TRANSFORMATION |
---|
274 | * |
---|
275 | * Constraint (1) can be replaced by equality constraint: |
---|
276 | * |
---|
277 | * sum a[p,j] x[j] + s = L[p], (2) |
---|
278 | * j |
---|
279 | * |
---|
280 | * where |
---|
281 | * |
---|
282 | * 0 <= s (<= U[p] - L[p]) (3) |
---|
283 | * |
---|
284 | * is a non-negative slack variable. |
---|
285 | * |
---|
286 | * Since in the primal system there appears column s having the only |
---|
287 | * non-zero coefficient in row p, in the dual system there appears a |
---|
288 | * new row: |
---|
289 | * |
---|
290 | * (+1) pi[p] + lambda = 0, (4) |
---|
291 | * |
---|
292 | * where (+1) is coefficient of column s in row p, pi[p] is multiplier |
---|
293 | * of row p, lambda is multiplier of column q, 0 is coefficient of |
---|
294 | * column s in the objective row. |
---|
295 | * |
---|
296 | * RECOVERING BASIC SOLUTION |
---|
297 | * |
---|
298 | * Status of row p in solution to the original problem is determined |
---|
299 | * by its status and status of column q in solution to the transformed |
---|
300 | * problem as follows: |
---|
301 | * |
---|
302 | * +--------------------------------------+------------------+ |
---|
303 | * | Transformed problem | Original problem | |
---|
304 | * +-----------------+--------------------+------------------+ |
---|
305 | * | Status of row p | Status of column s | Status of row p | |
---|
306 | * +-----------------+--------------------+------------------+ |
---|
307 | * | GLP_BS | GLP_BS | N/A | |
---|
308 | * | GLP_BS | GLP_NL | GLP_BS | |
---|
309 | * | GLP_BS | GLP_NU | GLP_BS | |
---|
310 | * | GLP_NS | GLP_BS | GLP_BS | |
---|
311 | * | GLP_NS | GLP_NL | GLP_NU | |
---|
312 | * | GLP_NS | GLP_NU | GLP_NL | |
---|
313 | * +-----------------+--------------------+------------------+ |
---|
314 | * |
---|
315 | * Value of row multiplier pi[p] in solution to the original problem |
---|
316 | * is the same as in solution to the transformed problem. |
---|
317 | * |
---|
318 | * 1. In solution to the transformed problem row p and column q cannot |
---|
319 | * be basic at the same time; otherwise the basis matrix would have |
---|
320 | * two linear dependent columns: unity column of auxiliary variable |
---|
321 | * of row p and unity column of variable s. |
---|
322 | * |
---|
323 | * 2. Though in the transformed problem row p is equality constraint, |
---|
324 | * it may be basic due to primal degeneracy. |
---|
325 | * |
---|
326 | * RECOVERING INTERIOR-POINT SOLUTION |
---|
327 | * |
---|
328 | * Value of row multiplier pi[p] in solution to the original problem |
---|
329 | * is the same as in solution to the transformed problem. |
---|
330 | * |
---|
331 | * RECOVERING MIP SOLUTION |
---|
332 | * |
---|
333 | * None needed. */ |
---|
334 | |
---|
335 | static int rcv_leq_row(NPP *npp, void *info); |
---|
336 | |
---|
337 | void npp_leq_row(NPP *npp, NPPROW *p) |
---|
338 | { /* process row of 'not greater than' type */ |
---|
339 | struct ineq_row *info; |
---|
340 | NPPCOL *s; |
---|
341 | /* the row must have upper bound */ |
---|
342 | xassert(p->ub != +DBL_MAX); |
---|
343 | xassert(p->lb < p->ub); |
---|
344 | /* create column for slack variable */ |
---|
345 | s = npp_add_col(npp); |
---|
346 | s->lb = 0.0; |
---|
347 | s->ub = (p->lb == -DBL_MAX ? +DBL_MAX : p->ub - p->lb); |
---|
348 | /* and add it to the transformed problem */ |
---|
349 | npp_add_aij(npp, p, s, +1.0); |
---|
350 | /* create transformation stack entry */ |
---|
351 | info = npp_push_tse(npp, |
---|
352 | rcv_leq_row, sizeof(struct ineq_row)); |
---|
353 | info->p = p->i; |
---|
354 | info->s = s->j; |
---|
355 | /* replace the row by equality constraint */ |
---|
356 | p->lb = p->ub; |
---|
357 | return; |
---|
358 | } |
---|
359 | |
---|
360 | static int rcv_leq_row(NPP *npp, void *_info) |
---|
361 | { /* recover row of 'not greater than' type */ |
---|
362 | struct ineq_row *info = _info; |
---|
363 | if (npp->sol == GLP_SOL) |
---|
364 | { if (npp->r_stat[info->p] == GLP_BS) |
---|
365 | { if (npp->c_stat[info->s] == GLP_BS) |
---|
366 | { npp_error(); |
---|
367 | return 1; |
---|
368 | } |
---|
369 | else if (npp->c_stat[info->s] == GLP_NL || |
---|
370 | npp->c_stat[info->s] == GLP_NU) |
---|
371 | npp->r_stat[info->p] = GLP_BS; |
---|
372 | else |
---|
373 | { npp_error(); |
---|
374 | return 1; |
---|
375 | } |
---|
376 | } |
---|
377 | else if (npp->r_stat[info->p] == GLP_NS) |
---|
378 | { if (npp->c_stat[info->s] == GLP_BS) |
---|
379 | npp->r_stat[info->p] = GLP_BS; |
---|
380 | else if (npp->c_stat[info->s] == GLP_NL) |
---|
381 | npp->r_stat[info->p] = GLP_NU; |
---|
382 | else if (npp->c_stat[info->s] == GLP_NU) |
---|
383 | npp->r_stat[info->p] = GLP_NL; |
---|
384 | else |
---|
385 | { npp_error(); |
---|
386 | return 1; |
---|
387 | } |
---|
388 | } |
---|
389 | else |
---|
390 | { npp_error(); |
---|
391 | return 1; |
---|
392 | } |
---|
393 | } |
---|
394 | return 0; |
---|
395 | } |
---|
396 | |
---|
397 | /*********************************************************************** |
---|
398 | * NAME |
---|
399 | * |
---|
400 | * npp_free_col - process free (unbounded) column |
---|
401 | * |
---|
402 | * SYNOPSIS |
---|
403 | * |
---|
404 | * #include "glpnpp.h" |
---|
405 | * void npp_free_col(NPP *npp, NPPCOL *q); |
---|
406 | * |
---|
407 | * DESCRIPTION |
---|
408 | * |
---|
409 | * The routine npp_free_col processes column q, which is free (i.e. has |
---|
410 | * no finite bounds): |
---|
411 | * |
---|
412 | * -oo < x[q] < +oo. (1) |
---|
413 | * |
---|
414 | * PROBLEM TRANSFORMATION |
---|
415 | * |
---|
416 | * Free (unbounded) variable can be replaced by the difference of two |
---|
417 | * non-negative variables: |
---|
418 | * |
---|
419 | * x[q] = s' - s'', s', s'' >= 0. (2) |
---|
420 | * |
---|
421 | * Assuming that in the transformed problem x[q] becomes s', |
---|
422 | * transformation (2) causes new column s'' to appear, which differs |
---|
423 | * from column s' only in the sign of coefficients in constraint and |
---|
424 | * objective rows. Thus, if in the dual system the following row |
---|
425 | * corresponds to column s': |
---|
426 | * |
---|
427 | * sum a[i,q] pi[i] + lambda' = c[q], (3) |
---|
428 | * i |
---|
429 | * |
---|
430 | * the row which corresponds to column s'' is the following: |
---|
431 | * |
---|
432 | * sum (-a[i,q]) pi[i] + lambda'' = -c[q]. (4) |
---|
433 | * i |
---|
434 | * |
---|
435 | * Then from (3) and (4) it follows that: |
---|
436 | * |
---|
437 | * lambda' + lambda'' = 0 => lambda' = lmabda'' = 0, (5) |
---|
438 | * |
---|
439 | * where lambda' and lambda'' are multipliers for columns s' and s'', |
---|
440 | * resp. |
---|
441 | * |
---|
442 | * RECOVERING BASIC SOLUTION |
---|
443 | * |
---|
444 | * With respect to (5) status of column q in solution to the original |
---|
445 | * problem is determined by statuses of columns s' and s'' in solution |
---|
446 | * to the transformed problem as follows: |
---|
447 | * |
---|
448 | * +--------------------------------------+------------------+ |
---|
449 | * | Transformed problem | Original problem | |
---|
450 | * +------------------+-------------------+------------------+ |
---|
451 | * | Status of col s' | Status of col s'' | Status of col q | |
---|
452 | * +------------------+-------------------+------------------+ |
---|
453 | * | GLP_BS | GLP_BS | N/A | |
---|
454 | * | GLP_BS | GLP_NL | GLP_BS | |
---|
455 | * | GLP_NL | GLP_BS | GLP_BS | |
---|
456 | * | GLP_NL | GLP_NL | GLP_NF | |
---|
457 | * +------------------+-------------------+------------------+ |
---|
458 | * |
---|
459 | * Value of column q is computed with formula (2). |
---|
460 | * |
---|
461 | * 1. In solution to the transformed problem columns s' and s'' cannot |
---|
462 | * be basic at the same time, because they differ only in the sign, |
---|
463 | * hence, are linear dependent. |
---|
464 | * |
---|
465 | * 2. Though column q is free, it can be non-basic due to dual |
---|
466 | * degeneracy. |
---|
467 | * |
---|
468 | * 3. If column q is integral, columns s' and s'' are also integral. |
---|
469 | * |
---|
470 | * RECOVERING INTERIOR-POINT SOLUTION |
---|
471 | * |
---|
472 | * Value of column q is computed with formula (2). |
---|
473 | * |
---|
474 | * RECOVERING MIP SOLUTION |
---|
475 | * |
---|
476 | * Value of column q is computed with formula (2). */ |
---|
477 | |
---|
478 | struct free_col |
---|
479 | { /* free (unbounded) column */ |
---|
480 | int q; |
---|
481 | /* column reference number for variables x[q] and s' */ |
---|
482 | int s; |
---|
483 | /* column reference number for variable s'' */ |
---|
484 | }; |
---|
485 | |
---|
486 | static int rcv_free_col(NPP *npp, void *info); |
---|
487 | |
---|
488 | void npp_free_col(NPP *npp, NPPCOL *q) |
---|
489 | { /* process free (unbounded) column */ |
---|
490 | struct free_col *info; |
---|
491 | NPPCOL *s; |
---|
492 | NPPAIJ *aij; |
---|
493 | /* the column must be free */ |
---|
494 | xassert(q->lb == -DBL_MAX && q->ub == +DBL_MAX); |
---|
495 | /* variable x[q] becomes s' */ |
---|
496 | q->lb = 0.0, q->ub = +DBL_MAX; |
---|
497 | /* create variable s'' */ |
---|
498 | s = npp_add_col(npp); |
---|
499 | s->is_int = q->is_int; |
---|
500 | s->lb = 0.0, s->ub = +DBL_MAX; |
---|
501 | /* duplicate objective coefficient */ |
---|
502 | s->coef = -q->coef; |
---|
503 | /* duplicate column of the constraint matrix */ |
---|
504 | for (aij = q->ptr; aij != NULL; aij = aij->c_next) |
---|
505 | npp_add_aij(npp, aij->row, s, -aij->val); |
---|
506 | /* create transformation stack entry */ |
---|
507 | info = npp_push_tse(npp, |
---|
508 | rcv_free_col, sizeof(struct free_col)); |
---|
509 | info->q = q->j; |
---|
510 | info->s = s->j; |
---|
511 | return; |
---|
512 | } |
---|
513 | |
---|
514 | static int rcv_free_col(NPP *npp, void *_info) |
---|
515 | { /* recover free (unbounded) column */ |
---|
516 | struct free_col *info = _info; |
---|
517 | if (npp->sol == GLP_SOL) |
---|
518 | { if (npp->c_stat[info->q] == GLP_BS) |
---|
519 | { if (npp->c_stat[info->s] == GLP_BS) |
---|
520 | { npp_error(); |
---|
521 | return 1; |
---|
522 | } |
---|
523 | else if (npp->c_stat[info->s] == GLP_NL) |
---|
524 | npp->c_stat[info->q] = GLP_BS; |
---|
525 | else |
---|
526 | { npp_error(); |
---|
527 | return -1; |
---|
528 | } |
---|
529 | } |
---|
530 | else if (npp->c_stat[info->q] == GLP_NL) |
---|
531 | { if (npp->c_stat[info->s] == GLP_BS) |
---|
532 | npp->c_stat[info->q] = GLP_BS; |
---|
533 | else if (npp->c_stat[info->s] == GLP_NL) |
---|
534 | npp->c_stat[info->q] = GLP_NF; |
---|
535 | else |
---|
536 | { npp_error(); |
---|
537 | return -1; |
---|
538 | } |
---|
539 | } |
---|
540 | else |
---|
541 | { npp_error(); |
---|
542 | return -1; |
---|
543 | } |
---|
544 | } |
---|
545 | /* compute value of x[q] with formula (2) */ |
---|
546 | npp->c_value[info->q] -= npp->c_value[info->s]; |
---|
547 | return 0; |
---|
548 | } |
---|
549 | |
---|
550 | /*********************************************************************** |
---|
551 | * NAME |
---|
552 | * |
---|
553 | * npp_lbnd_col - process column with (non-zero) lower bound |
---|
554 | * |
---|
555 | * SYNOPSIS |
---|
556 | * |
---|
557 | * #include "glpnpp.h" |
---|
558 | * void npp_lbnd_col(NPP *npp, NPPCOL *q); |
---|
559 | * |
---|
560 | * DESCRIPTION |
---|
561 | * |
---|
562 | * The routine npp_lbnd_col processes column q, which has (non-zero) |
---|
563 | * lower bound: |
---|
564 | * |
---|
565 | * l[q] <= x[q] (<= u[q]), (1) |
---|
566 | * |
---|
567 | * where l[q] < u[q], and upper bound may not exist (u[q] = +oo). |
---|
568 | * |
---|
569 | * PROBLEM TRANSFORMATION |
---|
570 | * |
---|
571 | * Column q can be replaced as follows: |
---|
572 | * |
---|
573 | * x[q] = l[q] + s, (2) |
---|
574 | * |
---|
575 | * where |
---|
576 | * |
---|
577 | * 0 <= s (<= u[q] - l[q]) (3) |
---|
578 | * |
---|
579 | * is a non-negative variable. |
---|
580 | * |
---|
581 | * Substituting x[q] from (2) into the objective row, we have: |
---|
582 | * |
---|
583 | * z = sum c[j] x[j] + c0 = |
---|
584 | * j |
---|
585 | * |
---|
586 | * = sum c[j] x[j] + c[q] x[q] + c0 = |
---|
587 | * j!=q |
---|
588 | * |
---|
589 | * = sum c[j] x[j] + c[q] (l[q] + s) + c0 = |
---|
590 | * j!=q |
---|
591 | * |
---|
592 | * = sum c[j] x[j] + c[q] s + c~0, |
---|
593 | * |
---|
594 | * where |
---|
595 | * |
---|
596 | * c~0 = c0 + c[q] l[q] (4) |
---|
597 | * |
---|
598 | * is the constant term of the objective in the transformed problem. |
---|
599 | * Similarly, substituting x[q] into constraint row i, we have: |
---|
600 | * |
---|
601 | * L[i] <= sum a[i,j] x[j] <= U[i] ==> |
---|
602 | * j |
---|
603 | * |
---|
604 | * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> |
---|
605 | * j!=q |
---|
606 | * |
---|
607 | * L[i] <= sum a[i,j] x[j] + a[i,q] (l[q] + s) <= U[i] ==> |
---|
608 | * j!=q |
---|
609 | * |
---|
610 | * L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i], |
---|
611 | * j!=q |
---|
612 | * |
---|
613 | * where |
---|
614 | * |
---|
615 | * L~[i] = L[i] - a[i,q] l[q], U~[i] = U[i] - a[i,q] l[q] (5) |
---|
616 | * |
---|
617 | * are lower and upper bounds of row i in the transformed problem, |
---|
618 | * resp. |
---|
619 | * |
---|
620 | * Transformation (2) does not affect the dual system. |
---|
621 | * |
---|
622 | * RECOVERING BASIC SOLUTION |
---|
623 | * |
---|
624 | * Status of column q in solution to the original problem is the same |
---|
625 | * as in solution to the transformed problem (GLP_BS, GLP_NL or GLP_NU). |
---|
626 | * Value of column q is computed with formula (2). |
---|
627 | * |
---|
628 | * RECOVERING INTERIOR-POINT SOLUTION |
---|
629 | * |
---|
630 | * Value of column q is computed with formula (2). |
---|
631 | * |
---|
632 | * RECOVERING MIP SOLUTION |
---|
633 | * |
---|
634 | * Value of column q is computed with formula (2). */ |
---|
635 | |
---|
636 | struct bnd_col |
---|
637 | { /* bounded column */ |
---|
638 | int q; |
---|
639 | /* column reference number for variables x[q] and s */ |
---|
640 | double bnd; |
---|
641 | /* lower/upper bound l[q] or u[q] */ |
---|
642 | }; |
---|
643 | |
---|
644 | static int rcv_lbnd_col(NPP *npp, void *info); |
---|
645 | |
---|
646 | void npp_lbnd_col(NPP *npp, NPPCOL *q) |
---|
647 | { /* process column with (non-zero) lower bound */ |
---|
648 | struct bnd_col *info; |
---|
649 | NPPROW *i; |
---|
650 | NPPAIJ *aij; |
---|
651 | /* the column must have non-zero lower bound */ |
---|
652 | xassert(q->lb != 0.0); |
---|
653 | xassert(q->lb != -DBL_MAX); |
---|
654 | xassert(q->lb < q->ub); |
---|
655 | /* create transformation stack entry */ |
---|
656 | info = npp_push_tse(npp, |
---|
657 | rcv_lbnd_col, sizeof(struct bnd_col)); |
---|
658 | info->q = q->j; |
---|
659 | info->bnd = q->lb; |
---|
660 | /* substitute x[q] into objective row */ |
---|
661 | npp->c0 += q->coef * q->lb; |
---|
662 | /* substitute x[q] into constraint rows */ |
---|
663 | for (aij = q->ptr; aij != NULL; aij = aij->c_next) |
---|
664 | { i = aij->row; |
---|
665 | if (i->lb == i->ub) |
---|
666 | i->ub = (i->lb -= aij->val * q->lb); |
---|
667 | else |
---|
668 | { if (i->lb != -DBL_MAX) |
---|
669 | i->lb -= aij->val * q->lb; |
---|
670 | if (i->ub != +DBL_MAX) |
---|
671 | i->ub -= aij->val * q->lb; |
---|
672 | } |
---|
673 | } |
---|
674 | /* column x[q] becomes column s */ |
---|
675 | if (q->ub != +DBL_MAX) |
---|
676 | q->ub -= q->lb; |
---|
677 | q->lb = 0.0; |
---|
678 | return; |
---|
679 | } |
---|
680 | |
---|
681 | static int rcv_lbnd_col(NPP *npp, void *_info) |
---|
682 | { /* recover column with (non-zero) lower bound */ |
---|
683 | struct bnd_col *info = _info; |
---|
684 | if (npp->sol == GLP_SOL) |
---|
685 | { if (npp->c_stat[info->q] == GLP_BS || |
---|
686 | npp->c_stat[info->q] == GLP_NL || |
---|
687 | npp->c_stat[info->q] == GLP_NU) |
---|
688 | npp->c_stat[info->q] = npp->c_stat[info->q]; |
---|
689 | else |
---|
690 | { npp_error(); |
---|
691 | return 1; |
---|
692 | } |
---|
693 | } |
---|
694 | /* compute value of x[q] with formula (2) */ |
---|
695 | npp->c_value[info->q] = info->bnd + npp->c_value[info->q]; |
---|
696 | return 0; |
---|
697 | } |
---|
698 | |
---|
699 | /*********************************************************************** |
---|
700 | * NAME |
---|
701 | * |
---|
702 | * npp_ubnd_col - process column with upper bound |
---|
703 | * |
---|
704 | * SYNOPSIS |
---|
705 | * |
---|
706 | * #include "glpnpp.h" |
---|
707 | * void npp_ubnd_col(NPP *npp, NPPCOL *q); |
---|
708 | * |
---|
709 | * DESCRIPTION |
---|
710 | * |
---|
711 | * The routine npp_ubnd_col processes column q, which has upper bound: |
---|
712 | * |
---|
713 | * (l[q] <=) x[q] <= u[q], (1) |
---|
714 | * |
---|
715 | * where l[q] < u[q], and lower bound may not exist (l[q] = -oo). |
---|
716 | * |
---|
717 | * PROBLEM TRANSFORMATION |
---|
718 | * |
---|
719 | * Column q can be replaced as follows: |
---|
720 | * |
---|
721 | * x[q] = u[q] - s, (2) |
---|
722 | * |
---|
723 | * where |
---|
724 | * |
---|
725 | * 0 <= s (<= u[q] - l[q]) (3) |
---|
726 | * |
---|
727 | * is a non-negative variable. |
---|
728 | * |
---|
729 | * Substituting x[q] from (2) into the objective row, we have: |
---|
730 | * |
---|
731 | * z = sum c[j] x[j] + c0 = |
---|
732 | * j |
---|
733 | * |
---|
734 | * = sum c[j] x[j] + c[q] x[q] + c0 = |
---|
735 | * j!=q |
---|
736 | * |
---|
737 | * = sum c[j] x[j] + c[q] (u[q] - s) + c0 = |
---|
738 | * j!=q |
---|
739 | * |
---|
740 | * = sum c[j] x[j] - c[q] s + c~0, |
---|
741 | * |
---|
742 | * where |
---|
743 | * |
---|
744 | * c~0 = c0 + c[q] u[q] (4) |
---|
745 | * |
---|
746 | * is the constant term of the objective in the transformed problem. |
---|
747 | * Similarly, substituting x[q] into constraint row i, we have: |
---|
748 | * |
---|
749 | * L[i] <= sum a[i,j] x[j] <= U[i] ==> |
---|
750 | * j |
---|
751 | * |
---|
752 | * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> |
---|
753 | * j!=q |
---|
754 | * |
---|
755 | * L[i] <= sum a[i,j] x[j] + a[i,q] (u[q] - s) <= U[i] ==> |
---|
756 | * j!=q |
---|
757 | * |
---|
758 | * L~[i] <= sum a[i,j] x[j] - a[i,q] s <= U~[i], |
---|
759 | * j!=q |
---|
760 | * |
---|
761 | * where |
---|
762 | * |
---|
763 | * L~[i] = L[i] - a[i,q] u[q], U~[i] = U[i] - a[i,q] u[q] (5) |
---|
764 | * |
---|
765 | * are lower and upper bounds of row i in the transformed problem, |
---|
766 | * resp. |
---|
767 | * |
---|
768 | * Note that in the transformed problem coefficients c[q] and a[i,q] |
---|
769 | * change their sign. Thus, the row of the dual system corresponding to |
---|
770 | * column q: |
---|
771 | * |
---|
772 | * sum a[i,q] pi[i] + lambda[q] = c[q] (6) |
---|
773 | * i |
---|
774 | * |
---|
775 | * in the transformed problem becomes the following: |
---|
776 | * |
---|
777 | * sum (-a[i,q]) pi[i] + lambda[s] = -c[q]. (7) |
---|
778 | * i |
---|
779 | * |
---|
780 | * Therefore: |
---|
781 | * |
---|
782 | * lambda[q] = - lambda[s], (8) |
---|
783 | * |
---|
784 | * where lambda[q] is multiplier for column q, lambda[s] is multiplier |
---|
785 | * for column s. |
---|
786 | * |
---|
787 | * RECOVERING BASIC SOLUTION |
---|
788 | * |
---|
789 | * With respect to (8) status of column q in solution to the original |
---|
790 | * problem is determined by status of column s in solution to the |
---|
791 | * transformed problem as follows: |
---|
792 | * |
---|
793 | * +-----------------------+--------------------+ |
---|
794 | * | Status of column s | Status of column q | |
---|
795 | * | (transformed problem) | (original problem) | |
---|
796 | * +-----------------------+--------------------+ |
---|
797 | * | GLP_BS | GLP_BS | |
---|
798 | * | GLP_NL | GLP_NU | |
---|
799 | * | GLP_NU | GLP_NL | |
---|
800 | * +-----------------------+--------------------+ |
---|
801 | * |
---|
802 | * Value of column q is computed with formula (2). |
---|
803 | * |
---|
804 | * RECOVERING INTERIOR-POINT SOLUTION |
---|
805 | * |
---|
806 | * Value of column q is computed with formula (2). |
---|
807 | * |
---|
808 | * RECOVERING MIP SOLUTION |
---|
809 | * |
---|
810 | * Value of column q is computed with formula (2). */ |
---|
811 | |
---|
812 | static int rcv_ubnd_col(NPP *npp, void *info); |
---|
813 | |
---|
814 | void npp_ubnd_col(NPP *npp, NPPCOL *q) |
---|
815 | { /* process column with upper bound */ |
---|
816 | struct bnd_col *info; |
---|
817 | NPPROW *i; |
---|
818 | NPPAIJ *aij; |
---|
819 | /* the column must have upper bound */ |
---|
820 | xassert(q->ub != +DBL_MAX); |
---|
821 | xassert(q->lb < q->ub); |
---|
822 | /* create transformation stack entry */ |
---|
823 | info = npp_push_tse(npp, |
---|
824 | rcv_ubnd_col, sizeof(struct bnd_col)); |
---|
825 | info->q = q->j; |
---|
826 | info->bnd = q->ub; |
---|
827 | /* substitute x[q] into objective row */ |
---|
828 | npp->c0 += q->coef * q->ub; |
---|
829 | q->coef = -q->coef; |
---|
830 | /* substitute x[q] into constraint rows */ |
---|
831 | for (aij = q->ptr; aij != NULL; aij = aij->c_next) |
---|
832 | { i = aij->row; |
---|
833 | if (i->lb == i->ub) |
---|
834 | i->ub = (i->lb -= aij->val * q->ub); |
---|
835 | else |
---|
836 | { if (i->lb != -DBL_MAX) |
---|
837 | i->lb -= aij->val * q->ub; |
---|
838 | if (i->ub != +DBL_MAX) |
---|
839 | i->ub -= aij->val * q->ub; |
---|
840 | } |
---|
841 | aij->val = -aij->val; |
---|
842 | } |
---|
843 | /* column x[q] becomes column s */ |
---|
844 | if (q->lb != -DBL_MAX) |
---|
845 | q->ub -= q->lb; |
---|
846 | else |
---|
847 | q->ub = +DBL_MAX; |
---|
848 | q->lb = 0.0; |
---|
849 | return; |
---|
850 | } |
---|
851 | |
---|
852 | static int rcv_ubnd_col(NPP *npp, void *_info) |
---|
853 | { /* recover column with upper bound */ |
---|
854 | struct bnd_col *info = _info; |
---|
855 | if (npp->sol == GLP_BS) |
---|
856 | { if (npp->c_stat[info->q] == GLP_BS) |
---|
857 | npp->c_stat[info->q] = GLP_BS; |
---|
858 | else if (npp->c_stat[info->q] == GLP_NL) |
---|
859 | npp->c_stat[info->q] = GLP_NU; |
---|
860 | else if (npp->c_stat[info->q] == GLP_NU) |
---|
861 | npp->c_stat[info->q] = GLP_NL; |
---|
862 | else |
---|
863 | { npp_error(); |
---|
864 | return 1; |
---|
865 | } |
---|
866 | } |
---|
867 | /* compute value of x[q] with formula (2) */ |
---|
868 | npp->c_value[info->q] = info->bnd - npp->c_value[info->q]; |
---|
869 | return 0; |
---|
870 | } |
---|
871 | |
---|
872 | /*********************************************************************** |
---|
873 | * NAME |
---|
874 | * |
---|
875 | * npp_dbnd_col - process non-negative column with upper bound |
---|
876 | * |
---|
877 | * SYNOPSIS |
---|
878 | * |
---|
879 | * #include "glpnpp.h" |
---|
880 | * void npp_dbnd_col(NPP *npp, NPPCOL *q); |
---|
881 | * |
---|
882 | * DESCRIPTION |
---|
883 | * |
---|
884 | * The routine npp_dbnd_col processes column q, which is non-negative |
---|
885 | * and has upper bound: |
---|
886 | * |
---|
887 | * 0 <= x[q] <= u[q], (1) |
---|
888 | * |
---|
889 | * where u[q] > 0. |
---|
890 | * |
---|
891 | * PROBLEM TRANSFORMATION |
---|
892 | * |
---|
893 | * Upper bound of column q can be replaced by the following equality |
---|
894 | * constraint: |
---|
895 | * |
---|
896 | * x[q] + s = u[q], (2) |
---|
897 | * |
---|
898 | * where s >= 0 is a non-negative complement variable. |
---|
899 | * |
---|
900 | * Since in the primal system along with new row (2) there appears a |
---|
901 | * new column s having the only non-zero coefficient in this row, in |
---|
902 | * the dual system there appears a new row: |
---|
903 | * |
---|
904 | * (+1)pi + lambda[s] = 0, (3) |
---|
905 | * |
---|
906 | * where (+1) is coefficient at column s in row (2), pi is multiplier |
---|
907 | * for row (2), lambda[s] is multiplier for column s, 0 is coefficient |
---|
908 | * at column s in the objective row. |
---|
909 | * |
---|
910 | * RECOVERING BASIC SOLUTION |
---|
911 | * |
---|
912 | * Status of column q in solution to the original problem is determined |
---|
913 | * by its status and status of column s in solution to the transformed |
---|
914 | * problem as follows: |
---|
915 | * |
---|
916 | * +-----------------------------------+------------------+ |
---|
917 | * | Transformed problem | Original problem | |
---|
918 | * +-----------------+-----------------+------------------+ |
---|
919 | * | Status of col q | Status of col s | Status of col q | |
---|
920 | * +-----------------+-----------------+------------------+ |
---|
921 | * | GLP_BS | GLP_BS | GLP_BS | |
---|
922 | * | GLP_BS | GLP_NL | GLP_NU | |
---|
923 | * | GLP_NL | GLP_BS | GLP_NL | |
---|
924 | * | GLP_NL | GLP_NL | GLP_NL (*) | |
---|
925 | * +-----------------+-----------------+------------------+ |
---|
926 | * |
---|
927 | * Value of column q in solution to the original problem is the same as |
---|
928 | * in solution to the transformed problem. |
---|
929 | * |
---|
930 | * 1. Formally, in solution to the transformed problem columns q and s |
---|
931 | * cannot be non-basic at the same time, since the constraint (2) |
---|
932 | * would be violated. However, if u[q] is close to zero, violation |
---|
933 | * may be less than a working precision even if both columns q and s |
---|
934 | * are non-basic. In this degenerate case row (2) can be only basic, |
---|
935 | * i.e. non-active constraint (otherwise corresponding row of the |
---|
936 | * basis matrix would be zero). This allows to pivot out auxiliary |
---|
937 | * variable and pivot in column s, in which case the row becomes |
---|
938 | * active while column s becomes basic. |
---|
939 | * |
---|
940 | * 2. If column q is integral, column s is also integral. |
---|
941 | * |
---|
942 | * RECOVERING INTERIOR-POINT SOLUTION |
---|
943 | * |
---|
944 | * Value of column q in solution to the original problem is the same as |
---|
945 | * in solution to the transformed problem. |
---|
946 | * |
---|
947 | * RECOVERING MIP SOLUTION |
---|
948 | * |
---|
949 | * Value of column q in solution to the original problem is the same as |
---|
950 | * in solution to the transformed problem. */ |
---|
951 | |
---|
952 | struct dbnd_col |
---|
953 | { /* double-bounded column */ |
---|
954 | int q; |
---|
955 | /* column reference number for variable x[q] */ |
---|
956 | int s; |
---|
957 | /* column reference number for complement variable s */ |
---|
958 | }; |
---|
959 | |
---|
960 | static int rcv_dbnd_col(NPP *npp, void *info); |
---|
961 | |
---|
962 | void npp_dbnd_col(NPP *npp, NPPCOL *q) |
---|
963 | { /* process non-negative column with upper bound */ |
---|
964 | struct dbnd_col *info; |
---|
965 | NPPROW *p; |
---|
966 | NPPCOL *s; |
---|
967 | /* the column must be non-negative with upper bound */ |
---|
968 | xassert(q->lb == 0.0); |
---|
969 | xassert(q->ub > 0.0); |
---|
970 | xassert(q->ub != +DBL_MAX); |
---|
971 | /* create variable s */ |
---|
972 | s = npp_add_col(npp); |
---|
973 | s->is_int = q->is_int; |
---|
974 | s->lb = 0.0, s->ub = +DBL_MAX; |
---|
975 | /* create equality constraint (2) */ |
---|
976 | p = npp_add_row(npp); |
---|
977 | p->lb = p->ub = q->ub; |
---|
978 | npp_add_aij(npp, p, q, +1.0); |
---|
979 | npp_add_aij(npp, p, s, +1.0); |
---|
980 | /* create transformation stack entry */ |
---|
981 | info = npp_push_tse(npp, |
---|
982 | rcv_dbnd_col, sizeof(struct dbnd_col)); |
---|
983 | info->q = q->j; |
---|
984 | info->s = s->j; |
---|
985 | /* remove upper bound of x[q] */ |
---|
986 | q->ub = +DBL_MAX; |
---|
987 | return; |
---|
988 | } |
---|
989 | |
---|
990 | static int rcv_dbnd_col(NPP *npp, void *_info) |
---|
991 | { /* recover non-negative column with upper bound */ |
---|
992 | struct dbnd_col *info = _info; |
---|
993 | if (npp->sol == GLP_BS) |
---|
994 | { if (npp->c_stat[info->q] == GLP_BS) |
---|
995 | { if (npp->c_stat[info->s] == GLP_BS) |
---|
996 | npp->c_stat[info->q] = GLP_BS; |
---|
997 | else if (npp->c_stat[info->s] == GLP_NL) |
---|
998 | npp->c_stat[info->q] = GLP_NU; |
---|
999 | else |
---|
1000 | { npp_error(); |
---|
1001 | return 1; |
---|
1002 | } |
---|
1003 | } |
---|
1004 | else if (npp->c_stat[info->q] == GLP_NL) |
---|
1005 | { if (npp->c_stat[info->s] == GLP_BS || |
---|
1006 | npp->c_stat[info->s] == GLP_NL) |
---|
1007 | npp->c_stat[info->q] = GLP_NL; |
---|
1008 | else |
---|
1009 | { npp_error(); |
---|
1010 | return 1; |
---|
1011 | } |
---|
1012 | } |
---|
1013 | else |
---|
1014 | { npp_error(); |
---|
1015 | return 1; |
---|
1016 | } |
---|
1017 | } |
---|
1018 | return 0; |
---|
1019 | } |
---|
1020 | |
---|
1021 | /*********************************************************************** |
---|
1022 | * NAME |
---|
1023 | * |
---|
1024 | * npp_fixed_col - process fixed column |
---|
1025 | * |
---|
1026 | * SYNOPSIS |
---|
1027 | * |
---|
1028 | * #include "glpnpp.h" |
---|
1029 | * void npp_fixed_col(NPP *npp, NPPCOL *q); |
---|
1030 | * |
---|
1031 | * DESCRIPTION |
---|
1032 | * |
---|
1033 | * The routine npp_fixed_col processes column q, which is fixed: |
---|
1034 | * |
---|
1035 | * x[q] = s[q], (1) |
---|
1036 | * |
---|
1037 | * where s[q] is a fixed column value. |
---|
1038 | * |
---|
1039 | * PROBLEM TRANSFORMATION |
---|
1040 | * |
---|
1041 | * The value of a fixed column can be substituted into the objective |
---|
1042 | * and constraint rows that allows removing the column from the problem. |
---|
1043 | * |
---|
1044 | * Substituting x[q] = s[q] into the objective row, we have: |
---|
1045 | * |
---|
1046 | * z = sum c[j] x[j] + c0 = |
---|
1047 | * j |
---|
1048 | * |
---|
1049 | * = sum c[j] x[j] + c[q] x[q] + c0 = |
---|
1050 | * j!=q |
---|
1051 | * |
---|
1052 | * = sum c[j] x[j] + c[q] s[q] + c0 = |
---|
1053 | * j!=q |
---|
1054 | * |
---|
1055 | * = sum c[j] x[j] + c~0, |
---|
1056 | * j!=q |
---|
1057 | * |
---|
1058 | * where |
---|
1059 | * |
---|
1060 | * c~0 = c0 + c[q] s[q] (2) |
---|
1061 | * |
---|
1062 | * is the constant term of the objective in the transformed problem. |
---|
1063 | * Similarly, substituting x[q] = s[q] into constraint row i, we have: |
---|
1064 | * |
---|
1065 | * L[i] <= sum a[i,j] x[j] <= U[i] ==> |
---|
1066 | * j |
---|
1067 | * |
---|
1068 | * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==> |
---|
1069 | * j!=q |
---|
1070 | * |
---|
1071 | * L[i] <= sum a[i,j] x[j] + a[i,q] s[q] <= U[i] ==> |
---|
1072 | * j!=q |
---|
1073 | * |
---|
1074 | * L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i], |
---|
1075 | * j!=q |
---|
1076 | * |
---|
1077 | * where |
---|
1078 | * |
---|
1079 | * L~[i] = L[i] - a[i,q] s[q], U~[i] = U[i] - a[i,q] s[q] (3) |
---|
1080 | * |
---|
1081 | * are lower and upper bounds of row i in the transformed problem, |
---|
1082 | * resp. |
---|
1083 | * |
---|
1084 | * RECOVERING BASIC SOLUTION |
---|
1085 | * |
---|
1086 | * Column q is assigned status GLP_NS and its value is assigned s[q]. |
---|
1087 | * |
---|
1088 | * RECOVERING INTERIOR-POINT SOLUTION |
---|
1089 | * |
---|
1090 | * Value of column q is assigned s[q]. |
---|
1091 | * |
---|
1092 | * RECOVERING MIP SOLUTION |
---|
1093 | * |
---|
1094 | * Value of column q is assigned s[q]. */ |
---|
1095 | |
---|
1096 | struct fixed_col |
---|
1097 | { /* fixed column */ |
---|
1098 | int q; |
---|
1099 | /* column reference number for variable x[q] */ |
---|
1100 | double s; |
---|
1101 | /* value, at which x[q] is fixed */ |
---|
1102 | }; |
---|
1103 | |
---|
1104 | static int rcv_fixed_col(NPP *npp, void *info); |
---|
1105 | |
---|
1106 | void npp_fixed_col(NPP *npp, NPPCOL *q) |
---|
1107 | { /* process fixed column */ |
---|
1108 | struct fixed_col *info; |
---|
1109 | NPPROW *i; |
---|
1110 | NPPAIJ *aij; |
---|
1111 | /* the column must be fixed */ |
---|
1112 | xassert(q->lb == q->ub); |
---|
1113 | /* create transformation stack entry */ |
---|
1114 | info = npp_push_tse(npp, |
---|
1115 | rcv_fixed_col, sizeof(struct fixed_col)); |
---|
1116 | info->q = q->j; |
---|
1117 | info->s = q->lb; |
---|
1118 | /* substitute x[q] = s[q] into objective row */ |
---|
1119 | npp->c0 += q->coef * q->lb; |
---|
1120 | /* substitute x[q] = s[q] into constraint rows */ |
---|
1121 | for (aij = q->ptr; aij != NULL; aij = aij->c_next) |
---|
1122 | { i = aij->row; |
---|
1123 | if (i->lb == i->ub) |
---|
1124 | i->ub = (i->lb -= aij->val * q->lb); |
---|
1125 | else |
---|
1126 | { if (i->lb != -DBL_MAX) |
---|
1127 | i->lb -= aij->val * q->lb; |
---|
1128 | if (i->ub != +DBL_MAX) |
---|
1129 | i->ub -= aij->val * q->lb; |
---|
1130 | } |
---|
1131 | } |
---|
1132 | /* remove the column from the problem */ |
---|
1133 | npp_del_col(npp, q); |
---|
1134 | return; |
---|
1135 | } |
---|
1136 | |
---|
1137 | static int rcv_fixed_col(NPP *npp, void *_info) |
---|
1138 | { /* recover fixed column */ |
---|
1139 | struct fixed_col *info = _info; |
---|
1140 | if (npp->sol == GLP_SOL) |
---|
1141 | npp->c_stat[info->q] = GLP_NS; |
---|
1142 | npp->c_value[info->q] = info->s; |
---|
1143 | return 0; |
---|
1144 | } |
---|
1145 | |
---|
1146 | /*********************************************************************** |
---|
1147 | * NAME |
---|
1148 | * |
---|
1149 | * npp_make_equality - process row with almost identical bounds |
---|
1150 | * |
---|
1151 | * SYNOPSIS |
---|
1152 | * |
---|
1153 | * #include "glpnpp.h" |
---|
1154 | * int npp_make_equality(NPP *npp, NPPROW *p); |
---|
1155 | * |
---|
1156 | * DESCRIPTION |
---|
1157 | * |
---|
1158 | * The routine npp_make_equality processes row p: |
---|
1159 | * |
---|
1160 | * L[p] <= sum a[p,j] x[j] <= U[p], (1) |
---|
1161 | * j |
---|
1162 | * |
---|
1163 | * where -oo < L[p] < U[p] < +oo, i.e. which is double-sided inequality |
---|
1164 | * constraint. |
---|
1165 | * |
---|
1166 | * RETURNS |
---|
1167 | * |
---|
1168 | * 0 - row bounds have not been changed; |
---|
1169 | * |
---|
1170 | * 1 - row has been replaced by equality constraint. |
---|
1171 | * |
---|
1172 | * PROBLEM TRANSFORMATION |
---|
1173 | * |
---|
1174 | * If bounds of row (1) are very close to each other: |
---|
1175 | * |
---|
1176 | * U[p] - L[p] <= eps, (2) |
---|
1177 | * |
---|
1178 | * where eps is an absolute tolerance for row value, the row can be |
---|
1179 | * replaced by the following almost equivalent equiality constraint: |
---|
1180 | * |
---|
1181 | * sum a[p,j] x[j] = b, (3) |
---|
1182 | * j |
---|
1183 | * |
---|
1184 | * where b = (L[p] + U[p]) / 2. If the right-hand side in (3) happens |
---|
1185 | * to be very close to its nearest integer: |
---|
1186 | * |
---|
1187 | * |b - floor(b + 0.5)| <= eps, (4) |
---|
1188 | * |
---|
1189 | * it is reasonable to use this nearest integer as the right-hand side. |
---|
1190 | * |
---|
1191 | * RECOVERING BASIC SOLUTION |
---|
1192 | * |
---|
1193 | * Status of row p in solution to the original problem is determined |
---|
1194 | * by its status and the sign of its multiplier pi[p] in solution to |
---|
1195 | * the transformed problem as follows: |
---|
1196 | * |
---|
1197 | * +-----------------------+---------+--------------------+ |
---|
1198 | * | Status of row p | Sign of | Status of row p | |
---|
1199 | * | (transformed problem) | pi[p] | (original problem) | |
---|
1200 | * +-----------------------+---------+--------------------+ |
---|
1201 | * | GLP_BS | + / - | GLP_BS | |
---|
1202 | * | GLP_NS | + | GLP_NL | |
---|
1203 | * | GLP_NS | - | GLP_NU | |
---|
1204 | * +-----------------------+---------+--------------------+ |
---|
1205 | * |
---|
1206 | * Value of row multiplier pi[p] in solution to the original problem is |
---|
1207 | * the same as in solution to the transformed problem. |
---|
1208 | * |
---|
1209 | * RECOVERING INTERIOR POINT SOLUTION |
---|
1210 | * |
---|
1211 | * Value of row multiplier pi[p] in solution to the original problem is |
---|
1212 | * the same as in solution to the transformed problem. |
---|
1213 | * |
---|
1214 | * RECOVERING MIP SOLUTION |
---|
1215 | * |
---|
1216 | * None needed. */ |
---|
1217 | |
---|
1218 | struct make_equality |
---|
1219 | { /* row with almost identical bounds */ |
---|
1220 | int p; |
---|
1221 | /* row reference number */ |
---|
1222 | }; |
---|
1223 | |
---|
1224 | static int rcv_make_equality(NPP *npp, void *info); |
---|
1225 | |
---|
1226 | int npp_make_equality(NPP *npp, NPPROW *p) |
---|
1227 | { /* process row with almost identical bounds */ |
---|
1228 | struct make_equality *info; |
---|
1229 | double b, eps, nint; |
---|
1230 | /* the row must be double-sided inequality */ |
---|
1231 | xassert(p->lb != -DBL_MAX); |
---|
1232 | xassert(p->ub != +DBL_MAX); |
---|
1233 | xassert(p->lb < p->ub); |
---|
1234 | /* check row bounds */ |
---|
1235 | eps = 1e-9 + 1e-12 * fabs(p->lb); |
---|
1236 | if (p->ub - p->lb > eps) return 0; |
---|
1237 | /* row bounds are very close to each other */ |
---|
1238 | /* create transformation stack entry */ |
---|
1239 | info = npp_push_tse(npp, |
---|
1240 | rcv_make_equality, sizeof(struct make_equality)); |
---|
1241 | info->p = p->i; |
---|
1242 | /* compute right-hand side */ |
---|
1243 | b = 0.5 * (p->ub + p->lb); |
---|
1244 | nint = floor(b + 0.5); |
---|
1245 | if (fabs(b - nint) <= eps) b = nint; |
---|
1246 | /* replace row p by almost equivalent equality constraint */ |
---|
1247 | p->lb = p->ub = b; |
---|
1248 | return 1; |
---|
1249 | } |
---|
1250 | |
---|
1251 | int rcv_make_equality(NPP *npp, void *_info) |
---|
1252 | { /* recover row with almost identical bounds */ |
---|
1253 | struct make_equality *info = _info; |
---|
1254 | if (npp->sol == GLP_SOL) |
---|
1255 | { if (npp->r_stat[info->p] == GLP_BS) |
---|
1256 | npp->r_stat[info->p] = GLP_BS; |
---|
1257 | else if (npp->r_stat[info->p] == GLP_NS) |
---|
1258 | { if (npp->r_pi[info->p] >= 0.0) |
---|
1259 | npp->r_stat[info->p] = GLP_NL; |
---|
1260 | else |
---|
1261 | npp->r_stat[info->p] = GLP_NU; |
---|
1262 | } |
---|
1263 | else |
---|
1264 | { npp_error(); |
---|
1265 | return 1; |
---|
1266 | } |
---|
1267 | } |
---|
1268 | return 0; |
---|
1269 | } |
---|
1270 | |
---|
1271 | /*********************************************************************** |
---|
1272 | * NAME |
---|
1273 | * |
---|
1274 | * npp_make_fixed - process column with almost identical bounds |
---|
1275 | * |
---|
1276 | * SYNOPSIS |
---|
1277 | * |
---|
1278 | * #include "glpnpp.h" |
---|
1279 | * int npp_make_fixed(NPP *npp, NPPCOL *q); |
---|
1280 | * |
---|
1281 | * DESCRIPTION |
---|
1282 | * |
---|
1283 | * The routine npp_make_fixed processes column q: |
---|
1284 | * |
---|
1285 | * l[q] <= x[q] <= u[q], (1) |
---|
1286 | * |
---|
1287 | * where -oo < l[q] < u[q] < +oo, i.e. which has both lower and upper |
---|
1288 | * bounds. |
---|
1289 | * |
---|
1290 | * RETURNS |
---|
1291 | * |
---|
1292 | * 0 - column bounds have not been changed; |
---|
1293 | * |
---|
1294 | * 1 - column has been fixed. |
---|
1295 | * |
---|
1296 | * PROBLEM TRANSFORMATION |
---|
1297 | * |
---|
1298 | * If bounds of column (1) are very close to each other: |
---|
1299 | * |
---|
1300 | * u[q] - l[q] <= eps, (2) |
---|
1301 | * |
---|
1302 | * where eps is an absolute tolerance for column value, the column can |
---|
1303 | * be fixed: |
---|
1304 | * |
---|
1305 | * x[q] = s[q], (3) |
---|
1306 | * |
---|
1307 | * where s[q] = (l[q] + u[q]) / 2. And if the fixed column value s[q] |
---|
1308 | * happens to be very close to its nearest integer: |
---|
1309 | * |
---|
1310 | * |s[q] - floor(s[q] + 0.5)| <= eps, (4) |
---|
1311 | * |
---|
1312 | * it is reasonable to use this nearest integer as the fixed value. |
---|
1313 | * |
---|
1314 | * RECOVERING BASIC SOLUTION |
---|
1315 | * |
---|
1316 | * In the dual system of the original (as well as transformed) problem |
---|
1317 | * column q corresponds to the following row: |
---|
1318 | * |
---|
1319 | * sum a[i,q] pi[i] + lambda[q] = c[q]. (5) |
---|
1320 | * i |
---|
1321 | * |
---|
1322 | * Since multipliers pi[i] are known for all rows from solution to the |
---|
1323 | * transformed problem, formula (5) allows computing value of multiplier |
---|
1324 | * (reduced cost) for column q: |
---|
1325 | * |
---|
1326 | * lambda[q] = c[q] - sum a[i,q] pi[i]. (6) |
---|
1327 | * i |
---|
1328 | * |
---|
1329 | * Status of column q in solution to the original problem is determined |
---|
1330 | * by its status and the sign of its multiplier lambda[q] in solution to |
---|
1331 | * the transformed problem as follows: |
---|
1332 | * |
---|
1333 | * +-----------------------+-----------+--------------------+ |
---|
1334 | * | Status of column q | Sign of | Status of column q | |
---|
1335 | * | (transformed problem) | lambda[q] | (original problem) | |
---|
1336 | * +-----------------------+-----------+--------------------+ |
---|
1337 | * | GLP_BS | + / - | GLP_BS | |
---|
1338 | * | GLP_NS | + | GLP_NL | |
---|
1339 | * | GLP_NS | - | GLP_NU | |
---|
1340 | * +-----------------------+-----------+--------------------+ |
---|
1341 | * |
---|
1342 | * Value of column q in solution to the original problem is the same as |
---|
1343 | * in solution to the transformed problem. |
---|
1344 | * |
---|
1345 | * RECOVERING INTERIOR POINT SOLUTION |
---|
1346 | * |
---|
1347 | * Value of column q in solution to the original problem is the same as |
---|
1348 | * in solution to the transformed problem. |
---|
1349 | * |
---|
1350 | * RECOVERING MIP SOLUTION |
---|
1351 | * |
---|
1352 | * None needed. */ |
---|
1353 | |
---|
1354 | struct make_fixed |
---|
1355 | { /* column with almost identical bounds */ |
---|
1356 | int q; |
---|
1357 | /* column reference number */ |
---|
1358 | double c; |
---|
1359 | /* objective coefficient at x[q] */ |
---|
1360 | NPPLFE *ptr; |
---|
1361 | /* list of non-zero coefficients a[i,q] */ |
---|
1362 | }; |
---|
1363 | |
---|
1364 | static int rcv_make_fixed(NPP *npp, void *info); |
---|
1365 | |
---|
1366 | int npp_make_fixed(NPP *npp, NPPCOL *q) |
---|
1367 | { /* process column with almost identical bounds */ |
---|
1368 | struct make_fixed *info; |
---|
1369 | NPPAIJ *aij; |
---|
1370 | NPPLFE *lfe; |
---|
1371 | double s, eps, nint; |
---|
1372 | /* the column must be double-bounded */ |
---|
1373 | xassert(q->lb != -DBL_MAX); |
---|
1374 | xassert(q->ub != +DBL_MAX); |
---|
1375 | xassert(q->lb < q->ub); |
---|
1376 | /* check column bounds */ |
---|
1377 | eps = 1e-9 + 1e-12 * fabs(q->lb); |
---|
1378 | if (q->ub - q->lb > eps) return 0; |
---|
1379 | /* column bounds are very close to each other */ |
---|
1380 | /* create transformation stack entry */ |
---|
1381 | info = npp_push_tse(npp, |
---|
1382 | rcv_make_fixed, sizeof(struct make_fixed)); |
---|
1383 | info->q = q->j; |
---|
1384 | info->c = q->coef; |
---|
1385 | info->ptr = NULL; |
---|
1386 | /* save column coefficients a[i,q] (needed for basic solution |
---|
1387 | only) */ |
---|
1388 | if (npp->sol == GLP_SOL) |
---|
1389 | { for (aij = q->ptr; aij != NULL; aij = aij->c_next) |
---|
1390 | { lfe = dmp_get_atom(npp->stack, sizeof(NPPLFE)); |
---|
1391 | lfe->ref = aij->row->i; |
---|
1392 | lfe->val = aij->val; |
---|
1393 | lfe->next = info->ptr; |
---|
1394 | info->ptr = lfe; |
---|
1395 | } |
---|
1396 | } |
---|
1397 | /* compute column fixed value */ |
---|
1398 | s = 0.5 * (q->ub + q->lb); |
---|
1399 | nint = floor(s + 0.5); |
---|
1400 | if (fabs(s - nint) <= eps) s = nint; |
---|
1401 | /* make column q fixed */ |
---|
1402 | q->lb = q->ub = s; |
---|
1403 | return 1; |
---|
1404 | } |
---|
1405 | |
---|
1406 | static int rcv_make_fixed(NPP *npp, void *_info) |
---|
1407 | { /* recover column with almost identical bounds */ |
---|
1408 | struct make_fixed *info = _info; |
---|
1409 | NPPLFE *lfe; |
---|
1410 | double lambda; |
---|
1411 | if (npp->sol == GLP_SOL) |
---|
1412 | { if (npp->c_stat[info->q] == GLP_BS) |
---|
1413 | npp->c_stat[info->q] = GLP_BS; |
---|
1414 | else if (npp->c_stat[info->q] == GLP_NS) |
---|
1415 | { /* compute multiplier for column q with formula (6) */ |
---|
1416 | lambda = info->c; |
---|
1417 | for (lfe = info->ptr; lfe != NULL; lfe = lfe->next) |
---|
1418 | lambda -= lfe->val * npp->r_pi[lfe->ref]; |
---|
1419 | /* assign status to non-basic column */ |
---|
1420 | if (lambda >= 0.0) |
---|
1421 | npp->c_stat[info->q] = GLP_NL; |
---|
1422 | else |
---|
1423 | npp->c_stat[info->q] = GLP_NU; |
---|
1424 | } |
---|
1425 | else |
---|
1426 | { npp_error(); |
---|
1427 | return 1; |
---|
1428 | } |
---|
1429 | } |
---|
1430 | return 0; |
---|
1431 | } |
---|
1432 | |
---|
1433 | /* eof */ |
---|