solve.c

生成直角Steiner树的程序包

      printf("%% col_prim:%d, reduced cost: %g\n",
	      colnr, (double)dpiv);
    else
      printf("%% col_prim: no negative reduced costs found, optimality!\n");
  }
  if(colnr == 0) {
    Doiter   = FALSE;
    DoInvert = FALSE;
    Status   = OPTIMAL;
  }
  return(colnr);
} /* colprim */

static int rowprim(lprec *lp,
		   int colnr,
		   REAL *theta,
		   REAL *pcol)
{
  int  i, row_nr;
  REAL f, quot;

  row_nr = 0;
  (*theta) = lp->infinite;
  for(i = 1; i <= lp->rows; i++) {
    f = pcol[i];
    if(f != 0) {
      if(my_abs(f) < Trej) {
	debug_print(lp, "pivot %g rejected, too small (limit %g)\n",
		    (double)f, (double)Trej);
      }
      else { /* pivot alright */
	quot = 2 * lp->infinite;
	if(f > 0)
	  quot = lp->rhs[i] / (REAL) f;
	else if(lp->upbo[lp->bas[i]] < lp->infinite)
	  quot = (lp->rhs[i] - lp->upbo[lp->bas[i]]) / (REAL) f;
	my_round(quot, lp->epsel);
	if(quot < (*theta)) {
	  (*theta) = quot;
	  row_nr = i;
	}
      }
    }
  }
  if(row_nr == 0)
    for(i = 1; i <= lp->rows; i++) {
      f = pcol[i];
      if(f != 0) {
	quot = 2 * lp->infinite;
	if(f > 0)
	  quot = lp->rhs[i] / (REAL) f;
	else
	  if(lp->upbo[lp->bas[i]] < lp->infinite)
	    quot = (lp->rhs[i] - lp->upbo[lp->bas[i]]) / (REAL) f;
	my_round(quot, lp->epsel);
	if(quot < (*theta)) {
	  (*theta) = quot;
	  row_nr = i;
	}
      }
    }

  if((*theta) < 0) {
    printf("%% Warning: Numerical instability, qout = %g\n",
	    (double)(*theta));
    printf("%% pcol[%d] = %18g, rhs[%d] = %18g , upbo = %g\n",
	    row_nr, (double)f, row_nr, (double)lp->rhs[row_nr],
	    (double)lp->upbo[lp->bas[row_nr]]);
  }
  if(row_nr == 0) {
    if(lp->upbo[colnr] == lp->infinite) {
      Doiter   = FALSE;
      DoInvert = FALSE;
      Status   = UNBOUNDED;
    }
    else {
      i = 1;
      while(pcol[i] >= 0 && i <= lp->rows)
	i++;
      if(i > lp->rows) { /* empty column with upperbound! */
	lp->lower[colnr] = FALSE;
	lp->rhs[0] += lp->upbo[colnr]*pcol[0];
	Doiter = FALSE;
	DoInvert = FALSE;
      }
      else if(pcol[i]<0) {
	row_nr = i;
      }
    }
  }
  if(row_nr > 0)
    Doiter = TRUE;
  if(lp->trace)
    printf("%% row_prim:%d, pivot element:%18g\n", row_nr,
	    (double)pcol[row_nr]);

  return(row_nr);
} /* rowprim */

static int rowdual(lprec *lp)
{
  int   i, row_nr;
  REAL  f, g, minrhs;
  short artifs;

  row_nr = 0;
  minrhs = -lp->epsb;
  i = 0;
  artifs = FALSE;
  while(i < lp->rows && !artifs) {
    i++;
    f = lp->upbo[lp->bas[i]];
    if(f == 0 && (lp->rhs[i] != 0)) {
      artifs = TRUE;
      row_nr = i;
    }
    else {
      if(lp->rhs[i] < f - lp->rhs[i])
	g = lp->rhs[i];
      else
	g = f - lp->rhs[i];
      if(g < minrhs) {
	minrhs = g;
	row_nr = i;
      }
    }
  }

  if(lp->trace) {
    if(row_nr > 0) {
      printf("%% row_dual:%d, rhs of selected row:           %18g\n",
	      row_nr, (double)lp->rhs[row_nr]);
      if(lp->upbo[lp->bas[row_nr]] < lp->infinite)
	printf("%%\t\tupper bound of basis variable:    %18g\n",
		(double)lp->upbo[lp->bas[row_nr]]);
    }
    else
      printf("%% row_dual: no infeasibilities found\n");
  }

  return(row_nr);
} /* rowdual */

static int coldual(lprec *lp,
		   int row_nr,
		   short minit,
		   REAL *prow,
		   REAL *drow)
{
  int  i, j, k, r, varnr, *rowp, row, colnr;
  REAL theta, quot, pivot, d, f, g, *valuep, value;

  Doiter = FALSE;
  if(!minit) {
    for(i = 0; i <= lp->rows; i++) {
      prow[i] = 0;
      drow[i] = 0;
    }

    drow[0] = 1;
    prow[row_nr] = 1;

    for(i = lp->eta_size; i >= 1; i--) {
      d = 0;
      f = 0;
      k = lp->eta_col_end[i] - 1;
      r = lp->eta_row_nr[k];
      j = lp->eta_col_end[i - 1];

      /* this is one of the loops where the program consumes a lot of CPU
	 time */
      /* let's help the compiler by doing some pointer arithmetic instead
	 of array indexing */
      for(rowp = lp->eta_row_nr + j, valuep = lp->eta_value + j;
	  j <= k;
	  j++, rowp++, valuep++) {
	f += prow[*rowp] * *valuep;
	d += drow[*rowp] * *valuep;
      }

      my_round(f, lp->epsel);
      prow[r] = f;
      my_round(d, lp->epsd);
      drow[r] = d;
    }

    for(i = 1; i <= lp->columns; i++) {
      varnr = lp->rows + i;
      if(!lp->basis[varnr]) {
	matrec *matentry;

	d = - Extrad * drow[0];
	f = 0;
	k = lp->col_end[i];
	j = lp->col_end[i - 1];

	/* this is one of the loops where the program consumes a lot
	   of cpu time */
	/* let's help the compiler with pointer arithmetic instead
	   of array indexing */
	for(matentry = lp->mat + j;
	    j < k;
	    j++, matentry++) {
	  row = (*matentry).row_nr;
	  value = (*matentry).value;
	  d += drow[row] * value;
	  f += prow[row] * value;
	}

	my_round(f, lp->epsel);
	prow[varnr] = f;
	my_round(d, lp->epsd);
	drow[varnr] = d;
      }
    }
  }

  if(lp->rhs[row_nr] > lp->upbo[lp->bas[row_nr]])
    g = -1;
  else
    g = 1;

  pivot = 0;
  colnr = 0;
  theta = lp->infinite;

  for(i = 1; i <= lp->sum; i++) {
    if(lp->lower[i])
      d = prow[i] * g;
    else
      d = -prow[i] * g;

    if((d < 0) && (!lp->basis[i]) && (lp->upbo[i] > 0)) {
      if(lp->lower[i])
	quot = -drow[i] / (REAL) d;
      else
	quot = drow[i] / (REAL) d;
      if(quot < theta) {
	theta = quot;
	pivot = d;
	colnr = i;
      }
      else if((quot == theta) && (my_abs(d) > my_abs(pivot))) {
	pivot = d;
	colnr = i;
      }
    }
  }

  if(lp->trace)
    printf("%% col_dual:%d, pivot element:  %18g\n", colnr,
	    (double)prow[colnr]);

  if(colnr > 0)
    Doiter = TRUE;

  return(colnr);
} /* coldual */

static void iteration(lprec *lp,
		      int row_nr,
		      int varin,
		      REAL *theta,
		      REAL up,
		      short *minit,
		      short *low,
		      short primal,
                      REAL *pcol)
{
  int  i, k, varout;
  REAL f;
  REAL pivot;

  lp->iter++;

  if(((*minit) = (*theta) > (up + lp->epsb))) {
    (*theta) = up;
    (*low) = !(*low);
  }

  k = lp->eta_col_end[lp->eta_size + 1];
  pivot = lp->eta_value[k - 1];

  for(i = lp->eta_col_end[lp->eta_size]; i < k; i++) {
    f = lp->rhs[lp->eta_row_nr[i]] - (*theta) * lp->eta_value[i];
    my_round(f, lp->epsb);
    lp->rhs[lp->eta_row_nr[i]] = f;
  }

  if(!(*minit)) {
    lp->rhs[row_nr] = (*theta);
    varout = lp->bas[row_nr];
    lp->bas[row_nr] = varin;
    lp->basis[varout] = FALSE;
    lp->basis[varin] = TRUE;

    if(primal && pivot < 0)
      lp->lower[varout] = FALSE;

    if(!(*low) && up < lp->infinite) {
      (*low) = TRUE;
      lp->rhs[row_nr] = up - lp->rhs[row_nr];
      for(i = lp->eta_col_end[lp->eta_size]; i < k; i++)
	lp->eta_value[i] = -lp->eta_value[i];
    }

    addetacol(lp);
    lp->num_inv++;
  }

  if(lp->trace) {
    printf("%% Theta = %g ", (double)(*theta));
    if((*minit)) {
      if(!lp->lower[varin])
	printf("%% Iteration: %d, variable %d changed from 0 to its upper bound of %g\n",
		lp->iter, varin, (double)lp->upbo[varin]);
      else
	printf("%% Iteration: %d, variable %d changed its upper bound of %g to 0\n",
		lp->iter, varin, (double)lp->upbo[varin]);
    }
    else
      printf("%% Iteration: %d, variable %d entered basis at: %g\n",
	      lp->iter, varin, (double)lp->rhs[row_nr]);
    if(!primal) {
      f = 0;
      for(i = 1; i <= lp->rows; i++)
	if(lp->rhs[i] < 0)
	  f -= lp->rhs[i];
	else
	  if(lp->rhs[i] > lp->upbo[lp->bas[i]])
	    f += lp->rhs[i] - lp->upbo[lp->bas[i]];
      printf("%% feasibility gap of this basis: %g\n",
	      (double)f);
    }
    else
      printf("%% objective function value of this feasible basis: %g\n",
	      (double)lp->rhs[0]);
  }
} /* iteration */


