1 // Copyright (C) 2002, International Business Machines
2 // Corporation and others. All Rights Reserved.
3 // This code is licensed under the terms of the Eclipse Public License (EPL).
4
5 #include <cstdlib>
6 #include <cstdio>
7 #include <cmath>
8 #include <cfloat>
9 #include <cassert>
10 #include <iostream>
11 //#define CGL_DEBUG 1
12 #ifdef NDEBUG
13 //#undef NDEBUG
14 #endif
15 #include "CoinPragma.hpp"
16 #include "CoinHelperFunctions.hpp"
17 #include "CoinPackedVector.hpp"
18 #include "CoinPackedMatrix.hpp"
19 #include "CoinIndexedVector.hpp"
20 #include "OsiRowCutDebugger.hpp"
21 #ifndef USE_CGL_RATIONAL
22 #define USE_CGL_RATIONAL 0
23 #endif
24 // -1 no rational, 0 as before, 1 use code from CglGmi
25 #if USE_CGL_RATIONAL>0
26 #if USE_CGL_RATIONAL<=10
27 #error "If USE_CGL_RATIONAL>0 must be at least 10 (maybe not more than 1000)"
28 #endif
29 #include "CoinRational.hpp"
30 #endif
31 #define COIN_HAS_CLP_GOMORY
32 #ifdef COIN_HAS_CLP_GOMORY
33 #include "OsiClpSolverInterface.hpp"
34 #endif
35 #include "CoinFactorization.hpp"
36 #undef CLP_OSL
37 #if COIN_BIG_INDEX==0
38 #define CLP_OSL 1
39 #if CLP_OSL!=1&&CLP_OSL!=3
40 #undef CLP_OSL
41 #else
42 #include "CoinOslFactorization.hpp"
43 #endif
44 #endif
45 #include "CoinWarmStartBasis.hpp"
46 #include "CglGomory.hpp"
47 #include "CoinFinite.hpp"
48 #ifdef CGL_DEBUG_GOMORY
49 int gomory_try=CGL_DEBUG_GOMORY;
50 #endif
51 //-------------------------------------------------------------------
52 // Generate Gomory cuts
53 //-------------------------------------------------------------------
generateCuts(const OsiSolverInterface & si,OsiCuts & cs,const CglTreeInfo info)54 void CglGomory::generateCuts(const OsiSolverInterface & si, OsiCuts & cs,
55 const CglTreeInfo info)
56 {
57 #ifdef CGL_DEBUG_GOMORY
58 gomory_try++;
59 #endif
60 // Get basic problem information
61 int numberColumns=si.getNumCols();
62
63 // get integer variables and basis
64 char * intVar = new char[numberColumns];
65 int i;
66 CoinWarmStart * warmstart = si.getWarmStart();
67 CoinWarmStartBasis* warm =
68 dynamic_cast<CoinWarmStartBasis*>(warmstart);
69 const double * colUpper = si.getColUpper();
70 const double * colLower = si.getColLower();
71 //#define CLP_INVESTIGATE2
72 #ifndef CLP_INVESTIGATE2
73 if ((info.options&16)!=0)
74 #endif
75 printf("%d %d %d\n",info.inTree,info.options,info.pass);
76 for (i=0;i<numberColumns;i++) {
77 if (si.isInteger(i)) {
78 if (colUpper[i]>colLower[i]+0.5) {
79 if (fabs(colUpper[i]-1.0)<1.0e-12&&fabs(colLower[i])<1.0e-12) {
80 intVar[i]=1; //0-1
81 } else if (colLower[i]>=0.0) {
82 intVar[i] = 2; // other
83 } else {
84 // negative bounds - I am not sure works
85 intVar[i] = 3;
86 }
87 } else {
88 intVar[i] = 4;
89 }
90 } else {
91 intVar[i]=0;
92 }
93 }
94 const OsiSolverInterface * useSolver=&si;
95 #ifdef COIN_HAS_CLP_GOMORY
96 double * objective = NULL;
97 OsiClpSolverInterface * clpSolver = dynamic_cast<OsiClpSolverInterface *>(originalSolver_);
98 int numberOriginalRows = -1;
99 if (clpSolver) {
100 useSolver = originalSolver_;
101 assert (gomoryType_);
102 // check simplex is plausible
103 if (!clpSolver->getNumRows()||numberColumns!=clpSolver->getNumCols()) {
104 delete originalSolver_;
105 originalSolver_=si.clone();
106 clpSolver = dynamic_cast<OsiClpSolverInterface *>(originalSolver_);
107 assert (clpSolver);
108 useSolver = originalSolver_;
109 }
110 ClpSimplex * simplex = clpSolver->getModelPtr();
111 numberOriginalRows = simplex->numberRows();
112 int numberRows = si.getNumRows();
113 assert (numberOriginalRows<=numberRows);
114 // only do if different (unless type 2x)
115 int gomoryType = gomoryType_%10;
116 int whenToDo = gomoryType_/10;
117 if (whenToDo==2 ||(numberRows>numberOriginalRows && whenToDo==1
118 && (info.options&512)==0) ||
119 ((info.options&1024)!=0 && (info.options&512)==0
120 && numberTimesStalled_<3)) {
121 // bounds
122 memcpy(simplex->columnLower(),colLower,numberColumns*sizeof(double));
123 memcpy(simplex->columnUpper(),colUpper,numberColumns*sizeof(double));
124 double * obj = simplex->objective();
125 objective = CoinCopyOfArray(obj,numberColumns);
126 const double * pi = si.getRowPrice();
127 const CoinPackedMatrix * rowCopy = si.getMatrixByRow();
128 const int * column = rowCopy->getIndices();
129 const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
130 const int * rowLength = rowCopy->getVectorLengths();
131 const double * rowElements = rowCopy->getElements();
132 const double * rowLower = si.getRowLower();
133 const double * rowUpper = si.getRowUpper();
134 int numberCopy;
135 int numberAdd;
136 double * rowLower2 = NULL;
137 double * rowUpper2 = NULL;
138 int * column2 = NULL;
139 CoinBigIndex * rowStart2 = NULL;
140 double * rowElements2 = NULL;
141 char * copy = new char [numberRows-numberOriginalRows];
142 memset(copy,0,numberRows-numberOriginalRows);
143 if (gomoryType==2) {
144 numberCopy=0;
145 numberAdd=0;
146 for (int iRow=numberOriginalRows;iRow<numberRows;iRow++) {
147 bool simple = true;
148 for (CoinBigIndex k=rowStart[iRow];
149 k<rowStart[iRow]+rowLength[iRow];k++) {
150 double value = rowElements[k];
151 if (value!=floor(value+0.5)) {
152 simple=false;
153 break;
154 }
155 }
156 if (simple) {
157 numberCopy++;
158 numberAdd+=rowLength[iRow];
159 copy[iRow-numberOriginalRows]=1;
160 }
161 }
162 if (numberCopy) {
163 //printf("Using %d rows out of %d\n",numberCopy,numberRows-numberOriginalRows);
164 rowLower2 = new double [numberCopy];
165 rowUpper2 = new double [numberCopy];
166 rowStart2 = new CoinBigIndex [numberCopy+1];
167 rowStart2[0]=0;
168 column2 = new int [numberAdd];
169 rowElements2 = new double [numberAdd];
170 }
171 }
172 numberCopy=0;
173 numberAdd=0;
174 const double * rowSolution = si.getRowActivity();
175 double offset=0.0;
176 for (int iRow=numberOriginalRows;iRow<numberRows;iRow++) {
177 if (!copy[iRow-numberOriginalRows]) {
178 double value = pi[iRow];
179 offset += rowSolution[iRow]*value;
180 for (CoinBigIndex k=rowStart[iRow];
181 k<rowStart[iRow]+rowLength[iRow];k++) {
182 int iColumn=column[k];
183 obj[iColumn] -= value*rowElements[k];
184 }
185 } else {
186 rowLower2[numberCopy]=rowLower[iRow];
187 rowUpper2[numberCopy]=rowUpper[iRow];
188 for (CoinBigIndex k=rowStart[iRow];
189 k<rowStart[iRow]+rowLength[iRow];k++) {
190 column2[numberAdd]=column[k];
191 rowElements2[numberAdd++]=rowElements[k];
192 }
193 numberCopy++;
194 rowStart2[numberCopy]=numberAdd;
195 }
196 }
197 #if 0
198 CoinThreadRandom randomNumberGenerator;
199 const double * solution = si.getColSolution();
200 for (int i=0;i<numberColumns;i++) {
201 if (intVar[i]==1) {
202 double randomNumber = randomNumberGenerator.randomDouble();
203 //randomNumber = 0.001*floor(randomNumber*1000.0);
204 if (solution[i]>0.5)
205 obj[i] -= randomNumber*0.001*fabs(obj[i]);
206 else
207 obj[i] += randomNumber*0.001*fabs(obj[i]);
208 }
209 }
210 #endif
211 if (numberCopy) {
212 clpSolver->addRows(numberCopy,
213 rowStart2,column2,rowElements2,
214 rowLower2,rowUpper2);
215 delete [] rowLower2 ;
216 delete [] rowUpper2 ;
217 delete [] column2 ;
218 delete [] rowStart2 ;
219 delete [] rowElements2 ;
220 }
221 delete [] copy;
222 memcpy(simplex->primalColumnSolution(),si.getColSolution(),
223 numberColumns*sizeof(double));
224 warm->resize(numberOriginalRows,numberColumns);
225 clpSolver->setBasis(*warm);
226 delete warm;
227 simplex->setDualObjectiveLimit(COIN_DBL_MAX);
228 simplex->setLogLevel(0);
229 simplex->primal(1);
230 // check basis
231 int numberTotal=simplex->numberRows()+simplex->numberColumns();
232 int superbasic=0;
233 for (int i=0;i<numberTotal;i++) {
234 if (simplex->getStatus(i)==ClpSimplex::superBasic)
235 superbasic++;
236 }
237 if (superbasic) {
238 //printf("%d superbasic!\n",superbasic);
239 simplex->dual();
240 superbasic=0;
241 for (int i=0;i<numberTotal;i++) {
242 if (simplex->getStatus(i)==ClpSimplex::superBasic)
243 superbasic++;
244 }
245 assert (!superbasic);
246 }
247 //printf("Trying - %d its status %d objs %g %g - with offset %g\n",
248 // simplex->numberIterations(),simplex->status(),
249 // simplex->objectiveValue(),si.getObjValue(),simplex->objectiveValue()+offset);
250 //simplex->setLogLevel(0);
251 warm=simplex->getBasis();
252 warmstart=warm;
253 if (simplex->status()) {
254 //printf("BAD status %d\n",simplex->status());
255 //clpSolver->writeMps("clp");
256 //si.writeMps("si");
257 delete [] objective;
258 objective=NULL;
259 useSolver=&si;
260 }
261 } else {
262 // don't do
263 delete warmstart;
264 warmstart=NULL;
265 if ((info.options&1024)==0)
266 numberTimesStalled_=0;
267 }
268 }
269 #endif
270 #ifdef CGL_DEBUG
271 const OsiRowCutDebugger * debugger = si.getRowCutDebugger();
272 if (debugger&&!debugger->onOptimalPath(si))
273 debugger = NULL;
274 #else
275 const OsiRowCutDebugger * debugger = NULL;
276 #endif
277 int numberRowCutsBefore = cs.sizeRowCuts();
278
279 if (warmstart)
280 generateCuts(debugger, cs, *useSolver->getMatrixByCol(),
281 *useSolver->getMatrixByRow(),
282 useSolver->getColSolution(),
283 useSolver->getColLower(), useSolver->getColUpper(),
284 useSolver->getRowLower(), useSolver->getRowUpper(),
285 intVar,warm,info);
286 #ifdef COIN_HAS_CLP_GOMORY
287 if (objective) {
288 ClpSimplex * simplex = clpSolver->getModelPtr();
289 memcpy(simplex->objective(),objective,numberColumns*sizeof(double));
290 delete [] objective;
291 // take out locally useless cuts
292 const double * solution = si.getColSolution();
293 double primalTolerance = 1.0e-7;
294 int numberRowCutsAfter = cs.sizeRowCuts();
295 for (int k = numberRowCutsAfter - 1; k >= numberRowCutsBefore; k--) {
296 const OsiRowCut * thisCut = cs.rowCutPtr(k) ;
297 double sum = 0.0;
298 int n = thisCut->row().getNumElements();
299 const int * column = thisCut->row().getIndices();
300 const double * element = thisCut->row().getElements();
301 assert (n);
302 for (int i = 0; i < n; i++) {
303 double value = element[i];
304 sum += value * solution[column[i]];
305 }
306 if (sum > thisCut->ub() + primalTolerance) {
307 sum = sum - thisCut->ub();
308 } else if (sum < thisCut->lb() - primalTolerance) {
309 sum = thisCut->lb() - sum;
310 } else {
311 sum = 0.0;
312 }
313 if (!sum) {
314 // take out
315 cs.eraseRowCut(k);
316 }
317 }
318 #ifdef CLP_INVESTIGATE2
319 printf("OR %p pass %d inTree %c - %d cuts (but %d deleted)\n",
320 originalSolver_,info.pass,info.inTree?'Y':'N',
321 numberRowCutsAfter-numberRowCutsBefore,
322 numberRowCutsAfter-cs.sizeRowCuts());
323 #endif
324 }
325 #endif
326
327 delete warmstart;
328 delete [] intVar;
329 if ((!info.inTree&&((info.options&4)==4||((info.options&8)&&!info.pass)))
330 ||(info.options&16)!=0) {
331 int limit = maximumLengthOfCutInTree();
332 int numberRowCutsAfter = cs.sizeRowCuts();
333 for (int i=numberRowCutsBefore;i<numberRowCutsAfter;i++) {
334 int length = cs.rowCutPtr(i)->row().getNumElements();
335 if (length<=limit)
336 cs.rowCutPtr(i)->setGloballyValid();
337 }
338 }
339 if ((gomoryType_%10)==2) {
340 // back to original
341 assert(clpSolver);
342 int numberRows = clpSolver->getNumRows();
343 if (numberRows>numberOriginalRows) {
344 int numberDelete = numberRows-numberOriginalRows;
345 int * delRow = new int [numberDelete];
346 for (int i=0;i<numberDelete;i++)
347 delRow[i]=i+numberOriginalRows;
348 clpSolver->deleteRows(numberDelete,delRow);
349 delete [] delRow;
350 }
351 }
352 }
353
354 // Returns value - floor but allowing for small errors
above_integer(double value)355 inline double above_integer(double value) {
356 double value2=floor(value);
357 double value3=floor(value+0.5);
358 if (fabs(value3-value)<1.0e-9*(fabs(value3)+1.0))
359 return 0.0;
360 return value-value2;
361 }
362 #if USE_CGL_RATIONAL <= 0
363 //-------------------------------------------------------------------
364 // Returns the greatest common denominator of two
365 // positive integers, a and b, found using Euclid's algorithm
366 //-------------------------------------------------------------------
gcd(long long int a,long long int b)367 static long long int gcd(long long int a, long long int b)
368 {
369 long long int remainder = -1;
370 #if CGL_DEBUG>1
371 printf("gcd of %ld and %ld\n",a,b);
372 int nLoop=0;
373 #endif
374 // make sure a<=b (will always remain so)
375 if(a > b) {
376 // Swap a and b
377 long long int temp = a;
378 a = b;
379 b = temp;
380 }
381 // if zero then gcd is nonzero (zero may occur in rhs of packed)
382 if (!a) {
383 if (b) {
384 return b;
385 } else {
386 printf("**** gcd given two zeros!!\n");
387 abort();
388 }
389 }
390 while (remainder) {
391
392 #if CGL_DEBUG>1
393 nLoop++;
394 if (nLoop>50) {
395 abort();
396 return -1;
397 }
398 #endif
399 remainder = b % a;
400 b = a;
401 a = remainder;
402 }
403 #if CGL_DEBUG>1
404 printf("=> %d\n",b);
405 #endif
406 return b;
407 }
408 #endif
409 #define GOMORY_LONG
410 #ifdef GOMORY_LONG
411 #define SMALL_VALUE1 1.0e-14
412 #define GOMORY_INT long long int
413 #else
414 #define SMALL_VALUE1 1.0e-10
415 #define GOMORY_INT int
416 #endif
417 #if USE_CGL_RATIONAL>0
computeGcd(long a,long b)418 static long computeGcd(long a, long b) {
419 // This is the standard Euclidean algorithm for gcd
420 long remainder = 1;
421 // Make sure a<=b (will always remain so)
422 if (a > b) {
423 // Swap a and b
424 long temp = a;
425 a = b;
426 b = temp;
427 }
428 // If zero then gcd is nonzero
429 if (!a) {
430 if (b) {
431 return b;
432 }
433 else {
434 printf("### WARNING: CglGMI::computeGcd() given two zeroes!\n");
435 exit(1);
436 }
437 }
438 while (remainder) {
439 remainder = b % a;
440 b = a;
441 a = remainder;
442 }
443 return b;
444 } /* computeGcd */
scaleCutIntegral(double * cutElem,int * cutIndex,int cutNz,double & cutRhs,double maxdelta)445 static bool scaleCutIntegral(double* cutElem, int* cutIndex, int cutNz,
446 double& cutRhs, double maxdelta) {
447 long gcd, lcm;
448 double maxscale = 1000;
449 long maxdnom = USE_CGL_RATIONAL;
450 //long numerator = 0, denominator = 0;
451 // Initialize gcd and lcm
452 CoinRational r = CoinRational(cutRhs, maxdelta, maxdnom);
453 if (r.getNumerator() != 0){
454 gcd = labs(r.getNumerator());
455 lcm = r.getDenominator();
456 }
457 else{
458 #if defined GMI_TRACE_CLEAN
459 printf("Cannot compute rational number, scaling procedure aborted\n");
460 #endif
461 return false;
462 }
463 for (int i = 0; i < cutNz; ++i) {
464 CoinRational r = CoinRational(cutElem[i], maxdelta, maxdnom);
465 if (r.getNumerator() != 0){
466 gcd = computeGcd(gcd, r.getNumerator());
467 lcm *= r.getDenominator()/(computeGcd(lcm,r.getDenominator()));
468 }
469 else{
470 #if defined GMI_TRACE_CLEAN
471 printf("Cannot compute rational number, scaling procedure aborted\n");
472 #endif
473 return false;
474 }
475 }
476 double scale = ((double)lcm)/((double)gcd);
477 if (fabs(scale) > maxscale) {
478 #if defined GMI_TRACE_CLEAN
479 printf("Scaling factor too large, scaling procedure aborted\n");
480 #endif
481 return false;
482 }
483 scale = fabs(scale);
484 // Looks like we have a good scaling factor; scale and return;
485 for (int i = 0; i < cutNz; ++i) {
486 double value = cutElem[i]*scale;
487 cutElem[i] = floor(value+0.5);
488 assert (fabs(cutElem[i]-value)<1.0e-9);
489 }
490 {
491 double value = cutRhs*scale;
492 cutRhs = floor(value+0.5);
493 assert (fabs(cutRhs-value)<1.0e-9);
494 }
495 return true;
496 } /* scaleCutIntegral */
497 #elif USE_CGL_RATIONAL == 0
498 //-------------------------------------------------------------------
499 // Returns the nearest rational with denominator < maxDenominator
500 //-------------------------------------------------------------------
501 typedef struct {
502 int numerator;
503 int denominator;
504 } Rational;
505 typedef struct {
506 GOMORY_INT numerator;
507 GOMORY_INT denominator;
508 } RationalLong;
nearestRational(double value,int maxDenominator)509 inline Rational nearestRational(double value, int maxDenominator)
510 {
511 RationalLong tryThis;
512 RationalLong tryA;
513 RationalLong tryB;
514 double integerPart;
515
516 #if CGL_DEBUG>1
517 printf("Rational of %g is ",value);
518 #endif
519 int nLoop=0;
520
521 tryA.numerator=0;
522 tryA.denominator=1;
523 tryB.numerator=1;
524 tryB.denominator=0;
525 Rational returnRational;
526
527 if (fabs(value)<SMALL_VALUE1) {
528 returnRational.numerator=static_cast<int>(tryA.numerator);
529 returnRational.denominator=static_cast<int>(tryA.denominator);
530 return returnRational;
531 }
532 integerPart = floor(value);
533 value -= integerPart;
534 tryThis.numerator = tryB.numerator* static_cast<GOMORY_INT> (integerPart) + tryA.numerator;
535 tryThis.denominator = tryB.denominator* static_cast<GOMORY_INT> (integerPart) + tryA.denominator;
536 tryA = tryB;
537 tryB = tryThis;
538
539 while (value>SMALL_VALUE1 && tryB.denominator <=maxDenominator) {
540 nLoop++;
541 if (nLoop>50) {
542 Rational bad;
543 bad.numerator=-1;
544 bad.denominator=-1;
545 #if CGL_DEBUG>1
546 printf(" *** bad rational\n");
547 #endif
548 return bad;
549 }
550 value = 1.0/value;
551 integerPart = floor(value+SMALL_VALUE1);
552 value -= integerPart;
553 tryThis.numerator = tryB.numerator* static_cast<GOMORY_INT> (integerPart) + tryA.numerator;
554 tryThis.denominator = tryB.denominator* static_cast<GOMORY_INT>(integerPart) + tryA.denominator;
555 tryA = tryB;
556 tryB = tryThis;
557 }
558 if (tryB.denominator <= maxDenominator) {
559 #if CGL_DEBUG>1
560 printf("%lld/%lld\n",tryB.numerator,tryB.denominator);
561 #endif
562 if (tryB.denominator<COIN_INT_MAX) {
563 returnRational.numerator=static_cast<int>(tryB.numerator);
564 returnRational.denominator=static_cast<int>(tryB.denominator);
565 } else {
566 returnRational.numerator=-1;
567 returnRational.denominator=-1;
568 #if CGL_DEBUG>1
569 printf(" *** bad rational\n");
570 #endif
571 }
572 } else {
573 #if CGL_DEBUG>1
574 printf("%lld/%lld\n",tryA.numerator,tryA.denominator);
575 #endif
576 if (tryA.denominator<COIN_INT_MAX) {
577 returnRational.numerator=static_cast<int>(tryA.numerator);
578 returnRational.denominator=static_cast<int>(tryA.denominator);
579 } else {
580 returnRational.numerator=-1;
581 returnRational.denominator=-1;
582 #if CGL_DEBUG>1
583 printf(" *** bad rational\n");
584 #endif
585 }
586 }
587 return returnRational;
588 }
589 #endif
590 // Does actual work - returns number of cuts
591 int
generateCuts(const OsiRowCutDebugger * debugger,OsiCuts & cs,const CoinPackedMatrix & columnCopy,const CoinPackedMatrix & rowCopy,const double * colsol,const double * colLower,const double * colUpper,const double * rowLower,const double * rowUpper,const char * intVar,const CoinWarmStartBasis * warm,const CglTreeInfo info)592 CglGomory::generateCuts(
593 #ifdef CGL_DEBUG
594 const OsiRowCutDebugger * debugger,
595 #else
596 const OsiRowCutDebugger * ,
597 #endif
598 OsiCuts & cs,
599 const CoinPackedMatrix & columnCopy,
600 const CoinPackedMatrix & rowCopy,
601 const double * colsol,
602 const double * colLower, const double * colUpper,
603 const double * rowLower, const double * rowUpper,
604 const char * intVar,
605 const CoinWarmStartBasis* warm,
606 const CglTreeInfo info)
607 {
608 int infoOptions=info.options;
609 bool globalCuts = (infoOptions&16)!=0;
610 double testFixed = (!globalCuts) ? 1.0e-8 : -1.0;
611 // get what to look at
612 double away = info.inTree ? away_ : CoinMin(away_,awayAtRoot_);
613 int numberRows=columnCopy.getNumRows();
614 int numberColumns=columnCopy.getNumCols();
615 CoinBigIndex numberElements=columnCopy.getNumElements();
616 // Allow bigger length on initial matrix (if special setting)
617 //if (limit==512&&!info.inTree&&!info.pass)
618 //limit=1024;
619 // Start of code to create a factorization from warm start (A) ====
620 // check factorization is okay
621 CoinFactorization factorization;
622 #ifdef CLP_OSL
623 CoinOslFactorization * factorization2=NULL;
624 if (alternateFactorization_) {
625 factorization2 = new CoinOslFactorization();
626 }
627 #endif
628 // We can either set increasing rows so ...IsBasic gives pivot row
629 // or we can just increment iBasic one by one
630 // for now let ...iBasic give pivot row
631 int status=-100;
632 // probably could use pivotVariables from OsiSimplexModel
633 int * rowIsBasic = new int[numberRows];
634 int * columnIsBasic = new int[numberColumns];
635 int i;
636 int numberBasic=0;
637 for (i=0;i<numberRows;i++) {
638 if (warm->getArtifStatus(i) == CoinWarmStartBasis::basic) {
639 rowIsBasic[i]=1;
640 numberBasic++;
641 } else {
642 rowIsBasic[i]=-1;
643 }
644 }
645 for (i=0;i<numberColumns;i++) {
646 if (warm->getStructStatus(i) == CoinWarmStartBasis::basic) {
647 columnIsBasic[i]=1;
648 numberBasic++;
649 } else {
650 columnIsBasic[i]=-1;
651 }
652 }
653 //returns 0 -okay, -1 singular, -2 too many in basis, -99 memory */
654 while (status<-98) {
655 #ifdef CLP_OSL
656 if (!alternateFactorization_) {
657 #endif
658 status=factorization.factorize(columnCopy,
659 rowIsBasic, columnIsBasic);
660 if (status==-99) factorization.areaFactor(factorization.areaFactor() * 2.0);
661 #ifdef CLP_OSL
662 } else {
663 double areaFactor=1.0;
664 status=factorization2->factorize(columnCopy,
665 rowIsBasic, columnIsBasic,areaFactor);
666 if (status==-99) areaFactor *= 2.0;
667 }
668 #endif
669 }
670 if (status) {
671 #ifdef COIN_DEVELOP
672 std::cout<<"Bad factorization of basis - status "<<status<<std::endl;
673 #endif
674 delete [] rowIsBasic;
675 delete [] columnIsBasic;
676 return -1;
677 }
678 // End of creation of factorization (A) ====
679
680 #ifdef CLP_OSL
681 double relaxation = !alternateFactorization_ ? factorization.conditionNumber() :
682 factorization2->conditionNumber();
683 #else
684 double relaxation = factorization.conditionNumber();
685 #endif
686 // if very small be a bit more careful
687 if (relaxation<1.0e-10)
688 relaxation=1.0/sqrt(relaxation);
689 #ifdef COIN_DEVELOP_z
690 if (relaxation>1.0e49)
691 printf("condition %g\n",relaxation);
692 #endif
693 //printf("condition %g %g\n",relaxation,conditionNumberMultiplier_);
694 relaxation *= conditionNumberMultiplier_;
695 double bounds[2]={-COIN_DBL_MAX,0.0};
696 int iColumn,iRow;
697
698 const int * column = rowCopy.getIndices();
699 const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
700 const int * rowLength = rowCopy.getVectorLengths();
701 const double * rowElements = rowCopy.getElements();
702 const int * row = columnCopy.getIndices();
703 const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
704 const int * columnLength = columnCopy.getVectorLengths();
705 const double * columnElements = columnCopy.getElements();
706
707 // we need to do book-keeping for variables at ub
708 double tolerance = 1.0e-7;
709 bool * swap= new bool [numberColumns];
710 for (iColumn=0;iColumn<numberColumns;iColumn++) {
711 if (columnIsBasic[iColumn]<0&&
712 colUpper[iColumn]-colsol[iColumn]<=tolerance) {
713 swap[iColumn]=true;
714 } else {
715 swap[iColumn]=false;
716 }
717 }
718
719 // get row activities (could use solver but lets do here )
720 double * rowActivity = new double [numberRows];
721 CoinFillN(rowActivity,numberRows,0.0);
722 for (iColumn=0;iColumn<numberColumns;iColumn++) {
723 double value = colsol[iColumn];
724 CoinBigIndex k;
725 for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];k++) {
726 iRow = row[k];
727 rowActivity[iRow] += columnElements[k]*value;
728 }
729 }
730 /* we need to mark rows:
731 0) equality
732 1) slack at lb (activity at ub)
733 2) slack at ub (activity at lb)
734 and 4 bit is set if slack must be integer
735
736 */
737 int * rowType = new int[numberRows];
738 for (iRow=0;iRow<numberRows;iRow++) {
739 if (rowIsBasic[iRow]<0&&rowUpper[iRow]>rowLower[iRow]+1.0e-7) {
740 int type=0;
741 double rhs=0.0;
742 if (rowActivity[iRow]>=rowUpper[iRow]-1.0e-7) {
743 type=1;
744 rhs=rowUpper[iRow];
745 } else if (rowActivity[iRow]<=rowLower[iRow]+1.0e-7) {
746 type=2;
747 rhs=rowLower[iRow];
748 } else {
749 // probably large rhs
750 if (rowActivity[iRow]-rowLower[iRow]<
751 rowUpper[iRow]-rowActivity[iRow])
752 rowType[iRow]=2;
753 else
754 rowType[iRow]=1;
755 #ifdef CGL_DEBUG
756 assert (CoinMin(rowUpper[iRow]-rowActivity[iRow],
757 rowActivity[iRow]-rowUpper[iRow])<1.0e-5);
758 //abort();
759 continue;
760 #else
761 continue;
762 #endif
763 }
764 if (above_integer(rhs)<1.0e-10) {
765 // could be integer slack
766 bool allInteger=true;
767 CoinBigIndex k;
768 for (k=rowStart[iRow];
769 k<rowStart[iRow]+rowLength[iRow];k++) {
770 int iColumn=column[k];
771 if (!intVar[iColumn]||above_integer(rowElements[k])>1.0e-10) {
772 // not integer slacks
773 allInteger=false;
774 break;
775 }
776 }
777 if (allInteger) {
778 type |= 4;
779 }
780 }
781 rowType[iRow]=type;
782 } else {
783 // row is equality or basic
784 rowType[iRow]=0;
785 }
786 }
787
788 // Start of code to create work arrays for factorization (B) ====
789 // two vectors for updating (one is work)
790 CoinIndexedVector work;
791 CoinIndexedVector array;
792 // make sure large enough
793 work.reserve(numberRows);
794 array.reserve(numberRows);
795 int * arrayRows = array.getIndices();
796 double * arrayElements = array.denseVector();
797 // End of code to create work arrays (B) ====
798
799 int numberAdded=0;
800 // we also need somewhere to accumulate cut
801 CoinIndexedVector cutVector;
802 cutVector.reserve(numberColumns+1);
803 int * cutIndex = cutVector.getIndices();
804 double * cutElement = cutVector.denseVector();
805 // and for packed form (as not necessarily in order)
806 // also space for sort
807 bool doSorted = (infoOptions&256)!=0;
808 int lengthArray = static_cast<int>(numberColumns+1+((numberColumns+1)*sizeof(int))/sizeof(double));
809 if (doSorted)
810 lengthArray+=numberColumns;
811 double * packed = new double[lengthArray];
812 double * sort = packed+numberColumns+1;
813 int * which = reinterpret_cast<int *>(doSorted ? (sort+numberColumns): (sort));
814 double tolerance1=1.0e-6;
815 double tolerance2=0.9;
816 #if USE_CGL_RATIONAL <= 0
817 double tolerance3=1.0e-10;
818 #endif
819 double tolerance6=1.0e-6;
820 double tolerance9=1.0e-4;
821 #define MORE_GOMORY_CUTS 1
822 #ifdef CLP_INVESTIGATE2
823 int saveLimit = info.inTree ? 50 : 1000;
824 #else
825 #if MORE_GOMORY_CUTS==2||MORE_GOMORY_CUTS==3
826 int saveLimit;
827 #endif
828 #endif
829 // get limit on length of cut
830 int limit = 0;
831 if (!limit_)
832 dynamicLimitInTree_ = CoinMax(50,numberColumns/40);
833 if (!info.inTree) {
834 limit = limitAtRoot_;
835 if (!info.pass) {
836 tolerance1=1.0;
837 tolerance2=1.0e-2;
838 #if USE_CGL_RATIONAL <= 0
839 //tolerance3=1.0e-6;
840 #endif
841 tolerance6=1.0e-7;
842 tolerance9=1.0e-5;
843 if (!limit)
844 limit=numberColumns;
845 } else {
846 if((infoOptions&32)==0/*&&numberTimesStalled_<3*/) {
847 if (!limit) {
848 if(numberElements>8*numberColumns)
849 limit=numberColumns;
850 else
851 limit = CoinMax(1000,numberRows/4);
852 }
853 } else {
854 limit=numberColumns;
855 numberTimesStalled_++;
856 }
857 }
858 } else {
859 limit = limit_;
860 if (!limit) {
861 if (!info.pass)
862 limit = dynamicLimitInTree_;
863 else
864 limit=50;
865 }
866 }
867 // If big - allow for rows
868 if (limit>=numberColumns)
869 limit += numberRows;
870 #ifdef CLP_INVESTIGATE2
871 if (limit>saveLimit&&!info.inTree&&(infoOptions&512)==0)
872 printf("Gomory limit changed from %d to %d, inTree %c, pass %d, r %d,c %d,e %d\n",
873 saveLimit,limit,info.inTree ? 'Y' : 'N',info.pass,
874 numberRows,numberColumns,numberElements);
875 #endif
876 int nCandidates=0;
877 for (iColumn=0;iColumn<numberColumns;iColumn++) {
878 // This returns pivot row for columns or -1 if not basic (C) ====
879 int iBasic=columnIsBasic[iColumn];
880 if (iBasic>=0&&intVar[iColumn]) {
881 double reducedValue=above_integer(colsol[iColumn]);
882 //printf("col %d bas %d val %.18g\n",iColumn,iBasic,colsol[iColumn]);
883 if(intVar[iColumn]&&reducedValue<1.0-away&&reducedValue>away) {
884 if (doSorted)
885 sort[nCandidates]=fabs(0.5-reducedValue);
886 which[nCandidates++]=iColumn;
887 }
888 }
889 }
890 CoinBigIndex nTotalEls=COIN_INT_MAX;
891 if (doSorted) {
892 CoinSort_2(sort,sort+nCandidates,which);
893 CoinBigIndex nElsNow = columnCopy.getNumElements();
894 CoinBigIndex nAdd;
895 CoinBigIndex nAdd2;
896 CoinBigIndex nReasonable;
897 int depth=info.level;
898 if (depth<2) {
899 nAdd=10000;
900 if (info.pass>0)
901 nAdd = CoinMin(nAdd,nElsNow+2*numberRows);
902 nAdd2 = 5*numberColumns;
903 nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
904 if (!depth&&!info.pass) {
905 // allow more
906 nAdd += nElsNow/2;
907 nAdd2 += nElsNow/2;
908 nReasonable += nElsNow/2;
909 limit=numberRows+numberColumns;
910 }
911 } else {
912 nAdd = 200;
913 nAdd2 = 2*numberColumns;
914 nReasonable = CoinMax(nAdd2,nElsNow/8+nAdd);
915 }
916 nTotalEls=nReasonable;
917 }
918 #ifdef MORE_GOMORY_CUTS
919 CoinBigIndex saveTotalEls=nTotalEls;
920 #endif
921 #if MORE_GOMORY_CUTS==2||MORE_GOMORY_CUTS==3
922 saveLimit=limit;
923 if (doSorted)
924 limit=numberRows+numberColumns;
925 OsiCuts longCuts;
926 #endif
927 #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
928 OsiCuts secondaryCuts;
929 #endif
930 for (int kColumn=0;kColumn<nCandidates;kColumn++) {
931 if (nTotalEls<=0)
932 break; // Got enough
933 iColumn=which[kColumn];
934 double reducedValue=above_integer(colsol[iColumn]);;
935 // This returns pivot row for columns or -1 if not basic (C) ====
936 int iBasic=columnIsBasic[iColumn];
937 double ratio=reducedValue/(1.0-reducedValue);
938 if (iBasic>=0) {
939 // Debug code below computes tableau column of basic ====
940 int j;
941 #ifdef CGL_DEBUG
942 {
943 // put column into array
944 array.setVector(columnLength[iColumn],row+columnStart[iColumn],
945 columnElements+columnStart[iColumn]);
946 // get column in tableau
947 #ifdef CLP_OSL
948 if (!alternateFactorization_)
949 #endif
950 factorization.updateColumn ( &work, &array );
951 #ifdef CLP_OSL
952 else
953 factorization2->updateColumn ( &work, &array );
954 #endif
955 int nn=0;
956 int numberInArray=array.getNumElements();
957 for (j=0;j<numberInArray;j++) {
958 int indexValue=arrayRows[j];
959 double value=arrayElements[indexValue];
960 if (fabs(value)>1.0e-5) {
961 assert (fabs(value-1.0)<1.0e-7);
962 assert (indexValue==iBasic);
963 nn++;
964 }
965 }
966 assert (nn==1);
967 array.clear();
968 //work.checkClear();
969 }
970 #endif
971 array.clear();
972 assert(intVar[iColumn]&&reducedValue<1.0-away&&reducedValue>away);
973 {
974 #ifdef CGL_DEBUG
975 //cutVector.checkClear();
976 #endif
977 // get row of tableau
978 double one =1.0;
979 array.setVector(1,&iBasic,&one);
980 int numberNonInteger=0;
981 //Code below computes tableau row ====
982 // get pi
983 #ifdef CLP_OSL
984 if (!alternateFactorization_)
985 #endif
986 factorization.updateColumnTranspose ( &work, &array );
987 #ifdef CLP_OSL
988 else
989 factorization2->updateColumnTranspose ( &work, &array );
990 #endif
991 int numberInArray=array.getNumElements();
992 #ifdef CGL_DEBUG
993 // check pivot on iColumn
994 {
995 double value=0.0;
996 int k;
997 // add in row of tableau
998 for (k=columnStart[iColumn];
999 k<columnStart[iColumn]+columnLength[iColumn];k++) {
1000 iRow = row[k];
1001 value += columnElements[k]*arrayElements[iRow];
1002 }
1003 // should be 1
1004 assert (fabs(value-1.0) < 1.0e-7);
1005 }
1006 #endif
1007 double largestFactor=0.0;
1008 for (j=0;j<numberInArray;j++) {
1009 int indexValue=arrayRows[j];
1010 double value=arrayElements[indexValue];
1011 largestFactor = CoinMax(largestFactor,fabs(value));
1012 }
1013 //reducedValue=colsol[iColumn];
1014 // coding from pg 130 of Wolsey
1015 // adjustment to rhs
1016 double rhs=0.0;
1017 int number=0;
1018 // number of terms (so includes modifications and cancellations)
1019 int numberCoefficients=0;
1020 #ifdef CGL_DEBUG_GOMORY
1021 if (!gomory_try)
1022 printf("start for basic column %d\n",iColumn);
1023 #endif
1024 // columns
1025 for (j=0;j<numberColumns;j++) {
1026 if (columnIsBasic[j]<0&&colUpper[j]>colLower[j]+testFixed) {
1027 double value=0.0;
1028 CoinBigIndex k;
1029 // add in row of tableau
1030 for (k=columnStart[j];k<columnStart[j]+columnLength[j];k++) {
1031 iRow = row[k];
1032 double value2 = columnElements[k]*arrayElements[iRow];
1033 largestFactor = CoinMax(largestFactor,fabs(value2));
1034 value += value2;
1035 }
1036 #ifdef CGL_DEBUG_GOMORY
1037 if (!gomory_try&&value)
1038 printf("col %d alpha %g colsol %g swap %c bounds %g %g\n",
1039 j,value,colsol[j],swap[j] ? 'Y' : 'N',
1040 colLower[j],colUpper[j]);
1041 #endif
1042 // value is entry in tableau row end (C) ====
1043 if (fabs(value)<1.0e-16) {
1044 // small value
1045 continue;
1046 } else {
1047 // left in to stop over compilation?
1048 //if (iColumn==-52) printf("for basic %d, column %d has alpha %g, colsol %g\n",
1049 // iColumn,j,value,colsol[j]);
1050 #if CGL_DEBUG>1
1051 if (iColumn==52) printf("for basic %d, column %d has alpha %g, colsol %g\n",
1052 iColumn,j,value,colsol[j]);
1053 #endif
1054 // deal with bounds
1055 if (swap[j]) {
1056 //reducedValue -= value*colUpper[j];
1057 // negate
1058 value = - value;
1059 } else {
1060 //reducedValue -= value*colLower[j];
1061 }
1062 #if CGL_DEBUG>1
1063 if (iColumn==52) printf("%d value %g reduced %g int %d rhs %g swap %d\n",
1064 j,value,reducedValue,intVar[j],rhs,swap[j]);
1065 #endif
1066 double coefficient;
1067 if (intVar[j]) {
1068 // integer
1069 coefficient = above_integer(value);
1070 if (coefficient > reducedValue) {
1071 coefficient = ratio * (1.0-coefficient);
1072 }
1073 } else {
1074 // continuous
1075 numberNonInteger++;
1076 if (value > 0.0) {
1077 coefficient = value;
1078 } else {
1079 //??? sign wrong in book
1080 coefficient = -ratio*value;
1081 }
1082 }
1083 if (swap[j]) {
1084 // negate
1085 coefficient = - coefficient;
1086 rhs += colUpper[j]*coefficient;
1087 } else {
1088 rhs += colLower[j]*coefficient;
1089 }
1090 if (fabs(coefficient)>= COIN_INDEXED_TINY_ELEMENT) {
1091 cutElement[j] = coefficient;
1092 numberCoefficients++;
1093 cutIndex[number++]=j;
1094 // If too many - break from loop
1095 if (number>limit)
1096 break;
1097 }
1098 }
1099 } else {
1100 // basic
1101 continue;
1102 }
1103 }
1104 cutVector.setNumElements(number);
1105 // If too many - just clear vector and skip
1106 if (number>limit) {
1107 cutVector.clear();
1108 continue;
1109 }
1110 //check will be cut
1111 //reducedValue=above_integer(reducedValue);
1112 rhs += reducedValue;
1113 double violation = reducedValue;
1114 #ifdef CGL_DEBUG
1115 std::cout<<"cut has violation of "<<violation
1116 <<" value "<<colsol[iColumn]<<std::endl;
1117 #endif
1118 // now do slacks part
1119 for (j=0;j<numberInArray;j++) {
1120 iRow=arrayRows[j];
1121 double value = arrayElements[iRow];
1122 int type=rowType[iRow];
1123 if (type&&fabs(value)>=1.0e-16) {
1124 if ((type&1)==0) {
1125 // negate to get correct coefficient
1126 value = - value;
1127 }
1128 double coefficient;
1129 if ((type&4)!=0) {
1130 // integer
1131 coefficient = above_integer(value);
1132 if (coefficient > reducedValue) {
1133 coefficient = ratio * (1.0-coefficient);
1134 }
1135 } else {
1136 numberNonInteger++;
1137 // continuous
1138 if (value > 0.0) {
1139 coefficient = value;
1140 } else {
1141 coefficient = -ratio*value;
1142 }
1143 }
1144 if ((type&1)!=0) {
1145 // slack at ub - treat as +1.0
1146 rhs -= coefficient*rowUpper[iRow];
1147 } else {
1148 // negate yet again ?
1149 coefficient = - coefficient;
1150 rhs -= coefficient*rowLower[iRow];
1151 }
1152 CoinBigIndex k;
1153 for (k=rowStart[iRow];
1154 k<rowStart[iRow]+rowLength[iRow];k++) {
1155 int jColumn=column[k];
1156 double value=rowElements[k];
1157 double oldValue=cutElement[jColumn];
1158 cutVector.quickAdd(jColumn,-coefficient*value);
1159 numberCoefficients++;
1160 if (!intVar[jColumn]&&!oldValue)
1161 numberNonInteger++;
1162 }
1163 }
1164 }
1165 //check again and pack down
1166 // also change signs
1167 // also zero cutElement
1168 double sum=0.0;
1169 rhs = - rhs;
1170 int n = cutVector.getNumElements();
1171 //assert (n); can be 0!
1172 // If too many - just clear vector and skip
1173 if (n>limit) {
1174 cutVector.clear();
1175 continue;
1176 }
1177 #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
1178 double violation2=violation;
1179 #endif
1180 number=0;
1181 numberNonInteger=0;
1182 double sumCoefficients=0.0;
1183 for (j=0;j<n;j++) {
1184 int jColumn =cutIndex[j];
1185 double value=-cutElement[jColumn];
1186 sumCoefficients += fabs(value);
1187 if (fabs(colsol[jColumn])>10.0)
1188 sumCoefficients += 2.0*fabs(value);
1189 cutElement[jColumn]=0.0;
1190 if (fabs(value)>1.0e-8) {
1191 sum+=value*colsol[jColumn];
1192 packed[number]=value;
1193 cutIndex[number++]=jColumn;
1194 if (!intVar[jColumn])
1195 numberNonInteger++;
1196 } else {
1197 #define LARGE_BOUND 1.0e20
1198 // small - adjust rhs if rhs reasonable
1199 if (value>0.0&&colLower[jColumn]>-LARGE_BOUND) {
1200 rhs -= value*colLower[jColumn];
1201 #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
1202 // weaken violation
1203 violation2 -= fabs(value*(colsol[jColumn]-colLower[jColumn]));
1204 #endif
1205 } else if (value<0.0&&colUpper[jColumn]<LARGE_BOUND) {
1206 rhs -= value*colUpper[jColumn];
1207 #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
1208 // weaken violation
1209 violation2 -= fabs(value*(colsol[jColumn]-colUpper[jColumn]));
1210 #endif
1211 } else if (fabs(value)>1.0e-13) {
1212 // take anyway
1213 sum+=value*colsol[jColumn];
1214 packed[number]=value;
1215 cutIndex[number++]=jColumn;
1216 if (!intVar[jColumn])
1217 numberNonInteger++;
1218 }
1219 }
1220 }
1221 // Final test on number
1222 //if (number>limit)
1223 //continue;
1224 // say zeroed out
1225 cutVector.setNumElements(0);
1226 bool accurate2=false;
1227 double difference=fabs((sum-rhs)-violation);
1228 double useTolerance;
1229 if (tolerance1>0.99) {
1230 // use absolute
1231 useTolerance = tolerance;
1232 } else {
1233 double rhs2=CoinMax(fabs(rhs),10.0);
1234 useTolerance=rhs2*0.1*tolerance1;
1235 }
1236 bool accurate = (difference<useTolerance);
1237 #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
1238 double difference2=fabs((sum-rhs)-violation2);
1239 #if MORE_GOMORY_CUTS==1
1240 if (difference2<useTolerance&&doSorted)
1241 #else
1242 if (difference2<useTolerance&&doSorted&&number<saveLimit)
1243 #endif
1244 accurate2=true;
1245 #endif
1246 if (sum >rhs+tolerance2*away&&
1247 (accurate||accurate2)) {
1248 //#ifdef CGL_DEBUG
1249 #ifdef CGL_DEBUG
1250 #if CGL_DEBUG<=1
1251 if (number<=10) {
1252 #endif
1253 std::cout<<"col "<<iColumn;
1254 for (j=0;j<number;j++) {
1255 std::cout<<" ["<<cutIndex[j]<<","<<packed[j]<<"]";
1256 }
1257 std::cout<<" <= "<<rhs<<std::endl;
1258 #if CGL_DEBUG<=1
1259 }
1260 #endif
1261 #endif
1262 bool cleanedCut=numberNonInteger>0;
1263 bool dontRelax = false;
1264 #define TRY7 1
1265 #if TRY7==2
1266 double rhsBeforeRelax=rhs;
1267 #endif
1268 if ((!numberNonInteger||number<6)&&number&&USE_CGL_RATIONAL>=0) {
1269 // pretend not integer
1270 numberNonInteger=0;
1271 #if USE_CGL_RATIONAL==0
1272 #ifdef CGL_DEBUG
1273 assert (sizeof(Rational)==sizeof(double));
1274 #endif
1275 Rational * cleaned = reinterpret_cast<Rational *> (cutElement);
1276 int * xInt = reinterpret_cast<int *> (cutElement);
1277 // cut should have an integer slack so try and simplify
1278 // add in rhs and put in cutElements (remember to zero later)
1279 cutIndex[number]=numberColumns+1;
1280 packed[number]=rhs;
1281 int numberNonSmall=0;
1282 long long int lcm = 1;
1283
1284 for (j=0;j<number+1;j++) {
1285 double value=above_integer(fabs(packed[j]));
1286 if (fabs(value)<tolerance3) {
1287 // very close to integer
1288 continue;
1289 } else {
1290 numberNonSmall++;
1291 }
1292
1293 cleaned[j]=nearestRational(value,100000);
1294 if (cleaned[j].denominator<0) {
1295 // bad rational
1296 lcm=-1;
1297 break;
1298 }
1299 long long int thisGcd = gcd(lcm,cleaned[j].denominator);
1300 // may need to check for overflow
1301 lcm /= thisGcd;
1302 if (static_cast<double>(lcm)*cleaned[j].denominator<
1303 1.0e18) {
1304 lcm *= cleaned[j].denominator;
1305 } else {
1306 lcm=-1;
1307 break;
1308 }
1309 }
1310 if (lcm>0&&numberNonSmall) {
1311 double multiplier = static_cast<double>(lcm);
1312 cleanedCut=true;
1313 int nOverflow = 0;
1314 for (j=0; j<number+1;j++) {
1315 double value = fabs(packed[j]);
1316 double dxInt = value*multiplier;
1317 xInt[j]= static_cast<int> (dxInt+0.5);
1318 #if CGL_DEBUG>1
1319 printf("%g => %g \n",value,dxInt);
1320 #endif
1321 if (dxInt>1.0e9||fabs(dxInt-xInt[j])> 1.0e-8) {
1322 nOverflow++;
1323 break;
1324 }
1325 }
1326
1327 if (nOverflow){
1328 #ifdef CGL_DEBUG
1329 printf("Gomory Scaling: Warning: Overflow detected \n");
1330 #endif
1331 numberNonInteger=-1;
1332 } else {
1333
1334 // find greatest common divisor of the elements
1335 j=0;
1336 while (!xInt[j])
1337 j++; // skip zeros
1338 long long int thisGcd = gcd(xInt[j],xInt[j+1]);
1339 j++;
1340 for (;j<number+1;j++) {
1341 if (xInt[j])
1342 thisGcd = gcd(thisGcd,xInt[j]);
1343 }
1344 #if 0
1345 // Check nothing too illegal - FIX this
1346 for (j=0;j<number+1;j++) {
1347 double old = lcm*packed[j];
1348 int newOne;
1349 if (old>0.0)
1350 newOne=xInt[j]/thisGcd;
1351 else
1352 newOne=-xInt[j]/thisGcd;
1353 if (fabs(((double) newOne)-old)>
1354 1.0e-10*(fabs(newOne)+1.0)) {
1355 // say no good - first see if happens
1356 printf("Fix this test 456 - just skip\n");
1357 abort();
1358 }
1359 }
1360 #endif
1361 #if CGL_DEBUG>1
1362 printf("The gcd of xInt is %i\n",thisGcd);
1363 #endif
1364
1365 dontRelax=true;
1366 // construct new cut by dividing through by gcd and
1367 double minMultiplier=1.0e100;
1368 double maxMultiplier=0.0;
1369 for (j=0; j<number+1;j++) {
1370 double old=packed[j];
1371 if (old>0.0) {
1372 packed[j]=static_cast<double>(xInt[j]/thisGcd);
1373 } else {
1374 packed[j]=-static_cast<double>(xInt[j]/thisGcd);
1375 }
1376 #if CGL_DEBUG>1
1377 printf("%g => %g \n",old,packed[j]);
1378 #endif
1379 if (packed[j]) {
1380 if (fabs(packed[j])>maxMultiplier*fabs(old))
1381 maxMultiplier = packed[j]/old;
1382 if (fabs(packed[j])<minMultiplier*fabs(old))
1383 minMultiplier = packed[j]/old;
1384 }
1385 }
1386 rhs = packed[number];
1387 #ifdef CGL_DEBUG
1388 printf("min, max multipliers - %g, %g\n",
1389 minMultiplier,maxMultiplier);
1390 #endif
1391 assert(maxMultiplier/minMultiplier>0.9999&&maxMultiplier/minMultiplier<1.0001);
1392 }
1393 } else if (!numberNonSmall) {
1394 dontRelax=true;
1395 }
1396 // erase cutElement
1397 CoinFillN(cutElement,number+1,0.0);
1398 #elif USE_CGL_RATIONAL>0
1399 // Use coding from CglGmi
1400 double maxdelta = 1.0e-12; //param.getEPS();
1401 if (0) {
1402 double total=0.0;
1403 for (j=0;j<number;j++) {
1404 int jColumn =cutIndex[j];
1405 double value=packed[j];
1406 std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn]
1407 <<") ";
1408 total += value*colsol[jColumn];
1409 }
1410 std::cout<<std::endl;
1411 printf("before %g <= %g <= %g\n",bounds[0],total,rhs);
1412 }
1413 dontRelax=scaleCutIntegral(packed,cutIndex,number,
1414 rhs,maxdelta);
1415 if (dontRelax&&false) {
1416 double total=0.0;
1417 for (j=0;j<number;j++) {
1418 int jColumn =cutIndex[j];
1419 double value=packed[j];
1420 std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn]
1421 <<") ";
1422 total += value*colsol[jColumn];
1423 }
1424 std::cout<<std::endl;
1425 printf("after %g <= %g <= %g\n",bounds[0],total,rhs);
1426 }
1427 #endif
1428 }
1429 if (!dontRelax) {
1430 // relax rhs a tiny bit
1431 //#define CGL_GOMORY_OLD_RELAX
1432 #ifndef CGL_GOMORY_OLD_RELAX
1433 #if 0
1434 double rhs2=rhs;
1435 rhs2 += 1.0e-8;
1436 // relax if lots of elements for mixed gomory
1437 if (number>=20) {
1438 rhs2 += 1.0e-7*(static_cast<double> (number/20));
1439 }
1440 #endif
1441 rhs += 1.0e-7;
1442 if (numberCoefficients>=10||true) {
1443 rhs += 1.0e-7*sumCoefficients+1.0e-8*numberCoefficients;
1444 }
1445 #if 0
1446 if (numberCoefficients>number*3)
1447 printf("old rhs %.18g new %.18g - n,nNon,nC,sumC %d,%d,%d %g\n",
1448 rhs2,rhs,number,numberNonInteger,numberCoefficients,
1449 sumCoefficients);
1450 #endif
1451 #else
1452 rhs += 1.0e-8;
1453 // relax if lots of elements for mixed gomory
1454 if (number>=20) {
1455 rhs += 1.0e-7*(static_cast<double> (number/20));
1456 }
1457 #endif
1458 }
1459 // Take off tiny elements
1460 // for first pass reject
1461 #ifndef CGL_GOMORY_TINY_ELEMENT
1462 #define CGL_GOMORY_TINY_ELEMENT 1.0e-12
1463 #endif
1464 {
1465 int i,number2=number;
1466 number=0;
1467 double largest=0.0;
1468 double smallest=1.0e30;
1469 for (i=0;i<number2;i++) {
1470 double value=fabs(packed[i]);
1471 if (value<CGL_GOMORY_TINY_ELEMENT) {
1472 int iColumn = cutIndex[i];
1473 if (colUpper[iColumn]-colLower[iColumn]<LARGE_BOUND) {
1474 if (intVar[iColumn])
1475 numberNonInteger--;
1476 // weaken cut
1477 if (packed[i]>0.0)
1478 rhs -= value*colLower[iColumn];
1479 else
1480 rhs += value*colUpper[iColumn];
1481 } else {
1482 // throw away
1483 number=limit+1;
1484 numberNonInteger=1;
1485 break;
1486 }
1487 } else {
1488 int iColumn = cutIndex[i];
1489 if (colUpper[iColumn]!=colLower[iColumn]||globalCuts) {
1490 value=fabs(value);
1491 largest=CoinMax(largest,value);
1492 smallest=CoinMin(smallest,value);
1493 cutIndex[number]=cutIndex[i];
1494 packed[number++]=packed[i];
1495 } else {
1496 // fixed so subtract out
1497 rhs -= packed[i]*colLower[iColumn];
1498 if (intVar[iColumn])
1499 numberNonInteger--;
1500 }
1501 }
1502 }
1503 if (largest>1.0e10*smallest||(number>20&&smallest<number*1.0e-6)||
1504 numberNonInteger<-10||!number) {
1505 //assert (number); // debug this - looks as if it could be
1506 number=limit+1; //reject
1507 numberNonInteger=1;
1508 } else if (largest>1.0e9*smallest) {
1509 #ifdef CLP_INVESTIGATE2
1510 printf("WOuld reject %g %g ratio %g\n",smallest,largest,
1511 smallest/largest);
1512 #endif
1513 #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
1514 accurate=false;
1515 #endif
1516 } else {
1517 #define PRINT_NUMBER 0
1518 #if PRINT_NUMBER
1519 if (number==PRINT_NUMBER) {
1520 printf("==========\n<= %.18g ",rhs);
1521 for (int i=0;i<PRINT_NUMBER;i++)
1522 printf("%.18g ",packed[i]);
1523 printf("\n");
1524 }
1525 #endif
1526 if (number>limit)
1527 continue;
1528 #if TRY7==1
1529 // Just scale - unless small integers
1530 bool dontScale=true;
1531 if (fabs(rhs-floor(rhs+0.5))>1.0e-14)
1532 dontScale=false;
1533 for (int i=0;i<number;i++) {
1534 if (fabs(packed[i]-floor(packed[i]+0.5))>1.0e-14||
1535 fabs(packed[i])>1.0e4)
1536 dontScale=false;
1537 }
1538 if (!dontScale) {
1539 double multiplier = 1.0/sqrt(largest*smallest);
1540 for (int i=0;i<number;i++)
1541 packed[i] *= multiplier;
1542 rhs *= multiplier;
1543 }
1544 #if PRINT_NUMBER
1545 if (number==PRINT_NUMBER) {
1546 printf("multiplier %g %g %g\n",
1547 multiplier,smallest,largest);
1548 }
1549 #endif
1550 #elif TRY7==2
1551 // Look at ratio
1552 double scaleFactor=fabs(rhsBeforeRelax);
1553 for (int i=0;i<number;i++) {
1554 double value=packed[i];
1555 double ratio = fabs(rhs/value);
1556 double nearest=floor(ratio+0.5);
1557 if (fabs(ratio-nearest)<1.0e-6 && nearest >= 1.0) {
1558 int iColumn=cutIndex[i];
1559 if (intVar[iColumn]) {
1560 if (colLower[iColumn]>=0.0) {
1561 if (value>0.0) {
1562 // better if smaller
1563 if (ratio>nearest) {
1564 packed[i]=scaleFactor/nearest;
1565 }
1566 } else {
1567 // better if larger
1568 if (ratio<nearest) {
1569 packed[i]=-scaleFactor/nearest;
1570 }
1571 }
1572 //if (value!=packed[i])
1573 //printf("column %d rhs %.18g element %g ratio %.18g - new element %.18g\n",
1574 // iColumn,rhs,value,ratio,packed[i]);
1575 //assert (fabs(value-packed[i])<1.0e-4);
1576 } else {
1577 //printf("column %d rhs %g element %g ratio %.18g - bounds %g,%g\n",
1578 // iColumn,rhs,value,ratio,colLower[iColumn],colUpper[iColumn]);
1579 }
1580 }
1581 }
1582 }
1583 #endif
1584 #if PRINT_NUMBER
1585 if (number==PRINT_NUMBER) {
1586 printf("after %.18g ",rhs);
1587 for (int i=0;i<PRINT_NUMBER;i++)
1588 printf("%.18g ",packed[i]);
1589 printf("\n");
1590 }
1591 #endif
1592 }
1593 }
1594 if (number<limit||!numberNonInteger) {
1595 bounds[1]=rhs;
1596 if (number>50&&numberNonInteger)
1597 bounds[1] = rhs+tolerance6+1.0e-8*fabs(rhs); // weaken
1598 #if GOMORY_RELAX_NUMBER
1599 else if (number>GOMORY_RELAX_NUMBER&&numberNonInteger>1)
1600 bounds[1] = rhs+tolerance6+1.0e-8*fabs(rhs); // weaken
1601 #endif
1602 // if close to integer - round up
1603 double nearest=floor(bounds[1]+0.5);
1604 if (bounds[1]<nearest&&bounds[1]>nearest-1.0e-4)
1605 bounds[1]=nearest;
1606 double test = CoinMin(largestFactor*largestFactorMultiplier_,
1607 relaxation);
1608 if (number>5&&numberNonInteger&&test>1.0e-20) {
1609 #ifdef CLP_INVESTIGATE2
1610 printf("relaxing rhs by %g - largestFactor was %g, rel %g\n",
1611 CoinMin(test*fabs(rhs),tolerance9),largestFactor,relaxation);
1612 #endif
1613 //bounds[1] = CoinMax(bounds[1],
1614 // rhs+CoinMin(test*fabs(rhs),tolerance9)); // weaken
1615 bounds[1] = bounds[1]+CoinMin(test*fabs(rhs),tolerance9); // weaken
1616 }
1617 #ifdef MORE_GOMORY_CUTS
1618 if (accurate) {
1619 #else
1620 {
1621 #endif
1622 // tidy
1623 if (!cleanedCut) {
1624 double range=0.0;
1625 for (int k=0;k<number;k++) {
1626 int iColumn=cutIndex[k];
1627 double thisRange=CoinMin(colUpper[iColumn]-colLower[iColumn],1000.0);
1628 range += thisRange;
1629 }
1630 // see if close to integer
1631 bool close=fabs(bounds[1]-floor(bounds[1]+0.5))*range<1.0e-6;
1632 if (close) {
1633 for (int k=0;k<number;k++) {
1634 if(fabs(packed[k]-floor(packed[k]+0.5))*range>1.0e-6) {
1635 close=false;
1636 break;
1637 }
1638 }
1639 if (close) {
1640 #ifdef PRINT_MORE
1641 bool printIt=false;
1642 for (int k=0;k<number;k++) {
1643 if(fabs(packed[k]-floor(packed[k]+0.5))>1.0e-12) {
1644 printIt=true;
1645 break;
1646 }
1647 }
1648 if (printIt) {
1649 printf("yy %.18g >= ",bounds[1]);
1650 for (int k=0;k<number;k++) {
1651 printf("(%d,%.18g) ",cutIndex[k],packed[k]);
1652 }
1653 printf("\n");
1654 }
1655 #endif
1656 bounds[1]=floor(bounds[1]+0.5);
1657 for (int k=0;k<number;k++) {
1658 packed[k]=floor(packed[k]+0.5);
1659 }
1660 }
1661 }
1662 }
1663 OsiRowCut rc;
1664 rc.setRow(number,cutIndex,packed,false);
1665 rc.setLb(bounds[0]);
1666 rc.setUb(bounds[1]);
1667 #if 0
1668 double total=0.0;
1669 for (j=0;j<number;j++) {
1670 int jColumn =cutIndex[j];
1671 double value=packed[j];
1672 std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn]
1673 <<") ";
1674 total += value*colsol[jColumn];
1675 }
1676 std::cout<<std::endl;
1677 printf("above CUT %g <= %g <= %g\n",bounds[0],total,bounds[1]);
1678 #endif
1679 #ifdef CGL_DEBUG
1680 if (debugger) {
1681 assert(!debugger->invalidCut(rc));
1682 if(debugger->invalidCut(rc))
1683 abort();
1684 }
1685 #endif
1686 #if MORE_GOMORY_CUTS<2
1687 nTotalEls -= number;
1688 cs.insertIfNotDuplicate(rc);
1689 #else
1690 if(number<saveLimit) {
1691 nTotalEls -= number;
1692 cs.insertIfNotDuplicate(rc);
1693 } else {
1694 longCuts.insertIfNotDuplicate(rc);
1695 }
1696 #endif
1697 //printf("nTot %d kCol %d iCol %d ibasic %d\n",
1698 // nTotalEls,kColumn,iColumn,iBasic);
1699 numberAdded++;
1700 #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
1701 } else if (accurate2) {
1702 OsiRowCut rc;
1703 rc.setRow(number,cutIndex,packed,false);
1704 rc.setLb(bounds[0]);
1705 rc.setUb(bounds[1]);
1706 secondaryCuts.insertIfNotDuplicate(rc);
1707 #endif
1708 }
1709 } else {
1710 #ifdef CGL_DEBUG
1711 std::cout<<"cut has "<<number<<" entries - skipped"
1712 <<std::endl;
1713 if (!number)
1714 std::cout<<"******* Empty cut - infeasible"<<std::endl;
1715 #endif
1716 }
1717 } else {
1718 // why dropped?
1719 #ifdef CGL_DEBUG
1720 std::cout<<"first violation "<<violation<<" now "
1721 <<sum-rhs<<" why?, rhs= "<<rhs<<std::endl;
1722
1723 for (j=0;j<number;j++) {
1724 int jColumn =cutIndex[j];
1725 double value=packed[j];
1726 std::cout<<"("<<jColumn<<","<<value<<","<<colsol[jColumn]
1727 <<") ";
1728 }
1729 std::cout<<std::endl;
1730 //abort();
1731 #endif
1732 }
1733 }
1734 } else {
1735 // not basic
1736 #if CGL_DEBUG>1
1737 {
1738 // put column into array
1739 array.setVector(columnLength[iColumn],row+columnStart[iColumn],
1740 columnElements+columnStart[iColumn]);
1741 // get column in tableau
1742 #ifdef CLP_OSL
1743 if (!alternateFactorization_)
1744 #endif
1745 factorization.updateColumn ( &work, &array );
1746 #ifdef CLP_OSL
1747 else
1748 factorization2->updateColumn ( &work, &array );
1749 #endif
1750 int numberInArray=array.getNumElements();
1751 printf("non-basic %d\n",iColumn);
1752 for (int j=0;j<numberInArray;j++) {
1753 int indexValue=arrayRows[j];
1754 double value=arrayElements[indexValue];
1755 if (fabs(value)>1.0e-6) {
1756 printf("%d %g\n",indexValue,value);
1757 }
1758 }
1759 }
1760 #endif
1761 }
1762 }
1763 #ifdef CLP_OSL
1764 delete factorization2;
1765 #endif
1766
1767 delete [] rowActivity;
1768 delete [] swap;
1769 delete [] rowType;
1770 delete [] packed;
1771 delete [] rowIsBasic;
1772 delete [] columnIsBasic;
1773 #ifdef MORE_GOMORY_CUTS
1774 #if MORE_GOMORY_CUTS==1
1775 int numberInaccurate = secondaryCuts.sizeRowCuts();
1776 #ifdef CLP_INVESTIGATE2
1777 int numberOrdinary = numberAdded-numberInaccurate;
1778 if (!info.inTree&&(infoOptions&512)==0)
1779 printf("Gomory has %d ordinary and %d less accurate cuts(%d els)\n",
1780 numberOrdinary,numberInaccurate,saveTotalEls-nTotalEls);
1781 #endif
1782 #elif MORE_GOMORY_CUTS==2
1783 int numberLong = longCuts.sizeRowCuts();
1784 #ifdef CLP_INVESTIGATE2
1785 int numberOrdinary = numberAdded-numberLong;
1786 if (!info.inTree&&(infoOptions&512)==0)
1787 printf("Gomory has %d ordinary and %d long cuts(%d els)\n",
1788 numberOrdinary,numberLong,saveTotalEls-nTotalEls);
1789 #endif
1790 #elif MORE_GOMORY_CUTS==3
1791 int numberLong = longCuts.sizeRowCuts();
1792 int numberInaccurate = secondaryCuts.sizeRowCuts();
1793 #ifdef CLP_INVESTIGATE2
1794 int numberOrdinary = numberAdded-numberLong-numberInaccurate;
1795 if (!info.inTree&&(infoOptions&512)==0)
1796 printf("Gomory has %d ordinary, %d long and %d less accurate cuts(%d els)\n",
1797 numberOrdinary,numberLong,numberInaccurate,saveTotalEls-nTotalEls);
1798 #endif
1799 #endif
1800 if (doSorted&&limit<numberColumns) {
1801 // Just half
1802 nTotalEls -= saveTotalEls/2;
1803 #if MORE_GOMORY_CUTS==2||MORE_GOMORY_CUTS==3
1804 while (nTotalEls>0) {
1805 for (int i=0;i<numberLong;i++) {
1806 nTotalEls -= longCuts.rowCutPtr(i)->row().getNumElements();
1807 cs.insertIfNotDuplicate(longCuts.rowCut(i));
1808 numberAdded ++;
1809 if (nTotalEls<=0)
1810 break;
1811 }
1812 break;
1813 }
1814 #endif
1815 #if MORE_GOMORY_CUTS==1||MORE_GOMORY_CUTS==3
1816 while (nTotalEls>0) {
1817 for (int i=0;i<numberInaccurate;i++) {
1818 nTotalEls -= secondaryCuts.rowCutPtr(i)->row().getNumElements();
1819 cs.insertIfNotDuplicate(secondaryCuts.rowCut(i));
1820 numberAdded ++;
1821 if (nTotalEls<=0)
1822 break;
1823 }
1824 break;
1825 }
1826 #endif
1827 }
1828 #else
1829 #ifdef CLP_INVESTIGATE2
1830 if (!info.inTree&&(infoOptions&512)==0)
1831 printf("Gomory added %d cuts(%d els)\n",numberAdded,saveTotalEls-nTotalEls);
1832 #endif
1833 #endif
1834 return numberAdded;
1835 }
1836 // Limit stuff
1837 void CglGomory::setLimit(int limit)
1838 {
1839 if (limit>=0)
1840 limit_=limit;
1841 }
1842 int CglGomory::getLimit() const
1843 {
1844 return limit_;
1845 }
1846 // Limit stuff at root
1847 void CglGomory::setLimitAtRoot(int limit)
1848 {
1849 if (limit>=0)
1850 limitAtRoot_=limit;
1851 }
1852 int CglGomory::getLimitAtRoot() const
1853 {
1854 return limitAtRoot_;
1855 }
1856 // Return maximum length of cut in tree
1857 int
1858 CglGomory::maximumLengthOfCutInTree() const
1859 {
1860 if (limit_)
1861 return limit_;
1862 else
1863 return dynamicLimitInTree_;
1864 }
1865
1866 // Away stuff
1867 void CglGomory::setAway(double value)
1868 {
1869 if (value>0.0&&value<=0.5)
1870 away_=value;
1871 }
1872 double CglGomory::getAway() const
1873 {
1874 return away_;
1875 }
1876
1877 // Away stuff at root
1878 void CglGomory::setAwayAtRoot(double value)
1879 {
1880 if (value>0.0&&value<=0.5)
1881 awayAtRoot_=value;
1882 }
1883 double CglGomory::getAwayAtRoot() const
1884 {
1885 return awayAtRoot_;
1886 }
1887
1888 // ConditionNumberMultiplier stuff
1889 void CglGomory::setConditionNumberMultiplier(double value)
1890 {
1891 if (value>=0.0)
1892 conditionNumberMultiplier_=value;
1893 }
1894 double CglGomory::getConditionNumberMultiplier() const
1895 {
1896 return conditionNumberMultiplier_;
1897 }
1898
1899 // LargestFactorMultiplier stuff
1900 void CglGomory::setLargestFactorMultiplier(double value)
1901 {
1902 if (value>=0.0)
1903 largestFactorMultiplier_=value;
1904 }
1905 double CglGomory::getLargestFactorMultiplier() const
1906 {
1907 return largestFactorMultiplier_;
1908 }
1909
1910 //-------------------------------------------------------------------
1911 // Default Constructor
1912 //-------------------------------------------------------------------
1913 CglGomory::CglGomory ()
1914 :
1915 CglCutGenerator(),
1916 away_(0.05),
1917 awayAtRoot_(0.05),
1918 conditionNumberMultiplier_(1.0e-18),
1919 largestFactorMultiplier_(1.0e-13),
1920 originalSolver_(NULL),
1921 limit_(50),
1922 limitAtRoot_(0),
1923 dynamicLimitInTree_(-1),
1924 numberTimesStalled_(0),
1925 alternateFactorization_(0),
1926 gomoryType_(0)
1927 {
1928
1929 }
1930
1931 //-------------------------------------------------------------------
1932 // Copy constructor
1933 //-------------------------------------------------------------------
1934 CglGomory::CglGomory (const CglGomory & source) :
1935 CglCutGenerator(source),
1936 away_(source.away_),
1937 awayAtRoot_(source.awayAtRoot_),
1938 conditionNumberMultiplier_(source.conditionNumberMultiplier_),
1939 largestFactorMultiplier_(source.largestFactorMultiplier_),
1940 originalSolver_(NULL),
1941 limit_(source.limit_),
1942 limitAtRoot_(source.limitAtRoot_),
1943 dynamicLimitInTree_(source.dynamicLimitInTree_),
1944 numberTimesStalled_(source.numberTimesStalled_),
1945 alternateFactorization_(source.alternateFactorization_),
1946 gomoryType_(source.gomoryType_)
1947 {
1948 if (source.originalSolver_)
1949 originalSolver_ = source.originalSolver_->clone();
1950 }
1951
1952 //-------------------------------------------------------------------
1953 // Clone
1954 //-------------------------------------------------------------------
1955 CglCutGenerator *
1956 CglGomory::clone() const
1957 {
1958 return new CglGomory(*this);
1959 }
1960
1961 //-------------------------------------------------------------------
1962 // Destructor
1963 //-------------------------------------------------------------------
1964 CglGomory::~CglGomory ()
1965 {
1966 delete originalSolver_;
1967 }
1968
1969 //----------------------------------------------------------------
1970 // Assignment operator
1971 //-------------------------------------------------------------------
1972 CglGomory &
1973 CglGomory::operator=(const CglGomory& rhs)
1974 {
1975 if (this != &rhs) {
1976 CglCutGenerator::operator=(rhs);
1977 away_=rhs.away_;
1978 awayAtRoot_=rhs.awayAtRoot_;
1979 conditionNumberMultiplier_ = rhs.conditionNumberMultiplier_;
1980 largestFactorMultiplier_ = rhs.largestFactorMultiplier_;
1981 limit_=rhs.limit_;
1982 limitAtRoot_=rhs.limitAtRoot_;
1983 dynamicLimitInTree_ = rhs.dynamicLimitInTree_;
1984 numberTimesStalled_ = rhs.numberTimesStalled_;
1985 alternateFactorization_=rhs.alternateFactorization_;
1986 gomoryType_ = rhs.gomoryType_;
1987 delete originalSolver_;
1988 if (rhs.originalSolver_)
1989 originalSolver_ = rhs.originalSolver_->clone();
1990 else
1991 originalSolver_=NULL;
1992 }
1993 return *this;
1994 }
1995 // This can be used to refresh any information
1996 void
1997 CglGomory::refreshSolver(OsiSolverInterface * solver)
1998 {
1999 int numberColumns=solver->getNumCols();
2000 const double * colUpper = solver->getColUpper();
2001 const double * colLower = solver->getColLower();
2002 canDoGlobalCuts_ = true;
2003 if (originalSolver_) {
2004 delete originalSolver_;
2005 originalSolver_ = solver->clone();
2006 }
2007 for (int i=0;i<numberColumns;i++) {
2008 if (solver->isInteger(i)) {
2009 if (colUpper[i]>colLower[i]+1.0) {
2010 canDoGlobalCuts_ = false;
2011 break;
2012 }
2013 }
2014 }
2015 }
2016 // Pass in a copy of original solver (clone it)
2017 void
2018 CglGomory::passInOriginalSolver(OsiSolverInterface * solver)
2019 {
2020 delete originalSolver_;
2021 if (solver) {
2022 if (!gomoryType_)
2023 gomoryType_=1;
2024 originalSolver_ = solver->clone();
2025 } else {
2026 gomoryType_=0;
2027 originalSolver_=NULL;
2028 }
2029 }
2030 // Does actual work - returns number of cuts
2031 int
2032 CglGomory::generateCuts( const OsiRowCutDebugger * debugger,
2033 OsiCuts & cs,
2034 const CoinPackedMatrix & columnCopy,
2035 const double * colsol,
2036 const double * colLower, const double * colUpper,
2037 const double * rowLower, const double * rowUpper,
2038 const char * intVar,
2039 const CoinWarmStartBasis* warm,
2040 const CglTreeInfo info)
2041 {
2042 CoinPackedMatrix rowCopy;
2043 rowCopy.reverseOrderedCopyOf(columnCopy);
2044 return generateCuts( debugger, cs, columnCopy, rowCopy,
2045 colsol, colLower, colUpper,
2046 rowLower, rowUpper, intVar, warm, info);
2047 }
2048 // Create C++ lines to get to current state
2049 std::string
2050 CglGomory::generateCpp( FILE * fp)
2051 {
2052 CglGomory other;
2053 fprintf(fp,"0#include \"CglGomory.hpp\"\n");
2054 fprintf(fp,"3 CglGomory gomory;\n");
2055 if (limit_!=other.limit_)
2056 fprintf(fp,"3 gomory.setLimit(%d);\n",limit_);
2057 else
2058 fprintf(fp,"4 gomory.setLimit(%d);\n",limit_);
2059 if (limitAtRoot_!=other.limitAtRoot_)
2060 fprintf(fp,"3 gomory.setLimitAtRoot(%d);\n",limitAtRoot_);
2061 else
2062 fprintf(fp,"4 gomory.setLimitAtRoot(%d);\n",limitAtRoot_);
2063 if (away_!=other.away_)
2064 fprintf(fp,"3 gomory.setAway(%g);\n",away_);
2065 else
2066 fprintf(fp,"4 gomory.setAway(%g);\n",away_);
2067 if (awayAtRoot_!=other.awayAtRoot_)
2068 fprintf(fp,"3 gomory.setAwayAtRoot(%g);\n",awayAtRoot_);
2069 else
2070 fprintf(fp,"4 gomory.setAwayAtRoot(%g);\n",awayAtRoot_);
2071 if (getAggressiveness()!=other.getAggressiveness())
2072 fprintf(fp,"3 gomory.setAggressiveness(%d);\n",getAggressiveness());
2073 else
2074 fprintf(fp,"4 gomory.setAggressiveness(%d);\n",getAggressiveness());
2075 return "gomory";
2076 }
2077