Context Navigation

glpios03.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f 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: 42.2 KB

Rev	Line
[1]	1	/* glpios03.c (branch-and-cut driver) */
	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	/***********************************************************************
	28	* show_progress - display current progress of the search
	29	*
	30	* This routine displays some information about current progress of the
	31	* search.
	32	*
	33	* The information includes:
	34	*
	35	* the current number of iterations performed by the simplex solver;
	36	*
	37	* the objective value for the best known integer feasible solution,
	38	* which is upper (minimization) or lower (maximization) global bound
	39	* for optimal solution of the original mip problem;
	40	*
	41	* the best local bound for active nodes, which is lower (minimization)
	42	* or upper (maximization) global bound for optimal solution of the
	43	* original mip problem;
	44	*
	45	* the relative mip gap, in percents;
	46	*
	47	* the number of open (active) subproblems;
	48	*
	49	* the number of completely explored subproblems, i.e. whose nodes have
	50	* been removed from the tree. */
	51
	52	static void show_progress(glp_tree *T, int bingo)
	53	{ int p;
	54	double temp;
	55	char best_mip[50], best_bound[50], *rho, rel_gap[50];
	56	/* format the best known integer feasible solution */
	57	if (T->mip->mip_stat == GLP_FEAS)
	58	sprintf(best_mip, "%17.9e", T->mip->mip_obj);
	59	else
	60	sprintf(best_mip, "%17s", "not found yet");
	61	/* determine reference number of an active subproblem whose local
	62	bound is best */
	63	p = ios_best_node(T);
	64	/* format the best bound */
	65	if (p == 0)
	66	sprintf(best_bound, "%17s", "tree is empty");
	67	else
	68	{ temp = T->slot[p].node->bound;
	69	if (temp == -DBL_MAX)
	70	sprintf(best_bound, "%17s", "-inf");
	71	else if (temp == +DBL_MAX)
	72	sprintf(best_bound, "%17s", "+inf");
	73	else
	74	sprintf(best_bound, "%17.9e", temp);
	75	}
	76	/* choose the relation sign between global bounds */
	77	if (T->mip->dir == GLP_MIN)
	78	rho = ">=";
	79	else if (T->mip->dir == GLP_MAX)
	80	rho = "<=";
	81	else
	82	xassert(T != T);
	83	/* format the relative mip gap */
	84	temp = ios_relative_gap(T);
	85	if (temp == 0.0)
	86	sprintf(rel_gap, " 0.0%%");
	87	else if (temp < 0.001)
	88	sprintf(rel_gap, "< 0.1%%");
	89	else if (temp <= 9.999)
	90	sprintf(rel_gap, "%5.1f%%", 100.0 * temp);
	91	else
	92	sprintf(rel_gap, "%6s", "");
	93	/* display progress of the search */
	94	xprintf("+%6d: %s %s %s %s %s (%d; %d)\n",
	95	T->mip->it_cnt, bingo ? ">>>>>" : "mip =", best_mip, rho,
	96	best_bound, rel_gap, T->a_cnt, T->t_cnt - T->n_cnt);
	97	T->tm_lag = xtime();
	98	return;
	99	}
	100
	101	/***********************************************************************
	102	* is_branch_hopeful - check if specified branch is hopeful
	103	*
	104	* This routine checks if the specified subproblem can have an integer
	105	* optimal solution which is better than the best known one.
	106	*
	107	* The check is based on comparison of the local objective bound stored
	108	* in the subproblem descriptor and the incumbent objective value which
	109	* is the global objective bound.
	110	*
	111	* If there is a chance that the specified subproblem can have a better
	112	* integer optimal solution, the routine returns non-zero. Otherwise, if
	113	* the corresponding branch can pruned, zero is returned. */
	114
	115	static int is_branch_hopeful(glp_tree *T, int p)
	116	{ xassert(1 <= p && p <= T->nslots);
	117	xassert(T->slot[p].node != NULL);
	118	return ios_is_hopeful(T, T->slot[p].node->bound);
	119	}
	120
	121	/***********************************************************************
	122	* check_integrality - check integrality of basic solution
	123	*
	124	* This routine checks if the basic solution of LP relaxation of the
	125	* current subproblem satisfies to integrality conditions, i.e. that all
	126	* variables of integer kind have integral primal values. (The solution
	127	* is assumed to be optimal.)
	128	*
	129	* For each variable of integer kind the routine computes the following
	130	* quantity:
	131	*
	132	* ii(x[j]) = min(x[j] - floor(x[j]), ceil(x[j]) - x[j]), (1)
	133	*
	134	* which is a measure of the integer infeasibility (non-integrality) of
	135	* x[j] (for example, ii(2.1) = 0.1, ii(3.7) = 0.3, ii(5.0) = 0). It is
	136	* understood that 0 <= ii(x[j]) <= 0.5, and variable x[j] is integer
	137	* feasible if ii(x[j]) = 0. However, due to floating-point arithmetic
	138	* the routine checks less restrictive condition:
	139	*
	140	* ii(x[j]) <= tol_int, (2)
	141	*
	142	* where tol_int is a given tolerance (small positive number) and marks
	143	* each variable which does not satisfy to (2) as integer infeasible by
	144	* setting its fractionality flag.
	145	*
	146	* In order to characterize integer infeasibility of the basic solution
	147	* in the whole the routine computes two parameters: ii_cnt, which is
	148	* the number of variables with the fractionality flag set, and ii_sum,
	149	* which is the sum of integer infeasibilities (1). */
	150
	151	static void check_integrality(glp_tree *T)
	152	{ glp_prob *mip = T->mip;
	153	int j, type, ii_cnt = 0;
	154	double lb, ub, x, temp1, temp2, ii_sum = 0.0;
	155	/* walk through the set of columns (structural variables) */
	156	for (j = 1; j <= mip->n; j++)
	157	{ GLPCOL *col = mip->col[j];
	158	T->non_int[j] = 0;
	159	/* if the column is not integer, skip it */
	160	if (col->kind != GLP_IV) continue;
	161	/* if the column is non-basic, it is integer feasible */
	162	if (col->stat != GLP_BS) continue;
	163	/* obtain the type and bounds of the column */
	164	type = col->type, lb = col->lb, ub = col->ub;
	165	/* obtain value of the column in optimal basic solution */
	166	x = col->prim;
	167	/* if the column's primal value is close to the lower bound,
	168	the column is integer feasible within given tolerance */
	169	if (type == GLP_LO \|\| type == GLP_DB \|\| type == GLP_FX)
	170	{ temp1 = lb - T->parm->tol_int;
	171	temp2 = lb + T->parm->tol_int;
	172	if (temp1 <= x && x <= temp2) continue;
	173	#if 0
	174	/* the lower bound must not be violated */
	175	xassert(x >= lb);
	176	#else
	177	if (x < lb) continue;
	178	#endif
	179	}
	180	/* if the column's primal value is close to the upper bound,
	181	the column is integer feasible within given tolerance */
	182	if (type == GLP_UP \|\| type == GLP_DB \|\| type == GLP_FX)
	183	{ temp1 = ub - T->parm->tol_int;
	184	temp2 = ub + T->parm->tol_int;
	185	if (temp1 <= x && x <= temp2) continue;
	186	#if 0
	187	/* the upper bound must not be violated */
	188	xassert(x <= ub);
	189	#else
	190	if (x > ub) continue;
	191	#endif
	192	}
	193	/* if the column's primal value is close to nearest integer,
	194	the column is integer feasible within given tolerance */
	195	temp1 = floor(x + 0.5) - T->parm->tol_int;
	196	temp2 = floor(x + 0.5) + T->parm->tol_int;
	197	if (temp1 <= x && x <= temp2) continue;
	198	/* otherwise the column is integer infeasible */
	199	T->non_int[j] = 1;
	200	/* increase the number of fractional-valued columns */
	201	ii_cnt++;
	202	/* compute the sum of integer infeasibilities */
	203	temp1 = x - floor(x);
	204	temp2 = ceil(x) - x;
	205	xassert(temp1 > 0.0 && temp2 > 0.0);
	206	ii_sum += (temp1 <= temp2 ? temp1 : temp2);
	207	}
	208	/* store ii_cnt and ii_sum to the current problem descriptor */
	209	xassert(T->curr != NULL);
	210	T->curr->ii_cnt = ii_cnt;
	211	T->curr->ii_sum = ii_sum;
	212	/* and also display these parameters */
	213	if (T->parm->msg_lev >= GLP_MSG_DBG)
	214	{ if (ii_cnt == 0)
	215	xprintf("There are no fractional columns\n");
	216	else if (ii_cnt == 1)
	217	xprintf("There is one fractional column, integer infeasibil"
	218	"ity is %.3e\n", ii_sum);
	219	else
	220	xprintf("There are %d fractional columns, integer infeasibi"
	221	"lity is %.3e\n", ii_cnt, ii_sum);
	222	}
	223	return;
	224	}
	225
	226	/***********************************************************************
	227	* record_solution - record better integer feasible solution
	228	*
	229	* This routine records optimal basic solution of LP relaxation of the
	230	* current subproblem, which being integer feasible is better than the
	231	* best known integer feasible solution. */
	232
	233	static void record_solution(glp_tree *T)
	234	{ glp_prob *mip = T->mip;
	235	int i, j;
	236	mip->mip_stat = GLP_FEAS;
	237	mip->mip_obj = mip->obj_val;
	238	for (i = 1; i <= mip->m; i++)
	239	{ GLPROW *row = mip->row[i];
	240	row->mipx = row->prim;
	241	}
	242	for (j = 1; j <= mip->n; j++)
	243	{ GLPCOL *col = mip->col[j];
	244	if (col->kind == GLP_CV)
	245	col->mipx = col->prim;
	246	else if (col->kind == GLP_IV)
	247	{ /* value of the integer column must be integral */
	248	col->mipx = floor(col->prim + 0.5);
	249	}
	250	else
	251	xassert(col != col);
	252	}
	253	T->sol_cnt++;
	254	return;
	255	}
	256
	257	/***********************************************************************
	258	* fix_by_red_cost - fix non-basic integer columns by reduced costs
	259	*
	260	* This routine fixes some non-basic integer columns if their reduced
	261	* costs indicate that increasing (decreasing) the column at least by
	262	* one involves the objective value becoming worse than the incumbent
	263	* objective value. */
	264
	265	static void fix_by_red_cost(glp_tree *T)
	266	{ glp_prob *mip = T->mip;
	267	int j, stat, fixed = 0;
	268	double obj, lb, ub, dj;
	269	/* the global bound must exist */
	270	xassert(T->mip->mip_stat == GLP_FEAS);
	271	/* basic solution of LP relaxation must be optimal */
	272	xassert(mip->pbs_stat == GLP_FEAS && mip->dbs_stat == GLP_FEAS);
	273	/* determine the objective function value */
	274	obj = mip->obj_val;
	275	/* walk through the column list */
	276	for (j = 1; j <= mip->n; j++)
	277	{ GLPCOL *col = mip->col[j];
	278	/* if the column is not integer, skip it */
	279	if (col->kind != GLP_IV) continue;
	280	/* obtain bounds of j-th column */
	281	lb = col->lb, ub = col->ub;
	282	/* and determine its status and reduced cost */
	283	stat = col->stat, dj = col->dual;
	284	/* analyze the reduced cost */
	285	switch (mip->dir)
	286	{ case GLP_MIN:
	287	/* minimization */
	288	if (stat == GLP_NL)
	289	{ /* j-th column is non-basic on its lower bound */
	290	if (dj < 0.0) dj = 0.0;
	291	if (obj + dj >= mip->mip_obj)
	292	glp_set_col_bnds(mip, j, GLP_FX, lb, lb), fixed++;
	293	}
	294	else if (stat == GLP_NU)
	295	{ /* j-th column is non-basic on its upper bound */
	296	if (dj > 0.0) dj = 0.0;
	297	if (obj - dj >= mip->mip_obj)
	298	glp_set_col_bnds(mip, j, GLP_FX, ub, ub), fixed++;
	299	}
	300	break;
	301	case GLP_MAX:
	302	/* maximization */
	303	if (stat == GLP_NL)
	304	{ /* j-th column is non-basic on its lower bound */
	305	if (dj > 0.0) dj = 0.0;
	306	if (obj + dj <= mip->mip_obj)
	307	glp_set_col_bnds(mip, j, GLP_FX, lb, lb), fixed++;
	308	}
	309	else if (stat == GLP_NU)
	310	{ /* j-th column is non-basic on its upper bound */
	311	if (dj < 0.0) dj = 0.0;
	312	if (obj - dj <= mip->mip_obj)
	313	glp_set_col_bnds(mip, j, GLP_FX, ub, ub), fixed++;
	314	}
	315	break;
	316	default:
	317	xassert(T != T);
	318	}
	319	}
	320	if (T->parm->msg_lev >= GLP_MSG_DBG)
	321	{ if (fixed == 0)
	322	/* nothing to say */;
	323	else if (fixed == 1)
	324	xprintf("One column has been fixed by reduced cost\n");
	325	else
	326	xprintf("%d columns have been fixed by reduced costs\n",
	327	fixed);
	328	}
	329	/* fixing non-basic columns on their current bounds does not
	330	change the basic solution */
	331	xassert(mip->pbs_stat == GLP_FEAS && mip->dbs_stat == GLP_FEAS);
	332	return;
	333	}
	334
	335	/***********************************************************************
	336	* branch_on - perform branching on specified variable
	337	*
	338	* This routine performs branching on j-th column (structural variable)
	339	* of the current subproblem. The specified column must be of integer
	340	* kind and must have a fractional value in optimal basic solution of
	341	* LP relaxation of the current subproblem (i.e. only columns for which
	342	* the flag non_int[j] is set are valid candidates to branch on).
	343	*
	344	* Let x be j-th structural variable, and beta be its primal fractional
	345	* value in the current basic solution. Branching on j-th variable is
	346	* dividing the current subproblem into two new subproblems, which are
	347	* identical to the current subproblem with the following exception: in
	348	* the first subproblem that begins the down-branch x has a new upper
	349	* bound x <= floor(beta), and in the second subproblem that begins the
	350	* up-branch x has a new lower bound x >= ceil(beta).
	351	*
	352	* Depending on estimation of local bounds for down- and up-branches
	353	* this routine returns the following:
	354	*
	355	* 0 - both branches have been created;
	356	* 1 - one branch is hopeless and has been pruned, so now the current
	357	* subproblem is other branch;
	358	* 2 - both branches are hopeless and have been pruned; new subproblem
	359	* selection is needed to continue the search. */
	360
	361	static int branch_on(glp_tree *T, int j, int next)
	362	{ glp_prob *mip = T->mip;
	363	IOSNPD *node;
	364	int m = mip->m;
	365	int n = mip->n;
	366	int type, dn_type, up_type, dn_bad, up_bad, p, ret, clone[1+2];
	367	double lb, ub, beta, new_ub, new_lb, dn_lp, up_lp, dn_bnd, up_bnd;
	368	/* determine bounds and value of x[j] in optimal solution to LP
	369	relaxation of the current subproblem */
	370	xassert(1 <= j && j <= n);
	371	type = mip->col[j]->type;
	372	lb = mip->col[j]->lb;
	373	ub = mip->col[j]->ub;
	374	beta = mip->col[j]->prim;
	375	/* determine new bounds of x[j] for down- and up-branches */
	376	new_ub = floor(beta);
	377	new_lb = ceil(beta);
	378	switch (type)
	379	{ case GLP_FR:
	380	dn_type = GLP_UP;
	381	up_type = GLP_LO;
	382	break;
	383	case GLP_LO:
	384	xassert(lb <= new_ub);
	385	dn_type = (lb == new_ub ? GLP_FX : GLP_DB);
	386	xassert(lb + 1.0 <= new_lb);
	387	up_type = GLP_LO;
	388	break;
	389	case GLP_UP:
	390	xassert(new_ub <= ub - 1.0);
	391	dn_type = GLP_UP;
	392	xassert(new_lb <= ub);
	393	up_type = (new_lb == ub ? GLP_FX : GLP_DB);
	394	break;
	395	case GLP_DB:
	396	xassert(lb <= new_ub && new_ub <= ub - 1.0);
	397	dn_type = (lb == new_ub ? GLP_FX : GLP_DB);
	398	xassert(lb + 1.0 <= new_lb && new_lb <= ub);
	399	up_type = (new_lb == ub ? GLP_FX : GLP_DB);
	400	break;
	401	default:
	402	xassert(type != type);
	403	}
	404	/* compute local bounds to LP relaxation for both branches */
	405	ios_eval_degrad(T, j, &dn_lp, &up_lp);
	406	/* and improve them by rounding */
	407	dn_bnd = ios_round_bound(T, dn_lp);
	408	up_bnd = ios_round_bound(T, up_lp);
	409	/* check local bounds for down- and up-branches */
	410	dn_bad = !ios_is_hopeful(T, dn_bnd);
	411	up_bad = !ios_is_hopeful(T, up_bnd);
	412	if (dn_bad && up_bad)
	413	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	414	xprintf("Both down- and up-branches are hopeless\n");
	415	ret = 2;
	416	goto done;
	417	}
	418	else if (up_bad)
	419	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	420	xprintf("Up-branch is hopeless\n");
	421	glp_set_col_bnds(mip, j, dn_type, lb, new_ub);
	422	T->curr->lp_obj = dn_lp;
	423	if (mip->dir == GLP_MIN)
	424	{ if (T->curr->bound < dn_bnd)
	425	T->curr->bound = dn_bnd;
	426	}
	427	else if (mip->dir == GLP_MAX)
	428	{ if (T->curr->bound > dn_bnd)
	429	T->curr->bound = dn_bnd;
	430	}
	431	else
	432	xassert(mip != mip);
	433	ret = 1;
	434	goto done;
	435	}
	436	else if (dn_bad)
	437	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	438	xprintf("Down-branch is hopeless\n");
	439	glp_set_col_bnds(mip, j, up_type, new_lb, ub);
	440	T->curr->lp_obj = up_lp;
	441	if (mip->dir == GLP_MIN)
	442	{ if (T->curr->bound < up_bnd)
	443	T->curr->bound = up_bnd;
	444	}
	445	else if (mip->dir == GLP_MAX)
	446	{ if (T->curr->bound > up_bnd)
	447	T->curr->bound = up_bnd;
	448	}
	449	else
	450	xassert(mip != mip);
	451	ret = 1;
	452	goto done;
	453	}
	454	/* both down- and up-branches seem to be hopeful */
	455	if (T->parm->msg_lev >= GLP_MSG_DBG)
	456	xprintf("Branching on column %d, primal value is %.9e\n",
	457	j, beta);
	458	/* determine the reference number of the current subproblem */
	459	xassert(T->curr != NULL);
	460	p = T->curr->p;
	461	T->curr->br_var = j;
	462	T->curr->br_val = beta;
	463	/* freeze the current subproblem */
	464	ios_freeze_node(T);
	465	/* create two clones of the current subproblem; the first clone
	466	begins the down-branch, the second one begins the up-branch */
	467	ios_clone_node(T, p, 2, clone);
	468	if (T->parm->msg_lev >= GLP_MSG_DBG)
	469	xprintf("Node %d begins down branch, node %d begins up branch "
	470	"\n", clone[1], clone[2]);
	471	/* set new upper bound of j-th column in the down-branch */
	472	node = T->slot[clone[1]].node;
	473	xassert(node != NULL);
	474	xassert(node->up != NULL);
	475	xassert(node->b_ptr == NULL);
	476	node->b_ptr = dmp_get_atom(T->pool, sizeof(IOSBND));
	477	node->b_ptr->k = m + j;
	478	node->b_ptr->type = (unsigned char)dn_type;
	479	node->b_ptr->lb = lb;
	480	node->b_ptr->ub = new_ub;
	481	node->b_ptr->next = NULL;
	482	node->lp_obj = dn_lp;
	483	if (mip->dir == GLP_MIN)
	484	{ if (node->bound < dn_bnd)
	485	node->bound = dn_bnd;
	486	}
	487	else if (mip->dir == GLP_MAX)
	488	{ if (node->bound > dn_bnd)
	489	node->bound = dn_bnd;
	490	}
	491	else
	492	xassert(mip != mip);
	493	/* set new lower bound of j-th column in the up-branch */
	494	node = T->slot[clone[2]].node;
	495	xassert(node != NULL);
	496	xassert(node->up != NULL);
	497	xassert(node->b_ptr == NULL);
	498	node->b_ptr = dmp_get_atom(T->pool, sizeof(IOSBND));
	499	node->b_ptr->k = m + j;
	500	node->b_ptr->type = (unsigned char)up_type;
	501	node->b_ptr->lb = new_lb;
	502	node->b_ptr->ub = ub;
	503	node->b_ptr->next = NULL;
	504	node->lp_obj = up_lp;
	505	if (mip->dir == GLP_MIN)
	506	{ if (node->bound < up_bnd)
	507	node->bound = up_bnd;
	508	}
	509	else if (mip->dir == GLP_MAX)
	510	{ if (node->bound > up_bnd)
	511	node->bound = up_bnd;
	512	}
	513	else
	514	xassert(mip != mip);
	515	/* suggest the subproblem to be solved next */
	516	xassert(T->child == 0);
	517	if (next == GLP_NO_BRNCH)
	518	T->child = 0;
	519	else if (next == GLP_DN_BRNCH)
	520	T->child = clone[1];
	521	else if (next == GLP_UP_BRNCH)
	522	T->child = clone[2];
	523	else
	524	xassert(next != next);
	525	ret = 0;
	526	done: return ret;
	527	}
	528
	529	/***********************************************************************
	530	* cleanup_the_tree - prune hopeless branches from the tree
	531	*
	532	* This routine walks through the active list and checks the local
	533	* bound for every active subproblem. If the local bound indicates that
	534	* the subproblem cannot have integer optimal solution better than the
	535	* incumbent objective value, the routine deletes such subproblem that,
	536	* in turn, involves pruning the corresponding branch of the tree. */
	537
	538	static void cleanup_the_tree(glp_tree *T)
	539	{ IOSNPD node, next_node;
	540	int count = 0;
	541	/* the global bound must exist */
	542	xassert(T->mip->mip_stat == GLP_FEAS);
	543	/* walk through the list of active subproblems */
	544	for (node = T->head; node != NULL; node = next_node)
	545	{ /* deleting some active problem node may involve deleting its
	546	parents recursively; however, all its parents being created
	547	before it are always precede it in the node list, so
	548	the next problem node is never affected by such deletion */
	549	next_node = node->next;
	550	/* if the branch is hopeless, prune it */
	551	if (!is_branch_hopeful(T, node->p))
	552	ios_delete_node(T, node->p), count++;
	553	}
	554	if (T->parm->msg_lev >= GLP_MSG_DBG)
	555	{ if (count == 1)
	556	xprintf("One hopeless branch has been pruned\n");
	557	else if (count > 1)
	558	xprintf("%d hopeless branches have been pruned\n", count);
	559	}
	560	return;
	561	}
	562
	563	/**********************************************************************/
	564
	565	static void generate_cuts(glp_tree *T)
	566	{ /* generate generic cuts with built-in generators */
	567	if (!(T->parm->mir_cuts == GLP_ON \|\|
	568	T->parm->gmi_cuts == GLP_ON \|\|
	569	T->parm->cov_cuts == GLP_ON \|\|
	570	T->parm->clq_cuts == GLP_ON)) goto done;
	571	#if 1 /* 20/IX-2008 */
	572	{ int i, max_cuts, added_cuts;
	573	max_cuts = T->n;
	574	if (max_cuts < 1000) max_cuts = 1000;
	575	added_cuts = 0;
	576	for (i = T->orig_m+1; i <= T->mip->m; i++)
	577	{ if (T->mip->row[i]->origin == GLP_RF_CUT)
	578	added_cuts++;
	579	}
	580	/* xprintf("added_cuts = %d\n", added_cuts); */
	581	if (added_cuts >= max_cuts) goto done;
	582	}
	583	#endif
	584	/* generate and add to POOL all cuts violated by x* */
	585	if (T->parm->gmi_cuts == GLP_ON)
	586	{ if (T->curr->changed < 5)
	587	ios_gmi_gen(T);
	588	}
	589	if (T->parm->mir_cuts == GLP_ON)
	590	{ xassert(T->mir_gen != NULL);
	591	ios_mir_gen(T, T->mir_gen);
	592	}
	593	if (T->parm->cov_cuts == GLP_ON)
	594	{ /* cover cuts works well along with mir cuts */
	595	/if (T->round <= 5)/
	596	ios_cov_gen(T);
	597	}
	598	if (T->parm->clq_cuts == GLP_ON)
	599	{ if (T->clq_gen != NULL)
	600	{ if (T->curr->level == 0 && T->curr->changed < 50 \|\|
	601	T->curr->level > 0 && T->curr->changed < 5)
	602	ios_clq_gen(T, T->clq_gen);
	603	}
	604	}
	605	done: return;
	606	}
	607
	608	/**********************************************************************/
	609
	610	static void remove_cuts(glp_tree *T)
	611	{ /* remove inactive cuts (some valueable globally valid cut might
	612	be saved in the global cut pool) */
	613	int i, cnt = 0, *num = NULL;
	614	xassert(T->curr != NULL);
	615	for (i = T->orig_m+1; i <= T->mip->m; i++)
	616	{ if (T->mip->row[i]->origin == GLP_RF_CUT &&
	617	T->mip->row[i]->level == T->curr->level &&
	618	T->mip->row[i]->stat == GLP_BS)
	619	{ if (num == NULL)
	620	num = xcalloc(1+T->mip->m, sizeof(int));
	621	num[++cnt] = i;
	622	}
	623	}
	624	if (cnt > 0)
	625	{ glp_del_rows(T->mip, cnt, num);
	626	#if 0
	627	xprintf("%d inactive cut(s) removed\n", cnt);
	628	#endif
	629	xfree(num);
	630	xassert(glp_factorize(T->mip) == 0);
	631	}
	632	return;
	633	}
	634
	635	/**********************************************************************/
	636
	637	static void display_cut_info(glp_tree *T)
	638	{ glp_prob *mip = T->mip;
	639	int i, gmi = 0, mir = 0, cov = 0, clq = 0, app = 0;
	640	for (i = mip->m; i > 0; i--)
	641	{ GLPROW *row;
	642	row = mip->row[i];
	643	/* if (row->level < T->curr->level) break; */
	644	if (row->origin == GLP_RF_CUT)
	645	{ if (row->klass == GLP_RF_GMI)
	646	gmi++;
	647	else if (row->klass == GLP_RF_MIR)
	648	mir++;
	649	else if (row->klass == GLP_RF_COV)
	650	cov++;
	651	else if (row->klass == GLP_RF_CLQ)
	652	clq++;
	653	else
	654	app++;
	655	}
	656	}
	657	xassert(T->curr != NULL);
	658	if (gmi + mir + cov + clq + app > 0)
	659	{ xprintf("Cuts on level %d:", T->curr->level);
	660	if (gmi > 0) xprintf(" gmi = %d;", gmi);
	661	if (mir > 0) xprintf(" mir = %d;", mir);
	662	if (cov > 0) xprintf(" cov = %d;", cov);
	663	if (clq > 0) xprintf(" clq = %d;", clq);
	664	if (app > 0) xprintf(" app = %d;", app);
	665	xprintf("\n");
	666	}
	667	return;
	668	}
	669
	670	/***********************************************************************
	671	* NAME
	672	*
	673	* ios_driver - branch-and-cut driver
	674	*
	675	* SYNOPSIS
	676	*
	677	* #include "glpios.h"
	678	* int ios_driver(glp_tree *T);
	679	*
	680	* DESCRIPTION
	681	*
	682	* The routine ios_driver is a branch-and-cut driver. It controls the
	683	* MIP solution process.
	684	*
	685	* RETURNS
	686	*
	687	* 0 The MIP problem instance has been successfully solved. This code
	688	* does not necessarily mean that the solver has found optimal
	689	* solution. It only means that the solution process was successful.
	690	*
	691	* GLP_EFAIL
	692	* The search was prematurely terminated due to the solver failure.
	693	*
	694	* GLP_EMIPGAP
	695	* The search was prematurely terminated, because the relative mip
	696	* gap tolerance has been reached.
	697	*
	698	* GLP_ETMLIM
	699	* The search was prematurely terminated, because the time limit has
	700	* been exceeded.
	701	*
	702	* GLP_ESTOP
	703	* The search was prematurely terminated by application. */
	704
	705	int ios_driver(glp_tree *T)
	706	{ int p, curr_p, p_stat, d_stat, ret;
	707	#if 1 /* carry out to glp_tree */
	708	int pred_p = 0;
	709	/* if the current subproblem has been just created due to
	710	branching, pred_p is the reference number of its parent
	711	subproblem, otherwise pred_p is zero */
	712	#endif
	713	glp_long ttt = T->tm_beg;
	714	#if 0
	715	((glp_iocp *)T->parm)->msg_lev = GLP_MSG_DBG;
	716	#endif
	717	/* on entry to the B&B driver it is assumed that the active list
	718	contains the only active (i.e. root) subproblem, which is the
	719	original MIP problem to be solved */
	720	loop: /* main loop starts here */
	721	/* at this point the current subproblem does not exist */
	722	xassert(T->curr == NULL);
	723	/* if the active list is empty, the search is finished */
	724	if (T->head == NULL)
	725	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	726	xprintf("Active list is empty!\n");
	727	xassert(dmp_in_use(T->pool).lo == 0);
	728	ret = 0;
	729	goto done;
	730	}
	731	/* select some active subproblem to continue the search */
	732	xassert(T->next_p == 0);
	733	/* let the application program select subproblem */
	734	if (T->parm->cb_func != NULL)
	735	{ xassert(T->reason == 0);
	736	T->reason = GLP_ISELECT;
	737	T->parm->cb_func(T, T->parm->cb_info);
	738	T->reason = 0;
	739	if (T->stop)
	740	{ ret = GLP_ESTOP;
	741	goto done;
	742	}
	743	}
	744	if (T->next_p != 0)
	745	{ /* the application program has selected something */
	746	;
	747	}
	748	else if (T->a_cnt == 1)
	749	{ /* the only active subproblem exists, so select it */
	750	xassert(T->head->next == NULL);
	751	T->next_p = T->head->p;
	752	}
	753	else if (T->child != 0)
	754	{ /* select one of branching childs suggested by the branching
	755	heuristic */
	756	T->next_p = T->child;
	757	}
	758	else
	759	{ /* select active subproblem as specified by the backtracking
	760	technique option */
	761	T->next_p = ios_choose_node(T);
	762	}
	763	/* the active subproblem just selected becomes current */
	764	ios_revive_node(T, T->next_p);
	765	T->next_p = T->child = 0;
	766	/* invalidate pred_p, if it is not the reference number of the
	767	parent of the current subproblem */
	768	if (T->curr->up != NULL && T->curr->up->p != pred_p) pred_p = 0;
	769	/* determine the reference number of the current subproblem */
	770	p = T->curr->p;
	771	if (T->parm->msg_lev >= GLP_MSG_DBG)
	772	{ xprintf("-----------------------------------------------------"
	773	"-------------------\n");
	774	xprintf("Processing node %d at level %d\n", p, T->curr->level);
	775	}
	776	/* if it is the root subproblem, initialize cut generators */
	777	if (p == 1)
	778	{ if (T->parm->gmi_cuts == GLP_ON)
	779	{ if (T->parm->msg_lev >= GLP_MSG_ALL)
	780	xprintf("Gomory's cuts enabled\n");
	781	}
	782	if (T->parm->mir_cuts == GLP_ON)
	783	{ if (T->parm->msg_lev >= GLP_MSG_ALL)
	784	xprintf("MIR cuts enabled\n");
	785	xassert(T->mir_gen == NULL);
	786	T->mir_gen = ios_mir_init(T);
	787	}
	788	if (T->parm->cov_cuts == GLP_ON)
	789	{ if (T->parm->msg_lev >= GLP_MSG_ALL)
	790	xprintf("Cover cuts enabled\n");
	791	}
	792	if (T->parm->clq_cuts == GLP_ON)
	793	{ xassert(T->clq_gen == NULL);
	794	if (T->parm->msg_lev >= GLP_MSG_ALL)
	795	xprintf("Clique cuts enabled\n");
	796	T->clq_gen = ios_clq_init(T);
	797	}
	798	}
	799	more: /* minor loop starts here */
	800	/* at this point the current subproblem needs either to be solved
	801	for the first time or re-optimized due to reformulation */
	802	/* display current progress of the search */
	803	if (T->parm->msg_lev >= GLP_MSG_DBG \|\|
	804	T->parm->msg_lev >= GLP_MSG_ON &&
	805	(double)(T->parm->out_frq - 1) <=
	806	1000.0 * xdifftime(xtime(), T->tm_lag))
	807	show_progress(T, 0);
	808	if (T->parm->msg_lev >= GLP_MSG_ALL &&
	809	xdifftime(xtime(), ttt) >= 60.0)
	810	{ glp_long total;
	811	glp_mem_usage(NULL, NULL, &total, NULL);
	812	xprintf("Time used: %.1f secs. Memory used: %.1f Mb.\n",
	813	xdifftime(xtime(), T->tm_beg), xltod(total) / 1048576.0);
	814	ttt = xtime();
	815	}
	816	/* check the mip gap */
	817	if (T->parm->mip_gap > 0.0 &&
	818	ios_relative_gap(T) <= T->parm->mip_gap)
	819	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	820	xprintf("Relative gap tolerance reached; search terminated "
	821	"\n");
	822	ret = GLP_EMIPGAP;
	823	goto done;
	824	}
	825	/* check if the time limit has been exhausted */
	826	if (T->parm->tm_lim < INT_MAX &&
	827	(double)(T->parm->tm_lim - 1) <=
	828	1000.0 * xdifftime(xtime(), T->tm_beg))
	829	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	830	xprintf("Time limit exhausted; search terminated\n");
	831	ret = GLP_ETMLIM;
	832	goto done;
	833	}
	834	/* let the application program preprocess the subproblem */
	835	if (T->parm->cb_func != NULL)
	836	{ xassert(T->reason == 0);
	837	T->reason = GLP_IPREPRO;
	838	T->parm->cb_func(T, T->parm->cb_info);
	839	T->reason = 0;
	840	if (T->stop)
	841	{ ret = GLP_ESTOP;
	842	goto done;
	843	}
	844	}
	845	/* perform basic preprocessing */
	846	if (T->parm->pp_tech == GLP_PP_NONE)
	847	;
	848	else if (T->parm->pp_tech == GLP_PP_ROOT)
	849	{ if (T->curr->level == 0)
	850	{ if (ios_preprocess_node(T, 100))
	851	goto fath;
	852	}
	853	}
	854	else if (T->parm->pp_tech == GLP_PP_ALL)
	855	{ if (ios_preprocess_node(T, T->curr->level == 0 ? 100 : 10))
	856	goto fath;
	857	}
	858	else
	859	xassert(T != T);
	860	/* preprocessing may improve the global bound */
	861	if (!is_branch_hopeful(T, p))
	862	{ xprintf("* not tested yet *\n");
	863	goto fath;
	864	}
	865	/* solve LP relaxation of the current subproblem */
	866	if (T->parm->msg_lev >= GLP_MSG_DBG)
	867	xprintf("Solving LP relaxation...\n");
	868	ret = ios_solve_node(T);
	869	if (!(ret == 0 \|\| ret == GLP_EOBJLL \|\| ret == GLP_EOBJUL))
	870	{ if (T->parm->msg_lev >= GLP_MSG_ERR)
	871	xprintf("ios_driver: unable to solve current LP relaxation;"
	872	" glp_simplex returned %d\n", ret);
	873	ret = GLP_EFAIL;
	874	goto done;
	875	}
	876	/* analyze status of the basic solution to LP relaxation found */
	877	p_stat = T->mip->pbs_stat;
	878	d_stat = T->mip->dbs_stat;
	879	if (p_stat == GLP_FEAS && d_stat == GLP_FEAS)
	880	{ /* LP relaxation has optimal solution */
	881	if (T->parm->msg_lev >= GLP_MSG_DBG)
	882	xprintf("Found optimal solution to LP relaxation\n");
	883	}
	884	else if (d_stat == GLP_NOFEAS)
	885	{ /* LP relaxation has no dual feasible solution */
	886	/* since the current subproblem cannot have a larger feasible
	887	region than its parent, there is something wrong */
	888	if (T->parm->msg_lev >= GLP_MSG_ERR)
	889	xprintf("ios_driver: current LP relaxation has no dual feas"
	890	"ible solution\n");
	891	ret = GLP_EFAIL;
	892	goto done;
	893	}
	894	else if (p_stat == GLP_INFEAS && d_stat == GLP_FEAS)
	895	{ /* LP relaxation has no primal solution which is better than
	896	the incumbent objective value */
	897	xassert(T->mip->mip_stat == GLP_FEAS);
	898	if (T->parm->msg_lev >= GLP_MSG_DBG)
	899	xprintf("LP relaxation has no solution better than incumben"
	900	"t objective value\n");
	901	/* prune the branch */
	902	goto fath;
	903	}
	904	else if (p_stat == GLP_NOFEAS)
	905	{ /* LP relaxation has no primal feasible solution */
	906	if (T->parm->msg_lev >= GLP_MSG_DBG)
	907	xprintf("LP relaxation has no feasible solution\n");
	908	/* prune the branch */
	909	goto fath;
	910	}
	911	else
	912	{ /* other cases cannot appear */
	913	xassert(T->mip != T->mip);
	914	}
	915	/* at this point basic solution to LP relaxation of the current
	916	subproblem is optimal */
	917	xassert(p_stat == GLP_FEAS && d_stat == GLP_FEAS);
	918	xassert(T->curr != NULL);
	919	T->curr->lp_obj = T->mip->obj_val;
	920	/* thus, it defines a local bound to integer optimal solution of
	921	the current subproblem */
	922	{ double bound = T->mip->obj_val;
	923	/* some local bound to the current subproblem could be already
	924	set before, so we should only improve it */
	925	bound = ios_round_bound(T, bound);
	926	if (T->mip->dir == GLP_MIN)
	927	{ if (T->curr->bound < bound)
	928	T->curr->bound = bound;
	929	}
	930	else if (T->mip->dir == GLP_MAX)
	931	{ if (T->curr->bound > bound)
	932	T->curr->bound = bound;
	933	}
	934	else
	935	xassert(T->mip != T->mip);
	936	if (T->parm->msg_lev >= GLP_MSG_DBG)
	937	xprintf("Local bound is %.9e\n", bound);
	938	}
	939	/* if the local bound indicates that integer optimal solution of
	940	the current subproblem cannot be better than the global bound,
	941	prune the branch */
	942	if (!is_branch_hopeful(T, p))
	943	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	944	xprintf("Current branch is hopeless and can be pruned\n");
	945	goto fath;
	946	}
	947	/* let the application program generate additional rows ("lazy"
	948	constraints) */
	949	xassert(T->reopt == 0);
	950	xassert(T->reinv == 0);
	951	if (T->parm->cb_func != NULL)
	952	{ xassert(T->reason == 0);
	953	T->reason = GLP_IROWGEN;
	954	T->parm->cb_func(T, T->parm->cb_info);
	955	T->reason = 0;
	956	if (T->stop)
	957	{ ret = GLP_ESTOP;
	958	goto done;
	959	}
	960	if (T->reopt)
	961	{ /* some rows were added; re-optimization is needed */
	962	T->reopt = T->reinv = 0;
	963	goto more;
	964	}
	965	if (T->reinv)
	966	{ /* no rows were added, however, some inactive rows were
	967	removed */
	968	T->reinv = 0;
	969	xassert(glp_factorize(T->mip) == 0);
	970	}
	971	}
	972	/* check if the basic solution is integer feasible */
	973	check_integrality(T);
	974	/* if the basic solution satisfies to all integrality conditions,
	975	it is a new, better integer feasible solution */
	976	if (T->curr->ii_cnt == 0)
	977	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	978	xprintf("New integer feasible solution found\n");
	979	if (T->parm->msg_lev >= GLP_MSG_ALL)
	980	display_cut_info(T);
	981	record_solution(T);
	982	if (T->parm->msg_lev >= GLP_MSG_ON)
	983	show_progress(T, 1);
	984	/* make the application program happy */
	985	if (T->parm->cb_func != NULL)
	986	{ xassert(T->reason == 0);
	987	T->reason = GLP_IBINGO;
	988	T->parm->cb_func(T, T->parm->cb_info);
	989	T->reason = 0;
	990	if (T->stop)
	991	{ ret = GLP_ESTOP;
	992	goto done;
	993	}
	994	}
	995	/* since the current subproblem has been fathomed, prune its
	996	branch */
	997	goto fath;
	998	}
	999	/* at this point basic solution to LP relaxation of the current
	1000	subproblem is optimal, but integer infeasible */
	1001	/* try to fix some non-basic structural variables of integer kind
	1002	on their current bounds due to reduced costs */
	1003	if (T->mip->mip_stat == GLP_FEAS)
	1004	fix_by_red_cost(T);
	1005	/* let the application program try to find some solution to the
	1006	original MIP with a primal heuristic */
	1007	if (T->parm->cb_func != NULL)
	1008	{ xassert(T->reason == 0);
	1009	T->reason = GLP_IHEUR;
	1010	T->parm->cb_func(T, T->parm->cb_info);
	1011	T->reason = 0;
	1012	if (T->stop)
	1013	{ ret = GLP_ESTOP;
	1014	goto done;
	1015	}
	1016	/* check if the current branch became hopeless */
	1017	if (!is_branch_hopeful(T, p))
	1018	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	1019	xprintf("Current branch became hopeless and can be prune"
	1020	"d\n");
	1021	goto fath;
	1022	}
	1023	}
	1024	/* try to find solution with the feasibility pump heuristic */
	1025	if (T->parm->fp_heur)
	1026	{ xassert(T->reason == 0);
	1027	T->reason = GLP_IHEUR;
	1028	ios_feas_pump(T);
	1029	T->reason = 0;
	1030	/* check if the current branch became hopeless */
	1031	if (!is_branch_hopeful(T, p))
	1032	{ if (T->parm->msg_lev >= GLP_MSG_DBG)
	1033	xprintf("Current branch became hopeless and can be prune"
	1034	"d\n");
	1035	goto fath;
	1036	}
	1037	}
	1038	/* it's time to generate cutting planes */
	1039	xassert(T->local != NULL);
	1040	xassert(T->local->size == 0);
	1041	/* let the application program generate some cuts; note that it
	1042	can add cuts either to the local cut pool or directly to the
	1043	current subproblem */
	1044	if (T->parm->cb_func != NULL)
	1045	{ xassert(T->reason == 0);
	1046	T->reason = GLP_ICUTGEN;
	1047	T->parm->cb_func(T, T->parm->cb_info);
	1048	T->reason = 0;
	1049	if (T->stop)
	1050	{ ret = GLP_ESTOP;
	1051	goto done;
	1052	}
	1053	}
	1054	/* try to generate generic cuts with built-in generators
	1055	(as suggested by Matteo Fischetti et al. the built-in cuts
	1056	are not generated at each branching node; an intense attempt
	1057	of generating new cuts is only made at the root node, and then
	1058	a moderate effort is spent after each backtracking step) */
	1059	if (T->curr->level == 0 \|\| pred_p == 0)
	1060	{ xassert(T->reason == 0);
	1061	T->reason = GLP_ICUTGEN;
	1062	generate_cuts(T);
	1063	T->reason = 0;
	1064	}
	1065	/* if the local cut pool is not empty, select useful cuts and add
	1066	them to the current subproblem */
	1067	if (T->local->size > 0)
	1068	{ xassert(T->reason == 0);
	1069	T->reason = GLP_ICUTGEN;
	1070	ios_process_cuts(T);
	1071	T->reason = 0;
	1072	}
	1073	/* clear the local cut pool */
	1074	ios_clear_pool(T, T->local);
	1075	/* perform re-optimization, if necessary */
	1076	if (T->reopt)
	1077	{ T->reopt = 0;
	1078	T->curr->changed++;
	1079	goto more;
	1080	}
	1081	/* no cuts were generated; remove inactive cuts */
	1082	remove_cuts(T);
	1083	if (T->parm->msg_lev >= GLP_MSG_ALL && T->curr->level == 0)
	1084	display_cut_info(T);
	1085	/* update history information used on pseudocost branching */
	1086	if (T->pcost != NULL) ios_pcost_update(T);
	1087	/* it's time to perform branching */
	1088	xassert(T->br_var == 0);
	1089	xassert(T->br_sel == 0);
	1090	/* let the application program choose variable to branch on */
	1091	if (T->parm->cb_func != NULL)
	1092	{ xassert(T->reason == 0);
	1093	xassert(T->br_var == 0);
	1094	xassert(T->br_sel == 0);
	1095	T->reason = GLP_IBRANCH;
	1096	T->parm->cb_func(T, T->parm->cb_info);
	1097	T->reason = 0;
	1098	if (T->stop)
	1099	{ ret = GLP_ESTOP;
	1100	goto done;
	1101	}
	1102	}
	1103	/* if nothing has been chosen, choose some variable as specified
	1104	by the branching technique option */
	1105	if (T->br_var == 0)
	1106	T->br_var = ios_choose_var(T, &T->br_sel);
	1107	/* perform actual branching */
	1108	curr_p = T->curr->p;
	1109	ret = branch_on(T, T->br_var, T->br_sel);
	1110	T->br_var = T->br_sel = 0;
	1111	if (ret == 0)
	1112	{ /* both branches have been created */
	1113	pred_p = curr_p;
	1114	goto loop;
	1115	}
	1116	else if (ret == 1)
	1117	{ /* one branch is hopeless and has been pruned, so now the
	1118	current subproblem is other branch */
	1119	/* the current subproblem should be considered as a new one,
	1120	since one bound of the branching variable was changed */
	1121	T->curr->solved = T->curr->changed = 0;
	1122	goto more;
	1123	}
	1124	else if (ret == 2)
	1125	{ /* both branches are hopeless and have been pruned; new
	1126	subproblem selection is needed to continue the search */
	1127	goto fath;
	1128	}
	1129	else
	1130	xassert(ret != ret);
	1131	fath: /* the current subproblem has been fathomed */
	1132	if (T->parm->msg_lev >= GLP_MSG_DBG)
	1133	xprintf("Node %d fathomed\n", p);
	1134	/* freeze the current subproblem */
	1135	ios_freeze_node(T);
	1136	/* and prune the corresponding branch of the tree */
	1137	ios_delete_node(T, p);
	1138	/* if a new integer feasible solution has just been found, other
	1139	branches may become hopeless and therefore must be pruned */
	1140	if (T->mip->mip_stat == GLP_FEAS) cleanup_the_tree(T);
	1141	/* new subproblem selection is needed due to backtracking */
	1142	pred_p = 0;
	1143	goto loop;
	1144	done: /* display progress of the search on exit from the solver */
	1145	if (T->parm->msg_lev >= GLP_MSG_ON)
	1146	show_progress(T, 0);
	1147	if (T->mir_gen != NULL)
	1148	ios_mir_term(T->mir_gen), T->mir_gen = NULL;
	1149	if (T->clq_gen != NULL)
	1150	ios_clq_term(T->clq_gen), T->clq_gen = NULL;
	1151	/* return to the calling program */
	1152	return ret;
	1153	}
	1154
	1155	/* eof */

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

Context Navigation

source: glpk-cmake/src/glpios03.c @ 1:c445c931472f

Download in other formats: