Context Navigation

glpios06.c @ 2:4c8956a7bdf4

Last change on this file since 2:4c8956a7bdf4 was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 14 years ago

Import glpk-4.45

Generated files and doc/notes are removed

File size: 47.8 KB

Rev	Line
[1]	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 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: glpk-cmake/src/glpios06.c @ 2:4c8956a7bdf4

Download in other formats: