solve.c

来自「生成直角Steiner树的程序包」· C语言 代码 · 共 2,233 行 · 第 1/4 页

  return(failure);
} /* milpsolve */


int solve(lprec *lp)
{
  int result, i;

  lp->total_iter  = 0;
  lp->max_level   = 1;
  lp->total_nodes = 0;

  if(isvalid(lp)) {
    if(lp->maximise && lp->obj_bound == lp->infinite)
      lp->best_solution[0] = -lp->infinite;
    else if(!lp->maximise && lp->obj_bound == -lp->infinite)
      lp->best_solution[0] = lp->infinite;
    else
      lp->best_solution[0] = lp->obj_bound;

    Level = 0;

    if(!lp->basis_valid) {
      for(i = 0; i <= lp->rows; i++) {
	lp->basis[i] = TRUE;
	lp->bas[i]   = i;
      }

      for(i = lp->rows + 1; i <= lp->sum; i++)
	lp->basis[i] = FALSE;

      for(i = 0; i <= lp->sum; i++)
	lp->lower[i] = TRUE;

      lp->basis_valid = TRUE;
    }

    lp->eta_valid = FALSE;
    Break_bb      = FALSE;
    result        = milpsolve(lp, lp->orig_upbo, lp->orig_lowbo, lp->basis,
			      lp->lower, lp->bas, FALSE);
    return(result);
  }

  /* if we get here, isvalid(lp) failed. I suggest we return FAILURE */
  printf("%% Error, the current LP seems to be invalid\n");
  return(FAILURE);
} /* solve */


/*
 * This routine saves the current basis of the given LP.
 */

	void
save_LP_basis (

lprec *			lp,		/* IN - LP to save basis for */
struct basis_save *	basp		/* OUT - saved basis info */
)
{
int		i;
int		j;
int		varnr;
int		rows;
int		sum;
double		theta;

	if (!lp->row_end_valid)
		set_row_end (lp);

	/* Re-initialize the various versions of things needed	*/
	/* to re-invert...					*/
	memcpy (lp->upbo,	lp->orig_upbo,	(lp->sum + 1) * sizeof (REAL));
	memcpy (lp->lowbo,	lp->orig_lowbo,	(lp->sum + 1) * sizeof (REAL));
	memcpy (lp->rh,		lp->orig_rh,	(lp->rows + 1) * sizeof (REAL));

	/* Process non-zero lower bounds... */
	for (i = 1; i <= lp->columns; i++) {
		varnr = lp->rows + i;
		theta = lp->lowbo[varnr];
		if (theta != 0.0) {
			if (lp->upbo[varnr] < lp->infinite) {
				lp->upbo[varnr] -= theta;
			}
			for (j = lp->col_end[i - 1]; j < lp->col_end[i]; j++) {
				lp->rh[lp->mat[j].row_nr]
					-= theta * lp->mat[j].value;
			}
		}
	}

	/* First things first -- reinvert so that the eta	*/
	/* vectors are as simple as possible before we save...	*/
	invert (lp);

	lp -> eta_valid = TRUE;

	rows	= lp -> rows;
	sum	= lp -> sum;

	/* Allocate buffers to save state into... */
	CALLOC (basp -> bas, rows + 1);
	CALLOC (basp -> basis, sum + 1);
	CALLOC (basp -> lower, sum + 1);
	CALLOC (basp -> rhs, rows + 1);
	CALLOC (basp -> drow, sum + 1);
	CALLOC (basp -> prow, sum + 1);
	CALLOC (basp -> Pcol, rows + 1);

	/* Save the basis state info... */
	basp -> eta_size = lp -> eta_size;
	memcpy (basp -> bas, lp -> bas, (rows + 1) * sizeof (int));
	memcpy (basp -> basis, lp -> basis, (sum + 1) * sizeof (short));
	memcpy (basp -> lower, lp -> lower, (sum + 1) * sizeof (short));
	memcpy (basp -> rhs, lp -> rhs, (rows + 1) * sizeof (REAL));
}

/*
 * Destroy the saved basis info...
 */

	void
destroy_LP_basis (

struct basis_save *	basp		/* IN - basis info to free up */
)
{
	free ((char *) (basp -> Pcol));
	free ((char *) (basp -> prow));
	free ((char *) (basp -> drow));
	free ((char *) (basp -> rhs));
	free ((char *) (basp -> lower));
	free ((char *) (basp -> basis));
	free ((char *) (basp -> bas));
}


