constrnt.c

生成直角Steiner树的程序包

	for (i = 0; i < nrows; i++) {
		row = pool -> lprows [i];
		rcp = &(pool -> rows [row]);
		matbeg [i] = 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 [i] = cp -> val;
		switch (var) {
		case RC_OP_LE:	ctype [i] = REL_LE;	break;
		case RC_OP_EQ:	ctype [i] = REL_EQ;	break;
		case RC_OP_GE:	ctype [i] = REL_GE;	break;
		default:
			fatal ("build_initial_formulation: Bug 1.");
			break;
		}
		rcp -> lprow = i;
	}
	matbeg [i] = nzi;
	if (nzi NE ncoeff) {
		fatal ("build_initial_formulation: Bug 2.");
	}

	if (nrows > pool -> hwmrow) {
		pool -> hwmrow = nrows;
	}
	if (nzi > pool -> hwmnz) {
		pool -> hwmnz = nzi;
	}

	add_rows (lp, 0, nrows, rhs, ctype, matbeg, matind, matval);

	free ((char *) matval);
	free ((char *) matind);
	free ((char *) matbeg);
	free ((char *) ctype);
	free ((char *) rhs);

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

	verify_pool (pool);

#if 0
	{ FILE *	fp;
		fp = fopen ("main.lp", "w");
		if (fp EQ NULL) {
			fatal ("main.lp -- cannot create");
		}
		write_LP (lp, fp);
		fclose (fp);
	}
#endif

	if (Use_Perturbations) {
		/* Turn on perturbations to deal with degeneracy... */
		lp -> anti_degen = TRUE;
	}
	if (Use_Scaling) {
		/* Turn on auto-scaling of the matrix... */
		auto_scale (lp);
	}

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

	return (lp);
}

#endif

/*
 * This routine solves the current LP relaxation over all constraints
 * currently residing in the constraint pool, regardless of how many
 * are actually in the current LP tableaux.  Each time it solves the
 * LP tableaux, it scans the entire constraint pool for violations.
 * Every violation that is found is appended to the tableaux and we
 * loop back to re-optimize the tableaux.  This procedure terminates
 * only when all constraints in the pool have been satisfied, or a
 * cutoff or infeasibility is encountered.
 *
 * Note also that this procedure NEVER deletes any constraints from
 * the tableaux, slack or otherwise.  Other code must do this.
 */

	int
