constrnt.c

生成直角Steiner树的程序包

struct rcoef *	p;
struct rcoef *	p2;
struct rcon *	rcp;
int *		hookp;
struct rblk *	blkp;
struct rblk *	blkp2;
int *		ip;
size_t		nbytes;

#define	_HASH(reg,value) \
	(reg) ^= (value); \
	(reg) = ((reg) < 0) ? ((reg) << 1) + 1 : ((reg) << 1);

	verify_pool (pool);

	/* Factor out the GCD of the row... */
	reduce_constraint (rp);

	/* Compute hash value and length of LHS... */
	hval = 0;
	len = 0;
	for (p = rp;; p++) {
		var = p -> var;
		if (var < RC_VAR_BASE) break;
		_HASH (hval, var);
		_HASH (hval, p -> val);
		++len;
	}
	hval %= CPOOL_HASH_SIZE;
	if (hval < 0) {
		hval += CPOOL_HASH_SIZE;
	}

	if ((hval < 0) OR (hval >= CPOOL_HASH_SIZE)) {
		fatal ("add_constraint_to_pool: Bug 0.");
	}

	nbytes = (len + 1) * sizeof (*rp);

	hookp = &(pool -> hash [hval]);

	for (row = *hookp; row >= 0;) {
		rcp = &(pool -> rows [row]);
		if ((rcp -> len EQ len) AND
		    (memcmp (rcp -> coefs, rp, nbytes) EQ 0)) {
			/* Constraint already here! */
			return (FALSE);
		}
		row = rcp -> next;
	}

	/* Constraint is not present -- add it.  Start by copying the	*/
	/* coefficients...  If no room, grab another block.		*/
	blkp = pool -> blocks;
	if (blkp EQ NULL) {
		fatal ("add_constraint_to_pool: Bug 2.");
	}
	pool -> num_nz += len;
	++len;		/* include op/rhs in length now... */
	if (blkp -> nfree < len) {
		/* Note: the free space (if any) at the end of the	*/
		/* current block NEVER gets used.  This is by design,	*/
		/* since this results in better cache behavior while	*/
		/* scanning the pool for violations (hitting all of	*/
		/* the rcoefs in each rblk in sequence).		*/

		/* Grab same number as last time... */
		n = (blkp -> ptr - blkp -> base) + blkp -> nfree;
		if (n < len) {
			n = len;
		}
		blkp2 = NEW (struct rblk);
		blkp2 -> next	= blkp;
		blkp2 -> base	= NEWA (n, struct rcoef);
		blkp2 -> ptr	= blkp2 -> base;
		blkp2 -> nfree	= n;
		pool -> blocks = blkp2;
		blkp = blkp2;
	}
	p = blkp -> ptr;
	blkp -> ptr	+= len;
	blkp -> nfree	-= len;
	memcpy (p, rp, len * sizeof (*rp));

	/* Now grab a new row header... */
	row = (pool -> nrows)++;
	if (row >= pool -> maxrows) {
		/* Must grab more rows.  Double it...	*/
		n = 2 * pool -> maxrows;
		rcp = pool -> rows;
		pool -> rows = NEWA (n, struct rcon);
		pool -> maxrows = n;
		memcpy (pool -> rows, rcp, row * sizeof (*rcp));
		free ((char *) rcp);

		/* Increase size of the lprows array also... */
		ip = NEWA (n, int);
		memcpy (ip, pool -> lprows, row * sizeof (*ip));
		free ((char *) (pool -> lprows));
		pool -> lprows = ip;
	}
	rcp = &(pool -> rows [row]);
	rcp -> len	= len - 1;	/* op/rhs not part of length here... */
	rcp -> coefs	= p;
	rcp -> next	= *hookp;
	rcp -> lprow	= -1;
	rcp -> biter	= pool -> iter;	/* assume binding (or violated) now */
	rcp -> hval	= hval;
	rcp -> flags	= 0;
	rcp -> uid	= (pool -> uid)++;
	rcp -> refc	= 0;		/* no OTHER node references it! */
	*hookp = row;

	if (add_to_lp) {
		/* This row is pending addition to the LP tableaux. */
		mark_row_pending_to_LP (pool, row);
	}

	verify_pool (pool);

	return (TRUE);
}

/*
 * This routine reduces the given constraint row to lowest terms by
 * dividing by the GCD.
 */

	static
	void
reduce_constraint (

struct rcoef *		rp	/* IN - coefficient row */
)
{
int			j;
int			k;
int			rem;
int			com_factor;
struct rcoef *		p;

	/* Initial common factor is first coefficient. */
	com_factor = rp -> val;
	if (com_factor <= 0) {
		if (com_factor EQ 0) {
			fatal ("reduce_constraint: Bug 1.");
		}
		com_factor = - com_factor;
	}

	if (com_factor EQ 1) return;

	for (p = rp + 1; ; p++) {
		k = p -> val;
		if (k <= 0) {
			if (k EQ 0) {
				fatal ("reduce_constraint: Bug 2.");
			}
			k = -k;
		}
		/* Euclid's algorithm: computes GCD... */
		j = com_factor;
		while (j > 0) {
			rem = k % j;
			k = j;
			j = rem;
		}
		com_factor = k;
		if (com_factor EQ 1) return;
		if (p -> var < RC_VAR_BASE) break;
	}

	/* We have a row to reduce! */
	for (p = rp; ; p++) {
		if ((p -> val % com_factor) NE 0) {
			fatal ("reduce_constraint: Bug 3.");
		}
		p -> val /= com_factor;
		if (p -> var < RC_VAR_BASE) break;
	}
}

/*
 * This routine sets up the LP problem instance for the initial
 * constraints of the LP relaxation.
 */

#if CPLEX

	LP_t *
