lemon-project-template-glpk
comparison deps/glpk/src/glpnpp02.c @ 9:33de93886c88
Import GLPK 4.47
author | Alpar Juttner <alpar@cs.elte.hu> |
---|---|
date | Sun, 06 Nov 2011 20:59:10 +0100 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:5c95e33a7511 |
---|---|
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, 2011 Andrew Makhorin, Department for Applied Informatics, | |
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved. | |
9 * E-mail: <mao@gnu.org>. | |
10 * | |
11 * GLPK is free software: you can redistribute it and/or modify it | |
12 * under the terms of the GNU General Public License as published by | |
13 * the Free Software Foundation, either version 3 of the License, or | |
14 * (at your option) any later version. | |
15 * | |
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT | |
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY | |
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public | |
19 * License for more details. | |
20 * | |
21 * You should have received a copy of the GNU General Public License | |
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>. | |
23 ***********************************************************************/ | |
24 | |
25 #include "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 */ |