solve_LP_over_constraint_pool (

struct bbinfo *		bbip		/* IN - branch and bound info */
)
{
int			i;
int			k;
int			status;
int			var;
int			ncols;
int			nrows;
struct rcon *		rcp;
struct rcoef *		cp;
double			sum;
LP_t *			lp;
struct bbnode *		nodep;
double *		x;
double *		dj;
struct cpool *		pool;
bool			any_violations;
bool			can_delete_slack;
int			pool_iteration;
double			slack;
double			prev_z;

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

	ncols	= GET_LP_NUM_COLS (lp);
	nrows	= GET_LP_NUM_ROWS (lp);

	if (nodep -> cpiter EQ pool -> uid) {
		/* nodep -> x is already the optimal solution	*/
		/* over this constraint pool.			*/
#if 1
		tracef ("%%	Constraint pool unchanged, skip LP solve.\n");
#endif
		/* Global solution info valid for this node... */

		/* Reallocate slack variables vector, if necessary... */
		if (bbip -> slack_size < pool -> nlprows) {
			if (bbip -> slack NE NULL) {
				free ((char *) (bbip -> slack));
			}
			bbip -> slack = NEWA (pool -> nlprows, double);
			bbip -> slack_size = pool -> nlprows;
		}

		for (i = 0; i < nrows; i++) {
			bbip -> slack [i] = 0.0;
		}
		return (BBLP_OPTIMAL);
	}

	x  = NEWA (2 * ncols, double);
	dj = x + ncols;

	pool_iteration = 0;

	for (;;) {
		prev_z = nodep -> z;

		verify_pool (bbip -> cpool);

		/* Reallocate slack variables vector, if necessary... */
		if (bbip -> slack_size < pool -> nlprows) {
			if (bbip -> slack NE NULL) {
				free ((char *) (bbip -> slack));
			}
			bbip -> slack = NEWA (pool -> nlprows, double);
			bbip -> slack_size = pool -> nlprows;
		}

		status = solve_single_LP (bbip, x, dj, pool_iteration);

		/* Record another "real" LP solved... */
		do {
			++(pool -> iter);
		} while (pool -> iter EQ -1);

		++pool_iteration;

		if (status NE BBLP_OPTIMAL) break;

		update_lp_solution_history (x, dj, bbip);

#if 0
		if (nodep -> z >= 1.0001 * prev_z) {
			/* Objective rose enough to go ahead	*/
			/* and delete slack rows.  This helps	*/
			/* keep memory usage down on VERY LARGE	*/
			/* problems...				*/
			delete_slack_rows_from_LP (bbip);
		}
#endif

		verify_pool (bbip -> cpool);

		/* Scan entire pool for violations... */
		rcp = &(pool -> rows [0]);
		any_violations = FALSE;
		for (i = 0; i < pool -> nrows; i++, rcp++) {
			slack = compute_slack_value (rcp -> coefs, x);
			if (slack > FUZZ) {
				/* Row is not binding, much less violated. */
				continue;
			}
			/* Consider this row to be binding now! */
			rcp -> biter = pool -> iter;
			if (rcp -> lprow >= 0) {
				/* Skip this row -- it is already in	*/
				/* the LP tableaux.			*/
				continue;
			}
			if (slack < -FUZZ) {
				/* Constraint "i" is not currently in	*/
				/* the LP tableaux, and is violated.	*/
				/* Add it.				*/
				mark_row_pending_to_LP (pool, i);
				any_violations = TRUE;
			}
		}

		/* Done if no violations were appended... */
		if (NOT any_violations) break;

		can_delete_slack = (nodep -> z >= prev_z + 0.0001 * fabs (prev_z));

		prune_pending_rows (bbip, can_delete_slack);

		/* Time to append these pool constraints to the	*/
		/* current LP tableaux...			*/
		add_pending_rows_to_LP (bbip);
	}

	free ((char *) x);

	if (status EQ BBLP_OPTIMAL) {
		/* Nodep -> x is optimal for the current version of the	*/
		/* constraint pool.  Skip re-solve if we re-enter this	*/
		/* routine with the same constraint pool.		*/
		nodep -> cpiter = pool -> uid;
	}
	else {
		/* Not optimal -- force re-solve next time so that	*/
		/* correct status is seen next time we're called.  Yes,	*/
		/* it would have been better if we also saved the LP	*/
		/* status in the node...				*/
		nodep -> cpiter = -1;
	}

	verify_pool (bbip -> cpool);

	return (status);
}

/*
 * This routine copies the current LP solution into the node's buffer,
 * and updates the branch heuristic values.
 */

	static
	void
update_lp_solution_history (

double *		srcx,		/* IN - source LP solution */
double *		dj,		/* IN - source reduced costs */
struct bbinfo *		bbip		/* IN - branch-and-bound info */
)
{
int			i;
int			j;
int			nedges;
int			dir;
int			dir2;
struct bbnode *		nodep;
double *		dstx;
double *		bheur;
double *		zlb;
double			lb;
double			z;

	nodep	= bbip -> node;
	nedges	= bbip -> cip -> num_edges;

	dstx	= nodep -> x;
	bheur	= nodep -> bheur;

	/* Update the branch heuristic value for each variable.		*/
	/* Variables that have been "stuck" at one value for a long	*/
	/* time receive low bheur values.  If fractional, these tend	*/
	/* to be good branch variables.					*/

	if ((nodep -> num EQ 0) AND (nodep -> iter EQ 0)) {
		/* First iteration: set initial values. */
		for (i = 0; i < nedges; i++) {
			dstx [i] = srcx [i];
			bheur [i] = 0;
		}
	}
	else {
		/* Susequent iterations: use time-decayed average. */
		for (i = 0; i < nedges; i++) {
			bheur [i] = 0.75 * bheur [i] + fabs (srcx [i] - dstx [i]);
			dstx [i] = srcx [i];
		}
	}

	/* Now update the Z lower bounds for each variable	*/
	/* using reduced costs.					*/

	zlb = nodep -> zlb;
	z   = nodep -> z;

	for (j = 0; j < nedges; j++) {
		lb = z + fabs (dj [j]);
		dir = (srcx [j] < 0.5);
		dir2 = 1 - dir;
		i = 2 * j;
		if (lb > zlb [i + dir]) {
			zlb [i + dir] = lb;
		}
		if (z > zlb [i + dir2]) {
			zlb [i + dir2] = z;
		}
	}
}

/*
 * This routine solves a single LP tableaux using LP_SOLVE.
 */

#ifdef LPSOLVE

	static
	int
solve_single_LP (

struct bbinfo *		bbip,		/* IN - branch and bound info */
double *		x,		/* OUT - LP solution variables */
double *		dj,		/* OUT - LP reduced costs */
int			pool_iteration	/* IN - pool iteration number */
)
{
int			i;
int			status;
double			z;
LP_t *			lp;
double *		slack;
int			nslack;
double *		djbuf;
struct cinfo *		cip;

	(void) pool_iteration;

	verify_pool (bbip -> cpool);

	cip	= bbip -> cip;
	lp	= bbip -> lp;

#if 0
	/* Debug code to dump each LP instance attempted... */
	{ char		buf [64];
		sprintf (buf,
			 "Node.%03d.%03d.%03d.lp",
			 bbip -> node -> num,
			 bbip -> node -> iter + 1,
			 pool_iteration);
		dump_lp (lp, buf);
	}
#endif

	/* Solve the current LP instance... */
	status = solve (lp);

	/* Get current LP solution... */
	z = lp -> best_solution [0];
	for (i = 0; i < cip -> num_edges; i++) {
		x [i] = lp -> best_solution [lp -> rows + i + 1];
	}

	bbip -> node -> z	= z;

	/* Get solution status into solver-independent form... */
	switch (status) {
	case OPTIMAL:		status = BBLP_OPTIMAL;		break;
	case MILP_FAIL:		status = BBLP_CUTOFF;		break;
	case INFEASIBLE:	status = BBLP_INFEASIBLE;	break;
	default:
		tracef ("solve status = %d\n", status);
		fatal ("solve_single_LP: Bug 1.");
	}

	/* Grab the reduced costs, for future reference. */
	djbuf = NEWA (lp -> sum + 1, double);
	get_reduced_costs (lp, djbuf);
	memcpy (dj,
		&djbuf [lp -> rows + 1],
		lp -> columns * sizeof (double));
	free ((char *) djbuf);

	/* Grab the values of the slack variables, for future reference. */
	slack = NEWA (lp -> rows + 1, double);
	get_slack_vars (lp, slack);
	memcpy (bbip -> slack, &slack [1], lp -> rows * sizeof (double));
	free ((char *) slack);

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

	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;
struct rcon *		rcp;
struct rcoef *		cp;
LP_t *			lp;
struct cpool *		pool;
double *		rhs;
short *			ctype;
int *			matbeg;
int *			matind;
double *		matval;

	verify_pool (bbip -> cpool);

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

	if (lp -> rows 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;
	}

	if (i2 > pool -> hwmrow) {
		pool -> hwmrow = i2;
	}
	if (lp -> non_zeros + ncoeff > pool -> hwmnz) {
		pool -> hwmnz = lp -> non_zeros + ncoeff;
	}

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

	/* Put the rows into the format that LP-solve 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:	ctype [j] = REL_LE;	break;
		case RC_OP_EQ:	ctype [j] = REL_EQ;	break;
		case RC_OP_GE:	ctype [j] = REL_GE;	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.");
	}

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

	add_rows (lp, 0, newrows, rhs, ctype, matbeg, matind, matval);

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

	verify_pool (bbip -> cpool);

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

#endif

/*
 * This routine solves a single LP tableaux using CPLEX.
 */

#ifdef CPLEX

	static
	int
solve_single_LP (

struct bbinfo *		bbip,		/* IN - branch and bound info */
double *		x,		/* OUT - LP solution variables */
double *		dj,		/* OUT - LP reduced costs */
int			pool_iteration	/* IN - pool iteration number */
)
{
int			i;
int			status;
double			z;