lemon-project-template-glpk

view 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
line source
1 /* glpnpp02.c */
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 ***********************************************************************/
25 #include "glpnpp.h"
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. */
69 struct free_row
70 { /* free (unbounded) row */
71 int p;
72 /* row reference number */
73 };
75 static int rcv_free_row(NPP *npp, void *info);
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 }
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 }
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. */
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 };
191 static int rcv_geq_row(NPP *npp, void *info);
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 }
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 }
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. */
335 static int rcv_leq_row(NPP *npp, void *info);
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 }
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 }
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). */
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 };
486 static int rcv_free_col(NPP *npp, void *info);
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 }
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 }
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). */
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 };
644 static int rcv_lbnd_col(NPP *npp, void *info);
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 }
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 }
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). */
812 static int rcv_ubnd_col(NPP *npp, void *info);
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 }
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 }
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. */
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 };
960 static int rcv_dbnd_col(NPP *npp, void *info);
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 }
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;
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;
1013 else
1014 { npp_error();
1015 return 1;
1018 return 0;
1021 /***********************************************************************
1022 * NAME
1024 * npp_fixed_col - process fixed column
1026 * SYNOPSIS
1028 * #include "glpnpp.h"
1029 * void npp_fixed_col(NPP *npp, NPPCOL *q);
1031 * DESCRIPTION
1033 * The routine npp_fixed_col processes column q, which is fixed:
1035 * x[q] = s[q], (1)
1037 * where s[q] is a fixed column value.
1039 * PROBLEM TRANSFORMATION
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.
1044 * Substituting x[q] = s[q] into the objective row, we have:
1046 * z = sum c[j] x[j] + c0 =
1047 * j
1049 * = sum c[j] x[j] + c[q] x[q] + c0 =
1050 * j!=q
1052 * = sum c[j] x[j] + c[q] s[q] + c0 =
1053 * j!=q
1055 * = sum c[j] x[j] + c~0,
1056 * j!=q
1058 * where
1060 * c~0 = c0 + c[q] s[q] (2)
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:
1065 * L[i] <= sum a[i,j] x[j] <= U[i] ==>
1066 * j
1068 * L[i] <= sum a[i,j] x[j] + a[i,q] x[q] <= U[i] ==>
1069 * j!=q
1071 * L[i] <= sum a[i,j] x[j] + a[i,q] s[q] <= U[i] ==>
1072 * j!=q
1074 * L~[i] <= sum a[i,j] x[j] + a[i,q] s <= U~[i],
1075 * j!=q
1077 * where
1079 * L~[i] = L[i] - a[i,q] s[q], U~[i] = U[i] - a[i,q] s[q] (3)
1081 * are lower and upper bounds of row i in the transformed problem,
1082 * resp.
1084 * RECOVERING BASIC SOLUTION
1086 * Column q is assigned status GLP_NS and its value is assigned s[q].
1088 * RECOVERING INTERIOR-POINT SOLUTION
1090 * Value of column q is assigned s[q].
1092 * RECOVERING MIP SOLUTION
1094 * Value of column q is assigned s[q]. */
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 };
1104 static int rcv_fixed_col(NPP *npp, void *info);
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;
1132 /* remove the column from the problem */
1133 npp_del_col(npp, q);
1134 return;
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;
1146 /***********************************************************************
1147 * NAME
1149 * npp_make_equality - process row with almost identical bounds
1151 * SYNOPSIS
1153 * #include "glpnpp.h"
1154 * int npp_make_equality(NPP *npp, NPPROW *p);
1156 * DESCRIPTION
1158 * The routine npp_make_equality processes row p:
1160 * L[p] <= sum a[p,j] x[j] <= U[p], (1)
1161 * j
1163 * where -oo < L[p] < U[p] < +oo, i.e. which is double-sided inequality
1164 * constraint.
1166 * RETURNS
1168 * 0 - row bounds have not been changed;
1170 * 1 - row has been replaced by equality constraint.
1172 * PROBLEM TRANSFORMATION
1174 * If bounds of row (1) are very close to each other:
1176 * U[p] - L[p] <= eps, (2)
1178 * where eps is an absolute tolerance for row value, the row can be
1179 * replaced by the following almost equivalent equiality constraint:
1181 * sum a[p,j] x[j] = b, (3)
1182 * j
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:
1187 * |b - floor(b + 0.5)| <= eps, (4)
1189 * it is reasonable to use this nearest integer as the right-hand side.
1191 * RECOVERING BASIC SOLUTION
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:
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 * +-----------------------+---------+--------------------+
1206 * Value of row multiplier pi[p] in solution to the original problem is
1207 * the same as in solution to the transformed problem.
1209 * RECOVERING INTERIOR POINT SOLUTION
1211 * Value of row multiplier pi[p] in solution to the original problem is
1212 * the same as in solution to the transformed problem.
1214 * RECOVERING MIP SOLUTION
1216 * None needed. */
1218 struct make_equality
1219 { /* row with almost identical bounds */
1220 int p;
1221 /* row reference number */
1222 };
1224 static int rcv_make_equality(NPP *npp, void *info);
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;
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;
1263 else
1264 { npp_error();
1265 return 1;
1268 return 0;
1271 /***********************************************************************
1272 * NAME
1274 * npp_make_fixed - process column with almost identical bounds
1276 * SYNOPSIS
1278 * #include "glpnpp.h"
1279 * int npp_make_fixed(NPP *npp, NPPCOL *q);
1281 * DESCRIPTION
1283 * The routine npp_make_fixed processes column q:
1285 * l[q] <= x[q] <= u[q], (1)
1287 * where -oo < l[q] < u[q] < +oo, i.e. which has both lower and upper
1288 * bounds.
1290 * RETURNS
1292 * 0 - column bounds have not been changed;
1294 * 1 - column has been fixed.
1296 * PROBLEM TRANSFORMATION
1298 * If bounds of column (1) are very close to each other:
1300 * u[q] - l[q] <= eps, (2)
1302 * where eps is an absolute tolerance for column value, the column can
1303 * be fixed:
1305 * x[q] = s[q], (3)
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:
1310 * |s[q] - floor(s[q] + 0.5)| <= eps, (4)
1312 * it is reasonable to use this nearest integer as the fixed value.
1314 * RECOVERING BASIC SOLUTION
1316 * In the dual system of the original (as well as transformed) problem
1317 * column q corresponds to the following row:
1319 * sum a[i,q] pi[i] + lambda[q] = c[q]. (5)
1320 * i
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:
1326 * lambda[q] = c[q] - sum a[i,q] pi[i]. (6)
1327 * i
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:
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 * +-----------------------+-----------+--------------------+
1342 * Value of column q in solution to the original problem is the same as
1343 * in solution to the transformed problem.
1345 * RECOVERING INTERIOR POINT SOLUTION
1347 * Value of column q in solution to the original problem is the same as
1348 * in solution to the transformed problem.
1350 * RECOVERING MIP SOLUTION
1352 * None needed. */
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 };
1364 static int rcv_make_fixed(NPP *npp, void *info);
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;
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;
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;
1425 else
1426 { npp_error();
1427 return 1;
1430 return 0;
1433 /* eof */