/*
 * This routine performs a quick "test run" of branching the given
 * variable in the given direction.  We start from the current basis
 * (which has been saved in basp, and is assumed optimal), and do at
 * most 50 or so iterations -- but NEVER reinverting, which would
 * clobber the saved eta vectors...  This means we can restore the
 * eta just by popping back to the initial eta size.
 */

	REAL
try_branch (

lprec *			lp,	/* IN - LP to test branching for */
int			var,	/* IN - variable to branch */
int			dir,	/* IN - direction to branch in, 0 or 1 */
REAL *			x,	/* OUT - LP solution found */
REAL			ival,	/* IN - value to yield if infeasible */
struct basis_save *	basp	/* IN - basis to restore when done */
)
{
int		i;
int		j;
int		varnr;
REAL		f;
REAL		theta;
short		primal;
REAL *		drow;
REAL *		prow;
REAL *		Pcol;
REAL		save_lb;
REAL		save_ub;
short		minit;
int		colnr;
int		row_nr;

	lp -> total_iter	= 0;
	lp -> max_level		= 1;
	lp -> total_nodes	= 0;

	if ((var < 1) ||
	    (var > lp -> columns) ||
	    !isvalid (lp) ||
	    (lp -> scaling_used) ||
	    !(lp -> basis_valid) ||
	    !(lp -> eta_valid) ||
	    !(lp -> basis [lp->rows + var])) {
		fprintf (stderr, "Bad call to try_branch...\n");
		abort ();
	}

	Level = 1;

	drow	= basp -> drow;
	prow	= basp -> prow;
	Pcol	= basp -> Pcol;

	lp -> iter = 0;
	minit = FALSE;
	Status = RUNNING;
	DoInvert = FALSE;
	Doiter = FALSE;

	/* We assume the caller starts us from an optimal basis, from	*/
	/* which we branch (up or down), creating an infeasibility.	*/
	/* Thus after branching we always use the Dual Simplex...	*/

	varnr = lp -> rows + var;
	save_lb = lp -> lowbo [varnr];
	save_ub = lp -> upbo [varnr];
	if (dir == 0) {
		/* Branching down -- fix upper bound and continue	*/
		/* optimizing from the same tableaux...			*/
		lp -> upbo [varnr] = 0.0;
	}
	else {
		/* Branching up -- adjust rhs.  The branching variable	*/
		/* MUST be basis!  (We verified this above.)  Therefore	*/
		/* there is only a single RHS value that we need to	*/
		/* adjust because this variables column is all zeros	*/
		/* except for a one in the row in which it is basic.	*/
		lp -> lowbo [varnr] = 1.0;
		lp -> upbo [varnr] = 0.0;
		for (i = 1; ; i++) {
			if (i > lp -> rows) {
				fprintf (stderr, "try_branch: Bug 1.\n");
				abort ();
			}
			if (lp -> bas [i] == varnr) {
				lp -> rhs [i] -= 1.0;
				break;
			}
		}
		/* Fix up objective function value also! */
		for (j = lp -> col_end [var - 1]; j < lp -> col_end [var]; j++) {
			if (lp -> mat [j].row_nr == 0) {
				lp -> rhs [0] -= lp -> mat [j].value;
				break;
			}
		}
	}

	primal = FALSE;

	/* Get most negative INITIAL objective coefficient into	*/
	/* Extrad...						*/
	drow [0] = 1;
	for (i = 1; i <= lp -> rows; i++) {
		drow [i] = 0;
	}
	Extrad = 0;
	for (i = 1; i <= lp -> columns; i++) {
		varnr = lp -> rows + i;
		drow [varnr] = 0;
		for (j = lp -> col_end [i - 1]; j < lp -> col_end [i]; j++) {
			if (drow [lp -> mat[j].row_nr] != 0) {
				drow [varnr] +=
					drow [lp -> mat [j].row_nr]
					* lp -> mat [j].value;
			}
		}
		if (drow [varnr] < Extrad) {
			Extrad = drow [varnr];
		}
	}
	if (lp -> trace) {
		printf ("%% Extrad = %f\n", Extrad);
	}

	row_nr = 0;
	colnr = 0;

	minit = FALSE;

	while (Status == RUNNING) {
		Doiter = FALSE;
		DoInvert = FALSE;

		if (!minit) {
			row_nr = rowdual (lp);
		}
		if (row_nr > 0 ) {
			colnr = coldual (lp, row_nr, minit, prow, drow);
			if (colnr > 0) {
				setpivcol (lp, colnr, Pcol);
				/* getting div by zero here ... MB */
				if (Pcol [row_nr] == 0) {
					printf ("%% An attempt was made to divide by zero (Pcol[%d])\n",
						 row_nr);
					printf ("%% This indicates numerical instability\n");
					Doiter = FALSE;
					Status = STOP_AT_INVERT;
				}
				else {
					condensecol (lp, row_nr, Pcol);
					f = lp -> rhs [row_nr]
						- lp -> upbo [lp -> bas [row_nr]];
					if (f > 0) {
						theta = f / (REAL) Pcol [row_nr];
						if (theta <= lp -> upbo [colnr] + lp -> epsb) {
							lp -> lower [lp -> bas [row_nr]]
								= !lp -> lower [lp -> bas [row_nr]];
						}
					}
					else { /* f <= 0 */
						theta = lp -> rhs [row_nr] / (REAL) Pcol [row_nr];
					}
				}
			}
			else {
				Status = INFEASIBLE;
			}
		}
		else {
			/* Assume optimal achieved... */
			Doiter   = FALSE;
			DoInvert = TRUE;
		}
		if (Doiter) {
			iteration (lp,
				   row_nr,
				   colnr,
				   &theta,
				   lp -> upbo [colnr],
				   &minit,
				   &(lp -> lower [colnr]),
				   primal,
				   Pcol);
		}
		if (lp -> num_inv >= lp -> max_num_inv) {
			DoInvert = TRUE;
		}
		if (DoInvert) {
			/* We DO NOT invert!  That would mess up the	*/
			/* eta vectors that we want to restore back to!	*/
			/* Instead we just stop here and we are done!	*/
			Status = STOP_AT_INVERT;
		}
	}

	lp -> total_iter += lp -> iter;

	construct_solution (lp);
	memcpy (x,
		&(lp -> solution [lp -> rows + 1]),
		lp -> columns * sizeof (REAL));

	/* Restore the LP state back to the saved basis... */
	varnr = lp -> rows + var;
	lp -> lowbo [varnr] = save_lb;
	lp -> upbo [varnr] = save_ub;
	lp -> eta_size = basp -> eta_size;

	memcpy (lp -> bas, basp -> bas, (lp -> rows + 1) * sizeof (int));
	memcpy (lp -> basis, basp -> basis, (lp -> sum + 1) * sizeof (short));
	memcpy (lp -> lower, basp -> lower, (lp -> sum + 1) * sizeof (short));
	memcpy (lp -> rhs, basp -> rhs, (lp -> rows + 1) * sizeof (REAL));

	if (Status == INFEASIBLE) {
		return (ival);
	}

	return (lp -> solution [0]);
}

