constrnt.c

生成直角Steiner树的程序包

/***********************************************************************

	File:	constrnt.c
	Rev:	a-2
	Date:	02/28/2001

	Copyright (c) 1996, 2001 by David M. Warme

************************************************************************

	Routines for handling constraints.

************************************************************************

	Modification Log:

	a-1:	11/18/96	warme
		: Created.  Split off from main file.
	a-2:	02/28/2001	warme
		: Numerous changes for 3.1 release.

************************************************************************/

#include "bb.h"
#include "config.h"
#include "constrnt.h"
#include "steiner.h"


/*
 * Global Routines
 */

bool		add_constraint_to_pool (struct cpool *,
					struct rcoef *,
					bool);
int		add_constraints (struct bbinfo *, struct constraint *);
void		add_pending_rows_to_LP (struct bbinfo *);
LP_t *		build_initial_formulation (struct cpool *,
					   bitmap_t *,
					   bitmap_t *,
					   struct cinfo *,
					   struct lpmem *);
void		debug_print_constraint (char *,
					char *,
					struct constraint *,
					double *,
					bitmap_t *,
					struct cinfo *);
void		delete_slack_rows_from_LP (struct bbinfo *);
void		destroy_initial_formulation (struct bbinfo *);
void		destroy_node_basis (struct bbnode *, struct bbinfo *);
void		initialize_constraint_pool (struct cpool *,
					    bitmap_t *,
					    bitmap_t *,
					    struct cinfo *);
bool		is_violation (struct rcoef *, double *);
void		mark_row_pending_to_LP (struct cpool *, int);
void		restore_node_basis (struct bbnode *, struct bbinfo *);
void		save_node_basis (struct bbnode *, struct bbinfo *);
int		solve_LP_over_constraint_pool (struct bbinfo *);
void		verify_pool (struct cpool *);

bool			Seed_Pool_With_2SECs = TRUE;
int			Target_Pool_Non_Zeros;

#if defined(CPLEX)
 int			min_cplex_rows;
 int			min_cplex_nzs;
#endif

#ifdef LPSOLVE
bool			Use_Perturbations = FALSE;
bool			Use_Scaling = FALSE;
#endif


/*
 * Local Routines
 */

static double		compute_slack_value (struct rcoef *, double *);
static void		garbage_collect_pool (struct cpool *, int, int);
static void		print_pool_memory_usage (struct cpool *);
static void		prune_pending_rows (struct bbinfo *, bool);
static void		reduce_constraint (struct rcoef *);
static struct rblk *	reverse_rblks (struct rblk *);
static int		solve_single_LP (struct bbinfo *,
					 double *,
					 double *,
					 int);
static void		sort_gc_candidates (int *, int32u *, int);
static bool		sprint_term (char *, bool, int, int);
static void		update_lp_solution_history (double *,
						    double *,
						    struct bbinfo *);

#if CPLEX
static void		reload_cplex_problem (struct bbinfo *);
#endif

#if LPSOLVE
static void		get_current_basis (LP_t *, int *, int *);
static void		set_current_basis (LP_t *, int *, int *);
#endif

/*
 * This routine initializes the given constraint pool and fills it with
 * the initial set of constraints:
 *
 *	- The total degree constraint.
 *	- One cutset constraint per terminal.
 *	- All two-terminal SECs.
 *	- All incompatibility constraints that aren't shadowed
 *	  by a two-terminal SEC.
 */

	void
