C++ CoinPackedVector::getElements方法代码示例

double LorentzCone::feasibility(int size, CoinPackedVector const & point) const {
  double * dense_point = new double[size]();
  int const * ind = point.getIndices();
  double const * val = point.getElements();
  int num_elem = point.getNumElements();
  for (int i=0; i<num_elem; ++i) {
    dense_point[ind[i]] = val[i];
  }
  double * par_point = new double[size_];
  for(int i=0; i<size_; ++i) {
    par_point[i] = dense_point[members_[i]];
  }
  delete[] dense_point;
  double * p = par_point;
  double term1;
  double term2;
  double feas;
  if (type()==LORENTZ) {
    term1 = p[0];
    term2 = std::inner_product(p+1, p+size_, p+1, 0.0);
    term2 = sqrt(term2);
  }
  else if (type()==RLORENTZ) {
    term1 = 2.0*p[0]*p[1];
    term2 = std::inner_product(p+2, p+size_, p+2, 0.0);
  }
  feas = term1-term2;
  delete[] par_point;
  return feas;
}

// ------------------------------------------------------------------------- //
void UtilPrintPackedVector(const CoinPackedVector& v,
                           ostream*                 os,
                           DecompApp*               app)
{
   (*os).precision(2);
   const int*     inds  = v.getIndices();
   const double* elems = v.getElements();
   const int      len   = v.getNumElements();

   for (int i = 0; i < len; i++) {
      if (!app) {
         (*os) << elems[i] << " x[" << inds[i] << "]  ";
      } else {
         (*os) << elems[i] << " ";
         app->printOriginalColumn(inds[i], os);
         (*os) << "  ";
      }

      if ((i + 1) % 5 == 0) {
         (*os) << "\n";
      }
   }

   (*os) << endl;
}

// take columns of matrix and add each to arrays for matrix under construction
void addSubMatr (int *start, int *len, int *ind, double *el,
		 CoinPackedMatrix &A, 
		 CoinPackedVector &v, 
		 int &cur,
		 int &ncols,
		 int dispM, 
		 int dispVec, 
		 int finalrow) {

  const int 
    *aLe  = A.getVectorLengths (),
    *aIn  = A.getIndices  (),
    *vIn  = v.getIndices  (),
     aCol = A.getMajorDim ();

  int vNum = v.getNumElements ();

  const double 
    *aEl = A.getElements (),
    *vEl = v.getElements ();

  // add each column
  for (int i=0; i<aCol; i++, len++) {

    *start++ = cur;
    *len     = *aLe++;

    // matrix entries
    for (int j = 0; j < *len; j++) {
      *ind++ = dispM + *aIn++;
      *el++  = *aEl++;
    }

    cur += *len;

    // check if there is a corresponding rhs
    if (vNum && (*vIn == i)) {

      ++*len;
      cur++;
      *ind++ = dispVec;
      *el++  = *vEl;

      vIn++;
      vEl++;
      --vNum;
    }

    // normalization entry
    ++*len;
    cur++;
    *ind++ = finalrow;
    *el++  = 1.;

    ++ncols;
  }

  *start = cur;
}

CoinPackedVector::CoinPackedVector(const CoinPackedVector & rhs) :
   CoinPackedVectorBase(),
   indices_(NULL),
   elements_(NULL),
   nElements_(0),
   origIndices_(NULL),
   capacity_(0)
{  
   gutsOfSetVector(rhs.getNumElements(), rhs.getIndices(), rhs.getElements(),
		   rhs.testForDuplicateIndex(), "copy constructor");
}

/// Returns the dot product of two CoinPackedVector objects whose elements are
/// doubles.  Use this version if the vectors are *not* guaranteed to be sorted.
inline double sparseDotProduct(const CoinPackedVectorBase& op1,
                        const CoinPackedVectorBase& op2){
  int len, i;
  double acc = 0.0;
  CoinPackedVector retVal;

  CoinPackedVector retval = op1*op2;
  len = retval.getNumElements();
  double * CParray = retval.getElements();

  for(i = 0; i < len; i++){
    acc += CParray[i];
  }
return acc;
}

