useVolume.cpp

// Copyright (C) 2003, International Business Machines
// Corporation and others.  All Rights Reserved.

#include "ClpSimplex.hpp"
#include "ClpFactorization.hpp"
#include "VolVolume.hpp"

//#############################################################################

class lpHook : public VOL_user_hooks {
private:
   lpHook(const lpHook&);
   lpHook& operator=(const lpHook&);
private:
   double  *colupper_;
   double  *collower_;
   double  *objcoeffs_;
   double  *rhs_;
   char    *sense_;

   CoinPackedMatrix rowMatrix_;
   CoinPackedMatrix colMatrix_;

public:
   lpHook(const double* clb, const double* cub, const double* obj,
          const double* rhs, const char* sense, const CoinPackedMatrix& mat);
   virtual ~lpHook();
   
public:
   // for all hooks: return value of -1 means that volume should quit
   virtual int compute_rc(const VOL_dvector& u, VOL_dvector& rc);

   virtual int solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
                                double& lcost, VOL_dvector& x, VOL_dvector& v,
                                double& pcost);
   virtual int heuristics(const VOL_problem& p, 
                          const VOL_dvector& x, double& heur_val) {
      return 0;
   }
};
 
//#############################################################################

lpHook::lpHook(const double* clb, const double* cub, const double* obj,
               const double* rhs, const char* sense,
               const CoinPackedMatrix& mat)
{
   const int colnum = mat.getNumCols();
   const int rownum = mat.getNumRows();

   colupper_ = new double[colnum];
   collower_ = new double[colnum];
   objcoeffs_ = new double[colnum];
   rhs_ = new double[rownum];
   sense_ = new char[rownum];

   std::copy(clb, clb + colnum, collower_);
   std::copy(cub, cub + colnum, colupper_);
   std::copy(obj, obj + colnum, objcoeffs_);
   std::copy(rhs, rhs + rownum, rhs_);
   std::copy(sense, sense + rownum, sense_);

   if (mat.isColOrdered()) {
      colMatrix_.copyOf(mat);
      rowMatrix_.reverseOrderedCopyOf(mat);
   } else {
      rowMatrix_.copyOf(mat);
      colMatrix_.reverseOrderedCopyOf(mat);
   }
}

//-----------------------------------------------------------------------------

lpHook::~lpHook()
{
   delete[] colupper_;
   delete[] collower_;
   delete[] objcoeffs_;
   delete[] rhs_;
   delete[] sense_;
}

//#############################################################################

int
lpHook::compute_rc(const VOL_dvector& u, VOL_dvector& rc)
{
   rowMatrix_.transposeTimes(u.v, rc.v);
   const int psize = rowMatrix_.getNumCols();

   for (int i = 0; i < psize; ++i)
      rc[i] = objcoeffs_[i] - rc[i];
   return 0;
}

//-----------------------------------------------------------------------------

int
lpHook::solve_subproblem(const VOL_dvector& dual, const VOL_dvector& rc,
                         double& lcost, VOL_dvector& x, VOL_dvector& v,
                         double& pcost)
{
   int i;
   const int psize = x.size();
   const int dsize = v.size();

   // compute the lagrangean solution corresponding to the reduced costs
   for (i = 0; i < psize; ++i)
      x[i] = (rc[i] >= 0.0) ? collower_[i] : colupper_[i];

   // compute the lagrangean value (rhs*dual + primal*rc)
   lcost = 0;
   for (i = 0; i < dsize; ++i)
      lcost += rhs_[i] * dual[i];
   for (i = 0; i < psize; ++i)
      lcost += x[i] * rc[i];

   // compute the rhs - lhs 
   colMatrix_.times(x.v, v.v);
   for (i = 0; i < dsize; ++i)
      v[i] = rhs_[i] - v[i];

   // compute the lagrangean primal objective
   pcost = 0;
   for (i = 0; i < psize; ++i)
      pcost += x[i] * objcoeffs_[i];

   return 0;
}

//#############################################################################

