|
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 */ |