Cutstk - Column generation for a cutting stock problem

/********************************************************
  Xpress-BCL C++ Example Problems
  ===============================

  file xbcutstk.cxx
  `````````````````
  Cutting stock problem, solved by column (= cutting 
  pattern) generation heuristic looping over the
  root node.

  (c) 2008 Fair Isaac Corporation
      author: S.Heipcke, 2001, rev. Mar. 2014
********************************************************/

#include 
#include 
#include "xprb_cpp.h"
#include "xprs.h"

using namespace std;
using namespace ::dashoptimization;

#define NWIDTHS 5
#define MAXWIDTH 94

#define EPS 1e-6
#define MAXCOL 10

/****DATA****/
double WIDTH[] = {17, 21, 22.5,  24, 29.5};  /* Possible widths */
int DEMAND[] =  {150, 96,   48, 108,  227};  /* Demand per width */
int PATTERNS[NWIDTHS][NWIDTHS];              /* (Basic) cutting patterns */

XPRBvar pat[NWIDTHS+MAXCOL];               /* Rolls per pattern */
XPRBctr dem[NWIDTHS];                      /* Demand constraints */
XPRBctr cobj;                              /* Objective function */

XPRBprob p("CutStock");                    /* Initialize a new problem in BCL */

double knapsack(int N, double *c, double *a, double R, int *d, int *xbest);

/***********************************************************************/

void modCutStock()
{
 int i,j;
 XPRBexpr le;

 for(j=0;j= DEMAND[i]);
 }
}

/**************************************************************************/
/*  Column generation loop at the top node:                               */
/*    solve the LP and save the basis                                     */
/*    get the solution values                                             */
/*    generate new column(s) (=cutting pattern)                           */
/*    load the modified problem and load the saved basis                  */
/**************************************************************************/
void solveCutStock()
{
 double objval;                  /* Objective value */
 int i,j;
 int starttime;
 int npatt, npass;               /* Counters for columns and passes */
 double solpat[NWIDTHS+MAXCOL];  /* Solution values for variables pat */
 double dualdem[NWIDTHS];        /* Dual values of demand constraints */
 XPRBbasis basis;
 double dw,z;
 int x[NWIDTHS];

 starttime=XPRB::getTime();
 npatt = NWIDTHS;

 for(npass=0;npass EPS)
    {
     dem[i] += x[i]*pat[npatt];
     if((int)ceil((double)DEMAND[i]/x[i]) > dw)
      dw = (int)ceil((double)DEMAND[i]/x[i]);
    }
   pat[npatt].setUB(dw);           /* Change the upper bound on the new var.*/  

   npatt++;

   p.loadMat();                    /* Reload the problem */
   p.loadBasis(basis);             /* Load the saved basis */
   basis.reset();                  /* No need to keep the basis any longer */
  }
 }

 p.mipOptimize("");                /* Solve the MIP */
  
 cout << "(" << (XPRB::getTime()-starttime)/1000.0 << " sec) Optimal solution: " << p.getObjVal() << " rolls, " << npatt << " patterns" << endl << "   "; 
 cout << "Rolls per pattern: ";
 for(i=0;i