int main (int argc, const char *argv[])
{
  ClpSimplex  model;
  int status;

  if (argc<2)
    status=model.readMps("../../Mps/Sample/p0033.mps",true);
  else
    status=model.readMps(argv[1],true);
  /*
    This driver uses volume algorithm
    then does dual - after adjusting costs
    then solves real problem
  */

  // do volume for a bit
  VOL_problem volprob;
  const CoinPackedMatrix* mat = model.matrix();
  const int psize = mat->getNumCols();
  const int dsize = mat->getNumRows();
  char * sense = new char[dsize];
  double * rhs = new double[dsize];
  const double * rowLower = model.rowLower();
  const double * rowUpper = model.rowUpper();
  // Set the lb/ub on the duals
  volprob.dsize = dsize;
  volprob.psize = psize;
  volprob.dual_lb.allocate(dsize);
  volprob.dual_ub.allocate(dsize);
  volprob.dsol.allocate(dsize);
  int i;
  for (i = 0; i < dsize; ++i) {
    if (rowUpper[i]==rowLower[i]) {
    // 'E':
      volprob.dual_lb[i] = -1.0e31;
      volprob.dual_ub[i] = 1.0e31;
      rhs[i]=rowUpper[i];
      sense[i]='E';
    } else if (rowLower[i]<-0.99e10&&rowUpper[i]<0.99e10) {
      // 'L':
      volprob.dual_lb[i] = -1.0e31;
      volprob.dual_ub[i] = 0.0;
      rhs[i]=rowUpper[i];
      sense[i]='L';
    } else if (rowLower[i]>-0.99e10&&rowUpper[i]>0.99e10) {
      // 'G':
      volprob.dual_lb[i] = 0.0;
      volprob.dual_ub[i] = 1.0e31;
      rhs[i]=rowLower[i];
      sense[i]='G';
    } else {
      printf("Volume Algorithm can't work if there is a non ELG row\n");
      abort();
    }
  }
  // Can't use read_param as private
  // anyway I want automatic use - so maybe this is problem
#if 0
  FILE* infile = fopen("parameters", "r");
  if (!infile) {
    printf("Failure to open parameter file\n");
  } else {
    volprob.read_params("parameters");
  }
#endif
#if 0
  // should save and restore bounds
  model.tightenPrimalBounds();
#else
  double * colUpper = model.columnUpper();
  for (i=0;i<psize;i++)
    colUpper[i]=1.0;
#endif
  lpHook myHook(model.getColLower(), model.getColUpper(),
                model.getObjCoefficients(),
                rhs, sense, *mat);
  // move duals
  double * pi = model.dualRowSolution();
  memcpy(volprob.dsol.v,pi,dsize*sizeof(double));
  volprob.solve(myHook,  false /* not warmstart */);
  // For now stop as not doing any good
  exit(77);
  // create objectives
  int numberRows = model.numberRows();
  int numberColumns = model.numberColumns();
  memcpy(pi,volprob.dsol.v,numberRows*sizeof(double));
#define MODIFYCOSTS
#ifdef MODIFYCOSTS
  double * saveObj = new double[numberColumns];
  memcpy(saveObj,model.objective(),numberColumns*sizeof(double));
  memcpy(model.dualColumnSolution(),model.objective(),
         numberColumns*sizeof(double));
  model.clpMatrix()->transposeTimes(-1.0,pi,model.dualColumnSolution());
  memcpy(model.objective(),model.dualColumnSolution(),
         numberColumns*sizeof(double));
  const double * rowsol = model.primalRowSolution();
  //const double * rowLower = model.rowLower();
  //const double * rowUpper = model.rowUpper();
  double offset=0.0;
  for (i=0;i<numberRows;i++) {
    offset += pi[i]*rowsol[i];
  }
  double value2;
  model.getDblParam(ClpObjOffset,value2);
  printf("Offset %g %g\n",offset,value2);
  model.setRowObjective(pi);
  // zero out pi
  memset(pi,0,numberRows*sizeof(double));
#endif
  // Could put some in basis - only partially tested
  model.allSlackBasis(); 
  model.factorization()->maximumPivots(1000);
  //model.setLogLevel(63);
  // solve
  model.dual(1);
  //model.primal(1);
#ifdef MODIFYCOSTS
  memcpy(model.objective(),saveObj,numberColumns*sizeof(double));
  // zero out pi
  memset(pi,0,numberRows*sizeof(double));
  model.setRowObjective(pi);
  delete [] saveObj;
  model.primal();
#endif
  
  return 0;
}    

---