Context Navigation

glpios06.c @ 10:5545663ca997

subpack-glpk

Last change on this file since 10:5545663ca997 was 9:33de93886c88, checked in by Alpar Juttner <alpar@…>, 13 years ago
Import GLPK 4.47
File size: 47.8 KB

Line
1	/* glpios06.c (MIR cut generator) */
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 "glpios.h"
26
27	#define _MIR_DEBUG 0
28
29	#define MAXAGGR 5
30	/* maximal number of rows which can be aggregated */
31
32	struct MIR
33	{ /* MIR cut generator working area */
34	/--------------------------------------------------------------/
35	/* global information valid for the root subproblem */
36	int m;
37	/* number of rows (in the root subproblem) */
38	int n;
39	/* number of columns */
40	char skip; / char skip[1+m]; */
41	/* skip[i], 1 <= i <= m, is a flag that means that row i should
42	not be used because (1) it is not suitable, or (2) because it
43	has been used in the aggregated constraint */
44	char isint; / char isint[1+m+n]; */
45	/* isint[k], 1 <= k <= m+n, is a flag that means that variable
46	x[k] is integer (otherwise, continuous) */
47	double lb; / double lb[1+m+n]; */
48	/* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means
49	that x[k] has no lower bound */
50	int vlb; / int vlb[1+m+n]; */
51	/* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable,
52	which defines variable lower bound x[k] >= lb[k] * x[k']; zero
53	means that x[k] has simple lower bound */
54	double ub; / double ub[1+m+n]; */
55	/* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means
56	that x[k] has no upper bound */
57	int vub; / int vub[1+m+n]; */
58	/* vub[k] = k', 1 <= k <= m+n, is the number of integer variable,
59	which defines variable upper bound x[k] <= ub[k] * x[k']; zero
60	means that x[k] has simple upper bound */
61	/--------------------------------------------------------------/
62	/* current (fractional) point to be separated */
63	double x; / double x[1+m+n]; */
64	/* x[k] is current value of auxiliary (1 <= k <= m) or structural
65	(m+1 <= k <= m+n) variable */
66	/--------------------------------------------------------------/
67	/* aggregated constraint sum a[k] * x[k] = b, which is a linear
68	combination of original constraints transformed to equalities
69	by introducing auxiliary variables */
70	int agg_cnt;
71	/* number of rows (original constraints) used to build aggregated
72	constraint, 1 <= agg_cnt <= MAXAGGR */
73	int agg_row; / int agg_row[1+MAXAGGR]; */
74	/* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build
75	aggregated constraint */
76	IOSVEC agg_vec; / IOSVEC agg_vec[1:m+n]; */
77	/* sparse vector of aggregated constraint coefficients, a[k] */
78	double agg_rhs;
79	/* right-hand side of the aggregated constraint, b */
80	/--------------------------------------------------------------/
81	/* bound substitution flags for modified constraint */
82	char subst; / char subst[1+m+n]; */
83	/* subst[k], 1 <= k <= m+n, is a bound substitution flag used for
84	variable x[k]:
85	'?' - x[k] is missing in modified constraint
86	'L' - x[k] = (lower bound) + x'[k]
87	'U' - x[k] = (upper bound) - x'[k] */
88	/--------------------------------------------------------------/
89	/* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0,
90	derived from aggregated constraint by substituting bounds;
91	note that due to substitution of variable bounds there may be
92	additional terms in the modified constraint */
93	IOSVEC mod_vec; / IOSVEC mod_vec[1:m+n]; */
94	/* sparse vector of modified constraint coefficients, a'[k] */
95	double mod_rhs;
96	/* right-hand side of the modified constraint, b' */
97	/--------------------------------------------------------------/
98	/* cutting plane sum alpha[k] * x[k] <= beta */
99	IOSVEC cut_vec; / IOSVEC cut_vec[1:m+n]; */
100	/* sparse vector of cutting plane coefficients, alpha[k] */
101	double cut_rhs;
102	/* right-hand size of the cutting plane, beta */
103	};
104
105	/***********************************************************************
106	* NAME
107	*
108	* ios_mir_init - initialize MIR cut generator
109	*
110	* SYNOPSIS
111	*
112	* #include "glpios.h"
113	* void ios_mir_init(glp_tree tree);
114	*
115	* DESCRIPTION
116	*
117	* The routine ios_mir_init initializes the MIR cut generator assuming
118	* that the current subproblem is the root subproblem.
119	*
120	* RETURNS
121	*
122	* The routine ios_mir_init returns a pointer to the MIR cut generator
123	* working area. */
124
125	static void set_row_attrib(glp_tree tree, struct MIR mir)
126	{ /* set global row attributes */
127	glp_prob *mip = tree->mip;
128	int m = mir->m;
129	int k;
130	for (k = 1; k <= m; k++)
131	{ GLPROW *row = mip->row[k];
132	mir->skip[k] = 0;
133	mir->isint[k] = 0;
134	switch (row->type)
135	{ case GLP_FR:
136	mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
137	case GLP_LO:
138	mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break;
139	case GLP_UP:
140	mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break;
141	case GLP_DB:
142	mir->lb[k] = row->lb, mir->ub[k] = row->ub; break;
143	case GLP_FX:
144	mir->lb[k] = mir->ub[k] = row->lb; break;
145	default:
146	xassert(row != row);
147	}
148	mir->vlb[k] = mir->vub[k] = 0;
149	}
150	return;
151	}
152
153	static void set_col_attrib(glp_tree tree, struct MIR mir)
154	{ /* set global column attributes */
155	glp_prob *mip = tree->mip;
156	int m = mir->m;
157	int n = mir->n;
158	int k;
159	for (k = m+1; k <= m+n; k++)
160	{ GLPCOL *col = mip->col[k-m];
161	switch (col->kind)
162	{ case GLP_CV:
163	mir->isint[k] = 0; break;
164	case GLP_IV:
165	mir->isint[k] = 1; break;
166	default:
167	xassert(col != col);
168	}
169	switch (col->type)
170	{ case GLP_FR:
171	mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
172	case GLP_LO:
173	mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break;
174	case GLP_UP:
175	mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break;
176	case GLP_DB:
177	mir->lb[k] = col->lb, mir->ub[k] = col->ub; break;
178	case GLP_FX:
179	mir->lb[k] = mir->ub[k] = col->lb; break;
180	default:
181	xassert(col != col);
182	}
183	mir->vlb[k] = mir->vub[k] = 0;
184	}
185	return;
186	}
187
188	static void set_var_bounds(glp_tree tree, struct MIR mir)
189	{ /* set variable bounds */
190	glp_prob *mip = tree->mip;
191	int m = mir->m;
192	GLPAIJ *aij;
193	int i, k1, k2;
194	double a1, a2;
195	for (i = 1; i <= m; i++)
196	{ /* we need the row to be '>= 0' or '<= 0' */
197	if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX \|\|
198	mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue;
199	/* take first term */
200	aij = mip->row[i]->ptr;
201	if (aij == NULL) continue;
202	k1 = m + aij->col->j, a1 = aij->val;
203	/* take second term */
204	aij = aij->r_next;
205	if (aij == NULL) continue;
206	k2 = m + aij->col->j, a2 = aij->val;
207	/* there must be only two terms */
208	if (aij->r_next != NULL) continue;
209	/* interchange terms, if needed */
210	if (!mir->isint[k1] && mir->isint[k2])
211	;
212	else if (mir->isint[k1] && !mir->isint[k2])
213	{ k2 = k1, a2 = a1;
214	k1 = m + aij->col->j, a1 = aij->val;
215	}
216	else
217	{ /* both terms are either continuous or integer */
218	continue;
219	}
220	/* x[k2] should be double-bounded */
221	if (mir->lb[k2] == -DBL_MAX \|\| mir->ub[k2] == +DBL_MAX \|\|
222	mir->lb[k2] == mir->ub[k2]) continue;
223	/* change signs, if necessary */
224	if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2;
225	/* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1
226	is continuous, x2 is integer */
227	if (a1 > 0.0)
228	{ /* x1 >= - (a2 / a1) * x2 */
229	if (mir->vlb[k1] == 0)
230	{ /* set variable lower bound for x1 */
231	mir->lb[k1] = - a2 / a1;
232	mir->vlb[k1] = k2;
233	/* the row should not be used */
234	mir->skip[i] = 1;
235	}
236	}
237	else /* a1 < 0.0 */
238	{ /* x1 <= - (a2 / a1) * x2 */
239	if (mir->vub[k1] == 0)
240	{ /* set variable upper bound for x1 */
241	mir->ub[k1] = - a2 / a1;
242	mir->vub[k1] = k2;
243	/* the row should not be used */
244	mir->skip[i] = 1;
245	}
246	}
247	}
248	return;
249	}
250
251	static void mark_useless_rows(glp_tree tree, struct MIR mir)
252	{ /* mark rows which should not be used */
253	glp_prob *mip = tree->mip;
254	int m = mir->m;
255	GLPAIJ *aij;
256	int i, k, nv;
257	for (i = 1; i <= m; i++)
258	{ /* free rows should not be used */
259	if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX)
260	{ mir->skip[i] = 1;
261	continue;
262	}
263	nv = 0;
264	for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
265	{ k = m + aij->col->j;
266	/* rows with free variables should not be used */
267	if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX)
268	{ mir->skip[i] = 1;
269	break;
270	}
271	/* rows with integer variables having infinite (lower or
272	upper) bound should not be used */
273	if (mir->isint[k] && mir->lb[k] == -DBL_MAX \|\|
274	mir->isint[k] && mir->ub[k] == +DBL_MAX)
275	{ mir->skip[i] = 1;
276	break;
277	}
278	/* count non-fixed variables */
279	if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 &&
280	mir->lb[k] == mir->ub[k])) nv++;
281	}
282	/* rows with all variables fixed should not be used */
283	if (nv == 0)
284	{ mir->skip[i] = 1;
285	continue;
286	}
287	}
288	return;
289	}
290
291	void ios_mir_init(glp_tree tree)
292	{ /* initialize MIR cut generator */
293	glp_prob *mip = tree->mip;
294	int m = mip->m;
295	int n = mip->n;
296	struct MIR *mir;
297	#if _MIR_DEBUG
298	xprintf("ios_mir_init: warning: debug mode enabled\n");
299	#endif
300	/* allocate working area */
301	mir = xmalloc(sizeof(struct MIR));
302	mir->m = m;
303	mir->n = n;
304	mir->skip = xcalloc(1+m, sizeof(char));
305	mir->isint = xcalloc(1+m+n, sizeof(char));
306	mir->lb = xcalloc(1+m+n, sizeof(double));
307	mir->vlb = xcalloc(1+m+n, sizeof(int));
308	mir->ub = xcalloc(1+m+n, sizeof(double));
309	mir->vub = xcalloc(1+m+n, sizeof(int));
310	mir->x = xcalloc(1+m+n, sizeof(double));
311	mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int));
312	mir->agg_vec = ios_create_vec(m+n);
313	mir->subst = xcalloc(1+m+n, sizeof(char));
314	mir->mod_vec = ios_create_vec(m+n);
315	mir->cut_vec = ios_create_vec(m+n);
316	/* set global row attributes */
317	set_row_attrib(tree, mir);
318	/* set global column attributes */
319	set_col_attrib(tree, mir);
320	/* set variable bounds */
321	set_var_bounds(tree, mir);
322	/* mark rows which should not be used */
323	mark_useless_rows(tree, mir);
324	return mir;
325	}
326
327	/***********************************************************************
328	* NAME
329	*
330	* ios_mir_gen - generate MIR cuts
331	*
332	* SYNOPSIS
333	*
334	* #include "glpios.h"
335	* void ios_mir_gen(glp_tree tree, void gen, IOSPOOL *pool);
336	*
337	* DESCRIPTION
338	*
339	* The routine ios_mir_gen generates MIR cuts for the current point and
340	* adds them to the cut pool. */
341
342	static void get_current_point(glp_tree tree, struct MIR mir)
343	{ /* obtain current point */
344	glp_prob *mip = tree->mip;
345	int m = mir->m;
346	int n = mir->n;
347	int k;
348	for (k = 1; k <= m; k++)
349	mir->x[k] = mip->row[k]->prim;
350	for (k = m+1; k <= m+n; k++)
351	mir->x[k] = mip->col[k-m]->prim;
352	return;
353	}
354
355	#if _MIR_DEBUG
356	static void check_current_point(struct MIR *mir)
357	{ /* check current point */
358	int m = mir->m;
359	int n = mir->n;
360	int k, kk;
361	double lb, ub, eps;
362	for (k = 1; k <= m+n; k++)
363	{ /* determine lower bound */
364	lb = mir->lb[k];
365	kk = mir->vlb[k];
366	if (kk != 0)
367	{ xassert(lb != -DBL_MAX);
368	xassert(!mir->isint[k]);
369	xassert(mir->isint[kk]);
370	lb *= mir->x[kk];
371	}
372	/* check lower bound */
373	if (lb != -DBL_MAX)
374	{ eps = 1e-6 * (1.0 + fabs(lb));
375	xassert(mir->x[k] >= lb - eps);
376	}
377	/* determine upper bound */
378	ub = mir->ub[k];
379	kk = mir->vub[k];
380	if (kk != 0)
381	{ xassert(ub != +DBL_MAX);
382	xassert(!mir->isint[k]);
383	xassert(mir->isint[kk]);
384	ub *= mir->x[kk];
385	}
386	/* check upper bound */
387	if (ub != +DBL_MAX)
388	{ eps = 1e-6 * (1.0 + fabs(ub));
389	xassert(mir->x[k] <= ub + eps);
390	}
391	}
392	return;
393	}
394	#endif
395
396	static void initial_agg_row(glp_tree tree, struct MIR mir, int i)
397	{ /* use original i-th row as initial aggregated constraint */
398	glp_prob *mip = tree->mip;
399	int m = mir->m;
400	GLPAIJ *aij;
401	xassert(1 <= i && i <= m);
402	xassert(!mir->skip[i]);
403	/* mark i-th row in order not to use it in the same aggregated
404	constraint */
405	mir->skip[i] = 2;
406	mir->agg_cnt = 1;
407	mir->agg_row[1] = i;
408	/* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary
409	variable of row i, x[m+j] are structural variables */
410	ios_clear_vec(mir->agg_vec);
411	ios_set_vj(mir->agg_vec, i, 1.0);
412	for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
413	ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val);
414	mir->agg_rhs = 0.0;
415	#if _MIR_DEBUG
416	ios_check_vec(mir->agg_vec);
417	#endif
418	return;
419	}
420
421	#if _MIR_DEBUG
422	static void check_agg_row(struct MIR *mir)
423	{ /* check aggregated constraint */
424	int m = mir->m;
425	int n = mir->n;
426	int j, k;
427	double r, big;
428	/* compute the residual r = sum a[k] * x[k] - b and determine
429	big = max(1, \|a[k]\|, \|b\|) */
430	r = 0.0, big = 1.0;
431	for (j = 1; j <= mir->agg_vec->nnz; j++)
432	{ k = mir->agg_vec->ind[j];
433	xassert(1 <= k && k <= m+n);
434	r += mir->agg_vec->val[j] * mir->x[k];
435	if (big < fabs(mir->agg_vec->val[j]))
436	big = fabs(mir->agg_vec->val[j]);
437	}
438	r -= mir->agg_rhs;
439	if (big < fabs(mir->agg_rhs))
440	big = fabs(mir->agg_rhs);
441	/* the residual must be close to zero */
442	xassert(fabs(r) <= 1e-6 * big);
443	return;
444	}
445	#endif
446
447	static void subst_fixed_vars(struct MIR *mir)
448	{ /* substitute fixed variables into aggregated constraint */
449	int m = mir->m;
450	int n = mir->n;
451	int j, k;
452	for (j = 1; j <= mir->agg_vec->nnz; j++)
453	{ k = mir->agg_vec->ind[j];
454	xassert(1 <= k && k <= m+n);
455	if (mir->vlb[k] == 0 && mir->vub[k] == 0 &&
456	mir->lb[k] == mir->ub[k])
457	{ /* x[k] is fixed */
458	mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k];
459	mir->agg_vec->val[j] = 0.0;
460	}
461	}
462	/* remove terms corresponding to fixed variables */
463	ios_clean_vec(mir->agg_vec, DBL_EPSILON);
464	#if _MIR_DEBUG
465	ios_check_vec(mir->agg_vec);
466	#endif
467	return;
468	}
469
470	static void bound_subst_heur(struct MIR *mir)
471	{ /* bound substitution heuristic */
472	int m = mir->m;
473	int n = mir->n;
474	int j, k, kk;
475	double d1, d2;
476	for (j = 1; j <= mir->agg_vec->nnz; j++)
477	{ k = mir->agg_vec->ind[j];
478	xassert(1 <= k && k <= m+n);
479	if (mir->isint[k]) continue; /* skip integer variable */
480	/* compute distance from x[k] to its lower bound */
481	kk = mir->vlb[k];
482	if (kk == 0)
483	{ if (mir->lb[k] == -DBL_MAX)
484	d1 = DBL_MAX;
485	else
486	d1 = mir->x[k] - mir->lb[k];
487	}
488	else
489	{ xassert(1 <= kk && kk <= m+n);
490	xassert(mir->isint[kk]);
491	xassert(mir->lb[k] != -DBL_MAX);
492	d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
493	}
494	/* compute distance from x[k] to its upper bound */
495	kk = mir->vub[k];
496	if (kk == 0)
497	{ if (mir->vub[k] == +DBL_MAX)
498	d2 = DBL_MAX;
499	else
500	d2 = mir->ub[k] - mir->x[k];
501	}
502	else
503	{ xassert(1 <= kk && kk <= m+n);
504	xassert(mir->isint[kk]);
505	xassert(mir->ub[k] != +DBL_MAX);
506	d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
507	}
508	/* x[k] cannot be free */
509	xassert(d1 != DBL_MAX \|\| d2 != DBL_MAX);
510	/* choose the bound which is closer to x[k] */
511	xassert(mir->subst[k] == '?');
512	if (d1 <= d2)
513	mir->subst[k] = 'L';
514	else
515	mir->subst[k] = 'U';
516	}
517	return;
518	}
519
520	static void build_mod_row(struct MIR *mir)
521	{ /* substitute bounds and build modified constraint */
522	int m = mir->m;
523	int n = mir->n;
524	int j, jj, k, kk;
525	/* initially modified constraint is aggregated constraint */
526	ios_copy_vec(mir->mod_vec, mir->agg_vec);
527	mir->mod_rhs = mir->agg_rhs;
528	#if _MIR_DEBUG
529	ios_check_vec(mir->mod_vec);
530	#endif
531	/* substitute bounds for continuous variables; note that due to
532	substitution of variable bounds additional terms may appear in
533	modified constraint */
534	for (j = mir->mod_vec->nnz; j >= 1; j--)
535	{ k = mir->mod_vec->ind[j];
536	xassert(1 <= k && k <= m+n);
537	if (mir->isint[k]) continue; /* skip integer variable */
538	if (mir->subst[k] == 'L')
539	{ /* x[k] = (lower bound) + x'[k] */
540	xassert(mir->lb[k] != -DBL_MAX);
541	kk = mir->vlb[k];
542	if (kk == 0)
543	{ /* x[k] = lb[k] + x'[k] */
544	mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
545	}
546	else
547	{ /* x[k] = lb[k] * x[kk] + x'[k] */
548	xassert(mir->isint[kk]);
549	jj = mir->mod_vec->pos[kk];
550	if (jj == 0)
551	{ ios_set_vj(mir->mod_vec, kk, 1.0);
552	jj = mir->mod_vec->pos[kk];
553	mir->mod_vec->val[jj] = 0.0;
554	}
555	mir->mod_vec->val[jj] +=
556	mir->mod_vec->val[j] * mir->lb[k];
557	}
558	}
559	else if (mir->subst[k] == 'U')
560	{ /* x[k] = (upper bound) - x'[k] */
561	xassert(mir->ub[k] != +DBL_MAX);
562	kk = mir->vub[k];
563	if (kk == 0)
564	{ /* x[k] = ub[k] - x'[k] */
565	mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
566	}
567	else
568	{ /* x[k] = ub[k] * x[kk] - x'[k] */
569	xassert(mir->isint[kk]);
570	jj = mir->mod_vec->pos[kk];
571	if (jj == 0)
572	{ ios_set_vj(mir->mod_vec, kk, 1.0);
573	jj = mir->mod_vec->pos[kk];
574	mir->mod_vec->val[jj] = 0.0;
575	}
576	mir->mod_vec->val[jj] +=
577	mir->mod_vec->val[j] * mir->ub[k];
578	}
579	mir->mod_vec->val[j] = - mir->mod_vec->val[j];
580	}
581	else
582	xassert(k != k);
583	}
584	#if _MIR_DEBUG
585	ios_check_vec(mir->mod_vec);
586	#endif
587	/* substitute bounds for integer variables */
588	for (j = 1; j <= mir->mod_vec->nnz; j++)
589	{ k = mir->mod_vec->ind[j];
590	xassert(1 <= k && k <= m+n);
591	if (!mir->isint[k]) continue; /* skip continuous variable */
592	xassert(mir->subst[k] == '?');
593	xassert(mir->vlb[k] == 0 && mir->vub[k] == 0);
594	xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX);
595	if (fabs(mir->lb[k]) <= fabs(mir->ub[k]))
596	{ /* x[k] = lb[k] + x'[k] */
597	mir->subst[k] = 'L';
598	mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
599	}
600	else
601	{ /* x[k] = ub[k] - x'[k] */
602	mir->subst[k] = 'U';
603	mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
604	mir->mod_vec->val[j] = - mir->mod_vec->val[j];
605	}
606	}
607	#if _MIR_DEBUG
608	ios_check_vec(mir->mod_vec);
609	#endif
610	return;
611	}
612
613	#if _MIR_DEBUG
614	static void check_mod_row(struct MIR *mir)
615	{ /* check modified constraint */
616	int m = mir->m;
617	int n = mir->n;
618	int j, k, kk;
619	double r, big, x;
620	/* compute the residual r = sum a'[k] * x'[k] - b' and determine
621	big = max(1, \|a[k]\|, \|b\|) */
622	r = 0.0, big = 1.0;
623	for (j = 1; j <= mir->mod_vec->nnz; j++)
624	{ k = mir->mod_vec->ind[j];
625	xassert(1 <= k && k <= m+n);
626	if (mir->subst[k] == 'L')
627	{ /* x'[k] = x[k] - (lower bound) */
628	xassert(mir->lb[k] != -DBL_MAX);
629	kk = mir->vlb[k];
630	if (kk == 0)
631	x = mir->x[k] - mir->lb[k];
632	else
633	x = mir->x[k] - mir->lb[k] * mir->x[kk];
634	}
635	else if (mir->subst[k] == 'U')
636	{ /* x'[k] = (upper bound) - x[k] */
637	xassert(mir->ub[k] != +DBL_MAX);
638	kk = mir->vub[k];
639	if (kk == 0)
640	x = mir->ub[k] - mir->x[k];
641	else
642	x = mir->ub[k] * mir->x[kk] - mir->x[k];
643	}
644	else
645	xassert(k != k);
646	r += mir->mod_vec->val[j] * x;
647	if (big < fabs(mir->mod_vec->val[j]))
648	big = fabs(mir->mod_vec->val[j]);
649	}
650	r -= mir->mod_rhs;
651	if (big < fabs(mir->mod_rhs))
652	big = fabs(mir->mod_rhs);
653	/* the residual must be close to zero */
654	xassert(fabs(r) <= 1e-6 * big);
655	return;
656	}
657	#endif
658
659	/***********************************************************************
660	* mir_ineq - construct MIR inequality
661	*
662	* Given the single constraint mixed integer set
663	*
664	* \|N\|
665	* X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s},
666	* + + j in N
667	*
668	* this routine constructs the mixed integer rounding (MIR) inequality
669	*
670	* sum alpha[j] * x[j] <= beta + gamma * s,
671	* j in N
672	*
673	* which is valid for X.
674	*
675	* If the MIR inequality has been successfully constructed, the routine
676	* returns zero. Otherwise, if b is close to nearest integer, there may
677	* be numeric difficulties due to big coefficients; so in this case the
678	* routine returns non-zero. */
679
680	static int mir_ineq(const int n, const double a[], const double b,
681	double alpha[], double beta, double gamma)
682	{ int j;
683	double f, t;
684	if (fabs(b - floor(b + .5)) < 0.01)
685	return 1;
686	f = b - floor(b);
687	for (j = 1; j <= n; j++)
688	{ t = (a[j] - floor(a[j])) - f;
689	if (t <= 0.0)
690	alpha[j] = floor(a[j]);
691	else
692	alpha[j] = floor(a[j]) + t / (1.0 - f);
693	}
694	*beta = floor(b);
695	*gamma = 1.0 / (1.0 - f);
696	return 0;
697	}
698
699	/***********************************************************************
700	* cmir_ineq - construct c-MIR inequality
701	*
702	* Given the mixed knapsack set
703	*
704	* MK \|N\|
705	* X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s,
706	* + + j in N
707	*
708	* x[j] <= u[j]},
709	*
710	* a subset C of variables to be complemented, and a divisor delta > 0,
711	* this routine constructs the complemented MIR (c-MIR) inequality
712	*
713	* sum alpha[j] * x[j] <= beta + gamma * s,
714	* j in N
715	* MK
716	* which is valid for X .
717	*
718	* If the c-MIR inequality has been successfully constructed, the
719	* routine returns zero. Otherwise, if there is a risk of numerical
720	* difficulties due to big coefficients (see comments to the routine
721	* mir_ineq), the routine cmir_ineq returns non-zero. */
722
723	static int cmir_ineq(const int n, const double a[], const double b,
724	const double u[], const char cset[], const double delta,
725	double alpha[], double beta, double gamma)
726	{ int j;
727	double *aa, bb;
728	aa = alpha, bb = b;
729	for (j = 1; j <= n; j++)
730	{ aa[j] = a[j] / delta;
731	if (cset[j])
732	aa[j] = - aa[j], bb -= a[j] * u[j];
733	}
734	bb /= delta;
735	if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1;
736	for (j = 1; j <= n; j++)
737	{ if (cset[j])
738	alpha[j] = - alpha[j], beta += alpha[j] u[j];
739	}
740	*gamma /= delta;
741	return 0;
742	}
743
744	/***********************************************************************
745	* cmir_sep - c-MIR separation heuristic
746	*
747	* Given the mixed knapsack set
748	*
749	* MK \|N\|
750	* X = {(x,s) in Z x R : sum a[j] * x[j] <= b + s,
751	* + + j in N
752	*
753	* x[j] <= u[j]}
754	*
755	* * *
756	* and a fractional point (x , s ), this routine tries to construct
757	* c-MIR inequality
758	*
759	* sum alpha[j] * x[j] <= beta + gamma * s,
760	* j in N
761	* MK
762	* which is valid for X and has (desirably maximal) violation at the
763	* fractional point given. This is attained by choosing an appropriate
764	* set C of variables to be complemented and a divisor delta > 0, which
765	* together define corresponding c-MIR inequality.
766	*
767	* If a violated c-MIR inequality has been successfully constructed,
768	* the routine returns its violation:
769	*
770	* * *
771	* sum alpha[j] * x [j] - beta - gamma * s ,
772	* j in N
773	*
774	* which is positive. In case of failure the routine returns zero. */
775
776	struct vset { int j; double v; };
777
778	static int cmir_cmp(const void p1, const void p2)
779	{ const struct vset v1 = p1, v2 = p2;
780	if (v1->v < v2->v) return -1;
781	if (v1->v > v2->v) return +1;
782	return 0;
783	}
784
785	static double cmir_sep(const int n, const double a[], const double b,
786	const double u[], const double x[], const double s,
787	double alpha[], double beta, double gamma)
788	{ int fail, j, k, nv, v;
789	double delta, eps, d_try[1+3], r, r_best;
790	char *cset;
791	struct vset *vset;
792	/* allocate working arrays */
793	cset = xcalloc(1+n, sizeof(char));
794	vset = xcalloc(1+n, sizeof(struct vset));
795	/* choose initial C */
796	for (j = 1; j <= n; j++)
797	cset[j] = (char)(x[j] >= 0.5 * u[j]);
798	/* choose initial delta */
799	r_best = delta = 0.0;
800	for (j = 1; j <= n; j++)
801	{ xassert(a[j] != 0.0);
802	/* if x[j] is close to its bounds, skip it */
803	eps = 1e-9 * (1.0 + fabs(u[j]));
804	if (x[j] < eps \|\| x[j] > u[j] - eps) continue;
805	/* try delta = \|a[j]\| to construct c-MIR inequality */
806	fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta,
807	gamma);
808	if (fail) continue;
809	/* compute violation */
810	r = - (beta) - (gamma) * s;
811	for (k = 1; k <= n; k++) r += alpha[k] * x[k];
812	if (r_best < r) r_best = r, delta = fabs(a[j]);
813	}
814	if (r_best < 0.001) r_best = 0.0;
815	if (r_best == 0.0) goto done;
816	xassert(delta > 0.0);
817	/* try to increase violation by dividing delta by 2, 4, and 8,
818	respectively */
819	d_try[1] = delta / 2.0;
820	d_try[2] = delta / 4.0;
821	d_try[3] = delta / 8.0;
822	for (j = 1; j <= 3; j++)
823	{ /* construct c-MIR inequality */
824	fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,
825	gamma);
826	if (fail) continue;
827	/* compute violation */
828	r = - (beta) - (gamma) * s;
829	for (k = 1; k <= n; k++) r += alpha[k] * x[k];
830	if (r_best < r) r_best = r, delta = d_try[j];
831	}
832	/* build subset of variables lying strictly between their bounds
833	and order it by nondecreasing values of \|x[j] - u[j]/2\| */
834	nv = 0;
835	for (j = 1; j <= n; j++)
836	{ /* if x[j] is close to its bounds, skip it */
837	eps = 1e-9 * (1.0 + fabs(u[j]));
838	if (x[j] < eps \|\| x[j] > u[j] - eps) continue;
839	/* add x[j] to the subset */
840	nv++;
841	vset[nv].j = j;
842	vset[nv].v = fabs(x[j] - 0.5 * u[j]);
843	}
844	qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp);
845	/* try to increase violation by successively complementing each
846	variable in the subset */
847	for (v = 1; v <= nv; v++)
848	{ j = vset[v].j;
849	/* replace x[j] by its complement or vice versa */
850	cset[j] = (char)!cset[j];
851	/* construct c-MIR inequality */
852	fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
853	/* restore the variable */
854	cset[j] = (char)!cset[j];
855	/* do not replace the variable in case of failure */
856	if (fail) continue;
857	/* compute violation */
858	r = - (beta) - (gamma) * s;
859	for (k = 1; k <= n; k++) r += alpha[k] * x[k];
860	if (r_best < r) r_best = r, cset[j] = (char)!cset[j];
861	}
862	/* construct the best c-MIR inequality chosen */
863	fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
864	xassert(!fail);
865	done: /* free working arrays */
866	xfree(cset);
867	xfree(vset);
868	/* return to the calling routine */
869	return r_best;
870	}
871
872	static double generate(struct MIR *mir)
873	{ /* try to generate violated c-MIR cut for modified constraint */
874	int m = mir->m;
875	int n = mir->n;
876	int j, k, kk, nint;
877	double s, u, x, *alpha, r_best = 0.0, b, beta, gamma;
878	ios_copy_vec(mir->cut_vec, mir->mod_vec);
879	mir->cut_rhs = mir->mod_rhs;
880	/* remove small terms, which can appear due to substitution of
881	variable bounds */
882	ios_clean_vec(mir->cut_vec, DBL_EPSILON);
883	#if _MIR_DEBUG
884	ios_check_vec(mir->cut_vec);
885	#endif
886	/* remove positive continuous terms to obtain MK relaxation */
887	for (j = 1; j <= mir->cut_vec->nnz; j++)
888	{ k = mir->cut_vec->ind[j];
889	xassert(1 <= k && k <= m+n);
890	if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0)
891	mir->cut_vec->val[j] = 0.0;
892	}
893	ios_clean_vec(mir->cut_vec, 0.0);
894	#if _MIR_DEBUG
895	ios_check_vec(mir->cut_vec);
896	#endif
897	/* move integer terms to the beginning of the sparse vector and
898	determine the number of integer variables */
899	nint = 0;
900	for (j = 1; j <= mir->cut_vec->nnz; j++)
901	{ k = mir->cut_vec->ind[j];
902	xassert(1 <= k && k <= m+n);
903	if (mir->isint[k])
904	{ double temp;
905	nint++;
906	/* interchange elements [nint] and [j] */
907	kk = mir->cut_vec->ind[nint];
908	mir->cut_vec->pos[k] = nint;
909	mir->cut_vec->pos[kk] = j;
910	mir->cut_vec->ind[nint] = k;
911	mir->cut_vec->ind[j] = kk;
912	temp = mir->cut_vec->val[nint];
913	mir->cut_vec->val[nint] = mir->cut_vec->val[j];
914	mir->cut_vec->val[j] = temp;
915	}
916	}
917	#if _MIR_DEBUG
918	ios_check_vec(mir->cut_vec);
919	#endif
920	/* if there is no integer variable, nothing to generate */
921	if (nint == 0) goto done;
922	/* allocate working arrays */
923	u = xcalloc(1+nint, sizeof(double));
924	x = xcalloc(1+nint, sizeof(double));
925	alpha = xcalloc(1+nint, sizeof(double));
926	/* determine u and x */
927	for (j = 1; j <= nint; j++)
928	{ k = mir->cut_vec->ind[j];
929	xassert(m+1 <= k && k <= m+n);
930	xassert(mir->isint[k]);
931	u[j] = mir->ub[k] - mir->lb[k];
932	xassert(u[j] >= 1.0);
933	if (mir->subst[k] == 'L')
934	x[j] = mir->x[k] - mir->lb[k];
935	else if (mir->subst[k] == 'U')
936	x[j] = mir->ub[k] - mir->x[k];
937	else
938	xassert(k != k);
939	xassert(x[j] >= -0.001);
940	if (x[j] < 0.0) x[j] = 0.0;
941	}
942	/* compute s = - sum of continuous terms */
943	s = 0.0;
944	for (j = nint+1; j <= mir->cut_vec->nnz; j++)
945	{ double x;
946	k = mir->cut_vec->ind[j];
947	xassert(1 <= k && k <= m+n);
948	/* must be continuous */
949	xassert(!mir->isint[k]);
950	if (mir->subst[k] == 'L')
951	{ xassert(mir->lb[k] != -DBL_MAX);
952	kk = mir->vlb[k];
953	if (kk == 0)
954	x = mir->x[k] - mir->lb[k];
955	else
956	x = mir->x[k] - mir->lb[k] * mir->x[kk];
957	}
958	else if (mir->subst[k] == 'U')
959	{ xassert(mir->ub[k] != +DBL_MAX);
960	kk = mir->vub[k];
961	if (kk == 0)
962	x = mir->ub[k] - mir->x[k];
963	else
964	x = mir->ub[k] * mir->x[kk] - mir->x[k];
965	}
966	else
967	xassert(k != k);
968	xassert(x >= -0.001);
969	if (x < 0.0) x = 0.0;
970	s -= mir->cut_vec->val[j] * x;
971	}
972	xassert(s >= 0.0);
973	/* apply heuristic to obtain most violated c-MIR inequality */
974	b = mir->cut_rhs;
975	r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,
976	&beta, &gamma);
977	if (r_best == 0.0) goto skip;
978	xassert(r_best > 0.0);
979	/* convert to raw cut */
980	/* sum alpha[j] * x[j] <= beta + gamma * s */
981	for (j = 1; j <= nint; j++)
982	mir->cut_vec->val[j] = alpha[j];
983	for (j = nint+1; j <= mir->cut_vec->nnz; j++)
984	{ k = mir->cut_vec->ind[j];
985	if (k <= m+n) mir->cut_vec->val[j] *= gamma;
986	}
987	mir->cut_rhs = beta;
988	#if _MIR_DEBUG
989	ios_check_vec(mir->cut_vec);
990	#endif
991	skip: /* free working arrays */
992	xfree(u);
993	xfree(x);
994	xfree(alpha);
995	done: return r_best;
996	}
997
998	#if _MIR_DEBUG
999	static void check_raw_cut(struct MIR *mir, double r_best)
1000	{ /* check raw cut before back bound substitution */
1001	int m = mir->m;
1002	int n = mir->n;
1003	int j, k, kk;
1004	double r, big, x;
1005	/* compute the residual r = sum a[k] * x[k] - b and determine
1006	big = max(1, \|a[k]\|, \|b\|) */
1007	r = 0.0, big = 1.0;
1008	for (j = 1; j <= mir->cut_vec->nnz; j++)
1009	{ k = mir->cut_vec->ind[j];
1010	xassert(1 <= k && k <= m+n);
1011	if (mir->subst[k] == 'L')
1012	{ xassert(mir->lb[k] != -DBL_MAX);
1013	kk = mir->vlb[k];
1014	if (kk == 0)
1015	x = mir->x[k] - mir->lb[k];
1016	else
1017	x = mir->x[k] - mir->lb[k] * mir->x[kk];
1018	}
1019	else if (mir->subst[k] == 'U')
1020	{ xassert(mir->ub[k] != +DBL_MAX);
1021	kk = mir->vub[k];
1022	if (kk == 0)
1023	x = mir->ub[k] - mir->x[k];
1024	else
1025	x = mir->ub[k] * mir->x[kk] - mir->x[k];
1026	}
1027	else
1028	xassert(k != k);
1029	r += mir->cut_vec->val[j] * x;
1030	if (big < fabs(mir->cut_vec->val[j]))
1031	big = fabs(mir->cut_vec->val[j]);
1032	}
1033	r -= mir->cut_rhs;
1034	if (big < fabs(mir->cut_rhs))
1035	big = fabs(mir->cut_rhs);
1036	/* the residual must be close to r_best */
1037	xassert(fabs(r - r_best) <= 1e-6 * big);
1038	return;
1039	}
1040	#endif
1041
1042	static void back_subst(struct MIR *mir)
1043	{ /* back substitution of original bounds */
1044	int m = mir->m;
1045	int n = mir->n;
1046	int j, jj, k, kk;
1047	/* at first, restore bounds of integer variables (because on
1048	restoring variable bounds of continuous variables we need
1049	original, not shifted, bounds of integer variables) */
1050	for (j = 1; j <= mir->cut_vec->nnz; j++)
1051	{ k = mir->cut_vec->ind[j];
1052	xassert(1 <= k && k <= m+n);
1053	if (!mir->isint[k]) continue; /* skip continuous */
1054	if (mir->subst[k] == 'L')
1055	{ /* x'[k] = x[k] - lb[k] */
1056	xassert(mir->lb[k] != -DBL_MAX);
1057	xassert(mir->vlb[k] == 0);
1058	mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
1059	}
1060	else if (mir->subst[k] == 'U')
1061	{ /* x'[k] = ub[k] - x[k] */
1062	xassert(mir->ub[k] != +DBL_MAX);
1063	xassert(mir->vub[k] == 0);
1064	mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
1065	mir->cut_vec->val[j] = - mir->cut_vec->val[j];
1066	}
1067	else
1068	xassert(k != k);
1069	}
1070	/* now restore bounds of continuous variables */
1071	for (j = 1; j <= mir->cut_vec->nnz; j++)
1072	{ k = mir->cut_vec->ind[j];
1073	xassert(1 <= k && k <= m+n);
1074	if (mir->isint[k]) continue; /* skip integer */
1075	if (mir->subst[k] == 'L')
1076	{ /* x'[k] = x[k] - (lower bound) */
1077	xassert(mir->lb[k] != -DBL_MAX);
1078	kk = mir->vlb[k];
1079	if (kk == 0)
1080	{ /* x'[k] = x[k] - lb[k] */
1081	mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
1082	}
1083	else
1084	{ /* x'[k] = x[k] - lb[k] * x[kk] */
1085	jj = mir->cut_vec->pos[kk];
1086	#if 0
1087	xassert(jj != 0);
1088	#else
1089	if (jj == 0)
1090	{ ios_set_vj(mir->cut_vec, kk, 1.0);
1091	jj = mir->cut_vec->pos[kk];
1092	xassert(jj != 0);
1093	mir->cut_vec->val[jj] = 0.0;
1094	}
1095	#endif
1096	mir->cut_vec->val[jj] -= mir->cut_vec->val[j] *
1097	mir->lb[k];
1098	}
1099	}
1100	else if (mir->subst[k] == 'U')
1101	{ /* x'[k] = (upper bound) - x[k] */
1102	xassert(mir->ub[k] != +DBL_MAX);
1103	kk = mir->vub[k];
1104	if (kk == 0)
1105	{ /* x'[k] = ub[k] - x[k] */
1106	mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
1107	}
1108	else
1109	{ /* x'[k] = ub[k] * x[kk] - x[k] */
1110	jj = mir->cut_vec->pos[kk];
1111	if (jj == 0)
1112	{ ios_set_vj(mir->cut_vec, kk, 1.0);
1113	jj = mir->cut_vec->pos[kk];
1114	xassert(jj != 0);
1115	mir->cut_vec->val[jj] = 0.0;
1116	}
1117	mir->cut_vec->val[jj] += mir->cut_vec->val[j] *
1118	mir->ub[k];
1119	}
1120	mir->cut_vec->val[j] = - mir->cut_vec->val[j];
1121	}
1122	else
1123	xassert(k != k);
1124	}
1125	#if _MIR_DEBUG
1126	ios_check_vec(mir->cut_vec);
1127	#endif
1128	return;
1129	}
1130
1131	#if _MIR_DEBUG
1132	static void check_cut_row(struct MIR *mir, double r_best)
1133	{ /* check the cut after back bound substitution or elimination of
1134	auxiliary variables */
1135	int m = mir->m;
1136	int n = mir->n;
1137	int j, k;
1138	double r, big;
1139	/* compute the residual r = sum a[k] * x[k] - b and determine
1140	big = max(1, \|a[k]\|, \|b\|) */
1141	r = 0.0, big = 1.0;
1142	for (j = 1; j <= mir->cut_vec->nnz; j++)
1143	{ k = mir->cut_vec->ind[j];
1144	xassert(1 <= k && k <= m+n);
1145	r += mir->cut_vec->val[j] * mir->x[k];
1146	if (big < fabs(mir->cut_vec->val[j]))
1147	big = fabs(mir->cut_vec->val[j]);
1148	}
1149	r -= mir->cut_rhs;
1150	if (big < fabs(mir->cut_rhs))
1151	big = fabs(mir->cut_rhs);
1152	/* the residual must be close to r_best */
1153	xassert(fabs(r - r_best) <= 1e-6 * big);
1154	return;
1155	}
1156	#endif
1157
1158	static void subst_aux_vars(glp_tree tree, struct MIR mir)
1159	{ /* final substitution to eliminate auxiliary variables */
1160	glp_prob *mip = tree->mip;
1161	int m = mir->m;
1162	int n = mir->n;
1163	GLPAIJ *aij;
1164	int j, k, kk, jj;
1165	for (j = mir->cut_vec->nnz; j >= 1; j--)
1166	{ k = mir->cut_vec->ind[j];
1167	xassert(1 <= k && k <= m+n);
1168	if (k > m) continue; /* skip structurals */
1169	for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next)
1170	{ kk = m + aij->col->j; /* structural */
1171	jj = mir->cut_vec->pos[kk];
1172	if (jj == 0)
1173	{ ios_set_vj(mir->cut_vec, kk, 1.0);
1174	jj = mir->cut_vec->pos[kk];
1175	mir->cut_vec->val[jj] = 0.0;
1176	}
1177	mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val;
1178	}
1179	mir->cut_vec->val[j] = 0.0;
1180	}
1181	ios_clean_vec(mir->cut_vec, 0.0);
1182	return;
1183	}
1184
1185	static void add_cut(glp_tree tree, struct MIR mir)
1186	{ /* add constructed cut inequality to the cut pool */
1187	int m = mir->m;
1188	int n = mir->n;
1189	int j, k, len;
1190	int *ind = xcalloc(1+n, sizeof(int));
1191	double *val = xcalloc(1+n, sizeof(double));
1192	len = 0;
1193	for (j = mir->cut_vec->nnz; j >= 1; j--)
1194	{ k = mir->cut_vec->ind[j];
1195	xassert(m+1 <= k && k <= m+n);
1196	len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j];
1197	}
1198	#if 0
1199	ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,
1200	mir->cut_rhs);
1201	#else
1202	glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
1203	mir->cut_rhs);
1204	#endif
1205	xfree(ind);
1206	xfree(val);
1207	return;
1208	}
1209
1210	static int aggregate_row(glp_tree tree, struct MIR mir)
1211	{ /* try to aggregate another row */
1212	glp_prob *mip = tree->mip;
1213	int m = mir->m;
1214	int n = mir->n;
1215	GLPAIJ *aij;
1216	IOSVEC *v;
1217	int ii, j, jj, k, kk, kappa = 0, ret = 0;
1218	double d1, d2, d, d_max = 0.0;
1219	/* choose appropriate structural variable in the aggregated row
1220	to be substituted */
1221	for (j = 1; j <= mir->agg_vec->nnz; j++)
1222	{ k = mir->agg_vec->ind[j];
1223	xassert(1 <= k && k <= m+n);
1224	if (k <= m) continue; /* skip auxiliary var */
1225	if (mir->isint[k]) continue; /* skip integer var */
1226	if (fabs(mir->agg_vec->val[j]) < 0.001) continue;
1227	/* compute distance from x[k] to its lower bound */
1228	kk = mir->vlb[k];
1229	if (kk == 0)
1230	{ if (mir->lb[k] == -DBL_MAX)
1231	d1 = DBL_MAX;
1232	else
1233	d1 = mir->x[k] - mir->lb[k];
1234	}
1235	else
1236	{ xassert(1 <= kk && kk <= m+n);
1237	xassert(mir->isint[kk]);
1238	xassert(mir->lb[k] != -DBL_MAX);
1239	d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
1240	}
1241	/* compute distance from x[k] to its upper bound */
1242	kk = mir->vub[k];
1243	if (kk == 0)
1244	{ if (mir->vub[k] == +DBL_MAX)
1245	d2 = DBL_MAX;
1246	else
1247	d2 = mir->ub[k] - mir->x[k];
1248	}
1249	else
1250	{ xassert(1 <= kk && kk <= m+n);
1251	xassert(mir->isint[kk]);
1252	xassert(mir->ub[k] != +DBL_MAX);
1253	d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
1254	}
1255	/* x[k] cannot be free */
1256	xassert(d1 != DBL_MAX \|\| d2 != DBL_MAX);
1257	/* d = min(d1, d2) */
1258	d = (d1 <= d2 ? d1 : d2);
1259	xassert(d != DBL_MAX);
1260	/* should not be close to corresponding bound */
1261	if (d < 0.001) continue;
1262	if (d_max < d) d_max = d, kappa = k;
1263	}
1264	if (kappa == 0)
1265	{ /* nothing chosen */
1266	ret = 1;
1267	goto done;
1268	}
1269	/* x[kappa] has been chosen */
1270	xassert(m+1 <= kappa && kappa <= m+n);
1271	xassert(!mir->isint[kappa]);
1272	/* find another row, which have not been used yet, to eliminate
1273	x[kappa] from the aggregated row */
1274	for (ii = 1; ii <= m; ii++)
1275	{ if (mir->skip[ii]) continue;
1276	for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
1277	if (aij->col->j == kappa - m) break;
1278	if (aij != NULL && fabs(aij->val) >= 0.001) break;
1279	}
1280	if (ii > m)
1281	{ /* nothing found */
1282	ret = 2;
1283	goto done;
1284	}
1285	/* row ii has been found; include it in the aggregated list */
1286	mir->agg_cnt++;
1287	xassert(mir->agg_cnt <= MAXAGGR);
1288	mir->agg_row[mir->agg_cnt] = ii;
1289	mir->skip[ii] = 2;
1290	/* v := new row */
1291	v = ios_create_vec(m+n);
1292	ios_set_vj(v, ii, 1.0);
1293	for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
1294	ios_set_vj(v, m + aij->col->j, - aij->val);
1295	#if _MIR_DEBUG
1296	ios_check_vec(v);
1297	#endif
1298	/* perform gaussian elimination to remove x[kappa] */
1299	j = mir->agg_vec->pos[kappa];
1300	xassert(j != 0);
1301	jj = v->pos[kappa];
1302	xassert(jj != 0);
1303	ios_linear_comb(mir->agg_vec,
1304	- mir->agg_vec->val[j] / v->val[jj], v);
1305	ios_delete_vec(v);
1306	ios_set_vj(mir->agg_vec, kappa, 0.0);
1307	#if _MIR_DEBUG
1308	ios_check_vec(mir->agg_vec);
1309	#endif
1310	done: return ret;
1311	}
1312
1313	void ios_mir_gen(glp_tree tree, void gen)
1314	{ /* main routine to generate MIR cuts */
1315	glp_prob *mip = tree->mip;
1316	struct MIR *mir = gen;
1317	int m = mir->m;
1318	int n = mir->n;
1319	int i;
1320	double r_best;
1321	xassert(mip->m >= m);
1322	xassert(mip->n == n);
1323	/* obtain current point */
1324	get_current_point(tree, mir);
1325	#if _MIR_DEBUG
1326	/* check current point */
1327	check_current_point(mir);
1328	#endif
1329	/* reset bound substitution flags */
1330	memset(&mir->subst[1], '?', m+n);
1331	/* try to generate a set of violated MIR cuts */
1332	for (i = 1; i <= m; i++)
1333	{ if (mir->skip[i]) continue;
1334	/* use original i-th row as initial aggregated constraint */
1335	initial_agg_row(tree, mir, i);
1336	loop: ;
1337	#if _MIR_DEBUG
1338	/* check aggregated row */
1339	check_agg_row(mir);
1340	#endif
1341	/* substitute fixed variables into aggregated constraint */
1342	subst_fixed_vars(mir);
1343	#if _MIR_DEBUG
1344	/* check aggregated row */
1345	check_agg_row(mir);
1346	#endif
1347	#if _MIR_DEBUG
1348	/* check bound substitution flags */
1349	{ int k;
1350	for (k = 1; k <= m+n; k++)
1351	xassert(mir->subst[k] == '?');
1352	}
1353	#endif
1354	/* apply bound substitution heuristic */
1355	bound_subst_heur(mir);
1356	/* substitute bounds and build modified constraint */
1357	build_mod_row(mir);
1358	#if _MIR_DEBUG
1359	/* check modified row */
1360	check_mod_row(mir);
1361	#endif
1362	/* try to generate violated c-MIR cut for modified row */
1363	r_best = generate(mir);
1364	if (r_best > 0.0)
1365	{ /* success */
1366	#if _MIR_DEBUG
1367	/* check raw cut before back bound substitution */
1368	check_raw_cut(mir, r_best);
1369	#endif
1370	/* back substitution of original bounds */
1371	back_subst(mir);
1372	#if _MIR_DEBUG
1373	/* check the cut after back bound substitution */
1374	check_cut_row(mir, r_best);
1375	#endif
1376	/* final substitution to eliminate auxiliary variables */
1377	subst_aux_vars(tree, mir);
1378	#if _MIR_DEBUG
1379	/* check the cut after elimination of auxiliaries */
1380	check_cut_row(mir, r_best);
1381	#endif
1382	/* add constructed cut inequality to the cut pool */
1383	add_cut(tree, mir);
1384	}
1385	/* reset bound substitution flags */
1386	{ int j, k;
1387	for (j = 1; j <= mir->mod_vec->nnz; j++)
1388	{ k = mir->mod_vec->ind[j];
1389	xassert(1 <= k && k <= m+n);
1390	xassert(mir->subst[k] != '?');
1391	mir->subst[k] = '?';
1392	}
1393	}
1394	if (r_best == 0.0)
1395	{ /* failure */
1396	if (mir->agg_cnt < MAXAGGR)
1397	{ /* try to aggregate another row */
1398	if (aggregate_row(tree, mir) == 0) goto loop;
1399	}
1400	}
1401	/* unmark rows used in the aggregated constraint */
1402	{ int k, ii;
1403	for (k = 1; k <= mir->agg_cnt; k++)
1404	{ ii = mir->agg_row[k];
1405	xassert(1 <= ii && ii <= m);
1406	xassert(mir->skip[ii] == 2);
1407	mir->skip[ii] = 0;
1408	}
1409	}
1410	}
1411	return;
1412	}
1413
1414	/***********************************************************************
1415	* NAME
1416	*
1417	* ios_mir_term - terminate MIR cut generator
1418	*
1419	* SYNOPSIS
1420	*
1421	* #include "glpios.h"
1422	* void ios_mir_term(void *gen);
1423	*
1424	* DESCRIPTION
1425	*
1426	* The routine ios_mir_term deletes the MIR cut generator working area
1427	* freeing all the memory allocated to it. */
1428
1429	void ios_mir_term(void *gen)
1430	{ struct MIR *mir = gen;
1431	xfree(mir->skip);
1432	xfree(mir->isint);
1433	xfree(mir->lb);
1434	xfree(mir->vlb);
1435	xfree(mir->ub);
1436	xfree(mir->vub);
1437	xfree(mir->x);
1438	xfree(mir->agg_row);
1439	ios_delete_vec(mir->agg_vec);
1440	xfree(mir->subst);
1441	ios_delete_vec(mir->mod_vec);
1442	ios_delete_vec(mir->cut_vec);
1443	xfree(mir);
1444	return;
1445	}
1446
1447	/* eof */

Note: See TracBrowser for help on using the repository browser.

Context Navigation

source: lemon-project-template-glpk/deps/glpk/src/glpios06.c @ 10:5545663ca997

Download in other formats: