constrnt.c

生成直角Steiner树的程序包

LP_t *			lp;
double *		slack;
int			nrows;
int			ncols;
int			non_zeros;
int			nslack;
bool			scaling_disabled;
int			small;
int			big;
int			obj_scale;

	(void) pool_iteration;

	lp	= bbip -> lp;

	scaling_disabled = FALSE;

retry_lp:
	/* Solve the current LP instance... */
	status = _MYCPX_dualopt (lp);
	if (status NE 0) {
		tracef (" %%  WARNING dualopt: status = %d\n", status);
	}

	/* Get current LP solution... */
	i = _MYCPX_solution (lp,
			     &status,		/* solution status */
			     &z,		/* objective value */
			     x,			/* solution variables */
			     NULL,		/* IGNORE dual values */
			     bbip -> slack,	/* slack variables */
			     dj);		/* reduced costs */
	if (i NE 0) {
		fprintf (stderr, "err_code = %d\n", i);
		fatal ("solve_single_LP: Bug 1.");
	}

	obj_scale = bbip -> lpmem -> obj_scale;
	ncols	  = _MYCPX_getnumcols (lp);

	if (obj_scale NE 0) {
		/* Unscale CPLEX results. */
		z = ldexp (z, obj_scale);
		for (i = 0; i < ncols; i++) {
			dj [i] = ldexp (dj [i], obj_scale);
		}
	}

	bbip -> node -> z	= z;

	/* Get solution status into solver-independent form... */
	switch (status) {
	case CPX_OPTIMAL:
		status = BBLP_OPTIMAL;
		break;

	case CPX_INFEASIBLE:
	case CPX_UNBOUNDED:	/* (CPLEX sometimes gives this for an	*/
				/* infeasible problem.)			*/
		status = BBLP_INFEASIBLE;
		break;

	case CPX_OBJ_LIM:	/* Objective limit exceeded... */
		status = BBLP_CUTOFF;
		break;

	case CPX_OPTIMAL_INFEAS:
		/* This means that CPLEX scaled the problem, found an	*/
		/* optimal solution, unscaled the solution, but that	*/
		/* the unscaled solution no longer satisfied all of the	*/
		/* bound or row feasibility tolerances (i.e. the the	*/
		/* unscaled solution is no longer feasible).  We fix	*/
		/* This by turning off scaling and trying again.  Note	*/
		/* that this happens very rarely, but that CPLEX runs	*/
		/* much slower with scaling turned off, so we don't	*/
		/* want to leave scaling off if we can help it...	*/
		if (scaling_disabled) {
			/* CPLEX is never supposed to return this code	*/
			/* when scaling is diabled!			*/
			fatal ("solve_single_LP: Bug 2.");
		}

		tracef (" %% TURNING OFF SCALING...\n");

		if (_MYCPX_setscaind (-1, &small, &big) NE 0) {
			fatal ("solve_single_LP: Bug 3.");
		}

		/* Must reload the entire problem for this to take effect! */
		reload_cplex_problem (bbip);
		lp = bbip -> lp;

		scaling_disabled = TRUE;

		goto retry_lp;

	default:
		fprintf (stderr, "Unexpected status = %d\n", status);
		_MYCPX_lpwrite (lp, "core.lp");
		fatal ("solve_single_LP: Bug 4.");
		break;
	}

	if (scaling_disabled) {
		/* Must re-enable scaling, or we'll be really slow! */
		tracef (" %% TURNING ON SCALING...\n");
		if (_MYCPX_setscaind (0, &small, &big) NE 0) {
			fatal ("solve_single_LP: Bug 6.");
		}

		/* Must reload entire problem for this to take affect! */
		reload_cplex_problem (bbip);
		lp = bbip -> lp;
	}

	/* Print info about the LP tableaux we just solved... */
	nrows	  = _MYCPX_getnumrows (lp);
	ncols	  = _MYCPX_getnumcols (lp);
	non_zeros = _MYCPX_getnumnz (lp);
	slack = bbip -> slack;
	nslack = 0;
	for (i = 0; i < nrows; i++) {
		if (slack [i] > FUZZ) {
			++nslack;
		}
	}
	(void) tracef ("  %% @PL %d rows, %d cols, %d nonzeros,"
		       " %d slack, %d tight.\n",
		       nrows, ncols, non_zeros,
		       nslack, nrows - nslack);

	return (status);
}

#endif

/*
 * This routine appends "pool -> npend" new rows onto the end of the
 * LP tableaux from the constraint pool.  The constraint numbers of
 * the actual pool constraints to add reside at the end of the
 * "pool -> lprows []" array (starting with element pool -> nlprows).
 * An additional detail is that for each pool constraint we add, we
 * must record which row of the LP tableaux it now resides in.
 */

#ifdef CPLEX

	void
add_pending_rows_to_LP (

struct bbinfo *		bbip		/* IN - branch and bound info */
)
{
int			i;
int			j;
int			i1;
int			i2;
int			newrows;
int			ncoeff;
int			row;
int			nzi;
int			var;
int			row_space;
int			nz_space;
int			num_nz;
struct rcon *		rcp;
struct rcoef *		cp;
LP_t *			lp;
struct cpool *		pool;
double *		rhs;
char *			sense;
int *			matbeg;
int *			matind;
double *		matval;

	verify_pool (bbip -> cpool);

	lp	= bbip -> lp;
	pool	= bbip -> cpool;

	if (_MYCPX_getnumrows (lp) NE pool -> nlprows) {
		/* LP is out of sync with the pool... */
		fatal ("add_pending_rows_to_LP: Bug 1.");
	}

	/* Get number of rows and non-zeros to add to LP... */
	newrows = pool -> npend;

	if (newrows < 0) {
		fatal ("add_pending_rows_to_LP: Bug 2.");
	}

	if (newrows EQ 0) return;

	i1	= pool -> nlprows;
	i2	= i1 + newrows;

	/* Get number of rows and non-zeros to add to LP... */
	ncoeff = 0;
	for (i = i1; i < i2; i++) {
		row = pool -> lprows [i];
		rcp = &(pool -> rows [row]);
		if (rcp -> lprow NE -2) {
			/* Constraint not pending? */
			fatal ("add_pending_rows_to_LP: Bug 3.");
		}
		rcp -> lprow = i;
		ncoeff += rcp -> len;
	}


	tracef ("  %% @PAP adding %d rows, %d nz to LP\n", newrows, ncoeff);

	/* Check to see if the current CPLEX allocations are	*/
	/* sufficient.  If not, we must reallocate...		*/
	row_space = _MYCPX_getrowspace (lp);
	nz_space  = _MYCPX_getnzspace (lp);
	num_nz	  = _MYCPX_getnumnz (lp);

	/* Update high-water marks... */
	if (i2 > pool -> hwmrow) {
		pool -> hwmrow = i2;
	}
	if (num_nz + ncoeff > pool -> hwmnz) {
		pool -> hwmnz = num_nz + ncoeff;
	}

	if ((i2 > row_space) OR (num_nz + ncoeff > nz_space)) {
		/* We must reallocate!  We do this by throwing away the */
		/* old LP completely and building it again from scratch	*/
		/* using only the info available in the constraint	*/
		/* pool.  Hopefully this way we avoid poor memory	*/
		/* utilization due to fragmentation...			*/

		reload_cplex_problem (bbip);

		return;
	}

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

	/* Put the rows into the format that CPLEX wants them in... */
	nzi = 0;
	j = 0;
	for (i = i1; i < i2; i++, j++) {
		row = pool -> lprows [i];
		rcp = &(pool -> rows [row]);
		matbeg [j] = nzi;
		for (cp = rcp -> coefs; ; cp++) {
			var = cp -> var;
			if (var < RC_VAR_BASE) break;
			matind [nzi] = var - RC_VAR_BASE;
			matval [nzi] = cp -> val;
			++nzi;
		}
		rhs [j] = cp -> val;
		switch (var) {
		case RC_OP_LE:	sense [j] = 'L';	break;
		case RC_OP_EQ:	sense [j] = 'E';	break;
		case RC_OP_GE:	sense [j] = 'G';	break;
		default:
			fatal ("add_pending_rows_to_LP: Bug 4.");
			break;
		}
	}
	matbeg [j] = nzi;
	if (nzi NE ncoeff) {
		fatal ("add_pending_rows_to_LP: Bug 5.");
	}

	i = _MYCPX_addrows (lp,
			    0,
			    newrows,
			    ncoeff,
			    rhs,
			    sense,
			    matbeg,
			    matind,
			    matval,
			    NULL,
			    NULL);

	if (i NE 0) {
		fatal ("add_pending_rows_to_LP: Bug 6.");
	}

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

	verify_pool (bbip -> cpool);

	free ((char *) matval);
	free ((char *) matind);
	free ((char *) matbeg);
	free ((char *) sense);
	free ((char *) rhs);
}

#endif

/*
 * This routine frees the current CPLEX problem, and reallocates/rebuilds
 * it from the current constraint pool.  This routine works even if
 * there are constraints pending addition to the LP tableaux.
 */

#if CPLEX

	static
	void
reload_cplex_problem (

struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
int			j;
int			i1;
int			i2;
int			newrows;
int			row;
int			nedges;
struct rcon *		rcp;
struct rcoef *		cp;
LP_t *			lp;
struct cpool *		pool;
int *			cstat;
int *			rstat;
int *			b_index;
char *			b_lu;
double *		b_bd;

	lp	= bbip -> lp;
	pool	= bbip -> cpool;

	newrows	= pool -> npend;
	i1	= pool -> nlprows;
	i2	= i1 + newrows;

	tracef (" %% REALLOCATING CPLEX PROBLEM...\n");

	/* Save off the current basis, setting the new	*/
	/* rows to be basic...				*/
	cstat = NEWA (bbip -> cip -> num_edges, int);
	rstat = NEWA (i2, int);
	if (_MYCPX_getbase (lp, cstat, rstat) NE 0) {
		fatal ("reload_cplex_problem: Bug 1.");
	}
	for (i = i1; i < i2; i++) {
		/* Set slack variables for new rows to be basic... */
		rstat [i] = 1;
	}

	/* Free up the current LP... */
	destroy_initial_formulation (bbip);

	/* Make all LP rows be pending again... */
	for (i = 0; i < pool -> nlprows; i++) {
		row = pool -> lprows [i];
		rcp = &(pool -> rows [row]);
		if (rcp -> lprow < 0) {
			/* Not currently in LP? */
			fatal ("reload_cplex_problem: Bug 2.");
		}
		rcp -> lprow = -2;	/* is now pending... */
	}
	pool -> npend += pool -> nlprows;
	pool -> nlprows = 0;

	/* Build the initial formulation from scratch again... */
	lp = build_initial_formulation (pool,
					bbip -> tmap,
					bbip -> fset_mask,
					bbip -> cip,
					bbip -> lpmem);
	bbip -> lp = lp;

	/* The initial formulation bounds all variables	*/
	/* from 0 to 1.  We must restore the proper	*/
	/* bounds for all variables that have been	*/
	/* fixed to 0 or 1...				*/

	nedges = bbip -> cip -> num_edges;
	b_index = NEWA (2 * nedges, int);
	b_lu	= NEWA (2 * nedges, char);
	b_bd	= NEWA (2 * nedges, double);
	j = 0;
	for (i = 0; i < nedges; i++) {
		if (NOT BITON (bbip -> fixed, i)) continue;
		b_index [j]	= i;
		b_lu [j]	= 'L';
		b_index [j+1]	= i;
		b_lu [j+1]	= 'U';
		if (NOT BITON (bbip -> value, i)) {
			b_bd [j]	= 0.0;
			b_bd [j+1]	= 0.0;
		}
		else {
			b_bd [j]	= 1.0;
			b_bd [j+1]	= 1.0;
		}
		j += 2;
	}

	if (j > 0) {
		if (_MYCPX_chgbds (lp, j, b_index, b_lu, b_bd) NE 0) {
			fatal ("reload_cplex_problem: Bug 3.");
		}
	}

	free ((char *) b_bd);
	free ((char *) b_lu);
	free ((char *) b_index);

	/* Restore augmented basis... */
	if (_MYCPX_copybase (lp, cstat, rstat) NE 0) {
		fatal ("reload_cplex_problem: Bug 4.");
	}
	free ((char *) rstat);
	free ((char *) cstat);
}

#endif

/*
 * This routine marks a single row as "pending addition to the LP tableaux,"
 * assuming it is not already pending or in the LP.
 */

	void
mark_row_pending_to_LP (

struct cpool *		pool,		/* IN - constraint pool */
int			row		/* IN - row to mark pending */
)
{
int			i;
struct rcon *		rcp;

	if ((row < 0) OR (row >= pool -> nrows)) {
		fatal ("mark_row_pending_to_LP: Bug 1.");
	}
	rcp = &(pool -> rows [row]);
	if ((rcp -> lprow >= 0) OR (rcp -> lprow EQ -2)) {
		/* Row is already in the LP, or was previously	*/
		/* made pending...				*/
		return;
	}
	if (rcp -> lprow NE -1) {
		/* Pool constraint has bad state... */
		fatal ("mark_row_pending_to_LP: Bug 2.");
	}

	/* row is now pending... */
	rcp -> lprow = -2;

	i = pool -> nlprows + (pool -> npend)++;
	pool -> lprows [i] = row;
}

/*
 * This routine adds the given list of LOGICAL constraints to the pool
 * as physical constraints (actual coefficient rows).  Any duplicates
 * that might exist in this process are discarded.  Those that
 * represent violations are also appended to the LP tableaux.
 */

	int
add_constraints (

struct bbinfo *		bbip,		/* IN - branch and bound info */
struct constraint *	lcp		/* IN - list of logical constraints */
)
{
int			nrows;
int			ncoeffs;
struct constraint *	p;
struct cpool *		pool;
struct rcoef *		cp;
bitmap_t *		fset_mask;
double *		x;
struct cinfo *		cip;
bool			add_to_lp;
bool			any_violations;
bool			violation;
bool			newly_added;
int			num_con;

	verify_pool (bbip -> cpool);

	cip		= bbip -> cip;
	x		= bbip -> node -> x;
	pool		= bbip -> cpool;
	fset_mask	= bbip -> fset_mask;

	/* Compute the space requirements for these constraints.  Don't	*/
	/* bother trying to account for duplicates or hits on		*/
	/* constraints already in the pool -- assume they are all	*/
	/* unique.							*/
	ncoeffs = 0;
	nrows = 0;
	for (p = lcp; p NE NULL; p = p -> next) {
		cp = expand_constraint (p, pool -> cbuf, fset_mask, cip);
		ncoeffs += (cp - pool -> cbuf);
		++nrows;
	}

	if (ncoeffs > pool -> blocks -> nfree) {
		/* Must delete some not recently used rows... */
		garbage_collect_pool (pool, ncoeffs, nrows);
	}

	any_violations = FALSE;

	num_con = 0;

	while (lcp NE NULL) {
		expand_constraint (lcp, pool -> cbuf, fset_mask, cip);

		add_to_lp = FALSE;
		violation = is_violation (pool -> cbuf, x);
		if (violation) {
			/* Add-to-lp only does anything if this constraint */
			/* is brand new...				   */
			add_to_lp = TRUE;
			any_violations = TRUE;
		}
		newly_added = add_constraint_to_pool (pool,
						      pool -> cbuf,
						      add_to_lp);
		if (newly_added AND violation) {
			++num_con;
		}

		lcp = lcp -> next;
	}

	if (any_violations) {
		/* Prune back the pending rows so that only the	*/
		/* smallest constraints are added to the LP	*/
		/* the first time around.  The hope is that	*/
		/* this gives a more sparse LP that can be	*/
		/* solved more quickly, and that many of the	*/
		/* larger constraints (that we didn't add in)	*/
		/* will no longer be violated by the new	*/
		/* solution.  In other words, why bog down the	*/
		/* LP solver with lots of dense rows that will	*/
		/* become slack with the slightest perturbation	*/
		/* of the solution?				*/

#if 1
		prune_pending_rows (bbip, FALSE);
#endif

		/* Add the new violations to the LP tableaux... */
		add_pending_rows_to_LP (bbip);
	}

	print_pool_memory_usage (pool);

	return (num_con);
}