int lag_solve(lprec *lp, REAL start_bound, int num_iter, short verbose)
{
  int i, j, result, citer;
  short status, OrigFeas, AnyFeas, same_basis;
  REAL *OrigObj, *ModObj, *SubGrad, *BestFeasSol;
  REAL Zub, Zlb, Ztmp, pie;
  REAL rhsmod, Step, SqrsumSubGrad;
  int   *old_bas;
  short *old_lower;

  /* allocate mem */
  MALLOC(OrigObj, lp->columns + 1);
  CALLOC(ModObj, lp->columns + 1);
  CALLOC(SubGrad, lp->nr_lagrange);
  CALLOC(BestFeasSol, lp->sum + 1);
  MALLOCCPY(old_bas, lp->bas, lp->rows + 1);
  MALLOCCPY(old_lower, lp->lower, lp->sum + 1);

  get_row(lp, 0, OrigObj);

  pie = 2;

  if(lp->maximise) {
    Zub = DEF_INFINITE;
    Zlb = start_bound;
  }
  else {
    Zlb = -DEF_INFINITE;
    Zub = start_bound;
  }
  status   = RUNNING;
  Step     = 1;
  OrigFeas = FALSE;
  AnyFeas  = FALSE;
  citer    = 0;

  for(i = 0 ; i < lp->nr_lagrange; i++)
    lp->lambda[i] = 0;

  while(status == RUNNING) {
    citer++;

    for(i = 1; i <= lp->columns; i++) {
      ModObj[i] = OrigObj[i];
      for(j = 0; j < lp->nr_lagrange; j++) {
	if(lp->maximise)
	  ModObj[i] -= lp->lambda[j] * lp->lag_row[j][i];
	else
	  ModObj[i] += lp->lambda[j] * lp->lag_row[j][i];
      }
    }
    for(i = 1; i <= lp->columns; i++) {
      set_mat(lp, 0, i, ModObj[i]);
    }
    rhsmod = 0;
    for(i = 0; i < lp->nr_lagrange; i++)
      if(lp->maximise)
	rhsmod += lp->lambda[i] * lp->lag_rhs[i];
      else
	rhsmod -= lp->lambda[i] * lp->lag_rhs[i];

    if(verbose) {
      printf("%% Zub: %10g Zlb: %10g Step: %10g pie: %10g Feas %d\n",
	      (double)Zub, (double)Zlb, (double)Step, (double)pie, OrigFeas);
      for(i = 0; i < lp->nr_lagrange; i++)
	printf("%% %3d SubGrad %10g lambda %10g\n", i,
		(double)SubGrad[i], (double)lp->lambda[i]);
    }

    if(verbose && lp->sum < 20)
      print_lp(lp);

    result = solve(lp);

    if(verbose && lp->sum < 20) {
      print_solution(lp);
    }

    same_basis = TRUE;
    i = 1;
    while(same_basis && i < lp->rows) {
      same_basis = (old_bas[i] == lp->bas[i]);
      i++;
    }
    i = 1;
    while(same_basis && i < lp->sum) {
      same_basis=(old_lower[i] == lp->lower[i]);
      i++;
    }
    if(!same_basis) {
      memcpy(old_lower, lp->lower, (lp->sum+1) * sizeof(short));
      memcpy(old_bas, lp->bas, (lp->rows+1) * sizeof(int));
      pie *= 0.95;
    }

    if(verbose)
      printf("%% result: %d  same basis: %d\n", result, same_basis);

    if(result == UNBOUNDED) {
      for(i = 1; i <= lp->columns; i++)
	printf("%g ", (double)ModObj[i]);
      exit(EXIT_FAILURE);
    }

    if(result == FAILURE)
      status = FAILURE;

    if(result == INFEASIBLE)
      status = INFEASIBLE;

    SqrsumSubGrad = 0;
    for(i = 0; i < lp->nr_lagrange; i++) {
      SubGrad[i]= -lp->lag_rhs[i];
      for(j = 1; j <= lp->columns; j++)
	SubGrad[i] += lp->best_solution[lp->rows + j] * lp->lag_row[i][j];
      SqrsumSubGrad += SubGrad[i] * SubGrad[i];
    }

    OrigFeas = TRUE;
    for(i = 0; i < lp->nr_lagrange; i++)
      if(lp->lag_con_type[i]) {
	if(my_abs(SubGrad[i]) > lp->epsb)
	  OrigFeas = FALSE;
      }
      else if(SubGrad[i] > lp->epsb)
	OrigFeas = FALSE;

    if(OrigFeas) {
      AnyFeas = TRUE;
      Ztmp = 0;
      for(i = 1; i <= lp->columns; i++)
	Ztmp += lp->best_solution[lp->rows + i] * OrigObj[i];
      if((lp->maximise) && (Ztmp > Zlb)) {
	Zlb = Ztmp;
	for(i = 1; i <= lp->sum; i++)
	  BestFeasSol[i] = lp->best_solution[i];
	BestFeasSol[0] = Zlb;
	if(verbose)
	  printf("%% Best feasible solution: %g\n", (double)Zlb);
      }
      else if(Ztmp < Zub) {
	Zub = Ztmp;
	for(i = 1; i <= lp->sum; i++)
	  BestFeasSol[i] = lp->best_solution[i];
	BestFeasSol[0] = Zub;
	if(verbose)
	  printf("%% Best feasible solution: %g\n", (double)Zub);
      }
        }

    if(lp->maximise)
      Zub = my_min(Zub, rhsmod + lp->best_solution[0]);
    else
      Zlb = my_max(Zlb, rhsmod + lp->best_solution[0]);

    if(my_abs(Zub-Zlb)<0.001) {
      status = OPTIMAL;
    }
    Step = pie * ((1.05*Zub) - Zlb) / SqrsumSubGrad;

    for(i = 0; i < lp->nr_lagrange; i++) {
      lp->lambda[i] += Step * SubGrad[i];
      if(!lp->lag_con_type[i] && lp->lambda[i] < 0)
	lp->lambda[i] = 0;
    }

    if(citer == num_iter && status==RUNNING) {
      if(AnyFeas)
	status = FEAS_FOUND;
      else
	status = NO_FEAS_FOUND;
    }
  }

  for(i = 0; i <= lp->sum; i++)
    lp->best_solution[i] = BestFeasSol[i];

  for(i = 1; i <= lp->columns; i++)
    set_mat(lp, 0, i, OrigObj[i]);

  if(lp->maximise)
    lp->lag_bound = Zub;
  else
    lp->lag_bound = Zlb;
  free(BestFeasSol);
  free(SubGrad);
  free(OrigObj);
  free(ModObj);
  free(old_bas);
  free(old_lower);

  return(status);
}