bool isWiped (OsiCuts &cs) {

  if (cs.sizeColCuts () == 0)
  //(cs.sizeColCuts () != 1))
    return false;

  CoinPackedVector 
    lbs = cs.colCutPtr (cs.sizeColCuts () - 1) -> lbs (),
    ubs = cs.colCutPtr (cs.sizeColCuts () - 1) -> ubs ();

  return ((lbs.getNumElements ()  == 1)  &&
	  (ubs.getNumElements ()  == 1)  &&
	  (*(lbs.getIndices   ()) == 0)  &&
	  (*(lbs.getElements  ()) == 1.) &&
	  (*(ubs.getIndices   ()) == 0)  &&
	  (*(ubs.getElements  ()) == -1.));
}

/* Insert a row cut unless it is a duplicate (CoinRelFltEq)*/
void 
OsiCuts::insertIfNotDuplicate( OsiRowCut & rc , CoinRelFltEq treatAsSame)
{
  double newLb = rc.lb();
  double newUb = rc.ub();
  CoinPackedVector vector = rc.row();
  int numberElements =vector.getNumElements();
  int * newIndices = vector.getIndices();
  double * newElements = vector.getElements();
  CoinSort_2(newIndices,newIndices+numberElements,newElements);
  bool notDuplicate=true;
  int numberRowCuts = sizeRowCuts();
  for ( int i =0; i<numberRowCuts;i++) {
    const OsiRowCut * cutPtr = rowCutPtr(i);
    if (cutPtr->row().getNumElements()!=numberElements)
      continue;
    if (!treatAsSame(cutPtr->lb(),newLb))
      continue;
    if (!treatAsSame(cutPtr->ub(),newUb))
      continue;
    const CoinPackedVector * thisVector = &(cutPtr->row());
    const int * indices = thisVector->getIndices();
    const double * elements = thisVector->getElements();
    int j;
    for(j=0;j<numberElements;j++) {
      if (indices[j]!=newIndices[j])
	break;
      if (!treatAsSame(elements[j],newElements[j]))
	break;
    }
    if (j==numberElements) {
      notDuplicate=false;
      break;
    }
  }
  if (notDuplicate) {
    OsiRowCut * newCutPtr = new OsiRowCut();
    newCutPtr->setLb(newLb);
    newCutPtr->setUb(newUb);
    newCutPtr->setRow(vector);
    rowCutPtrs_.push_back(newCutPtr);
  }
}

void OsiIF::_row(int i, Row &r) const
{
	const CoinPackedMatrix* coinMatrix;
	CoinPackedVector coinVector;
	int coinNumEl;
	const int* coinIndices;
	const double* coinElements;

	coinMatrix = osiLP_->getMatrixByRow();
	coinVector = coinMatrix->getVector(i);
	coinNumEl = coinVector.getNumElements();
	coinIndices = coinVector.getIndices();
	coinElements = coinVector.getElements();

	r.clear();

	for (int j = 0; j < coinNumEl; j++)
		r.insert(coinIndices[j], coinElements[j]);
	r.sense(osi2csense(rowsense_[i]));
	r.rhs(_rhs(i));
}

//.........这里部分代码省略.........
    siP->setColSolution(mycs);
    const double *colsol = siP->getColSolution(); 
    int k;
    for (k=0; k<nCols; k++){
      xstar[k]=colsol[k];
      complement[k]=0;
    }
    
    // test deriveAKnapsack
    int rind = ( siP->getRowSense()[0] == 'N' ) ? 1 : 0;
    const CoinShallowPackedVector reqdBySunCC = siP->getMatrixByRow()->getVector(rind) ;
    int deriveaknap = kccg.deriveAKnapsack(*siP, cs, krow,b,complement,xstar,rind,reqdBySunCC);
    assert(deriveaknap ==1);
    assert(complement[0]==0);
    assert(complement[1]==1);
    assert(complement[2]==1);
    int inx[3] = {0,1,2};
    double el[3] = {161, 120, 68};
    CoinPackedVector r;
    r.setVector(3,inx,el);
    assert (krow == r);
    //assert (b == 183.0); ????? but x1 and x2 at 1 is valid 
    
    // test findGreedyCover 
    CoinPackedVector cover,remainder;