build_initial_formulation (

struct cpool *		pool,		/* IN - initial constraint pool */
bitmap_t *		vert_mask,	/* IN - set of valid vertices */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
struct cinfo *		cip,		/* IN - compatibility info */
struct lpmem *		lpmem		/* OUT - dynamically allocated mem */
)
{
int			i, j, k;
int			nedges;
int			nrows;
int			ncoeff;
int			row;
int			var;
int *			tmp;
struct rcon *		rcp;
struct rcoef *		cp;
LP_t *			lp;
int			macsz, marsz, maesz, matsz;
int			mac, mar, mae;
int			objsen;
double *		objx;
double *		rhsx;
char *			senx;
double *		bdl;
double *		bdu;
int *			matbeg;
int *			matcnt;
int *			matind;
double *		matval;
cpu_time_t		T0;
cpu_time_t		T1;
double			min_c, max_c, ci;
int			min_exp, max_exp;
int			obj_scale;
char			tbuf [32];

	T0 = get_cpu_time ();

	nedges = cip -> num_edges;

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

	/* Build the objective function... */
	objx = NEWA (macsz, double);
	for (i = 0; i < macsz; i++) {
		objx [i] = 0.0;
	}
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;
		objx [i] = (double) (cip -> cost [i]);
	}

	/* CPLEX does not behave well if the objective coefficients	*/
	/* have very large magnitudes.  (If so, we often get "unscaled	*/
	/* infeasibility" error codes.)  Therefore, we scale objx here	*/
	/* (by an exact power of two so that the mantissas remain	*/
	/* unchanged).  Determine a power of two that brings the objx	*/
	/* magnitudes into a reasonable range.				*/

	min_c	= DBL_MAX;
	max_c	= 0.0;
	for (i = 0; i < macsz; i++) {
		ci = fabs (objx [i]);
		if (ci EQ 0.0) continue;
		if (ci < min_c) {
			min_c = ci;
		}
		if (ci > max_c) {
			max_c = ci;
		}
	}

	(void) frexp (min_c, &min_exp);
	(void) frexp (max_c, &max_exp);
	obj_scale = (min_exp + max_exp) / 2;

	/* Remember scale factor so we can unscale results. */
	lpmem -> obj_scale = obj_scale;

	obj_scale = - obj_scale;

	for (i = 0; i < macsz; i++) {
		objx [i] = ldexp (objx [i], obj_scale);
	}

	objsen = _MYCPX_MIN;	/* Minimize */

	/* Build variable bound arrays... */
	bdl = NEWA (macsz, double);
	bdu = NEWA (macsz, double);
	for (i = 0; i < macsz; i++) {
		bdl [i] = 0.0;
		bdu [i] = 1.0;
	}

	mar	= pool -> npend;
	if (pool -> hwmrow EQ 0) {
		/* Initial allocation.  Allocate space sufficiently	*/
		/* large that we are unlikely to need to reallocate the	*/
		/* CPLEX problem buffers...				*/

		/* Start with the total number of non-zeros in the	*/
		/* entire constraint pool...				*/
		ncoeff = 0;
		nrows = pool -> nrows;
		for (i = 0; i < nrows; i++) {
			rcp = &(pool -> rows [i]);
			ncoeff += rcp -> len;
		}

		marsz	= 2 * nrows;
		matsz	= 4 * ncoeff;
	}
	else {
		/* Reallocating CPLEX problem.  We want a moderate rate	*/
		/* of growth, but must trade this off against the	*/
		/* frequency of reallocation.  We expand both the rows	*/
		/* and the non-zeros by 25% over the largest need seen	*/
		/* now or previously.					*/
		ncoeff = 0;
		for (i = 0; i < pool -> npend; i++) {
			row = pool -> lprows [i];
			rcp = &(pool -> rows [row]);
			ncoeff += rcp -> len;
		}
		if ((mar > pool -> hwmrow) OR (ncoeff > pool -> hwmnz)) {
			/* high-water marks should be updated before! */
			fatal ("build_initial_formulation: Bug Z.");
		}
		marsz = 5 * pool -> hwmrow / 4;
		matsz = 5 * pool -> hwmnz / 4;
	}

	if (marsz < min_cplex_rows) {
		marsz = min_cplex_rows;
	}
	if (matsz < min_cplex_nzs) {
		matsz = min_cplex_nzs;
	}

	tracef (" %% cpx allocation: %d rows, %d cols, %d nz\n",
		marsz, macsz, matsz);

	/* Allocate arrays for constraint matrix... */
	rhsx = NEWA (marsz, double);
	senx = NEWA (marsz, char);
	matbeg = NEWA (macsz, int);
	matcnt = NEWA (macsz, int);
	matind = NEWA (matsz, int);
	matval = NEWA (matsz, double);

	for (i = 0; i < marsz; i++) {
		rhsx [i] = 0.0;
	}
	for (i = 0; i < macsz; i++) {
		matbeg [i] = 0;
		matcnt [i] = 0;
	}
	for (i = 0; i < matsz; i++) {
		matind [i] = 0;
		matval [i] = 0.0;
	}

	/* Now go through each row k and compute the number of	*/
	/* non-zero coefficients for each variable used...	*/
	tmp = NEWA (macsz, int);
	for (i = 0; i < macsz; i++) {
		tmp [i] = 0;
	}
	for (i = 0; i < pool -> npend; i++) {
		row = pool -> lprows [i];
		rcp = &(pool -> rows [row]);
		for (cp = rcp -> coefs; ; cp++) {
			var = cp -> var;
			if (var < RC_VAR_BASE) break;
			++(tmp [var - RC_VAR_BASE]);
		}
	}

	/* CPLEX wants columns, not rows... */
	j = 0;
	for (i = 0; i < mac; i++) {
		k = tmp [i];
		matbeg [i] = j;
		tmp [i] = j;
		matcnt [i] = k;
		j += k;
	}
	if (j > pool -> hwmnz) {
		pool -> hwmnz = j;
	}
	if (mar > pool -> hwmrow) {
		pool -> hwmrow = mar;
	}
	for (i = 0; i < pool -> npend; i++) {
		row = pool -> lprows [i];
		rcp = &(pool -> rows [row]);
		for (cp = rcp -> coefs; ; cp++) {
			var = cp -> var;
			if (var < RC_VAR_BASE) break;
			j = tmp [var - RC_VAR_BASE];
			matind [j] = i;
			matval [j] = cp -> val;
			++(tmp [var - RC_VAR_BASE]);
		}
		switch (var) {
		case RC_OP_LE:	senx [i] = 'L';		break;
		case RC_OP_EQ:	senx [i] = 'E';		break;
		case RC_OP_GE:	senx [i] = 'G';		break;
		default:
			fatal ("build_initial_formulation: Bug 1.");
		}
		rhsx [i] = cp -> val;
		rcp -> lprow = i;
	}

	/* Verify consistency of what we generated... */
	for (i = 0; i < mac; i++) {
		if (tmp [i] NE matbeg [i] + matcnt [i]) {
			fprintf (stderr,
				 " %% i = %d, tmp = %d, matbeg = %d, matcnt = %d\n",
				 i, tmp [i], matbeg [i], matcnt [i]);
			fatal ("build_initial_formulation: Bug 2.");
		}
	}

	free ((char *) tmp);

	pool -> nlprows	= pool -> npend;
	pool -> npend	= 0;