static void primloop(lprec *lp)
{
  int	 i;
  REAL   theta, x;
  short  primal;
  REAL   *drow, *prow, *Pcol;
  short  minit;
  int    colnr, row_nr;

  if(lp->trace)
    printf("%% Entering primal algorithm\n");

  CALLOC(drow, lp->sum + 1);
  CALLOC(prow, lp->sum + 1);
  CALLOC(Pcol, lp->rows + 1);

  Status = RUNNING;
  primal = TRUE;
  DoInvert = FALSE;
  Doiter = FALSE;
  Extrad = 0;

  row_nr = 0;
  colnr = 0;

  minit = FALSE;

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

    colnr = colprim(lp, minit, drow);
    if(colnr > 0) {
      setpivcol(lp, colnr, Pcol);

      row_nr = rowprim(lp, colnr, &theta, Pcol);
      if(row_nr > 0)
	condensecol(lp, row_nr, Pcol);
    }

    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) {
      if (!invert(lp)) {
	Status = SINGULAR_BASIS;
	break;
      }
      /* Check whether we are still feasible or not... */
      for(i = 1; i <= lp->rows; i++) {
	x = lp->rhs[i];
	if ((x < -lp->epsb) || (x > lp->upbo[lp->bas[i]] + lp->epsb)) {
	  Status = LOST_PRIMAL_FEASIBILITY;
	  break;
	}
      }
    }
  }

  free(drow);
  free(prow);
  free(Pcol);
} /* primloop */


static void dualloop(lprec *lp)
{
  int    i, j;
  REAL   f, theta;
  short  primal;
  REAL   *drow, *prow, *Pcol;
  short  minit;
  int    colnr, row_nr;

  if(lp->trace)
    printf("%% Entering dual algorithm\n");

  CALLOC(drow, lp->sum + 1);
  CALLOC(prow, lp->sum + 1);
  CALLOC(Pcol, lp->rows + 1);

  Status = RUNNING;
  primal = FALSE;
  DoInvert = FALSE;
  Doiter = FALSE;

  /* Set Extrad to be the most negative of the objective coefficients.	*/
  /* We effectively subtract Extrad from every element of the objective	*/
  /* row, thereby making the entire objective row non-negative.  Note	*/
  /* that this forces dual feasibility!  Although it also alters the	*/
  /* objective function, we don't really care about that too much	*/
  /* because we only use the dual algorithm to obtain a primal feasible	*/
  /* solution that we can start the primal algorithm with.  Virtually	*/
  /* any non-zero objective function will work for this!		*/
  Extrad = 0;
  for(i = 1; i <= lp->columns; i++) {
    f = 0;
    for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
      if(lp->mat[j].row_nr == 0)
	f += lp->mat[j].value;

    if(f < Extrad)
      Extrad = f;
  }

  if(lp->trace)
    printf("%% Extrad = %g\n", (double)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. Catch it and try to recover */
	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");
#if 0
	  printf("%%\tLower[%d] = %d\n"
		 "%%\tprow[%d] = %g\n"
		 "%%\tdrow[%d] = %g\n"
		 "%%\tPcol[%d] = %g\n"
		 "%%\tRhs[%d] = %g\n"
		 "%%\tBas[%d] = %d\n"
		 "%%\tUpbo[%d] = %g\n",
		  colnr, Lower[colnr],
		  colnr, prow[colnr],
		  colnr, drow[colnr],
		  row_nr, Pcol[row_nr],
		  row_nr, Rhs[row_nr],
		  row_nr, Bas[row_nr],
		  Bas[row_nr], Upbo[Bas[row_nr]]);
#if 1
	  for (i = 1; i <= lp->sum; i++)
	    if (prow[i] != 0.0)
	      printf("%%\t\tprow[%d] = %g\n", i, prow[i]);
	  printf ("%%");
	  for (i = 1; i <= lp->sum; i++) {
	    printf(" %d/%g", Lower[i], drow[i]);
	    if ((i % 10) == 0)
	      printf("\n%%");
	  }
	  printf("\n");
#endif
#endif
	  Doiter = FALSE;
	  if(!JustInverted) {
	    printf("%% Trying to recover. Reinverting Eta\n");
	    DoInvert = TRUE;
	  }
	  else {
	    printf("%% Can't reinvert, failure\n");
	    dump_lp (lp, "binary.lp");
	    Status = FAILURE;
	  }
	}
	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 {
      Status   = SWITCH_TO_PRIMAL;
      Doiter   = FALSE;
      Extrad   = 0;
      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) {
      if(!invert(lp)) {
	Status = SINGULAR_BASIS;
	break;
      }
    }
  }

  free(drow);
  free(prow);
  free(Pcol);
}