initialize_constraint_pool (

struct cpool *		pool,		/* OUT - the pool to initialize */
bitmap_t *		vert_mask,	/* IN - set of valid vertices */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
struct cinfo *		cip		/* IN - compatibility info */
)
{
int			i, j, k;
int			t;
int			nterms;
int			nedges;
int			nmasks;
int			kmasks;
int			nvt;
int			ncols;
int			nrows;
int			ncoeff;
int			num_total_degree_rows;
int			num_total_degree_coeffs;
int			num_cutset_rows;
int			num_cutset_coeffs;
int			num_incompat_rows;
int			num_incompat_coeffs;
int			num_2sec_rows;
int			num_2sec_coeffs;
int			rowsize;
int			nzsize;
int			fs;
int			common;
int			incompat_threshold;
int *			vp1;
int *			vp2;
int *			vp3;
int *			vp4;
int *			ep1;
int *			ep2;
int *			counts;
int *			tlist;
bitmap_t *		tmask;
bitmap_t *		fsmask;
struct rcoef *		rp;
struct rblk *		blkp;
cpu_time_t		T0;
cpu_time_t		T1;
char			tbuf [32];

	T0 = get_cpu_time ();

	nterms = cip -> num_verts;
	nedges = cip -> num_edges;
	kmasks = cip -> num_vert_masks;
	nmasks = cip -> num_edge_masks;

	if (nedges + RC_VAR_BASE > USHRT_MAX) {
		tracef (" %% Too many FSTs or hyperedges!  Max is %d.\n",
			USHRT_MAX - RC_VAR_BASE);
		exit (1);
	}

	/* We know exactly how many columns (variables) we will */
	/* ever need.  We never add additional variables. */
	ncols = nedges;

	tmask = NEWA (kmasks, bitmap_t);

	num_2sec_rows	= 0;
	num_2sec_coeffs	= 0;
	if (Seed_Pool_With_2SECs) {
		/* Count the number of non-trivial 2SEC constraints,	*/
		/* and their non-zero coefficients.			*/

		tlist	= NEWA (nterms, int);
		counts	= NEWA (nterms, int);
		for (i = 0; i < nterms; i++) {
			counts [i] = 0;
		}

		for (i = 0; i < kmasks; i++) {
			tmask [i] = 0;
		}

		for (i = 0; i < nterms; i++) {
			if (NOT BITON (vert_mask, i)) continue;
			vp1 = tlist;
			ep1 = cip -> term_trees [i];
			ep2 = cip -> term_trees [i + 1];
			while (ep1 < ep2) {
				fs = *ep1++;
				if (NOT BITON (edge_mask, fs)) continue;
				vp3 = cip -> edge [fs];
				vp4 = cip -> edge [fs + 1];
				while (vp3 < vp4) {
					j = *vp3++;
					if (j <= i) continue;
					if (NOT BITON (vert_mask, j)) continue;
					++(counts [j]);
					if (BITON (tmask, j)) continue;
					SETBIT (tmask, j);
					*vp1++ = j;
				}
			}
			vp2 = vp1;
			vp1 = tlist;
			while (vp1 < vp2) {
				j = *vp1++;
				if (counts [j] >= 2) {
					/* S={i,j} is a non-trivial SEC... */
					++num_2sec_rows;
					num_2sec_coeffs += counts [j];
				}
				counts [j] = 0;
				CLRBIT (tmask, j);
			}
		}
		free ((char *) counts);
		free ((char *) tlist);
	}

	/* Compute the number of "unshadowed" incompatibility	*/
	/* constraints, and the number of coefficients in the	*/
	/* total degree constraint.				*/
	num_total_degree_rows	= 1;
	num_total_degree_coeffs	= 0;
	num_incompat_rows	= 0;
	/* (The incompatibility constraint for edges having 2 or more */
	/* terminals in common are shadowed by at least one of the 2SECs */
	/* in the pool -- set threshold to 1 or fewer in common.) */
	incompat_threshold = Seed_Pool_With_2SECs ? 1 : cip -> num_verts;
	memset (tmask, 0, kmasks * sizeof (*tmask));
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;
		++num_total_degree_coeffs;
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			j = *vp1++;
			SETBIT (tmask, j);
		}
		ep1 = cip -> inc_edges [i];
		ep2 = cip -> inc_edges [i + 1];
		while (ep1 < ep2) {
			j = *ep1++;
			if (j >= i) break;
			if (NOT BITON (edge_mask, j)) continue;
			vp1 = cip -> edge [j];
			vp2 = cip -> edge [j + 1];
			common = 0;
			while (vp1 < vp2) {
				k = *vp1++;
				if (BITON (tmask, k)) {
					++common;
				}
			}
			if (common > incompat_threshold) continue;

			++num_incompat_rows;
		}
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			j = *vp1++;
			CLRBIT (tmask, j);
		}
	}
	num_incompat_coeffs = 2 * num_incompat_rows;

	/* Compute the total number of valid terminals, and the	*/
	/* number of coefficients in the 1-terminal cutsets.	*/
	nvt = 0;
	num_cutset_rows		= 0;
	num_cutset_coeffs	= 0;
	for (i = 0; i < cip -> num_verts; i++) {
		if (NOT BITON (vert_mask, i)) continue;
		/* This is a valid terminal.  There will be one	*/
		/* cutset row for it.				*/
		++nvt;
		++num_cutset_rows;
		ep1 = cip -> term_trees [i];
		ep2 = cip -> term_trees [i + 1];
		while (ep1 < ep2) {
			k = *ep1++;
			if (BITON (edge_mask, k)) {
				++num_cutset_coeffs;
			}
		}
	}

	/* Tally the total number of rows and coefficients... */
	nrows	=   num_total_degree_rows	+ num_cutset_rows
		  + num_incompat_rows		+ num_2sec_rows;

	ncoeff	=   num_total_degree_coeffs	+ num_cutset_coeffs
		  + num_incompat_coeffs		+ num_2sec_coeffs;

	rowsize = 2 * nrows;		/* extra space for more rows... */
	nzsize = 4 * ncoeff;		/* extra space for more coefficients */

	blkp		= NEW (struct rblk);
	blkp -> next	= NULL;
	blkp -> base	= NEWA (nzsize, struct rcoef);
	blkp -> ptr	= blkp -> base;
	blkp -> nfree	= nzsize;

	pool -> uid	= 0;
	pool -> rows	= NEWA (rowsize, struct rcon);
	pool -> nrows	= 0;
	pool -> maxrows	= rowsize;
	pool -> num_nz	= 0;
	pool -> lprows	= NEWA (rowsize, int);
	pool -> nlprows	= 0;
	pool -> npend	= 0;
	pool -> blocks	= blkp;
	pool -> cbuf	= NEWA (nedges + 1, struct rcoef);
	pool -> iter	= 0;
	pool -> initrows = 0;
	pool -> nvars	= nedges;
	pool -> hwmrow	= 0;
	pool -> hwmnz	= 0;

	/* Empty all of the hash table buckets... */
	for (i = 0; i < CPOOL_HASH_SIZE; i++) {
		pool -> hash [i] = -1;
	}

	/* Now generate the row for the spanning constraint... */
	rp = pool -> cbuf;
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;
		rp -> var = i + RC_VAR_BASE;
		rp -> val = (cip -> edge_size [i] - 1);
		++rp;
	}
	rp -> var = RC_OP_EQ;
	rp -> val = nvt - 1;
	add_constraint_to_pool (pool, pool -> cbuf, TRUE);

	/* Now generate one cutset constraint per terminal... */
	for (i = 0; i < cip -> num_verts; i++) {
		if (NOT BITON (vert_mask, i)) continue;
		rp = pool -> cbuf;
		ep1 = cip -> term_trees [i];
		ep2 = cip -> term_trees [i + 1];
		while (ep1 < ep2) {
			k = *ep1++;
			if (NOT BITON (edge_mask, k)) continue;
			rp -> var = k + RC_VAR_BASE;
			rp -> val = 1;
			++rp;
		}
		rp -> var = RC_OP_GE;
		rp -> val = 1;
		add_constraint_to_pool (pool, pool -> cbuf, TRUE);
	}

	/* Now generate one constraint per incompatible pair... */
	memset (tmask, 0, kmasks * sizeof (*tmask));
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			j = *vp1++;
			SETBIT (tmask, j);
		}
		ep1 = cip -> inc_edges [i];
		ep2 = cip -> inc_edges [i + 1];
		while (ep1 < ep2) {
			j = *ep1++;
			if (j >= i) break;
			if (NOT BITON (edge_mask, j)) continue;
			vp1 = cip -> edge [j];
			vp2 = cip -> edge [j + 1];
			common = 0;
			while (vp1 < vp2) {
				k = *vp1++;
				if (BITON (tmask, k)) {
					++common;
				}
			}
			if (common > incompat_threshold) continue;

			rp = pool -> cbuf;
			rp [0].var = j + RC_VAR_BASE;
			rp [0].val = 1;
			rp [1].var = i + RC_VAR_BASE;
			rp [1].val = 1;
			rp [2].var = RC_OP_LE;
			rp [2].val = 1;
			add_constraint_to_pool (pool, pool -> cbuf, FALSE);
		}
		vp1 = cip -> edge [i];
		vp2 = cip -> edge [i + 1];
		while (vp1 < vp2) {
			j = *vp1++;
			CLRBIT (tmask, j);
		}
	}

	if (Seed_Pool_With_2SECs) {
		/* Now generate one constraint for each 2-SEC... */

		tlist	= NEWA (nterms, int);
		counts	= NEWA (nterms, int);
		fsmask	= NEWA (nmasks, bitmap_t);
		memset (counts, 0, nterms * sizeof (*counts));
		memset (fsmask, 0, nmasks * sizeof (*fsmask));
		memset (tmask, 0, kmasks * sizeof (*tmask));

		for (i = 0; i < nterms; i++) {
			if (NOT BITON (vert_mask, i)) continue;
			vp1 = tlist;
			ep1 = cip -> term_trees [i];
			ep2 = cip -> term_trees [i + 1];
			while (ep1 < ep2) {
				fs = *ep1++;
				if (NOT BITON (edge_mask, fs)) continue;
				SETBIT (fsmask, fs);
				vp3 = cip -> edge [fs];
				vp4 = cip -> edge [fs + 1];
				while (vp3 < vp4) {
					j = *vp3++;
					if (j <= i) continue;
					if (NOT BITON (vert_mask, j)) continue;
					++(counts [j]);
					if (BITON (tmask, j)) continue;
					SETBIT (tmask, j);
					*vp1++ = j;
				}
			}
			vp2 = vp1;
			vp1 = tlist;
			while (vp1 < vp2) {
				j = *vp1++;
				if (counts [j] < 2) continue;
				/* Generate 2SEC {i,j} */
				rp = pool -> cbuf;
				ep1 = cip -> term_trees [j];
				ep2 = cip -> term_trees [j + 1];
				while (ep1 < ep2) {
					fs = *ep1++;
					if (NOT BITON (fsmask, fs)) continue;

					rp -> var = fs + RC_VAR_BASE;
					rp -> val = 1;
					++rp;
				}
				if (rp < &(pool -> cbuf [2])) {
					fatal ("initialize_constraint_pool: Bug 1.");
				}
				rp -> var = RC_OP_LE;
				rp -> val = 1;
				add_constraint_to_pool (pool, pool -> cbuf, FALSE);
			}
			vp1 = tlist;
			while (vp1 < vp2) {
				j = *vp1++;
				counts [j] = 0;
				CLRBIT (tmask, j);
			}
			ep1 = cip -> term_trees [i];
			ep2 = cip -> term_trees [i + 1];
			while (ep1 < ep2) {
				fs = *ep1++;
				CLRBIT (fsmask, fs);
			}
		}
		free ((char *) fsmask);
		free ((char *) counts);
		free ((char *) tlist);
	}

	free ((char *) tmask);

	/* Note those rows that are initially present.  Note that we do	*/
	/* not necessarily separate these, so we had better not delete	*/
	/* them from the pool...					*/

	pool -> initrows = pool -> nrows;

	T1 = get_cpu_time ();
	convert_cpu_time (T1 - T0, tbuf);
	tracef (" %% initialize_constraint_pool: %s seconds.\n", tbuf);

	/* Note: Due to the fact that duplicate rows can be produced	*/
	/* (for example by two nearby terminals having the same cutset	*/
	/* constraints), the final number of rows in the pool may be	*/
	/* smaller than the number of logical constraints seeded...	*/

	tracef (" %% Constraint pool initialized with:\n");
	tracef (" %%	%d	Total degree rows	%d	coeffs.\n",
		num_total_degree_rows, num_total_degree_coeffs);
	tracef (" %%	%d	Cutset rows		%d	coeffs.\n",
		num_cutset_rows, num_cutset_coeffs);
	tracef (" %%	%d	Incompatibility rows	%d	coeffs.\n",
		num_incompat_rows, num_incompat_coeffs);
	tracef (" %%	%d	2-terminal SEC rows	%d	coeffs.\n",
		num_2sec_rows, num_2sec_coeffs);
	tracef (" %%	%d	Total rows in pool	%d	in LP\n",
		pool -> nrows, pool -> npend);

	print_pool_memory_usage (pool);
}

/*
 * This routine frees up the constraint pool.
 */

	void
free_constraint_pool (

struct cpool *		pool		/* IN - pool to add constraint to */
)
{
struct rblk *		blkp;
struct rblk *		tmp;

	free ((char *) (pool -> cbuf));
	free ((char *) (pool -> lprows));
	free ((char *) (pool -> rows));

	blkp = pool -> blocks;
	while (blkp NE NULL) {
		tmp = blkp -> next;
		free ((char *) (blkp -> base));
		free ((char *) blkp);
		blkp = tmp;
	}

	free ((char *) pool);
}

/*
 * This routine adds a single constraint to the pool, unless it is
 * already present.  We use the hash table to determine this.
 */

	bool
add_constraint_to_pool (

struct cpool *		pool,		/* IN - pool to add constraint to */
struct rcoef *		rp,		/* IN - raw constraint to add */
bool			add_to_lp	/* IN - add it to LP tableaux also? */
)
{
int		hval;
int		len;
int		var;
int		row;
int		n;