solve.c

利用c语言编写

#include <string.h>
#include <sys/types.h>
#include <sys/timeb.h>
#include "lpkit.h"
#include "lpglob.h"
#include "debug.h"

static short milpsolve(lprec *lp, REAL *upbo, REAL *lowbo, MYBOOL *sbasis, MYBOOL *slower, int *sbas, int recursive);

static double timenow()
{
#ifdef INTEGERTIME
  return((double)time(NULL));
#else
  struct timeb buf;

  ftime(&buf);
  return((double)buf.time+((double) buf.millitm)/1000.0);
#endif
}

static int yieldformessages(lprec *lp)
{
  double currenttime = timenow();
  if(lp->sectimeout>0 && ((currenttime-lp->timestart)-(REAL)lp->sectimeout>0))
    lp->spx_status = TIMEOUT;

  if (lp->abort != NULL)
    return(lp->abort(lp, lp->aborthandle));
  else
    return(0);
}

#if USED
/* Get data column stored in a particular eta column */
static int eta_column(lprec *lp, int column)
{
  if(column > lp->eta_size) {
    column = lp->eta_col_end[column] - 1;
    return(lp->eta_row_nr[column]);
  }
  else
    return(0);
}
#endif

static int ftran(lprec *lp, REAL *pcol, REAL roundzero)
{
  int  i, j, k, r, *rowp, ok = TRUE;
  LREAL theta, *vcol;
  REAL  *valuep;

  if (MALLOC(vcol, lp->rows + 1) == NULL) {
    lp->spx_status = OUT_OF_MEMORY;
    ok = FALSE;
  }
  else {
    for(i = 0; i <= lp->rows; i++)
      vcol[i] = pcol[i];

    for(i = 1; i <= lp->eta_size; i++) {
      k = lp->eta_col_end[i] - 1;
      r = lp->eta_row_nr[k];
      theta = vcol[r];
      if(theta != 0) {
	j = lp->eta_col_end[i - 1];

	/* CPU intensive loop, let's do pointer arithmetic */
	for(rowp = lp->eta_row_nr + j, valuep = lp->eta_value + j;
	    j < k;
	    j++, rowp++, valuep++)
	  vcol[*rowp] += theta * *valuep;

	vcol[r] *= lp->eta_value[k];
      }
    }

    for(i = 0; i <= lp->rows; i++)
      pcol[i] = (REAL) vcol[i];

    free(vcol);

    /* round small values to zero */
    if(roundzero != 0)
      for(i = 0; i <= lp->rows; i++)
	my_round(pcol[i], roundzero);
  }
  return(ok);
} /* ftran */


void btran(lprec *lp, REAL *row, REAL roundzero)
{
  int   i, j, k, *rowp;
  LREAL f;
  REAL  *valuep;

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

    for(rowp = lp->eta_row_nr + j, valuep = lp->eta_value + j;
	j <= k;
	j++, rowp++, valuep++)
      f += row[*rowp] * *valuep;

    if(roundzero != 0)
      my_round(f, roundzero);
    row[lp->eta_row_nr[k]] = (REAL) f;
  }
} /* btran */


MYBOOL isvalid(lprec *lp)
{
  int i, j, *rownum = NULL, *colnum;
  int *num = NULL, row_nr;

  /* Check consistency of bounds */
  for(i = 1; i <= lp->columns; i++)
    if(lp->orig_upbo[lp->rows + i] < lp->orig_lowbo[lp->rows + i]){
      report(lp, IMPORTANT, "Error: Column %d has inconsistent bounds (lowbo > upbo).", i);
      return(FALSE);
    }

  /* Check consistency of ranges */
  for(i = 1; i <= lp->rows; i++) 
    if(lp->orig_upbo[i] < lp->orig_lowbo[i]){
      report(lp, IMPORTANT, "Error: Row %d has inconsistent range (lowbo > upbo).", i);
      return(FALSE);
    }

  if(!lp->row_end_valid) {
    if ((CALLOC(num, lp->rows + 1) == NULL) ||
        (CALLOC(rownum, lp->rows + 1) == NULL)
       ) {
      FREE(num);
      FREE(rownum);
      return(FALSE);
    }
    for(i = 0; i < lp->non_zeros; i++)
      rownum[lp->mat[i].row_nr]++;

    lp->row_end[0] = 0;

    for(i = 1; i <= lp->rows; i++)
      lp->row_end[i] = lp->row_end[i - 1] + rownum[i];

    for(i = 1; i <= lp->columns; i++)
      for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++) {
	row_nr = lp->mat[j].row_nr;
	if(row_nr != 0) {
	  num[row_nr]++;
	  lp->col_no[lp->row_end[row_nr - 1] + num[row_nr]] = i;
	}
      }

    FREE(num);
    FREE(rownum);
    lp->row_end_valid = TRUE;
  }

  if(lp->valid)
    return(TRUE);

  rownum = colnum = NULL;
  if ((CALLOC(rownum, lp->rows + 1) == NULL) ||
      (CALLOC(colnum, lp->columns + 1) == NULL)
      ) {
    FREE(rownum);
    FREE(colnum);
    return(FALSE);
  }

  for(i = 1 ; i <= lp->columns; i++)
    for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++) {
      colnum[i]++;
      rownum[lp->mat[j].row_nr]++;
    }

  for(i = 1; i <= lp->columns; i++)
    if(colnum[i] == 0) {
      report(lp, NORMAL, "Warning: Variable %s not used in any constraints",
                get_col_name(lp, i));
    }
  free(rownum);
  free(colnum);
  lp->valid = TRUE;
  return(TRUE);
}

static int resize_eta(lprec *lp, int min_size)
{
  while(lp->eta_alloc <= min_size)
    lp->eta_alloc = (int)((double) lp->eta_alloc*RESIZEFACTOR);
  /* report(lp, FULL, "resizing eta to size %d", lp->eta_alloc); */
  if ((REALLOC(lp->eta_value, lp->eta_alloc + 1) == NULL) ||
      (REALLOC(lp->eta_row_nr, lp->eta_alloc + 1) == NULL)) {
    lp->spx_status = OUT_OF_MEMORY;
    return(FALSE);
  }
  else
    return(TRUE);
} /* resize_eta */

static int condensecol(lprec *lp,
			int row_nr,
			REAL *pcol)
{
  int i, elnr, min_size;

  elnr = lp->eta_col_end[lp->eta_size];

  min_size = elnr + lp->rows + 2;
  if(min_size >= lp->eta_alloc) /* maximum local growth of Eta */
    if (!resize_eta(lp, min_size))
      return(FALSE);

  for(i = 0; i <= lp->rows; i++)
    if(i != row_nr && pcol[i] != 0) {
      lp->eta_row_nr[elnr] = i;
      lp->eta_value[elnr] = pcol[i];
      elnr++;
    }

  lp->eta_row_nr[elnr] = row_nr;
  lp->eta_value[elnr] = pcol[row_nr];
  elnr++;
  lp->eta_col_end[lp->eta_size + 1] = elnr;
  return(TRUE);
} /* condensecol */


static void addetacol(lprec *lp, int colnr)
{
  int  i, j, k;
  LREAL theta;

  j = lp->eta_col_end[lp->eta_size];
  lp->eta_size++;

  k = lp->eta_col_end[lp->eta_size] - 1;
  theta = 1 / lp->eta_value[k];
  lp->eta_value[k] = (REAL) theta;
  for(i = j; i < k; i++)
    lp->eta_value[i] = (REAL) (-theta * lp->eta_value[i]);
  lp->justInverted = FALSE;
} /* addetacol */


static int setpivcol(lprec *lp, int varin, REAL *pcol)
{
  int  i, j, colnr, ok = TRUE;

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

  if(lp->lower[varin]) {
    if(varin > lp->rows) {
      colnr = varin - lp->rows;
      for(i = lp->col_end[colnr - 1]; i < lp->col_end[colnr]; i++)
	pcol[lp->mat[i].row_nr] = lp->mat[i].value;
      pcol[0] -= lp->Extrad;
    }
    else
      pcol[varin] = 1;
  }
  else { /* !lower */
    if(varin > lp->rows) {
      colnr = varin - lp->rows;
      for(i = lp->col_end[colnr - 1]; i < lp->col_end[colnr]; i++)
	pcol[lp->mat[i].row_nr] = -lp->mat[i].value;
      pcol[0] += lp->Extrad;
    }
    else
      pcol[varin] = -1;
  }

 /* Test if we should do the error-correction version */
  if((lp->improve & IMPROVE_FTRAN) && lp->num_inv) {
    REAL   *errors = NULL, sdp;
    matrec *matentry;
    int    k, ie;

    if ((MALLOCCPY(errors, pcol, lp->rows + 1) == NULL) || (!ftran(lp, pcol, lp->epsel))) {
      FREE(errors);
      lp->spx_status = OUT_OF_MEMORY;
      ok = FALSE;
    }
    else {
      for(j = 1; j <=lp->rows; j++) {
	colnr = lp->bas[j];
	sdp = pcol[j];
	if(colnr <= lp->rows)        /* A slack variable is in the basis */
	  errors[colnr] -= sdp;
	else {                       /* A normal variable is in the basis */
	  colnr -= lp->rows;
	  ie = lp->col_end[colnr];
	  i = lp->col_end[colnr - 1];
	  for(matentry = lp->mat + i; i < ie; i++, matentry++) {
	    k = (*matentry).row_nr;
	    errors[k] -= (*matentry).value*sdp;
	  }
	}
      }
      if (!ftran(lp, errors, lp->epsel))
        ok = FALSE;
      else {
	sdp = 0;
	for(j = 1; j <=lp->rows; j++)
	  if(fabs(errors[j])>sdp)
	    sdp = fabs(errors[j]);
    /*      sdp += pow(errors[j],2); */
    /*    sdp = sqrt(sdp/lp->rows); */
	if(sdp > lp->epsel) {
	  if(lp->debug)
	    report(lp, DETAILED, "Iterative FTRAN correction metric %g", sdp);
	  for(j = 1; j <=lp->rows; j++)
	    pcol[j] += errors[j];
	}
      }
      free(errors);
    }
  }
  else
    if (!ftran(lp, pcol, lp->epsel))
      ok = FALSE;

  return(ok);
} /* setpivcol */

static int minoriteration(lprec *lp, int colnr, int row_nr)
{
  int  i, j, k, wk, varin, varout, elnr;
  LREAL piv = 0, theta;

  varin = colnr + lp->rows;
  elnr = lp->eta_col_end[lp->eta_size];
  wk = elnr;
  lp->eta_size++;

  if(lp->Extrad != 0) {
    lp->eta_row_nr[elnr] = 0;
    lp->eta_value[elnr] = -lp->Extrad;
    elnr++;
    if(elnr >= lp->eta_alloc)
      if (!resize_eta(lp, elnr))
        return(FALSE);
  }

  /* Move pivot column data to Eta (but not the pivot row item yet) */
  for(j = lp->col_end[colnr - 1] ; j < lp->col_end[colnr]; j++) {
    k = lp->mat[j].row_nr;

    if(k == 0 && lp->Extrad != 0)
      lp->eta_value[lp->eta_col_end[lp->eta_size - 1]] += lp->mat[j].value;
    else if(k != row_nr) {
      lp->eta_row_nr[elnr] = k;
      lp->eta_value[elnr] = lp->mat[j].value;
      elnr++;
      if(elnr >= lp->eta_alloc)
	if (!resize_eta(lp, elnr))
	  return(FALSE);
    }
    else
      piv = lp->mat[j].value;
  }

  /* Now add the pivot row item to Eta */
  lp->eta_row_nr[elnr] = row_nr;
  lp->eta_value[elnr] = (REAL) (1 / piv);
  theta = lp->rhs[row_nr] / piv;
  lp->rhs[row_nr] = (REAL) theta;

  /* Update RHS for new pivot column */
  for(i = wk; i < elnr; i++)
    lp->rhs[lp->eta_row_nr[i]] = (REAL) (lp->rhs[lp->eta_row_nr[i]] - theta * lp->eta_value[i]);

  varout = lp->bas[row_nr];
  lp->bas[row_nr] = varin;
  lp->basis[varout] = FALSE;
  lp->basis[varin] = TRUE;

  /* Scale the pivoted column in Eta by the pivot value */
  for(i = wk; i < elnr; i++)
    lp->eta_value[i] = (REAL) (-lp->eta_value[i] / piv);

  lp->eta_col_end[lp->eta_size] = elnr + 1;
  return(TRUE);
} /* minoriteration */


static void rhsmincol(lprec *lp,
		      LREAL theta,
		      int row_nr,
		      int varin)
{
  int   i, j, k, varout;
  LREAL f;

  if(row_nr > lp->rows + 1) {
    if (lp->trace) {
      report(lp, IMPORTANT, "Error: rhsmincol called with row_nr: %d, rows: %d",
             row_nr, lp->rows);
      report(lp, IMPORTANT, "This indicates numerical instability");
    }
    lp->spx_status = FAILURE;
    return;
  }

  j = lp->eta_col_end[lp->eta_size];
  k = lp->eta_col_end[lp->eta_size + 1];
  for(i = j; 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]] = (REAL) f;
  }

  lp->rhs[row_nr] = theta;
  varout = lp->bas[row_nr];
  lp->bas[row_nr] = varin;
  lp->basis[varout] = FALSE;
  lp->basis[varin] = TRUE;
} /* rhsmincol */

#if USED
static void get_markowitz(lprec *lp, int *rownz, int *colnz, MYBOOL *frow, MYBOOL *fcol, int *rownr, int *colnr)
{
  int elmnr;
  int i, j, holdval, minval;
  unsigned int _MAXINT;

  _MAXINT =(unsigned int)  ~0;
  _MAXINT = _MAXINT / 2 - 1;
  minval = _MAXINT;
  (*rownr) = 0;
  (*colnr) = 0;

  for(j = 1; j <= lp->columns; j++) {
    if(!fcol[j] || !colnz[j]) continue;

   /* Make sure that we have a non-zero value to pivot on */
    for(elmnr = lp->col_end[j-1]; elmnr < lp->col_end[j]; elmnr++) {
      i = lp->mat[elmnr].row_nr;
      if(!frow[i] || !rownz[i]) continue;
      if(fabs(lp->mat[elmnr].value) < lp->epspivot) continue;

     /* Now compute the statistic */
      holdval = (rownz[i]-1)*(colnz[j]-1);
      if(holdval > 0 && holdval < minval) {
        minval = holdval;
        (*rownr) = i;
        (*colnr) = j;
      }
    }
  }
}
#endif

static MYBOOL invert(lprec *lp)
{
  LREAL  theta, hold;
  matrec *matentry;
  REAL   *pcol = NULL;
  MYBOOL   *fcol = NULL, *frow = NULL;
  int    *colnum = NULL, *rownum = NULL, *col = NULL, *row = NULL;
  int    k, kk, kkk, i, j, v, numit, varnr, rownr, colnr, varin;
  int	 singularities;
  short  spx_save;
  MYBOOL   Restart;

 /* Must save spx_status since it is used to carry information about
    the presence and handling of singular columns in the matrix */
  spx_save = lp->spx_status;
  lp->spx_status = RUNNING;
  if(yieldformessages(lp)!=0)
    lp->spx_status = USERABORT;

  if((lp->usermessage != NULL) && (lp->msgmask & MSG_INVERT))
    lp->usermessage(lp, lp->msghandle, MSG_INVERT);

  if(lp->spx_status != RUNNING)
    return(FALSE);
  lp->spx_status = spx_save;
  singularities = 0;

  if(lp->print_at_invert)
    report(lp, DETAILED, "Start Invert iter %d eta_size %d rhs[0] %g ",
	    lp->iter, lp->eta_size, (double) - lp->rhs[0]);

  if ((CALLOC(col, lp->rows + 1) == NULL) ||
      (CALLOC(row, lp->rows + 1) == NULL) ||
      (CALLOC(pcol, lp->rows + 1) == NULL) ||
      (CALLOC(frow, lp->rows + 1) == NULL) ||
      (CALLOC(fcol, lp->columns + 1) == NULL) ||
      (CALLOC(rownum, lp->rows + 1) == NULL) ||
      (CALLOC(colnum, lp->columns + 1) == NULL)
     ) {
    lp->spx_status = OUT_OF_MEMORY;