solve.c

利用c语言编写

    FREE(colnum);
    FREE(rownum);
    FREE(fcol);
    FREE(frow);
    FREE(pcol);
    FREE(row);
    FREE(col);
    return(FALSE);
  }

  lp->time_refactstart = timenow();    /* Time of start of current cyle */

  /* Initialize working basis indicators to all slacks ... */
  for(i = 0; i <= lp->rows; i++)
    frow[i] = TRUE;                   /* Row slack is in the basis */
  for(i = 0; i < lp->columns; i++)
    fcol[i] = FALSE;                  /* Column has not been pivoted in */

  /* ... then store state of pre-existing basis */
  for(i = 0; i <= lp->rows; i++)
    if(lp->bas[i] > lp->rows)
      fcol[lp->bas[i] - lp->rows] = TRUE;
    else
      frow[lp->bas[i]] = FALSE;

  /* Get row and column number entry counts for basic slacks (includes OF row) */
  for(i = 1; i <= lp->rows; i++) {
    if(frow[i])
      for(j = lp->row_end[i - 1] + 1; j <= lp->row_end[i]; j++) {
	v = lp->col_no[j];
	if(fcol[v]) {
	  colnum[v]++;
	  rownum[i]++;
	}
      }
  }

  /* Reset basis indicators to all slacks */
  for(i = 1; i <= lp->rows; i++) {
    lp->bas[i] = i;
    lp->basis[i] = TRUE;
  }
  for(i = 1; i <= lp->columns; i++)
    lp->basis[i + lp->rows] = FALSE;

  /* Save lower bound-adjusted RHS */
  for(i = 0; i <= lp->rows; i++)
    lp->rhs[i] = lp->rh[i];

  /* Adjust active RHS for state of variable */
  for(i = 1; i <= lp->columns; i++) {
    varnr = lp->rows + i;
    if(!lp->lower[varnr]) {
      theta = lp->upbo[varnr];
      k = lp->col_end[i];
      j = lp->col_end[i - 1];
      for(matentry = lp->mat + j; j < k; j++, matentry++) {
        v = (*matentry).row_nr;
        lp->rhs[v] -= theta * (*matentry).value;
      }
    }
  }

  /* Finally, adjust for row state if it is at its upper bound */
  for(i = 1; i <= lp->rows; i++)
    if(!lp->lower[i])
      lp->rhs[i] -= lp->upbo[i];

  /* Check timeout and user abort again */
  spx_save = lp->spx_status;
  lp->spx_status = RUNNING;
  if(yieldformessages(lp) != 0)
    lp->spx_status = USERABORT;
  else
    lp->spx_status = spx_save;

  k = 0;    /* Total number of rows pivoted */
  numit = 0;
  singularities = 0;
  kkk = 0;  /* Number of singleton rows pivoted in current iteration */
  if(!((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT) || (lp->spx_status == OUT_OF_MEMORY))) {
    /* Progress to the inversion process proper */
    lp->num_inv = 0;
    lp->num_refact++;
    lp->eta_size = 0;

    /* Loop reentry point */

    kk = 0;   /* Singleton pivot iteration counter */
    do {
      kk++;
      kkk = 0;  /* Number of singleton rows pivoted in current iteration */

      /* Loop over rows, hunting for row singletons */
      rownr = 0;
      v = 0;
      while(v < lp->rows) {
        rownr++;
        if(rownr > lp->rows)
          rownr = 1;

        v++;
        if(rownum[rownr] == 1)
          if(frow[rownr]) {
    	    v = 0;
            k++;
            kkk++;

            /* Find first column available to be pivoted */
            j = lp->row_end[rownr - 1] + 1;
    	    while(!(fcol[lp->col_no[j]]))
    	      j++;
    	    colnr = lp->col_no[j];

            /* Reduce item counts for the selected pivot column/row */
            colnum[colnr] = 0;
    	    fcol[colnr] = FALSE;
    	    for(j = lp->col_end[colnr - 1]; j < lp->col_end[colnr]; j++)
    	      if(frow[lp->mat[j].row_nr])
    		rownum[lp->mat[j].row_nr]--;
    	    frow[rownr] = FALSE;
            /* if(rownum[rownr]) */
            /*   rownum[rownr] = 0; */

            /* Perform the pivot */
            if (!minoriteration(lp, colnr, rownr))
	      break;
          }
      }

      if(!((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT) || (lp->spx_status == OUT_OF_MEMORY))) {
	/* Loop over columns, hunting for column singletons */
	colnr = 0;
	v = 0;
	if((k < lp->rows) && ((kk <= 0) || (kkk > 0)))
	  while(v < lp->columns) {
	    colnr++;
	    if(colnr > lp->columns)
	      colnr = 1;

	    v++;
	    if(colnum[colnr] == 1)
	      if(fcol[colnr]) {
      		v = 0;
		k++;
		kkk++;

		/* Find first available row to be pivoted */
		j = lp->col_end[colnr - 1];
      		while(!(frow[lp->mat[j].row_nr]))
      		  j++;
      		rownr = lp->mat[j].row_nr;

		/* Reduce item counts for the selected pivot column/row */
		rownum[rownr] = 0;
      		frow[rownr] = FALSE;
      		for(j = lp->row_end[rownr - 1] + 1; j <= lp->row_end[rownr]; j++)
      		  if(fcol[lp->col_no[j]])
      		    colnum[lp->col_no[j]]--;
      		fcol[colnr] = FALSE;
		/* if(colnum[colnr]) */
		/*   colnum[colnr] = 0; */

		/* Store pivot information */
		numit++;
      		col[numit - 1] = colnr;
      		row[numit - 1] = rownr;
	      }
	  }

	/* Check timeout and user abort again */
	spx_save = lp->spx_status;
	lp->spx_status = RUNNING;
	if(yieldformessages(lp)!=0)
	  lp->spx_status = USERABORT;
	else
	  lp->spx_status = spx_save;
      }

      Restart = FALSE;
      if(!((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT) || (lp->spx_status == OUT_OF_MEMORY))) {
        /* Check for more singletons, exhaust the supply - Added by KE */
        if((kkk > 0) && (k < lp->rows))
          Restart = TRUE;
        else if (k < lp->rows) {
#if 0
          for(j = 1; j <= lp->columns; j++) {
            get_markowitz(lp, rownum, colnum, frow, fcol, &rownr, &colnr);
            if(colnr) {
              fcol[colnr] = FALSE;
              /*    colnum[colnr] = 0; */
              if (!setpivcol(lp, lp->rows + colnr, pcol)) {
  		Restart = FALSE;
                break;
              }
  	      else {
                frow[rownr] = FALSE;
            /*    rownum[rownr]--; */
                condensecol(lp, rownr, pcol);
                theta = lp->rhs[rownr] / (LREAL)pcol[rownr];
                rhsmincol(lp, theta, rownr, lp->rows + colnr);
                addetacol(lp, colnr);
                k++;
                kkk++;
  		if(k >= lp->rows)
                  break;
  	      }
            }
          }
#endif
        }
      }
    } while (Restart);
  }
  if(!((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT) || (lp->spx_status == OUT_OF_MEMORY))) {
    /* Find pivots for remaining non-singleton cases where a column is in the basis */
    if(k < lp->rows) {
      for(j = 1; j <= lp->columns; j++) {
        colnr = j;
        if(fcol[colnr]) {
          fcol[colnr] = FALSE;
          if (!setpivcol(lp, lp->rows + colnr, pcol))
	     break;

          /* Find first coefficient available row to be pivoted; **** KE deleted!
             (original lp_solve version that brings a lot of numerical instability) */
         /*
          rownr = 1;
          while((rownr <= lp->rows) && (!(frow[rownr] && pcol[rownr])))
    	    rownr++;
         */

          /* Find largest coefficient available row to be pivoted; **** KE added!
             much better numerical stability, but requires more CPU processing */
          rownr = lp->rows + 1;
          hold = 0;
          for(i = 1; i <= lp->rows; i++) {
            if(frow[i] && fabs(pcol[i])>hold) {
              hold = fabs(pcol[i]);
              rownr = i;
            }
          }

          if(rownr > lp->rows) {
            /* This column is singular!  Just skip it, leaving one of the
               slack variables basic in its place... (Source: Geosteiner changes!) */
            report(lp, DETAILED, "--> Column %d is singular! Skipped.", colnr);
            singularities++;
          }
          else {
            frow[rownr] = FALSE;
            if (!condensecol(lp, rownr, pcol))
	      break;
            theta = lp->rhs[rownr] / (LREAL)pcol[rownr];
            rhsmincol(lp, (REAL) theta, rownr, lp->rows + colnr);
            addetacol(lp, colnr);
          }
          k++;
          kkk++;
          if(k >= lp->rows)
            break;
        }
      }

      if(!((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT) || (lp->spx_status == OUT_OF_MEMORY))) {
	/* Check timeout and user abort again */
	spx_save = lp->spx_status;
	lp->spx_status = RUNNING;
	if(yieldformessages(lp)!=0)
	  lp->spx_status = USERABORT;
	else
	  lp->spx_status = spx_save;
      }
    }
  }
  if(!((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT) || (lp->spx_status == OUT_OF_MEMORY))) {
    /* Perform pivoting of the row-column combinations stored above */
    for(i = numit - 1; i >= 0; i--) {
      colnr = col[i];
      rownr = row[i];
      varin = lp->rows + colnr;

      /* Move the constraint column to the dense pcol vector */
      for(j = 0; j <= lp->rows; j++)
        pcol[j] = 0;
      for(j = lp->col_end[colnr - 1]; j < lp->col_end[colnr]; j++)
        pcol[lp->mat[j].row_nr] = lp->mat[j].value;
      pcol[0] -= lp->Extrad;

      if (!condensecol(lp, rownr, pcol))
        break;
      theta = lp->rhs[rownr] / (LREAL)pcol[rownr];
      rhsmincol(lp, (REAL) theta, rownr, varin);
      addetacol(lp, colnr);
    }

    if(!((lp->spx_status == USERABORT) || (lp->spx_status == TIMEOUT) || (lp->spx_status == OUT_OF_MEMORY))) {
      /* Round net RHS values */
      for(i = 1; i <= lp->rows; i++)
	my_round(lp->rhs[i], lp->epsel);

      /* Do user reporting */
      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->print_at_invert)
	report(lp, DETAILED,
  		"End Invert                eta_size %d rhs[0] %g",
  		lp->eta_size, (double) - lp->rhs[0]);

      /* Set inversion completion status */
      lp->justInverted = TRUE;
      lp->doInvert = FALSE;
    }
  }
  free(rownum);
  free(col);
  free(row);
  free(pcol);
  free(frow);
  free(fcol);
  free(colnum);

  return((MYBOOL) (singularities <= 0));
} /* invert */

static int colprim(lprec *lp,
		     MYBOOL minit,
		     REAL   *drow)
{
  int  varnr, i, j, k, ie, ok = TRUE;
  int  colnr;
  LREAL f, dpiv, opiv;
  matrec *matentry;

  if(!minit) {
    for(i = 1; i <= lp->sum; i++)
      drow[i] = 0;
    drow[0] = 1;

    /* Test if we should do the error-correction version */
    if(FALSE && (lp->improve & IMPROVE_BTRAN) && lp->num_inv) {
      REAL  *errors;

      if (MALLOCCPY(errors, drow, lp->rows + 1) == NULL) {
        lp->spx_status = OUT_OF_MEMORY;
        ok = FALSE;
      }
      else {
	btran(lp, drow, lp->epsel);

	for(j = 1; j <=lp->rows; j++) {
	  colnr = lp->bas[j];
	  if(colnr <= lp->rows)        /* A slack variable is in the basis */
            f = drow[j];
	  else {                       /* A normal variable is in the basis */
	    colnr -= lp->rows;
	    f = 0;
	    ie = lp->col_end[colnr];
	    i = lp->col_end[colnr - 1];
	    for(matentry = lp->mat + i; i < ie; i++, matentry++) {
	      k = (*matentry).row_nr;
	      f += drow[k] * (*matentry).value;
	    }
          }
          errors[j] -= (REAL) f;
	}
	btran(lp, errors, lp->epsel);

	f = 0;
	for(j = 1; j <=lp->rows; j++)
  /*        f += pow(errors[j],2); */
	  if(fabs(errors[j])>f)
	    f = fabs(errors[j]);
  /*      f = sqrt(f/lp->rows); */
	if(f > lp->epsel) {
	  if(lp->debug)
	    report(lp, DETAILED, "Iterative BTRAN correction metric %g", f);
	  for(j = 1; j <=lp->rows; j++)
	    drow[j] += errors[j];
	}
	free(errors);
      }
    }
    else
      btran(lp, drow, lp->epsel);

    if (ok) {
      /* Continue */
      for(i = 1; i <= lp->columns; i++) {
	varnr = lp->rows + i;
	if(!lp->basis[varnr])
	  if(lp->upbo[varnr] > 0) {
	    f = 0;
	    ie = lp->col_end[i];
	    j = lp->col_end[i - 1];
	    for(matentry = lp->mat + j; j < ie; j++, matentry++) {
	      k = (*matentry).row_nr;
	      f += drow[k] * (*matentry).value;
	    }
	    drow[varnr] = (REAL) f;
	  }
      }
    }
    if (ok)
      for(i = 1; i <= lp->sum; i++)
        my_round(drow[i], lp->epsd);
  }

  if (ok) {
    /* Identify pivot column (fall-through is 'largest coefficient') */
    dpiv = lp->epsd;
    colnr = 0;
    opiv = 0;
    varnr = 0;

    for(i = lp->sum; i > 0; i--) {
      if(!lp->basis[i])
	if(lp->upbo[i] > 0) {

       /* Retrieve the reduced cost, skip if non-positive */
          if(lp->lower[i])
	    f = -drow[i];
	  else
	    f = drow[i];
          if(f < lp->epsd) continue;

          /* Save largest reduced cost */
	  if(f > dpiv) {
            dpiv = f;
            colnr = i;
          }

          /* Compute objective contribution */
	  if(lp->piv_rule == GREEDY_SELECT && i > lp->rows) {
            f = get_mat_raw(lp, 0, i - lp->rows);
   	    if(f < opiv) {
              opiv = f;
	      varnr = i;
	    }
	  }

 	  if(varnr > 0 && i <= lp->rows+1) {
	    colnr = varnr;
	    break;
	  }
          if(colnr > 0 && lp->piv_rule == FIRST_SELECT) break;
	}
    }

    if(lp->trace) {
      if(colnr>0)
	report(lp, NORMAL, "col_prim:%d, reduced cost: %g",
		colnr, (double)dpiv);
      else
	report(lp, NORMAL,
		"col_prim: no positive reduced costs found, optimality!\n");
    }
    if(colnr == 0) {
      lp->doIterate = FALSE;
      lp->doInvert = FALSE;
      lp->spx_status = OPTIMAL;
    }
  }
  else {
    colnr = -1;
    lp->doIterate = FALSE;
    lp->doInvert = FALSE;
    lp->spx_status = OUT_OF_MEMORY;
  }
  return(colnr);
} /* colprim */

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

  row_nr = 0;
  (*theta) = lp->infinite;
  savef = 0;
  quot = 0;

  for(i = 1; i <= lp->rows; i++) {
    f = pcol[i];
    if(f != 0) {
      if(my_abs(f) < lp->epspivot) {
        if(lp->trace)
  	  report(lp, FULL, "Pivot %g rejected, too small (limit %g)",
		 (double)f, (double)lp->epspivot);
      }
      else { /* pivot alright */
	if(f > 0)
	  quot = lp->rhs[i] / f;