solve.c

利用c语言编写

	else if(lp->upbo[lp->bas[i]] < lp->infinite)
	  quot = (lp->rhs[i] - lp->upbo[lp->bas[i]]) / f;
        else {
	  savef = f;
          quot = 2 * lp->infinite;
	}
        my_round(quot, lp->epsel);
        if(quot < (*theta))
          if(quot >= 0)  /* Added by KE 19052002 */
        {
          (*theta) = quot;
          row_nr = i;
          if(lp->piv_rule == FIRST_SELECT) break;
        }
      }
    }
  }

  /* No pivot greater than epspivot was found; accept a smaller one */
  if(row_nr == 0)
    for(quot = 0, i = 1; i <= lp->rows; i++) {
      f = pcol[i];
      if(f != 0) {
	if(f > 0)
	  quot = lp->rhs[i] / f;
	else if(lp->upbo[lp->bas[i]] < lp->infinite)
	  quot = (lp->rhs[i] - lp->upbo[lp->bas[i]]) / f;
        else {
	  savef = f;
          quot = 2 * lp->infinite;
	}
        my_round(quot, lp->epsel);
        if(quot < (*theta))
          if(quot >= 0)    /* Added by KE 19052002 */
        {
          (*theta) = quot;
          row_nr = i;
          if(lp->piv_rule == FIRST_SELECT) break;
        }
      }
    }

  if(row_nr == 0) {
    if(lp->upbo[colnr] == lp->infinite) {
      lp->doIterate = FALSE;
      lp->doInvert = FALSE;
      lp->spx_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];
        lp->doIterate = FALSE;
        lp->doInvert = FALSE;
      }
      else /* if(pcol[i]<0) */ {
	row_nr = i;
      }
    }
  }
  else
/*  if((*theta) < 0) { */
    if((*theta) >= lp->infinite) {   /* Added by KE 19052002 */
      (*theta) = -1;                 /* Added by KE 19052002 */
      report(lp, SEVERE, "Warning: Numerical instability, quot = %g",
	      (double)(*theta));
      report(lp, IMPORTANT, "pcol[%d] = %18g, rhs[%d] = %18g , upbo = %g",
	      row_nr, (double)savef, row_nr, (double)lp->rhs[row_nr],
	      (double)lp->upbo[lp->bas[row_nr]]);
    }

  if(row_nr > 0)
    lp->doIterate = TRUE;
  if(lp->trace)
    report(lp, NORMAL, "row_prim:%d, pivot element:%18g", row_nr,
	    (double)pcol[row_nr]);

  return(row_nr);
} /* rowprim */

static int rowdual(lprec *lp)
{
  int   i;
  int   row_nr;
  RREAL f, g, minrhs;
  MYBOOL  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) {
      report(lp, NORMAL,
	      "row_dual:%d, rhs of selected row:           %18g",
	      row_nr, (double)lp->rhs[row_nr]);
      if(lp->upbo[lp->bas[row_nr]] < lp->infinite)
	report(lp, NORMAL,
		"\t\tupper bound of basis variable:    %18g",
		(double)lp->upbo[lp->bas[row_nr]]);
    }
    else
      report(lp, FULL, "row_dual: no infeasibilities found");
  }

  return(row_nr);
} /* rowdual */

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

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

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

    /* A double BTRAN equation solver process is implemented "in-line" below in
       order to save time and to implement different rounding for the two */
/*     btran(lp, drow, lp->epsd);
     btran(lp, prow, lp->epsel); */

    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] = (REAL) f;
      my_round(d, lp->epsd);
      drow[r] = (REAL) d;
    }

    /* Multiply solution vectors with matrix values */
    for(i = 1; i <= lp->columns; i++) {
      varnr = lp->rows + i;
      if(!lp->basis[varnr]) {
	matrec *matentry;

	d = - lp->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;
	  d += drow[row] * (*matentry).value;
	  f += prow[row] * (*matentry).value;
	}

	my_round(f, lp->epsel);
	prow[varnr] = (REAL) f;
	my_round(d, lp->epsd);
	drow[varnr] = (REAL) 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] / d;
      else
	quot = drow[i] / 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)
    report(lp, NORMAL, "coldual:%d, pivot element:  %18g", colnr,
	    (double)prow[colnr]);

  if(colnr > 0)
    lp->doIterate = TRUE;

  return(colnr);
} /* coldual */

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

  if(yieldformessages(lp)!=0)
    lp->spx_status = USERABORT;

  lp->iter++;
  if((lp->usermessage != NULL) && (lp->msgmask & MSG_ITERATION))
    lp->usermessage(lp, lp->msghandle, MSG_ITERATION);

  if(lp->spx_status != RUNNING)
    return;

  if(((*minit) = (MYBOOL) ((*theta) > (up + lp->epsb)))) {
    (*theta) = up;
    (*low) = (MYBOOL) !(*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++) {
    varout = lp->eta_row_nr[i];
    f = lp->rhs[varout] - (*theta) * lp->eta_value[i];
/*    my_round(f, lp->epsb); */ /* **** Could a large value here actually introduce errors? */
    my_round(f, lp->epsel);
    lp->rhs[varout] = (REAL) f;
  }

  if(!(*minit)) {
    lp->rhs[row_nr] = (REAL) (*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, 0);
    lp->num_inv++;
  }

  if(lp->trace) {
    report(lp, NORMAL, "Theta = %g", (double)(*theta));
    if((*minit)) {
      if(!lp->lower[varin])
	report(lp, NORMAL,
		"Iteration: %d, variable %d changed from 0 to its upper bound of %g",
		lp->iter, varin, (double)lp->upbo[varin]);
      else
	report(lp, NORMAL,
		"Iteration: %d, variable %d changed its upper bound of %g to 0",
		lp->iter, varin, (double)lp->upbo[varin]);
    }
    else
      report(lp, NORMAL,
	      "Iteration: %d, variable %d entered basis at: %g",
	      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]];
      report(lp, NORMAL, "Feasibility gap of this basis: %g",
	      (double)f);
    }
    else
      report(lp, NORMAL,
	      "Objective function value of this feasible basis: %g",
	      (double)lp->rhs[0]);
  }
} /* iteration */

static int primloop(lprec *lp, double *refacttime)
{
  int	 i, ok = TRUE;
  RREAL  theta, f;
  REAL   *drow = NULL, *prow = NULL, *pcol = NULL;
  MYBOOL   primal;
  MYBOOL   minit;
  int    colnr, row_nr;

  if(lp->trace)
    report(lp, DETAILED, "Entering primal algorithm");

  if ((CALLOC(drow, lp->sum + 1) == NULL) ||
      (CALLOC(prow, lp->sum + 1) == NULL) ||
      (CALLOC(pcol, lp->rows + 1) == NULL)
     ) {
    lp->spx_status = OUT_OF_MEMORY;
    ok = FALSE;
  }
  else {
    primal = TRUE;
    lp->spx_status = RUNNING;
    lp->doInvert = FALSE;
    lp->doIterate = FALSE;
    lp->Extrad = 0;

    row_nr = 0;
    colnr = 0;

    minit = FALSE;

    while(lp->spx_status == RUNNING) {
      lp->doIterate = FALSE;
      lp->doInvert = FALSE;

      if ((colnr = colprim(lp, minit, drow)) == -1) {  /* Solve BTRAN here */
        ok = FALSE;
	break;
      }
      if(colnr > 0) {
	if (!setpivcol(lp, colnr, pcol)) {        /* Solve FTRAN here */
	  ok = FALSE;
	  break;
	}

	row_nr = rowprim(lp, colnr, &theta, pcol);
	if(row_nr > 0)
	  if (!condensecol(lp, row_nr, pcol)) {
	    ok = FALSE;
	    break;
	  }
      }

      if(lp->doIterate) {
	iteration(lp, row_nr, colnr, &theta, lp->upbo[colnr], &minit,
		  &lp->lower[colnr], primal);
	if((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT))
	  break;
      }

      if(lp->num_inv > 0)
	f = (timenow()-lp->time_refactstart) / lp->num_inv;
      else
	f = 0;
      if((lp->num_inv >= lp->max_num_inv) ||
	 ((lp->num_inv > 1) && (f >= MINTIMEPIVOT) && (f >= (*refacttime))))
	lp->doInvert = TRUE;

      if(lp->doInvert) {
	if(lp->print_at_invert)
	  report(lp, DETAILED, "Inverting: Primal = %d", primal);
	i = invert(lp);
	if((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT) || (lp->spx_status == OUT_OF_MEMORY)) {
	  ok = FALSE;
	  break;
	}
	else if(!i) {
	  lp->spx_status = SINGULAR_BASIS;
	  break;
	}
	/* Check whether we are still feasible or not... */
	for(i = 1; i <= lp->rows; i++) {
	  f = lp->rhs[i];
	  if((f < -lp->epsb) || (f > lp->upbo[lp->bas[i]] + lp->epsb)) {
	    lp->spx_status = LOST_PRIMAL_FEASIBILITY;
	    break;
	  }
	}
      }
      else
	(*refacttime) = (REAL) f;

      if(yieldformessages(lp)!=0)
	lp->spx_status = USERABORT;

    }
  }

  FREE(drow);
  FREE(prow);
  FREE(pcol);
  return(ok);
} /* primloop */

static int dualloop(lprec *lp, double *refacttime)
{
  int    i, j, ok = TRUE;
  RREAL  f, theta;
  MYBOOL   primal;
  REAL   *drow = NULL, *prow = NULL, *pcol = NULL;
  MYBOOL   minit;
  int    colnr, row_nr;

  if(lp->trace)
    report(lp, DETAILED, "Entering dual algorithm");

  if ((CALLOC(drow, lp->sum + 1) == NULL) ||
      (CALLOC(prow, lp->sum + 1) == NULL) ||
      (CALLOC(pcol, lp->rows + 1) == NULL)
     ) {
    lp->spx_status = OUT_OF_MEMORY;
    FREE(pcol);
    ok = FALSE;
  }
  else {
    primal = FALSE;
    lp->spx_status = RUNNING;
    lp->doInvert = FALSE;
    lp->doIterate = 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!	          */
    lp->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;
	else       /* Since the A-matrix is sorted with the objective function */
	  break;   /* first we do not need to scan the entire matrix  - *** KE added */

      if(f < lp->Extrad)
	lp->Extrad = (REAL) f;
    }

    if(lp->trace)
      report(lp, DETAILED, "Extrad = %g", (double)lp->Extrad);

    row_nr = 0;
    colnr = 0;