solve.c

利用c语言编写

    minit = FALSE;

    while(lp->spx_status == RUNNING) {

      lp->doIterate = FALSE;
      lp->doInvert = FALSE;

      if(!minit)
	row_nr = rowdual(lp);

      if(row_nr > 0 ) {
	colnr = coldual(lp, row_nr, minit, prow, drow); /* What about BTRAN in here? */
	/* report(lp, NORMAL, "Dual-phase pivots (minit=%d) col=%d, row=%d", minit, colnr, row_nr); */

	if(colnr > 0) {
	  if (!setpivcol(lp, colnr, pcol)) {                   /* Solve FTRAN here */
	    ok = FALSE;
	    break;
	  }

	 /* getting div by zero here. Catch it and try to recover */
	  if(pcol[row_nr] == 0) {
	    report(lp, NORMAL, "An attempt was made to divide by zero (pcol[%d])", row_nr);
	    lp->doIterate = FALSE;
	    if(!lp->justInverted) {
	      report(lp, NORMAL, "Trying to recover. Reinverting Eta");
	      lp->doInvert = TRUE;
	    }
	    else {
	      report(lp, NORMAL, "Failed to recover. Can't reinvert");
	      lp->spx_status = FAILURE;
	    }
	  }
	  else {
	    int bnx;

	    if (!condensecol(lp, row_nr, pcol)) {
	      ok = FALSE;
	      break;
	    }
	    bnx = lp->bas[row_nr];
	    f = lp->rhs[row_nr] - lp->upbo[bnx];

	    if(f > 0) {
	      theta = f / (RREAL) pcol[row_nr];
	      if(theta <= lp->upbo[colnr] + lp->epsb)
		lp->lower[bnx] = (MYBOOL) !lp->lower[bnx];
	    }
	    else /* f <= 0 */
	      theta = lp->rhs[row_nr] / (RREAL) pcol[row_nr];
	  }
	}
	else
	  lp->spx_status = INFEASIBLE;
      }
      else {
	lp->spx_status = SWITCH_TO_PRIMAL;
	lp->doIterate = FALSE;
	lp->Extrad = 0;
	lp->doInvert = TRUE;
      }

      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) / (REAL)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;
	}
      }
      else
	(*refacttime) = (REAL) f;

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

    }
  }

  FREE(drow);
  FREE(prow);
  FREE(pcol);

  return(ok);
}

static short solvelp(lprec *lp)
{
  int    i, singular_count, lost_feas_count;
  MYBOOL   feasible;
  REAL	 x, refacttime;
  int    iter;

  lp->iter = 0;
  iter = 0;      /* Set to -1 to use once-through loop */

  refacttime = 0;
  singular_count = 0;
  lost_feas_count = 0;
  lp->spx_status = RUNNING;

  while(lp->spx_status == RUNNING) {

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

    /* Check whether we are feasible or infeasible. */
    if(iter >= 0)
      iter = lp->iter;
    feasible = TRUE;
    for(i = 1; i <= lp->rows; i++) {
      x = (REAL) lp->rhs[i];
      if((x < 0) || (x > lp->upbo[lp->bas[i]])) {
        feasible = FALSE;
        break;
      }
    }

   /* Now do the simplex magic */
    if(feasible) {
      if(lp->trace)
        report(lp, NORMAL, "Start at feasible basis");
      if (!primloop(lp, &refacttime))
        break;
    }
    else {
      if(lp->trace) {
        if(lost_feas_count > 0)
          report(lp, NORMAL, "Continuing at infeasible basis");
        else
          report(lp, NORMAL, "Start at infeasible basis");
      }
      if (!dualloop(lp, &refacttime))
        break;
      if(lp->spx_status == SWITCH_TO_PRIMAL)
        if (!primloop(lp, &refacttime))
	  break;
    }

    if(lp->spx_status == SINGULAR_BASIS) {
      singular_count++;
      if(singular_count >= DEF_MAXSINGULARITIES) {
        report(lp, NORMAL, "SINGULAR BASIS!  Too many singularities - aborting.");
        lp->spx_status = FAILURE;
        break;
      }
      report(lp, DETAILED, "SINGULAR BASIS!  Will attempt to recover.");
      lp->spx_status = RUNNING;
      /* Singular pivots are simply skipped by the inversion, leaving a row
         a row's slack var in the basis instead of the singular problem var
         This basis could be feasible or infeasible.  Check how to restart. */
    }
    else if((lp->spx_status == LOST_PRIMAL_FEASIBILITY) || ((iter > 0) && (iter == lp->iter))) {
      lost_feas_count++;
      if(lost_feas_count >= DEF_MAXSINGULARITIES) {
        report(lp, NORMAL, "LOST PRIMAL FEASIBILITY too many times, aborting.");
        lp->spx_status = FAILURE;
        break;
      }
      report(lp, DETAILED, "LOST PRIMAL FEASIBILITY!  Recovering.");
      lp->spx_status = RUNNING;
    }

  }

  lp->total_iter += lp->iter;

  return(lp->spx_status);
} /* solvelp */

static MYBOOL solution_is_int(lprec *lp, int i)
{
  REAL value;

  value = lp->solution[i];
  value = value - (REAL)floor((double)value);

  if(value < lp->epsilon) {
/*    lp->solution[i] = (REAL)floor((double)lp->solution[i]); */
    return(TRUE);
  }

  if(value > (1 - lp->epsilon)) {
/*    lp->solution[i] = (REAL)floor((double)lp->solution[i]+1); */
    return(TRUE);
  }

  return(FALSE);
} /* solution_is_int */


static void construct_solution(lprec *lp)
{
  int    i, j, basi;
  REAL   f;

  /* zero all results of rows */
  for (i = 1; i <= lp->rows; i++)
    lp->solution[i] = 0.0;
  lp->solution[0] = -lp->orig_rh[0];

  if(lp->scaling_used) {
    lp->solution[0] /= lp->scale[0];

    for(i = lp->rows + 1; i <= lp->sum; i++)
      lp->solution[i] = lp->lowbo[i] * lp->scale[i];

    for(i = 1; i <= lp->rows; i++) {
      basi = lp->bas[i];
      if(basi > lp->rows)
	lp->solution[basi] += (REAL) (lp->rhs[i] * lp->scale[basi]);
    }
    for(i = lp->rows + 1; i <= lp->sum; i++)
      if(!lp->basis[i] && !lp->lower[i])
	lp->solution[i] += lp->upbo[i] * lp->scale[i];

    for(j = 1; j <= lp->columns; j++) {
      f = lp->solution[lp->rows + j];
      if(f != 0)
	for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++) {
          basi = lp->mat[i].row_nr;
          lp->solution[basi] += (f / lp->scale[lp->rows+j])
	    * (lp->mat[i].value / lp->scale[basi]);
        }
    }
  }
  else { /* no scaling */
    for(i = lp->rows + 1; i <= lp->sum; i++)
      lp->solution[i] = lp->lowbo[i];

    for(i = 1; i <= lp->rows; i++) {
      basi = lp->bas[i];
      if(basi > lp->rows)
	lp->solution[basi] += (REAL) lp->rhs[i];
    }

    for(i = lp->rows + 1; i <= lp->sum; i++)
      if(!lp->basis[i] && !lp->lower[i])
	lp->solution[i] += lp->upbo[i];

    for(j = 1; j <= lp->columns; j++) {
      f = lp->solution[lp->rows + j];
      if(f != 0)
	for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
	  lp->solution[lp->mat[i].row_nr] += f * lp->mat[i].value;
    }

  }

  /* clean out near-zero slack values */
  for(i = 0; i <= lp->rows; i++) {
    if(my_abs(lp->solution[i]) < lp->epsb)
      lp->solution[i] = 0;
    else if(lp->ch_sign[i])
      lp->solution[i] = -lp->solution[i];
  }

} /* construct_solution */

static void transfer_solution(lprec *lp)
{
  int i;

  MEMCPY(lp->best_solution, lp->solution, lp->sum + 1);

  /* round integer solution values to actual integers */
  if(lp->scalemode & INTEGERSCALE)
    for(i = 1; i <= lp->columns; i++)
      if(is_int(lp, i))
        lp->best_solution[lp->rows + i] = floor(lp->best_solution[lp->rows + i] + 0.5);

}

static void calculate_duals(lprec *lp)
{
  int varnr, i, j;
  REAL scale0;
  LREAL f;

  if(!lp->eta_valid)
    return;

  /* initialize */
  lp->duals[0] = 1;
  for(i = 1; i <= lp->sum; i++)
    lp->duals[i] = 0;

  btran(lp, lp->duals, lp->epsel);

  for(i = 1; i <= lp->columns; i++) {
    varnr = lp->rows + i;
    if(!lp->basis[varnr])
      if((lp->upbo[varnr] > 0) && (lp->solution[varnr] > 0.0)) {
	f = 0;
	for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
	  f += (LREAL)lp->duals[lp->mat[j].row_nr] * (LREAL)lp->mat[j].value;
	lp->duals[varnr] = (REAL)f;
      }
  }

  /* the dual values are the reduced costs of the slacks */
  /* When the slack is at its upper bound, change the sign. */
  for(i = 1; i <= lp->rows; i++) {
    if(lp->basis[i])
      lp->duals[i] = 0;
    /* added a test if variable is different from 0 because sometime you get
       -0 and this is different from 0 on for example INTEL processors (ie 0
       != -0 on INTEL !) peno */
    else if((lp->ch_sign[0] == lp->ch_sign[i]) && lp->duals[i])
      lp->duals[i] = - lp->duals[i];
  }

  if (lp->scaling_used)
    scale0 = lp->scale[0];
  else
    scale0 = 1;
  for(i = 1; i <= lp->sum; i++) {
    if(lp->scaling_used) {
      lp->duals[i] /= scale0;
      if(i <= lp->rows)
        lp->duals[i] *= lp->scale[i];
      else
        lp->duals[i] /= lp->scale[i];
    }
    my_round(lp->duals[i], lp->epsd);
  }
} /* calculate_duals */

/* calculate sensitivity duals */
static int calculate_sensitivity_duals(lprec *lp)
 {
  int k,varnr, ok = TRUE;
  REAL *pcol,a,infinite,epsel,from,till;

  /* one column of the matrix */
  if (CALLOC(pcol, lp->rows + 1) == NULL) {
    lp->spx_status = OUT_OF_MEMORY;
    ok = FALSE;
  }
  else {
    infinite=lp->infinite;
    epsel=lp->epsel;
    for(varnr=1;varnr<=lp->sum;varnr++) {
     from=infinite;
     till=infinite;
     if ((!lp->basis[varnr]) && ((varnr<=lp->rows) || (lp->solution[varnr]>0.0))) {
      if (!setpivcol(lp,varnr,pcol)) { /* construct one column of the tableau */
       ok = FALSE;
       break;
      }
      for (k=1;k<=lp->rows;k++)   /* search for the rows(s) which first results in further iterations */
       if (my_abs(pcol[k])>epsel) {
	a=(REAL) (lp->rhs[k]/pcol[k]);
	if (lp->scaling_used) {
	 if (varnr<=lp->rows)
	  a/=lp->scale[varnr];
	 else
	  a*=lp->scale[varnr];
        }
	if ((a<=0.0) && (pcol[k]<0.0) && (-a<from)) from=-a;
	if ((a>=0.0) && (pcol[k]>0.0) && ( a<till)) till= a;
	if (lp->upbo[lp->bas[k]] < infinite) {
	 a=(REAL) ((lp->rhs[k]-lp->upbo[lp->bas[k]])/pcol[k]);
	 if (lp->scaling_used) {
	  if (varnr<=lp->rows)
	   a/=lp->scale[varnr];
	  else
	   a*=lp->scale[varnr];
         }
	 if ((a<=0.0) && (pcol[k]>0.0) && (-a<from)) from=-a;
	 if ((a>=0.0) && (pcol[k]<0.0) && ( a<till)) till= a;
	}
       }
      if (!lp->lower[varnr]) {
       a=from;
       from=till;
       till=a;
      }
      if ((varnr<=lp->rows) && (!lp->ch_sign[varnr])) {
       a=from;
       from=till;
       till=a;
      }
     }
     if (from!=infinite)
      lp->dualsfrom[varnr]=lp->solution[varnr]-from;
     else
      lp->dualsfrom[varnr]=-infinite;
     if (till!=infinite)
      lp->dualstill[varnr]=lp->solution[varnr]+till;
     else
      lp->dualstill[varnr]=infinite;
    }

    free(pcol);
  }
  return(ok);
 } /* calculate_sensitivity_duals */

static void setpivrow(lprec *lp,
		      int   row_nr,
		      REAL *prow)
 {
  int i,j,k,r,*rowp,row,varnr;
  REAL f,*valuep,value;

  for(i = 0; i <= lp->rows; i++)
    prow[i] = 0;

  prow[row_nr] = 1;

  for(i = lp->eta_size; i >= 1; i--) {
    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;
    }

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

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

      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;
	f += prow[row] * value;
      }

      my_round(f, lp->epsel);
      prow[varnr] = f;
    }
  }
} /* setpivrow */

/* calculate sensitivity objective function */
static int calculate_sensitivity_obj(lprec *lp)
 {
  int i,j,l,varnr,row_nr,ok = TRUE;
  REAL *OrigObj = NULL,*drow = NULL,*prow = NULL,f,a,min1,min2,infinite,epsel,from,till;

  /* objective function */
  if ((CALLOC(drow, lp->sum + 1) == NULL) ||
      (MALLOC(OrigObj, lp->columns + 1) == NULL) ||
      (CALLOC(prow,lp->sum + 1) == NULL)
     ) {
    lp->spx_status = OUT_OF_MEMORY;
    FREE(prow);
    FREE(OrigObj);
    FREE(drow);
    ok = FALSE;
  }
  else {
    for(i = 1; i <= lp->sum; i++)
      drow[i] = 0;