#if 0
    int findgreedy =  kccg.findGreedyCover( 0, krow, b, xstar, cover, remainder );
    assert( findgreedy == 1 );
    int coveri = cover.getNumElements();
    assert( cover.getNumElements() == 2);
    coveri = cover.getIndices()[0];
    assert( cover.getIndices()[0] == 0);
    assert( cover.getIndices()[1] == 1);
    assert( cover.getElements()[0] == 161.0);
    assert( cover.getElements()[1] == 120.0);
    assert( remainder.getNumElements() == 1);
    assert( remainder.getIndices()[0] == 2);
    assert( remainder.getElements()[0] == 68.0);

    // test liftCoverCut
    CoinPackedVector cut;
    double * rowupper = ekk_rowupper(model);
    double cutRhs = cover.getNumElements() - 1.0;
    kccg.liftCoverCut(b, krow.getNumElements(),
      cover, remainder,
      cut);
    assert ( cut.getNumElements() == 3 );
    assert ( cut.getIndices()[0] == 0 );
    assert ( cut.getIndices()[1] == 1 );
    assert ( cut.getIndices()[2] == 2 );
    assert( cut.getElements()[0] == 1 );
    assert( cut.getElements()[1] == 1 );
    assert( eq(cut.getElements()[2], 0.087719) );
    
    // test liftAndUncomplementAndAdd
    OsiCuts cuts;    
    kccg.liftAndUncomplementAndAdd(*siP.getRowUpper()[0],krow,b,complement,0,
      cover,remainder,cuts);   
    int sizerowcuts = cuts.sizeRowCuts();
    assert ( sizerowcuts== 1 );
    OsiRowCut testRowCut = cuts.rowCut(0);
    CoinPackedVector testRowPV = testRowCut.row(); 
    OsiRowCut sampleRowCut;
    const int sampleSize = 3;
    int sampleCols[sampleSize]={0,1,2};
    double sampleElems[sampleSize]={1.0,-1.0,-0.087719};

/**Clean cut 2, different algorithm. First check the dynamic of the cut if < maxRatio scale to a biggest coef of 1
   otherwise scale it so that biggest coeff is 1 and try removing tinys ( < 1/maxRatio) either succeed or fail */