#if 0
	_MYCPX_setadvind (1);		/* continue from previous basis. */
#endif

	lp = _MYCPX_loadlp ("root",
			    mac,
			    mar,
			    objsen,
			    objx,
			    rhsx,
			    senx,
			    matbeg,
			    matcnt,
			    matind,
			    matval,
			    bdl,
			    bdu,
			    NULL,
			    macsz,
			    marsz,
			    matsz);

	if (lp EQ NULL) {
		fatal ("build_initial_formulation: Bug 3.");
	}

	/* Remember addresses of each buffer for when we need to free them. */
	lpmem -> objx		= objx;
	lpmem -> rhsx		= rhsx;
	lpmem -> senx		= senx;
	lpmem -> matbeg		= matbeg;
	lpmem -> matcnt		= matcnt;
	lpmem -> matind		= matind;
	lpmem -> matval		= matval;
	lpmem -> bdl		= bdl;
	lpmem -> bdu		= bdu;

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

	return (lp);
}

#endif

/*
 * This routine sets up the LP problem instance for the initial
 * constraints of the LP relaxation.
 */

#if LPSOLVE

	LP_t *
build_initial_formulation (

struct cpool *		pool,		/* IN - initial constraint pool */
bitmap_t *		vert_mask,	/* IN - set of valid vertices */
bitmap_t *		edge_mask,	/* IN - set of valid hyperedges */
struct cinfo *		cip,		/* IN - compatibility info */
struct lpmem *		lpmem		/* OUT - dynamically allocated mem */
)
{
int			i;
int			j;
int			k;
int			nedges;
int			nrows;
int			ncols;
int			ncoeff;
int			nzi;
int			row;
int			var;
struct rcon *		rcp;
struct rcoef *		cp;
LP_t *			lp;
double *		rowvec;
double *		rhs;
short *			ctype;
int *			matbeg;
int *			matind;
double *		matval;
int *			ip1;
int *			ip2;
cpu_time_t		T0;
cpu_time_t		T1;
char			tbuf [32];

	T0 = get_cpu_time ();

	nedges = cip -> num_edges;

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

	/* Compute the total number of non-zeros in the INITIAL	*/
	/* constraints of the constraint pool...		*/
	ncoeff = 0;
	nrows = pool -> npend;
	for (i = 0; i < nrows; i++) {
		row = pool -> lprows [i];
		rcp = &(pool -> rows [row]);
		ncoeff += rcp -> len;
	}

	/* Make the initial LP... */
	lp = make_lp (0, ncols);

	/* All variables are 0-1 variables... */
	for (i = 1; i <= nedges; i++) {
		set_bounds (lp, i, 0.0, 1.0);
	}

	/* Minimize */
	set_minim (lp);

	/* Build the objective function... */
	rowvec = NEWA (ncols + 1, double);
	for (i = 0; i <= ncols; i++) {
		rowvec [i] = 0.0;
	}
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (edge_mask, i)) continue;
		rowvec [i + 1] = (double) (cip -> cost [i]);
	}
#if 1
	inc_mat_space (lp, ncols + 1);
#endif
	set_obj_fn (lp, rowvec);

	free ((char *) rowvec);

	/* Allocate arrays for setting the rows... */
	rhs	= NEWA (nrows, double);
	ctype	= NEWA (nrows, short);
	matbeg	= NEWA (nrows + 1, int);
	matind	= NEWA (ncoeff, int);
	matval	= NEWA (ncoeff, double);

	/* Put the rows into the format that LP-solve wants them in...	*/
	nzi = 0;