int
Validator::cleanCut2(OsiRowCut & aCut, const double * solCut, const OsiSolverInterface &si, const CglParam &/* par */,
                     const double * origColLower, const double * origColUpper) const
{
    /** Compute fill-in in si */
    int numcols = si.getNumCols();
    // int numrows = si.getNumRows();
    const double * colLower = (origColLower) ? origColLower : si.getColLower();
    const double * colUpper = (origColUpper) ? origColUpper : si.getColUpper();

    int maxNnz = static_cast<int> ( maxFillIn_ * static_cast<double> (numcols));

    double rhs = aCut.lb();
    assert (aCut.ub()> 1e50);

    CoinPackedVector *vec = const_cast<CoinPackedVector *>(&aCut.row());
    //  vec->sortIncrIndex();

    int * indices = vec->getIndices();
    double * elems = vec->getElements();
    int n = vec->getNumElements();
    if (n==0)
    {
        numRejected_[EmptyCut]++;
        return EmptyCut;
    }
    /** First compute violation if it is too small exit */
    double violation = aCut.violated(solCut);
    if (violation < minViolation_)
        return 1;

    /** Now relax get dynamic and remove tiny elements */
    int offset = 0;
    rhs -= 1e-10;
    double smallest = fabs(rhs);
    double biggest = smallest;
    double veryTiny = 1e-20;
    for (int i = 0 ; i < n ; i++)
    {
        double val = fabs(elems[i]);
        if (val > veryTiny)   //tiny should be very very small
        {
            smallest = std::min(val,smallest);
            biggest = std::max (val,biggest);
        }
    }

    if (biggest > 1e9)
    {
#ifdef DEBUG
        std::cout<<"Whaooo "<<biggest/smallest<<std::endl;
#endif
        numRejected_[BigDynamic]++;
        return BigDynamic;
    }

    //rescale the cut so that biggest is 1e1.
    double toBeBiggest = rhsScale_;
    rhs *= (toBeBiggest / biggest);
    toBeBiggest /= biggest;
    for (int i = 0 ; i < n ; i++)
    {
        elems[i] *= toBeBiggest;
    }


    if (biggest > maxRatio_ * smallest)   //we have to remove some small coefficients
    {
        double myTiny = biggest * toBeBiggest / maxRatio_;
        veryTiny *= toBeBiggest ;
        for (int i = 0 ; i < n ; i++)
        {
            double val = fabs(elems[i]);
            if (val < myTiny)
            {
                if (val< veryTiny)
                {
                    offset++;
                    continue;
                }
                int & iCol = indices[i];
                if (elems[i]>0. && colUpper[iCol] < 1000.)
                {
                    offset++;
                    rhs -= elems[i] * colUpper[iCol];
                    elems[i]=0;
                }
                else if (elems[i]<0. && colLower[iCol] > -1000.)
                {
                    offset++;
                    rhs -= elems[i] * colLower[iCol];
                    elems[i]=0.;
                }
                else
                {
                    numRejected_[SmallCoefficient]++;
                    return SmallCoefficient;
                }
//.........这里部分代码省略.........

void
CoinShallowPackedVectorUnitTest()
{
  CoinRelFltEq eq;
  int i;
  // Test default constructor
  {
    CoinShallowPackedVector r;
    assert( r.indices_==NULL );
    assert( r.elements_==NULL );
    assert( r.nElements_==0 );
  }

  // Test set and get methods
  const int ne = 4;
  int inx[ne] = { 1, 3, 4, 7 };
  double el[ne] = { 1.2, 3.4, 5.6, 7.8 };
  {
    CoinShallowPackedVector r;    
    assert( r.getNumElements()==0 );
    
    // Test setting/getting elements with int* & double* vectors
    r.setVector( ne, inx, el );
    assert( r.getNumElements()==ne );
    for ( i=0; i<ne; i++ ) {
      assert( r.getIndices()[i]  == inx[i] );
      assert( r.getElements()[i] == el[i]  );
    }
    assert ( r.getMaxIndex()==7 );
    assert ( r.getMinIndex()==1 );

    // try to clear it
    r.clear();
    assert( r.indices_==NULL );
    assert( r.elements_==NULL );
    assert( r.nElements_==0 );

    // Test setting/getting elements with indices out of order  
    const int ne2 = 5;
    int inx2[ne2] = { 2, 4, 8, 14, 3 };
    double el2[ne2] = { 2.2, 4.4, 6.6, 8.8, 3.3 };
 
    r.setVector(ne2,inx2,el2);
    
    assert( r.getNumElements()==ne2 );    
    for (i = 0; i < ne2; ++i) {
       assert( r.getIndices()[i]==inx2[i] );
       assert( r.getElements()[i]==el2[i] );
    }
    
    assert ( r.getMaxIndex()==14 );
    assert ( r.getMinIndex()==2 );
    // try to call it once more
    assert ( r.getMaxIndex()==14 );
    assert ( r.getMinIndex()==2 );

    CoinShallowPackedVector r1(ne2,inx2,el2);
    assert( r == r1 );

    // assignment operator
    r1.clear();
    r1 = r;
    assert( r == r1 );

    // assignment from packed vector
    CoinPackedVector pv1(ne2,inx2,el2);
    r1 = pv1;
    assert( r == r1 );

    // construction
    CoinShallowPackedVector r2(r1);
    assert( r2 == r );
    
    // construction from packed vector
    CoinShallowPackedVector r3(pv1);
    assert( r3 == r );

    // test duplicate indices
    {
      const int ne3 = 4;
      int inx3[ne3] = { 2, 4, 2, 3 };
      double el3[ne3] = { 2.2, 4.4, 8.8, 6.6 };
      r.setVector(ne3,inx3,el3, false);
      assert(r.testForDuplicateIndex() == false);
      bool errorThrown = false;
      try {
        r.setTestForDuplicateIndex(true);
      }
      catch (CoinError& e) {
        errorThrown = true;
      }
      assert( errorThrown );

      r.clear();
      errorThrown = false;
      try {
	 r.setVector(ne3,inx3,el3);
      }
      catch (CoinError& e) {
        errorThrown = true;
//.........这里部分代码省略.........

//--------------------------------------------------------------------------
// test EKKsolution methods.
void
CglProbingUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
# ifdef CGL_DEBUG
  int i ;	// define just once
# endif
  CoinRelFltEq eq(0.000001);

  // Test default constructor
  {
    CglProbing aGenerator;
  }
  
  // Test copy & assignment
  {
    CglProbing rhs;
    {
      CglProbing bGenerator;
      CglProbing cGenerator(bGenerator);
      rhs=bGenerator;
    }
  }

  {
    OsiCuts osicuts;
    CglProbing test1;
    OsiSolverInterface  * siP = baseSiP->clone();
    int nColCuts;
    int nRowCuts;
    
    std::string fn = mpsDir+"p0033";
    siP->readMps(fn.c_str(),"mps");
    siP->initialSolve();
    // just unsatisfied variables
    test1.generateCuts(*siP,osicuts);
    nColCuts = osicuts.sizeColCuts();
    nRowCuts = osicuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" probing cuts"<<std::endl;
    {
      std::cout<<"there are "<<nColCuts<<" probing column cuts"<<std::endl;

#ifdef CGL_DEBUG
      const double * lo = siP->getColLower();
      const double * up = siP->getColUpper();
      for (i=0; i<nColCuts; i++){
	OsiColCut ccut;
	CoinPackedVector cpv;
	ccut = osicuts.colCut(i);
	cpv = ccut.lbs();
	int n = cpv.getNumElements();
        int j;
	const int * indices = cpv.getIndices();
	double* elements = cpv.getElements();
	for (j=0;j<n;j++) {
	  int icol=indices[j];
	  if (elements[j]>lo[icol])
	    std::cout<<"Can increase lb on "<<icol<<" from "<<lo[icol]<<
	      " to "<<elements[j]<<std::endl;
	}
	cpv = ccut.ubs();
	n = cpv.getNumElements();
	indices = cpv.getIndices();
	elements = cpv.getElements();

	for (j=0;j<n;j++) {
	  int icol=indices[j];
	  if (elements[j]<up[icol])
	    std::cout<<"Can decrease ub on "<<icol<<" from "<<up[icol]<<
	      " to "<<elements[j]<<std::endl;
	}
      }

#endif

    }

#ifdef CGL_DEBUG
    for (i=0; i<nRowCuts; i++){
      OsiRowCut rcut;
      CoinPackedVector rpv;
      const double * colsol = siP->getColSolution();
      rcut = osicuts.rowCut(i);
      rpv = rcut.row();
      const int n = rpv.getNumElements();
      const int * indices = rpv.getIndices();
      double* elements = rpv.getElements();
      double sum2=0.0;
      int k=0;
      double lb=rcut.lb();
      double ub=rcut.ub();
      for (k=0; k<n; k++){
	int column=indices[k];
	sum2 += colsol[column]*elements[k];
      }
      if (sum2 >ub + 1.0e-7 ||sum2 < lb - 1.0e-7) {
	std::cout<<"Cut "<<i<<" lb "<<lb<<" solution "<<sum2<<" ub "<<ub<<std::endl;
//.........这里部分代码省略.........

//.........这里部分代码省略.........
            int start[2] = {0,4};
            int length[2] = {4,4};
            int rows[8] = {0,1,2,3,0,1,2,3};
            double elements[8] = {2.0,-2.0,7.0,-7.0,2.0,2.0,4.0,4.0};
            CoinPackedMatrix  columnCopy(true,4,2,8,elements,rows,start,length);

            double rowLower[4]={-COIN_DBL_MAX,-COIN_DBL_MAX,
                                -COIN_DBL_MAX,-COIN_DBL_MAX};
            double rowUpper[4]={3.,1.,8.,1.};
            double colLower[2]={0.0,0.0};
            double colUpper[2]={1.0,1.0};
            double obj[2]={-1,-1};
            int intVar[2]={0,1};

            OsiSolverInterface  * siP = si->clone();
            siP->loadProblem(columnCopy, colLower, colUpper, obj, rowLower, rowUpper);
            siP->setInteger(intVar,2);
            CglLandP test;
            test.setLogLevel(2);
            test.parameter().sepSpace = CglLandP::Full;
            siP->resolve();
            // Test generateCuts method
            {
                OsiCuts cuts;
                test.generateCuts(*siP,cuts);
                cuts.printCuts();
                assert(cuts.sizeRowCuts()==1);
                OsiRowCut aCut = cuts.rowCut(0);
                assert(eq(aCut.lb(), -.0714286));
                CoinPackedVector row = aCut.row();
                if (row.getNumElements() == 1)
                {
                    assert(row.getIndices()[0]==1);
                    assert(eq(row.getElements()[0], -4*.0714286));
                }
                else if (row.getNumElements() == 2)
                {
                    assert(row.getIndices()[0]==0);
                    assert(eq(row.getElements()[0], 0.));
                    assert(row.getIndices()[1]==1);
                    assert(eq(row.getElements()[1], -1));
                }
                OsiSolverInterface::ApplyCutsReturnCode rc = siP->applyCuts(cuts);

                siP->resolve();
            }
            if (0)
            {
                OsiCuts cuts;
                test.generateCuts(*siP,cuts);
                cuts.printCuts();
                assert(cuts.sizeRowCuts()==1);
                OsiRowCut aCut = cuts.rowCut(0);
                CoinPackedVector row = aCut.row();
                if (row.getNumElements() == 1)
                {
                    assert(row.getIndices()[0]==1);
                    assert(eq(row.getElements()[0], -1));
                }
                else if (row.getNumElements() == 2)
                {
                    assert(row.getIndices()[0]==0);
                    assert(eq(row.getElements()[0], 0.));
                    assert(row.getIndices()[1]==1);
                    assert(eq(row.getElements()[1], -1));
                }

//--------------------------------------------------------------------------
// ** At present this does not use any solver
void
CglGomoryUnitTest(
  const OsiSolverInterface * baseSiP,
  const std::string mpsDir )
{
  CoinRelFltEq eq(0.000001);

  // Test default constructor
  {
    CglGomory aGenerator;
    assert (aGenerator.getLimit()==50);
    assert (aGenerator.getAway()==0.05);
  }
  
  // Test copy & assignment etc
  {
    CglGomory rhs;
    {
      CglGomory bGenerator;
      bGenerator.setLimit(99);
      bGenerator.setAway(0.2);
      CglGomory cGenerator(bGenerator);
      rhs=bGenerator;
      assert (rhs.getLimit()==99);
      assert (rhs.getAway()==0.2);
    }
  }

  // Test explicit form - all integer (pg 125 Wolsey)
  if (1) {
    OsiCuts osicuts;
    CglGomory test1;
    int i;
    int nOldCuts=0,nRowCuts;
 
    // matrix data
    //deliberate hiccup of 2 between 0 and 1
    CoinBigIndex start[5]={0,4,7,8,9};
    int length[5]={2,3,1,1,1};
    int rows[11]={0,2,-1,-1,0,1,2,0,1,2};
    double elements[11]={7.0,2.0,1.0e10,1.0e10,-2.0,1.0,-2.0,1,1,1};
    CoinPackedMatrix matrix(true,3,5,8,elements,rows,start,length);
    
    // rim data (objective not used just yet)
    double rowLower[5]={14.0,3.0,3.0,1.0e10,1.0e10};
    double rowUpper[5]={14.0,3.0,3.0,-1.0e10,-1.0e10};
    double colLower[7]={0.0,0.0,0.0,0.0,0.0,0.0,0.0};
    double colUpper[7]={100.0,100.0,100.0,100.0,100.0,100.0,100.0};
  
    // integer
    char intVar[7]={2,2,2,2,2,2,2};

    // basis 1
    int rowBasis1[3]={-1,-1,-1};
    int colBasis1[5]={1,1,-1,-1,1};
    CoinWarmStartBasis warm;
    warm.setSize(5,3);
    for (i=0;i<3;i++) {
      if (rowBasis1[i]<0) {
	warm.setArtifStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setArtifStatus(i,CoinWarmStartBasis::basic);
      }
    }
    for (i=0;i<5;i++) {
      if (colBasis1[i]<0) {
	warm.setStructStatus(i,CoinWarmStartBasis::atLowerBound);
      } else {
	warm.setStructStatus(i,CoinWarmStartBasis::basic);
      }
    }

    // solution 1
    double colsol1[5]={20.0/7.0,3.0,0.0,0.0,23.0/7.0};
    test1.generateCuts(NULL, osicuts, matrix,
		       /*objective,*/ colsol1,
		 colLower, colUpper,
		 rowLower, rowUpper, intVar, &warm);
    nRowCuts = osicuts.sizeRowCuts();
    std::cout<<"There are "<<nRowCuts<<" gomory cuts"<<std::endl;
    assert (nRowCuts==2);
    // cuts always <=
    int testCut=0; // test first cut as stronger
    double rhs=-6.0;
    double testCut1[5]={0.0,0.0,-1.0,-2.0,0.0};
    double * cut = testCut1;
    double * colsol = colsol1;
    for (i=nOldCuts; i<nRowCuts; i++){
      OsiRowCut rcut;
      CoinPackedVector rpv;
      rcut = osicuts.rowCut(i);
      rpv = rcut.row();
      const int n = rpv.getNumElements();
      const int * indices = rpv.getIndices();
      double* elements = rpv.getElements();
      double sum2=0.0;
      int k=0;
      for (k=0; k<n; k++){
//.........这里部分代码省略.........

//.........这里部分代码省略.........
		if (cutrow.getNumElements() == 0) continue;

		/** is optimality cut? */
		bool isOptimalityCut = false;
		for (int j = nvars_ - naux_; j < nvars_; ++j)
		{
			if (cutrow.getMaxIndex() == j)
			{
				isOptimalityCut = true;
				break;
			}
		}

		double efficacy = rc->violated(vals) / cutrow.twoNorm();
		SCIP_Bool isEfficacious = efficacy > 1.e-6;

#define KK_TEST
#ifdef KK_TEST
		if (SCIPgetStage(scip) == SCIP_STAGE_INITSOLVE ||
			SCIPgetStage(scip) == SCIP_STAGE_SOLVING)
		{
			/** create empty row */
			SCIP_ROW * row = NULL;
			SCIP_CALL(SCIPcreateEmptyRowCons(scip, &row, conshdlr, "benders", rc->lb(), SCIPinfinity(scip),
					FALSE, /**< is row local? */
					FALSE, /**< is row modifiable? */
					FALSE  /**< is row removable? can this be TRUE? */));

			/** cache the row extension and only flush them if the cut gets added */
			SCIP_CALL(SCIPcacheRowExtensions(scip, row));

			/** collect all non-zero coefficients */
			for (int j = 0; j < cutrow.getNumElements(); ++j)
				SCIP_CALL(SCIPaddVarToRow(scip, row, vars_[cutrow.getIndices()[j]], cutrow.getElements()[j]));

			DSPdebugMessage("found Benders (%s) cut: act=%f, lhs=%f, norm=%f, eff=%f, min=%f, max=%f (range=%f)\n",
				isOptimalityCut ? "opti" : "feas",
				SCIPgetRowLPActivity(scip, row), SCIProwGetLhs(row), SCIProwGetNorm(row),
				SCIPgetCutEfficacy(scip, sol, row),
				SCIPgetRowMinCoef(scip, row), SCIPgetRowMaxCoef(scip, row),
				SCIPgetRowMaxCoef(scip, row)/SCIPgetRowMinCoef(scip, row));

			/** flush all changes before adding cut */
			SCIP_CALL(SCIPflushRowExtensions(scip, row));

			DSPdebugMessage("efficacy %e isEfficatious %d\n", efficacy, isEfficacious);

			if (isEfficacious)
			{
				if (where == from_scip_sepalp ||
					where == from_scip_sepasol ||
					where == from_scip_enfolp)
				{
					/** add cut */
					SCIP_Bool infeasible;
					SCIP_CALL(SCIPaddCut(scip, sol, row,
							FALSE, /**< force cut */
							&infeasible));

					if (infeasible)
						*result = SCIP_CUTOFF;
					else //if (*result != SCIP_CUTOFF)
						*result = SCIP_SEPARATED;
				}
				else
					*result = SCIP_INFEASIBLE;