1 /* $Id: AbcSimplex.cpp 2431 2019-03-15 15:56:51Z stefan $ */
2 // Copyright (C) 2000, International Business Machines
3 // Corporation and others, Copyright (C) 2012, FasterCoin. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6 #include "ClpConfig.h"
7
8 #include "CoinPragma.hpp"
9 #include <math.h>
10 //#define ABC_DEBUG 2
11
12 #if SLIM_CLP == 2
13 #define SLIM_NOIO
14 #endif
15 #include "CoinHelperFunctions.hpp"
16 #include "CoinFloatEqual.hpp"
17 #include "ClpSimplex.hpp"
18 #include "AbcSimplex.hpp"
19 #include "AbcSimplexDual.hpp"
20 #include "AbcSimplexFactorization.hpp"
21 #include "AbcNonLinearCost.hpp"
22 #include "CoinAbcCommon.hpp"
23 #include "AbcMatrix.hpp"
24 #include "CoinIndexedVector.hpp"
25 #include "AbcDualRowDantzig.hpp"
26 #include "AbcDualRowSteepest.hpp"
27 #include "AbcPrimalColumnDantzig.hpp"
28 #include "AbcPrimalColumnSteepest.hpp"
29 #include "ClpMessage.hpp"
30 #include "ClpEventHandler.hpp"
31 #include "ClpLinearObjective.hpp"
32 #include "CoinAbcHelperFunctions.hpp"
33 #include "CoinModel.hpp"
34 #include "CoinLpIO.hpp"
35 #include <cfloat>
36
37 #include <string>
38 #include <stdio.h>
39 #include <iostream>
40 //#############################################################################
AbcSimplex(bool emptyMessages)41 AbcSimplex::AbcSimplex(bool emptyMessages)
42 :
43
44 ClpSimplex(emptyMessages)
45 {
46 gutsOfDelete(0);
47 gutsOfInitialize(0, 0, true);
48 assert(maximumAbcNumberRows_ >= 0);
49 //printf("XX %x AbcSimplex constructor\n",this);
50 }
51
52 //-----------------------------------------------------------------------------
53
~AbcSimplex()54 AbcSimplex::~AbcSimplex()
55 {
56 //printf("XX %x AbcSimplex destructor\n",this);
57 gutsOfDelete(1);
58 }
59 // Copy constructor.
AbcSimplex(const AbcSimplex & rhs)60 AbcSimplex::AbcSimplex(const AbcSimplex &rhs)
61 : ClpSimplex(rhs)
62 {
63 //printf("XX %x AbcSimplex constructor from %x\n",this,&rhs);
64 gutsOfDelete(0);
65 gutsOfInitialize(numberRows_, numberColumns_, false);
66 gutsOfCopy(rhs);
67 assert(maximumAbcNumberRows_ >= 0);
68 }
69 #include "ClpDualRowSteepest.hpp"
70 #include "ClpPrimalColumnSteepest.hpp"
71 #include "ClpFactorization.hpp"
72 // Copy constructor from model
AbcSimplex(const ClpSimplex & rhs)73 AbcSimplex::AbcSimplex(const ClpSimplex &rhs)
74 : ClpSimplex(rhs)
75 {
76 //printf("XX %x AbcSimplex constructor from ClpSimplex\n",this);
77 gutsOfDelete(0);
78 gutsOfInitialize(numberRows_, numberColumns_, true);
79 gutsOfResize(numberRows_, numberColumns_);
80 // Set up row/column selection and factorization type
81 ClpDualRowSteepest *pivot = dynamic_cast< ClpDualRowSteepest * >(rhs.dualRowPivot());
82 if (pivot) {
83 AbcDualRowSteepest steep(pivot->mode());
84 setDualRowPivotAlgorithm(steep);
85 } else {
86 AbcDualRowDantzig dantzig;
87 setDualRowPivotAlgorithm(dantzig);
88 }
89 ClpPrimalColumnSteepest *pivotColumn = dynamic_cast< ClpPrimalColumnSteepest * >(rhs.primalColumnPivot());
90 if (pivotColumn) {
91 AbcPrimalColumnSteepest steep(pivotColumn->mode());
92 setPrimalColumnPivotAlgorithm(steep);
93 } else {
94 AbcPrimalColumnDantzig dantzig;
95 setPrimalColumnPivotAlgorithm(dantzig);
96 }
97 //if (rhs.factorization()->)
98 //factorization_->forceOtherFactorization();
99 factorization()->synchronize(rhs.factorization(), this);
100 //factorization_->setGoDenseThreshold(rhs.factorization()->goDenseThreshold());
101 //factorization_->setGoSmallThreshold(rhs.factorization()->goSmallThreshold());
102 translate(DO_SCALE_AND_MATRIX | DO_BASIS_AND_ORDER | DO_STATUS | DO_SOLUTION);
103 assert(maximumAbcNumberRows_ >= 0);
104 }
105 // Assignment operator. This copies the data
106 AbcSimplex &
operator =(const AbcSimplex & rhs)107 AbcSimplex::operator=(const AbcSimplex &rhs)
108 {
109 if (this != &rhs) {
110 gutsOfDelete(1);
111 ClpSimplex::operator=(rhs);
112 gutsOfCopy(rhs);
113 assert(maximumAbcNumberRows_ >= 0);
114 }
115 return *this;
116 }
117 // fills in perturbationSaved_ from start with 0.5+random
fillPerturbation(int start,int number)118 void AbcSimplex::fillPerturbation(int start, int number)
119 {
120 double *array = perturbationSaved_ + start;
121 for (int i = 0; i < number; i++)
122 array[i] = 0.5 + 0.5 * randomNumberGenerator_.randomDouble();
123 }
124 // Sets up all extra pointers
setupPointers(int maxRows,int maxColumns)125 void AbcSimplex::setupPointers(int maxRows, int maxColumns)
126 {
127 int numberTotal = maxRows + maxColumns;
128 scaleToExternal_ = scaleFromExternal_ + numberTotal;
129 tempArray_ = offset_ + numberTotal;
130 lowerSaved_ = abcLower_ + numberTotal;
131 upperSaved_ = abcUpper_ + numberTotal;
132 costSaved_ = abcCost_ + numberTotal;
133 solutionSaved_ = abcSolution_ + numberTotal;
134 djSaved_ = abcDj_ + numberTotal;
135 internalStatusSaved_ = internalStatus_ + numberTotal;
136 perturbationSaved_ = abcPerturbation_ + numberTotal;
137 offsetRhs_ = tempArray_ + numberTotal;
138 lowerBasic_ = lowerSaved_ + numberTotal;
139 upperBasic_ = upperSaved_ + numberTotal;
140 costBasic_ = costSaved_ + 2 * numberTotal;
141 solutionBasic_ = solutionSaved_ + numberTotal;
142 djBasic_ = djSaved_ + numberTotal;
143 perturbationBasic_ = perturbationSaved_ + numberTotal;
144 columnUseScale_ = scaleToExternal_ + maxRows;
145 inverseColumnUseScale_ = scaleFromExternal_ + maxRows;
146 }
147 // Copies all saved versions to working versions and may do something for perturbation
copyFromSaved(int which)148 void AbcSimplex::copyFromSaved(int which)
149 {
150 if ((which & 1) != 0)
151 CoinAbcMemcpy(abcSolution_, solutionSaved_, maximumNumberTotal_);
152 if ((which & 2) != 0)
153 CoinAbcMemcpy(abcCost_, costSaved_, maximumNumberTotal_);
154 if ((which & 4) != 0)
155 CoinAbcMemcpy(abcLower_, lowerSaved_, maximumNumberTotal_);
156 if ((which & 8) != 0)
157 CoinAbcMemcpy(abcUpper_, upperSaved_, maximumNumberTotal_);
158 if ((which & 16) != 0)
159 CoinAbcMemcpy(abcDj_, djSaved_, maximumNumberTotal_);
160 if ((which & 32) != 0) {
161 CoinAbcMemcpy(abcLower_, abcPerturbation_, numberTotal_);
162 CoinAbcMemcpy(abcUpper_, abcPerturbation_ + numberTotal_, numberTotal_);
163 }
164 }
gutsOfCopy(const AbcSimplex & rhs)165 void AbcSimplex::gutsOfCopy(const AbcSimplex &rhs)
166 {
167 #ifdef ABC_INHERIT
168 abcSimplex_ = this;
169 #endif
170 assert(numberRows_ == rhs.numberRows_);
171 assert(numberColumns_ == rhs.numberColumns_);
172 numberTotal_ = rhs.numberTotal_;
173 maximumNumberTotal_ = rhs.maximumNumberTotal_;
174 // special options here to be safe
175 specialOptions_ = rhs.specialOptions_;
176 //assert (maximumInternalRows_ >= numberRows_);
177 //assert (maximumInternalColumns_ >= numberColumns_);
178 // Copy all scalars
179 CoinAbcMemcpy(reinterpret_cast< int * >(&sumNonBasicCosts_),
180 reinterpret_cast< const int * >(&rhs.sumNonBasicCosts_),
181 (&swappedAlgorithm_ - reinterpret_cast< int * >(&sumNonBasicCosts_)) + 1);
182 // could add 2 so can go off end
183 int sizeArray = 2 * maximumNumberTotal_ + maximumAbcNumberRows_;
184 internalStatus_ = ClpCopyOfArray(rhs.internalStatus_, sizeArray + maximumNumberTotal_);
185 abcLower_ = ClpCopyOfArray(rhs.abcLower_, sizeArray);
186 abcUpper_ = ClpCopyOfArray(rhs.abcUpper_, sizeArray);
187 abcCost_ = ClpCopyOfArray(rhs.abcCost_, sizeArray + maximumNumberTotal_);
188 abcDj_ = ClpCopyOfArray(rhs.abcDj_, sizeArray);
189
190 abcSolution_ = ClpCopyOfArray(rhs.abcSolution_, sizeArray + maximumNumberTotal_);
191 abcPerturbation_ = ClpCopyOfArray(rhs.abcPerturbation_, sizeArray);
192 abcPivotVariable_ = ClpCopyOfArray(rhs.abcPivotVariable_, maximumAbcNumberRows_);
193 //fromExternal_ = ClpCopyOfArray(rhs.fromExternal_,sizeArray);
194 //toExternal_ = ClpCopyOfArray(rhs.toExternal_,sizeArray);
195 scaleFromExternal_ = ClpCopyOfArray(rhs.scaleFromExternal_, sizeArray);
196 offset_ = ClpCopyOfArray(rhs.offset_, sizeArray);
197 setupPointers(maximumAbcNumberRows_, maximumAbcNumberColumns_);
198 if (rhs.abcMatrix_)
199 abcMatrix_ = new AbcMatrix(*rhs.abcMatrix_);
200 else
201 abcMatrix_ = NULL;
202 for (int i = 0; i < ABC_NUMBER_USEFUL; i++) {
203 usefulArray_[i] = rhs.usefulArray_[i];
204 }
205 if (rhs.abcFactorization_) {
206 setFactorization(*rhs.abcFactorization_);
207 } else {
208 delete abcFactorization_;
209 abcFactorization_ = NULL;
210 }
211 #ifdef EARLY_FACTORIZE
212 delete abcEarlyFactorization_;
213 if (rhs.abcEarlyFactorization_) {
214 abcEarlyFactorization_ = new AbcSimplexFactorization(*rhs.abcEarlyFactorization_);
215 } else {
216 abcEarlyFactorization_ = NULL;
217 }
218 #endif
219 abcDualRowPivot_ = rhs.abcDualRowPivot_->clone(true);
220 abcDualRowPivot_->setModel(this);
221 abcPrimalColumnPivot_ = rhs.abcPrimalColumnPivot_->clone(true);
222 abcPrimalColumnPivot_->setModel(this);
223 if (rhs.abcBaseModel_) {
224 abcBaseModel_ = new AbcSimplex(*rhs.abcBaseModel_);
225 } else {
226 abcBaseModel_ = NULL;
227 }
228 if (rhs.clpModel_) {
229 clpModel_ = new ClpSimplex(*rhs.clpModel_);
230 } else {
231 clpModel_ = NULL;
232 }
233 abcProgress_ = rhs.abcProgress_;
234 solveType_ = rhs.solveType_;
235 }
236 // type == 0 nullify, 1 also delete
gutsOfDelete(int type)237 void AbcSimplex::gutsOfDelete(int type)
238 {
239 if (type) {
240 delete[] abcLower_;
241 delete[] abcUpper_;
242 delete[] abcCost_;
243 delete[] abcDj_;
244 delete[] abcSolution_;
245 //delete [] fromExternal_ ;
246 //delete [] toExternal_ ;
247 delete[] scaleFromExternal_;
248 //delete [] scaleToExternal_ ;
249 delete[] offset_;
250 delete[] internalStatus_;
251 delete[] abcPerturbation_;
252 delete abcMatrix_;
253 delete abcFactorization_;
254 #ifdef EARLY_FACTORIZE
255 delete abcEarlyFactorization_;
256 #endif
257 delete[] abcPivotVariable_;
258 delete abcDualRowPivot_;
259 delete abcPrimalColumnPivot_;
260 delete abcBaseModel_;
261 delete clpModel_;
262 delete abcNonLinearCost_;
263 }
264 CoinAbcMemset0(reinterpret_cast< char * >(&scaleToExternal_),
265 reinterpret_cast< char * >(&usefulArray_[0]) - reinterpret_cast< char * >(&scaleToExternal_));
266 }
267 template < class T >
newArray(T *,int size)268 T *newArray(T * /*nullArray*/, int size)
269 {
270 if (size) {
271 T *arrayNew = new T[size];
272 #ifndef NDEBUG
273 memset(arrayNew, 'A', (size) * sizeof(T));
274 #endif
275 return arrayNew;
276 } else {
277 return NULL;
278 }
279 }
280 template < class T >
resizeArray(T * array,int oldSize1,int oldSize2,int newSize1,int newSize2,int extra)281 T *resizeArray(T *array, int oldSize1, int oldSize2, int newSize1, int newSize2, int extra)
282 {
283 if ((array || !oldSize1) && (newSize1 != oldSize1 || newSize2 != oldSize2)) {
284 int newTotal = newSize1 + newSize2;
285 int oldTotal = oldSize1 + oldSize2;
286 T *arrayNew;
287 if (newSize1 > oldSize1 || newSize2 > oldSize2) {
288 arrayNew = new T[2 * newTotal + newSize1 + extra];
289 #ifndef NDEBUG
290 memset(arrayNew, 'A', (2 * newTotal + newSize1 + extra) * sizeof(T));
291 #endif
292 CoinAbcMemcpy(arrayNew, array, oldSize1);
293 CoinAbcMemcpy(arrayNew + newSize1, array + oldSize1, oldSize2);
294 // and second half
295 CoinAbcMemcpy(arrayNew + newTotal, array, oldSize1 + oldTotal);
296 CoinAbcMemcpy(arrayNew + newSize1 + newTotal, array + oldSize1 + oldTotal, oldSize2);
297 delete[] array;
298 } else {
299 arrayNew = array;
300 for (int i = 0; i < newSize2; i++)
301 array[newSize1 + i] = array[oldSize1 + i];
302 for (int i = 0; i < newSize1; i++)
303 array[newTotal + i] = array[oldTotal + i];
304 for (int i = 0; i < newSize2; i++)
305 array[newSize1 + newTotal + i] = array[oldSize1 + oldTotal + i];
306 }
307 return arrayNew;
308 } else {
309 return array;
310 }
311 }
312 // Initializes arrays
gutsOfInitialize(int numberRows,int numberColumns,bool doMore)313 void AbcSimplex::gutsOfInitialize(int numberRows, int numberColumns, bool doMore)
314 {
315 #ifdef ABC_INHERIT
316 abcSimplex_ = this;
317 #endif
318 // Zero all
319 CoinAbcMemset0(&sumNonBasicCosts_, (reinterpret_cast< double * >(&usefulArray_[0]) - &sumNonBasicCosts_));
320 zeroTolerance_ = 1.0e-13;
321 bestObjectiveValue_ = -COIN_DBL_MAX;
322 primalToleranceToGetOptimal_ = -1.0;
323 primalTolerance_ = 1.0e-6;
324 //dualTolerance_ = 1.0e-6;
325 alphaAccuracy_ = -1.0;
326 upperIn_ = -COIN_DBL_MAX;
327 lowerOut_ = -1;
328 valueOut_ = -1;
329 upperOut_ = -1;
330 dualOut_ = -1;
331 acceptablePivot_ = 1.0e-8;
332 //dualBound_=1.0e9;
333 sequenceIn_ = -1;
334 directionIn_ = -1;
335 sequenceOut_ = -1;
336 directionOut_ = -1;
337 pivotRow_ = -1;
338 lastGoodIteration_ = -100;
339 numberPrimalInfeasibilities_ = 100;
340 numberRefinements_ = 1;
341 changeMade_ = 1;
342 forceFactorization_ = -1;
343 if (perturbation_ < 50 || (perturbation_ > 60 && perturbation_ < 100))
344 perturbation_ = 50;
345 lastBadIteration_ = -999999;
346 lastFlaggedIteration_ = -999999; // doesn't seem to be used
347 firstFree_ = -1;
348 incomingInfeasibility_ = 1.0;
349 allowedInfeasibility_ = 10.0;
350 solveType_ = 1; // say simplex based life form
351 //specialOptions_|=65536;
352 //ClpModel::startPermanentArrays();
353 maximumInternalRows_ = 0;
354 maximumInternalColumns_ = 0;
355 numberRows_ = numberRows;
356 numberColumns_ = numberColumns;
357 numberTotal_ = numberRows_ + numberColumns_;
358 maximumAbcNumberRows_ = numberRows;
359 maximumAbcNumberColumns_ = numberColumns;
360 maximumNumberTotal_ = numberTotal_;
361 int sizeArray = 2 * maximumNumberTotal_ + maximumAbcNumberRows_;
362 if (doMore) {
363 // say Steepest pricing
364 abcDualRowPivot_ = new AbcDualRowSteepest();
365 abcPrimalColumnPivot_ = new AbcPrimalColumnSteepest();
366 internalStatus_ = newArray((unsigned char *)NULL,
367 sizeArray + maximumNumberTotal_);
368 abcLower_ = newArray((double *)NULL, sizeArray);
369 abcUpper_ = newArray((double *)NULL, sizeArray);
370 abcCost_ = newArray((double *)NULL, sizeArray + maximumNumberTotal_);
371 abcDj_ = newArray((double *)NULL, sizeArray);
372 abcSolution_ = newArray((double *)NULL, sizeArray + maximumNumberTotal_);
373 //fromExternal_ = newArray((int *)NULL,sizeArray);
374 //toExternal_ = newArray((int *)NULL,sizeArray);
375 scaleFromExternal_ = newArray((double *)NULL, sizeArray);
376 offset_ = newArray((double *)NULL, sizeArray);
377 abcPerturbation_ = newArray((double *)NULL, sizeArray);
378 abcPivotVariable_ = newArray((int *)NULL, maximumAbcNumberRows_);
379 // Fill perturbation array
380 setupPointers(maximumAbcNumberRows_, maximumAbcNumberColumns_);
381 fillPerturbation(0, maximumNumberTotal_);
382 }
383 // get an empty factorization so we can set tolerances etc
384 getEmptyFactorization();
385 for (int i = 0; i < ABC_NUMBER_USEFUL; i++)
386 usefulArray_[i].reserve(CoinMax(CoinMax(numberTotal_, maximumAbcNumberRows_ + 200), 2 * numberRows_));
387 //savedStatus_ = internalStatus_+maximumNumberTotal_;
388 //startPermanentArrays();
389 }
390 // resizes arrays
gutsOfResize(int numberRows,int numberColumns)391 void AbcSimplex::gutsOfResize(int numberRows, int numberColumns)
392 {
393 if (!abcSolution_) {
394 numberRows_ = 0;
395 numberColumns_ = 0;
396 numberTotal_ = 0;
397 maximumAbcNumberRows_ = 0;
398 maximumAbcNumberColumns_ = 0;
399 maximumNumberTotal_ = 0;
400 }
401 if (numberRows == numberRows_ && numberColumns == numberColumns_)
402 return;
403 // can do on state bit
404 int newSize1 = CoinMax(numberRows, maximumAbcNumberRows_);
405 if ((stateOfProblem_ & ADD_A_BIT) != 0 && numberRows > maximumAbcNumberRows_)
406 newSize1 = CoinMax(numberRows, maximumAbcNumberRows_ + CoinMin(100, numberRows_ / 10));
407 int newSize2 = CoinMax(numberColumns, maximumAbcNumberColumns_);
408 numberRows_ = numberRows;
409 numberColumns_ = numberColumns;
410 numberTotal_ = numberRows_ + numberColumns_;
411 //fromExternal_ = resizeArray(fromExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1);
412 //toExternal_ = resizeArray(toExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0);
413 scaleFromExternal_ = resizeArray(scaleFromExternal_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
414 //scaleToExternal_ = resizeArray(scaleToExternal_,maximumAbcNumberRows_,maximumAbcNumberColumns_,newSize1,newSize2,1,0);
415 internalStatus_ = resizeArray(internalStatus_, maximumAbcNumberRows_,
416 maximumAbcNumberColumns_,
417 newSize1, newSize2, numberTotal_);
418 abcLower_ = resizeArray(abcLower_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
419 abcUpper_ = resizeArray(abcUpper_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
420 abcCost_ = resizeArray(abcCost_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, numberTotal_);
421 abcDj_ = resizeArray(abcDj_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
422 abcSolution_ = resizeArray(abcSolution_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, numberTotal_);
423 abcPerturbation_ = resizeArray(abcPerturbation_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
424 offset_ = resizeArray(offset_, maximumAbcNumberRows_, maximumAbcNumberColumns_, newSize1, newSize2, 0);
425 setupPointers(newSize1, newSize2);
426 // Fill gaps in perturbation array
427 fillPerturbation(maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_);
428 fillPerturbation(newSize1 + maximumAbcNumberColumns_, newSize2 - maximumAbcNumberColumns_);
429 // Clean gap
430 //CoinIotaN(fromExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_);
431 //CoinIotaN(toExternal_+maximumAbcNumberRows_,newSize1-maximumAbcNumberRows_,maximumAbcNumberRows_);
432 CoinFillN(scaleFromExternal_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 1.0);
433 CoinFillN(scaleToExternal_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 1.0);
434 CoinFillN(internalStatusSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, static_cast< unsigned char >(1));
435 CoinFillN(lowerSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, -COIN_DBL_MAX);
436 CoinFillN(upperSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, COIN_DBL_MAX);
437 CoinFillN(costSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0);
438 CoinFillN(djSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0);
439 CoinFillN(solutionSaved_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0);
440 CoinFillN(offset_ + maximumAbcNumberRows_, newSize1 - maximumAbcNumberRows_, 0.0);
441 maximumAbcNumberRows_ = newSize1;
442 maximumAbcNumberColumns_ = newSize2;
443 maximumNumberTotal_ = newSize1 + newSize2;
444 delete[] abcPivotVariable_;
445 abcPivotVariable_ = new int[maximumAbcNumberRows_];
446 for (int i = 0; i < ABC_NUMBER_USEFUL; i++)
447 usefulArray_[i].reserve(CoinMax(numberTotal_, maximumAbcNumberRows_ + 200));
448 }
translate(int type)449 void AbcSimplex::translate(int type)
450 {
451 // clear bottom bits
452 stateOfProblem_ &= ~(VALUES_PASS - 1);
453 if ((type & DO_SCALE_AND_MATRIX) != 0) {
454 //stateOfProblem_ |= DO_SCALE_AND_MATRIX;
455 stateOfProblem_ |= DO_SCALE_AND_MATRIX;
456 delete abcMatrix_;
457 abcMatrix_ = new AbcMatrix(*matrix());
458 abcMatrix_->setModel(this);
459 abcMatrix_->scale(scalingFlag_ ? 0 : -1);
460 }
461 if ((type & DO_STATUS) != 0 && (type & DO_BASIS_AND_ORDER) == 0) {
462 // from Clp enum to Abc enum (and bound flip)
463 unsigned char lookupToAbcSlack[6] = { 4, 6, 0 /*1*/, 1 /*0*/, 5, 7 };
464 unsigned char *COIN_RESTRICT statusAbc = internalStatus_;
465 const unsigned char *COIN_RESTRICT statusClp = status_;
466 int i;
467 for (i = numberRows_ - 1; i >= 0; i--) {
468 unsigned char status = statusClp[i] & 7;
469 bool basicClp = status == 1;
470 bool basicAbc = (statusAbc[i] & 7) == 4;
471 if (basicClp == basicAbc)
472 statusAbc[i] = lookupToAbcSlack[status];
473 else
474 break;
475 }
476 if (!i) {
477 // from Clp enum to Abc enum
478 unsigned char lookupToAbc[6] = { 4, 6, 1, 0, 5, 7 };
479 statusAbc += maximumAbcNumberRows_;
480 statusClp += maximumAbcNumberRows_;
481 for (i = numberColumns_ - 1; i >= 0; i--) {
482 unsigned char status = statusClp[i] & 7;
483 bool basicClp = status == 1;
484 bool basicAbc = (statusAbc[i] & 7) == 4;
485 if (basicClp == basicAbc)
486 statusAbc[i] = lookupToAbc[status];
487 else
488 break;
489 }
490 if (i)
491 type |= DO_BASIS_AND_ORDER;
492 } else {
493 type |= DO_BASIS_AND_ORDER;
494 }
495 stateOfProblem_ |= DO_STATUS;
496 }
497 if ((type & DO_BASIS_AND_ORDER) != 0) {
498 permuteIn();
499 permuteBasis();
500 stateOfProblem_ |= DO_BASIS_AND_ORDER;
501 type &= ~DO_SOLUTION;
502 }
503 if ((type & DO_SOLUTION) != 0) {
504 permuteBasis();
505 stateOfProblem_ |= DO_SOLUTION;
506 }
507 if ((type & DO_JUST_BOUNDS) != 0) {
508 stateOfProblem_ |= DO_JUST_BOUNDS;
509 }
510 if (!type) {
511 // just move stuff down
512 CoinAbcMemcpy(abcLower_, abcLower_ + maximumNumberTotal_, numberTotal_);
513 CoinAbcMemcpy(abcUpper_, abcUpper_ + maximumNumberTotal_, numberTotal_);
514 CoinAbcMemcpy(abcCost_, abcCost_ + maximumNumberTotal_, numberTotal_);
515 }
516 }
517 #ifdef ABC_SPRINT
518 // Overwrite to create sub problem (just internal arrays) - save full stuff
519 AbcSimplex *
createSubProblem(int numberColumns,const int * whichColumn)520 AbcSimplex::createSubProblem(int numberColumns, const int *whichColumn)
521 {
522 int numberFullColumns = numberColumns_;
523 numberColumns_ = numberColumns;
524 AbcSimplex *subProblem = this;
525 AbcSimplex *fullProblem = reinterpret_cast< AbcSimplex * >(new char[sizeof(AbcSimplex)]);
526 memset(fullProblem, 0, sizeof(AbcSimplex));
527 fullProblem->maximumAbcNumberRows_ = maximumAbcNumberRows_;
528 fullProblem->numberColumns_ = numberFullColumns;
529 fullProblem->numberTotal_ = numberTotal_;
530 fullProblem->maximumNumberTotal_ = maximumNumberTotal_;
531 fullProblem->numberTotalWithoutFixed_ = numberTotalWithoutFixed_;
532 fullProblem->abcPrimalColumnPivot_ = abcPrimalColumnPivot_;
533 fullProblem->internalStatus_ = internalStatus_;
534 fullProblem->abcLower_ = abcLower_;
535 fullProblem->abcUpper_ = abcUpper_;
536 fullProblem->abcCost_ = abcCost_;
537 fullProblem->abcDj_ = abcDj_;
538 fullProblem->abcSolution_ = abcSolution_;
539 fullProblem->scaleFromExternal_ = scaleFromExternal_;
540 fullProblem->offset_ = offset_;
541 fullProblem->abcPerturbation_ = abcPerturbation_;
542 fullProblem->abcPivotVariable_ = abcPivotVariable_;
543 fullProblem->abcMatrix_ = abcMatrix_;
544 fullProblem->abcNonLinearCost_ = abcNonLinearCost_;
545 fullProblem->setupPointers(maximumAbcNumberRows_, numberFullColumns);
546 subProblem->numberTotal_ = maximumAbcNumberRows_ + numberColumns;
547 subProblem->maximumNumberTotal_ = maximumAbcNumberRows_ + numberColumns;
548 subProblem->numberTotalWithoutFixed_ = subProblem->numberTotal_;
549 int sizeArray = 2 * subProblem->maximumNumberTotal_ + maximumAbcNumberRows_;
550 subProblem->internalStatus_ = newArray((unsigned char *)NULL,
551 sizeArray + subProblem->maximumNumberTotal_);
552 subProblem->abcLower_ = newArray((double *)NULL, sizeArray);
553 subProblem->abcUpper_ = newArray((double *)NULL, sizeArray);
554 subProblem->abcCost_ = newArray((double *)NULL, sizeArray + subProblem->maximumNumberTotal_);
555 subProblem->abcDj_ = newArray((double *)NULL, sizeArray);
556 subProblem->abcSolution_ = newArray((double *)NULL, sizeArray + subProblem->maximumNumberTotal_);
557 //fromExternal_ = newArray((int *)NULL,sizeArray);
558 //toExternal_ = newArray((int *)NULL,sizeArray);
559 subProblem->scaleFromExternal_ = newArray((double *)NULL, sizeArray);
560 subProblem->offset_ = newArray((double *)NULL, sizeArray);
561 subProblem->abcPerturbation_ = newArray((double *)NULL, sizeArray);
562 subProblem->abcPivotVariable_ = newArray((int *)NULL, maximumAbcNumberRows_);
563 subProblem->setupPointers(maximumAbcNumberRows_, numberColumns);
564 // could use arrays - but for now be safe
565 int *backward = new int[numberFullColumns + numberRows_];
566 int *whichRow = backward + numberFullColumns;
567 for (int i = 0; i < numberFullColumns; i++)
568 backward[i] = -1;
569 for (int i = 0; i < numberColumns; i++)
570 backward[whichColumn[i]] = i + numberRows_;
571 for (int i = 0; i < numberRows_; i++) {
572 whichRow[i] = i;
573 int iPivot = fullProblem->abcPivotVariable_[i];
574 if (iPivot < numberRows_) {
575 subProblem->abcPivotVariable_[i] = iPivot;
576 } else {
577 assert(backward[iPivot - numberRows_] >= 0);
578 subProblem->abcPivotVariable_[i] = backward[iPivot - numberRows_];
579 }
580 subProblem->internalStatus_[i] = fullProblem->internalStatus_[i];
581 subProblem->abcLower_[i] = fullProblem->abcLower_[i];
582 subProblem->abcUpper_[i] = fullProblem->abcUpper_[i];
583 subProblem->abcCost_[i] = fullProblem->abcCost_[i];
584 subProblem->abcDj_[i] = fullProblem->abcDj_[i];
585 subProblem->abcSolution_[i] = fullProblem->abcSolution_[i];
586 subProblem->abcPerturbation_[i] = fullProblem->abcPerturbation_[i];
587 subProblem->internalStatusSaved_[i] = fullProblem->internalStatusSaved_[i];
588 subProblem->perturbationSaved_[i] = fullProblem->perturbationSaved_[i];
589 subProblem->lowerSaved_[i] = fullProblem->lowerSaved_[i];
590 subProblem->upperSaved_[i] = fullProblem->upperSaved_[i];
591 subProblem->costSaved_[i] = fullProblem->costSaved_[i];
592 subProblem->djSaved_[i] = fullProblem->djSaved_[i];
593 subProblem->solutionSaved_[i] = fullProblem->solutionSaved_[i];
594 subProblem->offset_[i] = fullProblem->offset_[i];
595 subProblem->lowerBasic_[i] = fullProblem->lowerBasic_[i];
596 subProblem->upperBasic_[i] = fullProblem->upperBasic_[i];
597 subProblem->costBasic_[i] = fullProblem->costBasic_[i];
598 subProblem->solutionBasic_[i] = fullProblem->solutionBasic_[i];
599 subProblem->djBasic_[i] = fullProblem->djBasic_[i];
600 }
601 for (int i = 0; i < numberColumns; i++) {
602 int k = whichColumn[i];
603 subProblem->internalStatus_[maximumAbcNumberRows_ + i] = fullProblem->internalStatus_[maximumAbcNumberRows_ + k];
604 subProblem->internalStatus_[maximumAbcNumberRows_ + i] = fullProblem->internalStatus_[maximumAbcNumberRows_ + k];
605 subProblem->abcLower_[maximumAbcNumberRows_ + i] = fullProblem->abcLower_[maximumAbcNumberRows_ + k];
606 subProblem->abcUpper_[maximumAbcNumberRows_ + i] = fullProblem->abcUpper_[maximumAbcNumberRows_ + k];
607 subProblem->abcCost_[maximumAbcNumberRows_ + i] = fullProblem->abcCost_[maximumAbcNumberRows_ + k];
608 subProblem->abcDj_[maximumAbcNumberRows_ + i] = fullProblem->abcDj_[maximumAbcNumberRows_ + k];
609 subProblem->abcSolution_[maximumAbcNumberRows_ + i] = fullProblem->abcSolution_[maximumAbcNumberRows_ + k];
610 subProblem->abcPerturbation_[maximumAbcNumberRows_ + i] = fullProblem->abcPerturbation_[maximumAbcNumberRows_ + k];
611 subProblem->internalStatusSaved_[maximumAbcNumberRows_ + i] = fullProblem->internalStatusSaved_[maximumAbcNumberRows_ + k];
612 subProblem->perturbationSaved_[maximumAbcNumberRows_ + i] = fullProblem->perturbationSaved_[maximumAbcNumberRows_ + k];
613 subProblem->lowerSaved_[maximumAbcNumberRows_ + i] = fullProblem->lowerSaved_[maximumAbcNumberRows_ + k];
614 subProblem->upperSaved_[maximumAbcNumberRows_ + i] = fullProblem->upperSaved_[maximumAbcNumberRows_ + k];
615 subProblem->costSaved_[maximumAbcNumberRows_ + i] = fullProblem->costSaved_[maximumAbcNumberRows_ + k];
616 subProblem->djSaved_[maximumAbcNumberRows_ + i] = fullProblem->djSaved_[maximumAbcNumberRows_ + k];
617 subProblem->solutionSaved_[maximumAbcNumberRows_ + i] = fullProblem->solutionSaved_[maximumAbcNumberRows_ + k];
618 subProblem->offset_[maximumAbcNumberRows_ + i] = fullProblem->offset_[maximumAbcNumberRows_ + k];
619 }
620 subProblem->abcNonLinearCost_ = new AbcNonLinearCost(subProblem);
621 subProblem->abcNonLinearCost_->checkInfeasibilities(0.0);
622 subProblem->abcMatrix_ = new AbcMatrix(*fullProblem->abcMatrix_, numberRows_, whichRow,
623 numberColumns, whichColumn);
624 subProblem->abcMatrix_->setModel(subProblem);
625 subProblem->abcMatrix_->rebalance();
626 subProblem->abcPrimalColumnPivot_ = new AbcPrimalColumnSteepest();
627 subProblem->abcPrimalColumnPivot_->saveWeights(subProblem, 2);
628 delete[] backward;
629 // swap
630 return fullProblem;
631 }
632 // Restore stuff from sub problem (and delete sub problem)
restoreFromSubProblem(AbcSimplex * fullProblem,const int * whichColumn)633 void AbcSimplex::restoreFromSubProblem(AbcSimplex *fullProblem, const int *whichColumn)
634 {
635 AbcSimplex *subProblem = this;
636 for (int i = 0; i < numberRows_; i++) {
637 int iPivot = subProblem->abcPivotVariable_[i];
638 if (iPivot < numberRows_) {
639 fullProblem->abcPivotVariable_[i] = iPivot;
640 } else {
641 fullProblem->abcPivotVariable_[i] = whichColumn[iPivot - numberRows_] + numberRows_;
642 }
643 fullProblem->internalStatus_[i] = subProblem->internalStatus_[i];
644 fullProblem->abcLower_[i] = subProblem->abcLower_[i];
645 fullProblem->abcUpper_[i] = subProblem->abcUpper_[i];
646 fullProblem->abcCost_[i] = subProblem->abcCost_[i];
647 fullProblem->abcDj_[i] = subProblem->abcDj_[i];
648 fullProblem->abcSolution_[i] = subProblem->abcSolution_[i];
649 fullProblem->abcPerturbation_[i] = subProblem->abcPerturbation_[i];
650 fullProblem->internalStatusSaved_[i] = subProblem->internalStatusSaved_[i];
651 fullProblem->perturbationSaved_[i] = subProblem->perturbationSaved_[i];
652 fullProblem->lowerSaved_[i] = subProblem->lowerSaved_[i];
653 fullProblem->upperSaved_[i] = subProblem->upperSaved_[i];
654 fullProblem->costSaved_[i] = subProblem->costSaved_[i];
655 fullProblem->djSaved_[i] = subProblem->djSaved_[i];
656 fullProblem->solutionSaved_[i] = subProblem->solutionSaved_[i];
657 fullProblem->offset_[i] = subProblem->offset_[i];
658 fullProblem->lowerBasic_[i] = subProblem->lowerBasic_[i];
659 fullProblem->upperBasic_[i] = subProblem->upperBasic_[i];
660 fullProblem->costBasic_[i] = subProblem->costBasic_[i];
661 fullProblem->solutionBasic_[i] = subProblem->solutionBasic_[i];
662 fullProblem->djBasic_[i] = subProblem->djBasic_[i];
663 }
664 int numberColumns = subProblem->numberColumns_;
665 for (int i = 0; i < numberColumns; i++) {
666 int k = whichColumn[i];
667 fullProblem->internalStatus_[maximumAbcNumberRows_ + k] = subProblem->internalStatus_[maximumAbcNumberRows_ + i];
668 fullProblem->abcLower_[maximumAbcNumberRows_ + k] = subProblem->abcLower_[maximumAbcNumberRows_ + i];
669 fullProblem->abcUpper_[maximumAbcNumberRows_ + k] = subProblem->abcUpper_[maximumAbcNumberRows_ + i];
670 fullProblem->abcCost_[maximumAbcNumberRows_ + k] = subProblem->abcCost_[maximumAbcNumberRows_ + i];
671 fullProblem->abcDj_[maximumAbcNumberRows_ + k] = subProblem->abcDj_[maximumAbcNumberRows_ + i];
672 fullProblem->abcSolution_[maximumAbcNumberRows_ + k] = subProblem->abcSolution_[maximumAbcNumberRows_ + i];
673 fullProblem->abcPerturbation_[maximumAbcNumberRows_ + k] = subProblem->abcPerturbation_[maximumAbcNumberRows_ + i];
674 fullProblem->internalStatusSaved_[maximumAbcNumberRows_ + k] = subProblem->internalStatusSaved_[maximumAbcNumberRows_ + i];
675 fullProblem->perturbationSaved_[maximumAbcNumberRows_ + k] = subProblem->perturbationSaved_[maximumAbcNumberRows_ + i];
676 fullProblem->lowerSaved_[maximumAbcNumberRows_ + k] = subProblem->lowerSaved_[maximumAbcNumberRows_ + i];
677 fullProblem->upperSaved_[maximumAbcNumberRows_ + k] = subProblem->upperSaved_[maximumAbcNumberRows_ + i];
678 fullProblem->costSaved_[maximumAbcNumberRows_ + k] = subProblem->costSaved_[maximumAbcNumberRows_ + i];
679 fullProblem->djSaved_[maximumAbcNumberRows_ + k] = subProblem->djSaved_[maximumAbcNumberRows_ + i];
680 fullProblem->solutionSaved_[maximumAbcNumberRows_ + k] = subProblem->solutionSaved_[maximumAbcNumberRows_ + i];
681 fullProblem->offset_[maximumAbcNumberRows_ + k] = subProblem->offset_[maximumAbcNumberRows_ + i];
682 }
683 delete[] subProblem->internalStatus_;
684 delete[] subProblem->abcPerturbation_;
685 delete subProblem->abcMatrix_;
686 delete[] subProblem->abcLower_;
687 delete[] subProblem->abcUpper_;
688 delete[] subProblem->abcCost_;
689 delete[] subProblem->abcSolution_;
690 delete[] subProblem->abcDj_;
691 delete subProblem->abcPrimalColumnPivot_;
692 delete[] subProblem->scaleFromExternal_;
693 delete[] subProblem->offset_;
694 delete[] subProblem->abcPivotVariable_;
695 delete[] subProblem->reversePivotVariable_;
696 delete subProblem->abcNonLinearCost_;
697 numberColumns_ = fullProblem->numberColumns_;
698 numberTotal_ = fullProblem->numberTotal_;
699 maximumNumberTotal_ = fullProblem->maximumNumberTotal_;
700 numberTotalWithoutFixed_ = fullProblem->numberTotalWithoutFixed_;
701 abcPrimalColumnPivot_ = fullProblem->abcPrimalColumnPivot_;
702 internalStatus_ = fullProblem->internalStatus_;
703 abcLower_ = fullProblem->abcLower_;
704 abcUpper_ = fullProblem->abcUpper_;
705 abcCost_ = fullProblem->abcCost_;
706 abcDj_ = fullProblem->abcDj_;
707 abcSolution_ = fullProblem->abcSolution_;
708 scaleFromExternal_ = fullProblem->scaleFromExternal_;
709 offset_ = fullProblem->offset_;
710 abcPerturbation_ = fullProblem->abcPerturbation_;
711 abcPivotVariable_ = fullProblem->abcPivotVariable_;
712 abcMatrix_ = fullProblem->abcMatrix_;
713 setupPointers(maximumAbcNumberRows_, numberColumns);
714 // ? redo nonlinearcost
715 abcNonLinearCost_ = fullProblem->abcNonLinearCost_;
716 abcNonLinearCost_->refresh();
717 delete[] reinterpret_cast< char * >(fullProblem);
718 }
719 #endif
720 /* Sets dual values pass djs using unscaled duals
721 type 1 - values pass
722 type 2 - just use as infeasibility weights
723 type 3 - as 2 but crash
724 */
setupDualValuesPass(const double * fakeDuals,const double * fakePrimals,int type)725 void AbcSimplex::setupDualValuesPass(const double *fakeDuals,
726 const double *fakePrimals,
727 int type)
728 {
729 // allslack basis
730 memset(internalStatus_, 6, numberRows_);
731 // temp
732 if (type == 1) {
733 bool allEqual = true;
734 for (int i = 0; i < numberRows_; i++) {
735 if (rowUpper_[i] > rowLower_[i]) {
736 allEqual = false;
737 break;
738 }
739 }
740 if (allEqual) {
741 // just modify costs
742 transposeTimes(-1.0, fakeDuals, objective());
743 return;
744 }
745 }
746 // compute unscaled djs
747 CoinAbcMemcpy(djSaved_ + maximumAbcNumberRows_, objective(), numberColumns_);
748 matrix_->transposeTimes(-1.0, fakeDuals, djSaved_ + maximumAbcNumberRows_);
749 // save fake solution
750 assert(solution_);
751 //solution_ = new double[numberTotal_];
752 CoinAbcMemset0(solution_, numberRows_);
753 CoinAbcMemcpy(solution_ + maximumAbcNumberRows_, fakePrimals, numberColumns_);
754 // adjust
755 for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++)
756 solution_[iSequence] -= offset_[iSequence];
757 matrix_->times(-1.0, solution_ + maximumAbcNumberRows_, solution_);
758 double direction = optimizationDirection_;
759 const double *COIN_RESTRICT rowScale = scaleFromExternal_;
760 const double *COIN_RESTRICT inverseRowScale = scaleToExternal_;
761 for (int iRow = 0; iRow < numberRows_; iRow++) {
762 djSaved_[iRow] = direction * fakeDuals[iRow] * inverseRowScale[iRow];
763 solution_[iRow] *= rowScale[iRow];
764 }
765 const double *COIN_RESTRICT columnScale = scaleToExternal_;
766 for (int iColumn = maximumAbcNumberRows_; iColumn < numberTotal_; iColumn++)
767 djSaved_[iColumn] *= direction * columnScale[iColumn];
768 // Set solution values
769 const double *lower = abcLower_ + maximumAbcNumberRows_;
770 const double *upper = abcUpper_ + maximumAbcNumberRows_;
771 double *solution = abcSolution_ + maximumAbcNumberRows_;
772 const double *djs = djSaved_ + maximumAbcNumberRows_;
773 const double *inverseColumnScale = scaleFromExternal_ + maximumAbcNumberRows_;
774 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
775 double thisValue = fakePrimals[iColumn];
776 double thisLower = columnLower_[iColumn];
777 double thisUpper = columnUpper_[iColumn];
778 double thisDj = djs[iColumn];
779 solution_[iColumn + maximumAbcNumberRows_] = solution_[iColumn + maximumAbcNumberRows_] * inverseColumnScale[iColumn];
780 if (thisLower > -1.0e30) {
781 if (thisUpper < 1.0e30) {
782 double gapUp = thisUpper - thisValue;
783 double gapDown = thisValue - thisLower;
784 bool goUp;
785 if (gapUp > gapDown && thisDj > -dualTolerance_) {
786 goUp = false;
787 } else if (gapUp < gapDown && thisDj < dualTolerance_) {
788 goUp = true;
789 } else {
790 // wild guess
791 double badUp;
792 double badDown;
793 if (gapUp > gapDown) {
794 badUp = gapUp * dualTolerance_;
795 badDown = -gapDown * thisDj;
796 } else {
797 badUp = gapUp * thisDj;
798 badDown = gapDown * dualTolerance_;
799 }
800 goUp = badDown > badUp;
801 }
802 if (goUp) {
803 solution[iColumn] = upper[iColumn];
804 setInternalStatus(iColumn + maximumAbcNumberRows_, atUpperBound);
805 setStatus(iColumn, ClpSimplex::atUpperBound);
806 } else {
807 solution[iColumn] = lower[iColumn];
808 setInternalStatus(iColumn + maximumAbcNumberRows_, atLowerBound);
809 setStatus(iColumn, ClpSimplex::atLowerBound);
810 }
811 } else {
812 solution[iColumn] = lower[iColumn];
813 setInternalStatus(iColumn + maximumAbcNumberRows_, atLowerBound);
814 setStatus(iColumn, ClpSimplex::atLowerBound);
815 }
816 } else if (thisUpper < 1.0e30) {
817 solution[iColumn] = upper[iColumn];
818 setInternalStatus(iColumn + maximumAbcNumberRows_, atUpperBound);
819 setStatus(iColumn, ClpSimplex::atUpperBound);
820 } else {
821 // free
822 solution[iColumn] = thisValue * inverseColumnScale[iColumn];
823 setInternalStatus(iColumn + maximumAbcNumberRows_, isFree);
824 setStatus(iColumn, ClpSimplex::isFree);
825 }
826 }
827 switch (type) {
828 case 1:
829 stateOfProblem_ |= VALUES_PASS;
830 break;
831 case 2:
832 stateOfProblem_ |= VALUES_PASS2;
833 delete[] solution_;
834 solution_ = NULL;
835 break;
836 case 3:
837 stateOfProblem_ |= VALUES_PASS2;
838 break;
839 }
840 }
841 //#############################################################################
computePrimals(CoinIndexedVector * arrayVector,CoinIndexedVector * previousVector)842 int AbcSimplex::computePrimals(CoinIndexedVector *arrayVector, CoinIndexedVector *previousVector)
843 {
844
845 arrayVector->clear();
846 previousVector->clear();
847 // accumulate non basic stuff
848 double *COIN_RESTRICT array = arrayVector->denseVector();
849 CoinAbcScatterZeroTo(abcSolution_, abcPivotVariable_, numberRows_);
850 abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, array);
851 #if 0
852 static int xxxxxx=0;
853 if (xxxxxx==0)
854 for (int i=0;i<numberRows_;i++)
855 printf("%d %.18g\n",i,array[i]);
856 #endif
857 //arrayVector->scan(0,numberRows_,zeroTolerance_);
858 // Ftran adjusted RHS and iterate to improve accuracy
859 double lastError = COIN_DBL_MAX;
860 CoinIndexedVector *thisVector = arrayVector;
861 CoinIndexedVector *lastVector = previousVector;
862 //if (arrayVector->getNumElements())
863 #if 0
864 double largest=0.0;
865 int iLargest=-1;
866 for (int i=0;i<numberRows_;i++) {
867 if (fabs(array[i])>largest) {
868 largest=array[i];
869 iLargest=i;
870 }
871 }
872 printf("largest %g at row %d\n",array[iLargest],iLargest);
873 #endif
874 abcFactorization_->updateFullColumn(*thisVector);
875 #if 0
876 largest=0.0;
877 iLargest=-1;
878 for (int i=0;i<numberRows_;i++) {
879 if (fabs(array[i])>largest) {
880 largest=array[i];
881 iLargest=i;
882 }
883 }
884 printf("after largest %g at row %d\n",array[iLargest],iLargest);
885 #endif
886 #if 0
887 if (xxxxxx==0)
888 for (int i=0;i<numberRows_;i++)
889 printf("xx %d %.19g\n",i,array[i]);
890 if (numberIterations_>300)
891 exit(0);
892 #endif
893 int numberRefinements = 0;
894 for (int iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
895 int numberIn = thisVector->getNumElements();
896 const int *COIN_RESTRICT indexIn = thisVector->getIndices();
897 const double *COIN_RESTRICT arrayIn = thisVector->denseVector();
898 CoinAbcScatterToList(arrayIn, abcSolution_, indexIn, abcPivotVariable_, numberIn);
899 // check Ax == b (for all)
900 abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, solutionBasic_);
901 #if 0
902 if (xxxxxx==0)
903 for (int i=0;i<numberTotal_;i++)
904 printf("sol %d %.19g\n",i,abcSolution_[i]);
905 if (xxxxxx==0)
906 for (int i=0;i<numberRows_;i++)
907 printf("basic %d %.19g\n",i,solutionBasic_[i]);
908 #endif
909 double multiplier = 131072.0;
910 largestPrimalError_ = CoinAbcMaximumAbsElementAndScale(solutionBasic_, multiplier, numberRows_);
911 double maxValue = 0.0;
912 for (int i = 0; i < numberRows_; i++) {
913 double value = fabs(solutionBasic_[i]);
914 if (value > maxValue) {
915 #if 0
916 if (xxxxxx==0)
917 printf("larger %.19gg at row %d\n",value,i);
918 maxValue=value;
919 #endif
920 }
921 }
922 if (largestPrimalError_ >= lastError) {
923 // restore
924 CoinIndexedVector *temp = thisVector;
925 thisVector = lastVector;
926 lastVector = temp;
927 //goodSolution = false;
928 break;
929 }
930 if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) {
931 // try and make better
932 numberRefinements++;
933 // save this
934 CoinIndexedVector *temp = thisVector;
935 thisVector = lastVector;
936 lastVector = temp;
937 int *COIN_RESTRICT indexOut = thisVector->getIndices();
938 int number = 0;
939 array = thisVector->denseVector();
940 thisVector->clear();
941 for (int iRow = 0; iRow < numberRows_; iRow++) {
942 double value = solutionBasic_[iRow];
943 if (value) {
944 array[iRow] = value;
945 indexOut[number++] = iRow;
946 solutionBasic_[iRow] = 0.0;
947 }
948 }
949 thisVector->setNumElements(number);
950 lastError = largestPrimalError_;
951 abcFactorization_->updateFullColumn(*thisVector);
952 double *previous = lastVector->denseVector();
953 number = 0;
954 multiplier = 1.0 / multiplier;
955 for (int iRow = 0; iRow < numberRows_; iRow++) {
956 double value = previous[iRow] + multiplier * array[iRow];
957 if (value) {
958 array[iRow] = value;
959 indexOut[number++] = iRow;
960 } else {
961 array[iRow] = 0.0;
962 }
963 }
964 thisVector->setNumElements(number);
965 } else {
966 break;
967 }
968 }
969
970 // solution as accurate as we are going to get
971 //if (!goodSolution) {
972 CoinAbcMemcpy(solutionBasic_, thisVector->denseVector(), numberRows_);
973 CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_);
974 arrayVector->clear();
975 previousVector->clear();
976 return numberRefinements;
977 }
978 // Moves basic stuff to basic area
moveToBasic(int which)979 void AbcSimplex::moveToBasic(int which)
980 {
981 if ((which & 1) != 0)
982 CoinAbcGatherFrom(abcSolution_, solutionBasic_, abcPivotVariable_, numberRows_);
983 if ((which & 2) != 0)
984 CoinAbcGatherFrom(abcCost_, costBasic_, abcPivotVariable_, numberRows_);
985 if ((which & 4) != 0)
986 CoinAbcGatherFrom(abcLower_, lowerBasic_, abcPivotVariable_, numberRows_);
987 if ((which & 8) != 0)
988 CoinAbcGatherFrom(abcUpper_, upperBasic_, abcPivotVariable_, numberRows_);
989 }
990 // now dual side
computeDuals(double * givenDjs,CoinIndexedVector * arrayVector,CoinIndexedVector * previousVector)991 int AbcSimplex::computeDuals(double *givenDjs, CoinIndexedVector *arrayVector, CoinIndexedVector *previousVector)
992 {
993 double *COIN_RESTRICT array = arrayVector->denseVector();
994 int *COIN_RESTRICT index = arrayVector->getIndices();
995 int number = 0;
996 if (!givenDjs) {
997 for (int iRow = 0; iRow < numberRows_; iRow++) {
998 double value = costBasic_[iRow];
999 if (value) {
1000 array[iRow] = value;
1001 index[number++] = iRow;
1002 }
1003 }
1004 } else {
1005 // dual values pass - djs may not be zero
1006 for (int iRow = 0; iRow < numberRows_; iRow++) {
1007 int iPivot = abcPivotVariable_[iRow];
1008 // make sure zero if done
1009 if (!pivoted(iPivot))
1010 givenDjs[iPivot] = 0.0;
1011 double value = abcCost_[iPivot] - givenDjs[iPivot];
1012 if (value) {
1013 array[iRow] = value;
1014 index[number++] = iRow;
1015 }
1016 }
1017 }
1018 arrayVector->setNumElements(number);
1019 // Btran basic costs and get as accurate as possible
1020 double lastError = COIN_DBL_MAX;
1021 CoinIndexedVector *thisVector = arrayVector;
1022 CoinIndexedVector *lastVector = previousVector;
1023 abcFactorization_->updateFullColumnTranspose(*thisVector);
1024
1025 int numberRefinements = 0;
1026 for (int iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
1027 // check basic reduced costs zero
1028 // put reduced cost of basic into djBasic_
1029 CoinAbcMemcpy(djBasic_, costBasic_, numberRows_);
1030 abcMatrix_->transposeTimesBasic(-1.0, thisVector->denseVector(), djBasic_);
1031 largestDualError_ = CoinAbcMaximumAbsElement(djBasic_, numberRows_);
1032 if (largestDualError_ >= lastError) {
1033 // restore
1034 CoinIndexedVector *temp = thisVector;
1035 thisVector = lastVector;
1036 lastVector = temp;
1037 break;
1038 }
1039 if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10
1040 && !givenDjs) {
1041 numberRefinements++;
1042 // try and make better
1043 // save this
1044 CoinIndexedVector *temp = thisVector;
1045 thisVector = lastVector;
1046 lastVector = temp;
1047 array = thisVector->denseVector();
1048 double multiplier = 131072.0;
1049 //array=djBasic_*multiplier
1050 CoinAbcMultiplyAdd(djBasic_, numberRows_, multiplier, array, 0.0);
1051 lastError = largestDualError_;
1052 abcFactorization_->updateFullColumnTranspose(*thisVector);
1053 #if ABC_DEBUG
1054 thisVector->checkClean();
1055 #endif
1056 multiplier = 1.0 / multiplier;
1057 double *COIN_RESTRICT previous = lastVector->denseVector();
1058 // array = array*multiplier+previous
1059 CoinAbcMultiplyAdd(previous, numberRows_, 1.0, array, multiplier);
1060 } else {
1061 break;
1062 }
1063 }
1064 // now look at dual solution
1065 CoinAbcMemcpy(dual_, thisVector->denseVector(), numberRows_);
1066 CoinAbcMemset0(thisVector->denseVector(), numberRows_);
1067 thisVector->setNumElements(0);
1068 if (numberRefinements) {
1069 CoinAbcMemset0(lastVector->denseVector(), numberRows_);
1070 lastVector->setNumElements(0);
1071 }
1072 CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_);
1073 abcMatrix_->transposeTimesNonBasic(-1.0, dual_, abcDj_);
1074 // If necessary - override results
1075 if (givenDjs) {
1076 // restore accurate duals
1077 CoinMemcpyN(abcDj_, numberTotal_, givenDjs);
1078 }
1079 //arrayVector->clear();
1080 //previousVector->clear();
1081 #if ABC_DEBUG
1082 arrayVector->checkClean();
1083 previousVector->checkClean();
1084 #endif
1085 return numberRefinements;
1086 }
1087
1088 /* Factorizes using current basis.
1089 solveType - 1 iterating, 0 initial, -1 external
1090 - 2 then iterating but can throw out of basis
1091 If 10 added then in primal values pass
1092 Return codes are as from AbcSimplexFactorization unless initial factorization
1093 when total number of singularities is returned.
1094 Special case is numberRows_+1 -> all slack basis.
1095 */
internalFactorize(int solveType)1096 int AbcSimplex::internalFactorize(int solveType)
1097 {
1098 assert(status_);
1099
1100 bool valuesPass = false;
1101 if (solveType >= 10) {
1102 valuesPass = true;
1103 solveType -= 10;
1104 }
1105 #if 0
1106 // Make sure everything is clean
1107 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1108 if(getInternalStatus(iSequence) == isFixed) {
1109 // double check fixed
1110 assert (abcUpper_[iSequence] == abcLower_[iSequence]);
1111 assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<100.0*primalTolerance_);
1112 } else if (getInternalStatus(iSequence) == isFree) {
1113 assert (abcUpper_[iSequence] == COIN_DBL_MAX && abcLower_[iSequence]==-COIN_DBL_MAX);
1114 } else if (getInternalStatus(iSequence) == atLowerBound) {
1115 assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<1000.0*primalTolerance_);
1116 } else if (getInternalStatus(iSequence) == atUpperBound) {
1117 assert (fabs(abcSolution_[iSequence]-abcUpper_[iSequence])<1000.0*primalTolerance_);
1118 } else if (getInternalStatus(iSequence) == superBasic) {
1119 assert (!valuesPass);
1120 }
1121 }
1122 #endif
1123 #if 0 //ndef NDEBUG
1124 // Make sure everything is clean
1125 double sumOutside=0.0;
1126 int numberOutside=0;
1127 //double sumOutsideLarge=0.0;
1128 int numberOutsideLarge=0;
1129 double sumInside=0.0;
1130 int numberInside=0;
1131 //double sumInsideLarge=0.0;
1132 int numberInsideLarge=0;
1133 char rowcol[] = {'R', 'C'};
1134 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1135 if(getInternalStatus(iSequence) == isFixed) {
1136 // double check fixed
1137 assert (abcUpper_[iSequence] == abcLower_[iSequence]);
1138 assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<primalTolerance_);
1139 } else if (getInternalStatus(iSequence) == isFree) {
1140 assert (abcUpper_[iSequence] == COIN_DBL_MAX && abcLower_[iSequence]==-COIN_DBL_MAX);
1141 } else if (getInternalStatus(iSequence) == atLowerBound) {
1142 assert (fabs(abcSolution_[iSequence]-abcLower_[iSequence])<1000.0*primalTolerance_);
1143 if (abcSolution_[iSequence]<abcLower_[iSequence]) {
1144 numberOutside++;
1145 #if ABC_NORMAL_DEBUG > 1
1146 if (handler_->logLevel()==63)
1147 printf("%c%d below by %g\n",rowcol[isColumn(iSequence)],sequenceWithin(iSequence),
1148 abcLower_[iSequence]-abcSolution_[iSequence]);
1149 #endif
1150 sumOutside-=abcSolution_[iSequence]-abcLower_[iSequence];
1151 if (abcSolution_[iSequence]<abcLower_[iSequence]-primalTolerance_)
1152 numberOutsideLarge++;
1153 } else if (abcSolution_[iSequence]>abcLower_[iSequence]) {
1154 numberInside++;
1155 sumInside+=abcSolution_[iSequence]-abcLower_[iSequence];
1156 if (abcSolution_[iSequence]>abcLower_[iSequence]+primalTolerance_)
1157 numberInsideLarge++;
1158 }
1159 } else if (getInternalStatus(iSequence) == atUpperBound) {
1160 assert (fabs(abcSolution_[iSequence]-abcUpper_[iSequence])<1000.0*primalTolerance_);
1161 if (abcSolution_[iSequence]>abcUpper_[iSequence]) {
1162 numberOutside++;
1163 #if ABC_NORMAL_DEBUG > 0
1164 if (handler_->logLevel()==63)
1165 printf("%c%d above by %g\n",rowcol[isColumn(iSequence)],sequenceWithin(iSequence),
1166 -(abcUpper_[iSequence]-abcSolution_[iSequence]));
1167 #endif
1168 sumOutside+=abcSolution_[iSequence]-abcUpper_[iSequence];
1169 if (abcSolution_[iSequence]>abcUpper_[iSequence]+primalTolerance_)
1170 numberOutsideLarge++;
1171 } else if (abcSolution_[iSequence]<abcUpper_[iSequence]) {
1172 numberInside++;
1173 sumInside-=abcSolution_[iSequence]-abcUpper_[iSequence];
1174 if (abcSolution_[iSequence]<abcUpper_[iSequence]-primalTolerance_)
1175 numberInsideLarge++;
1176 }
1177 } else if (getInternalStatus(iSequence) == superBasic) {
1178 //assert (!valuesPass);
1179 }
1180 }
1181 #if ABC_NORMAL_DEBUG > 0
1182 if (numberInside+numberOutside)
1183 printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n",
1184 numberOutside,sumOutside,numberOutsideLarge,
1185 numberInside,sumInside,numberInsideLarge);
1186 #endif
1187 #endif
1188 // *** replace below by cleanStatus
1189 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1190 AbcSimplex::Status status = getInternalStatus(iSequence);
1191 if (status != basic && status != isFixed && abcUpper_[iSequence] == abcLower_[iSequence])
1192 setInternalStatus(iSequence, isFixed);
1193 }
1194 if (numberIterations_ == baseIteration_ && !valuesPass) {
1195 double largeValue = this->largeValue();
1196 double *COIN_RESTRICT solution = abcSolution_;
1197 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1198 AbcSimplex::Status status = getInternalStatus(iSequence);
1199 if (status == superBasic) {
1200 double lower = abcLower_[iSequence];
1201 double upper = abcUpper_[iSequence];
1202 double value = solution[iSequence];
1203 AbcSimplex::Status thisStatus = isFree;
1204 if (lower > -largeValue || upper < largeValue) {
1205 if (lower != upper) {
1206 if (fabs(value - lower) < fabs(value - upper)) {
1207 thisStatus = AbcSimplex::atLowerBound;
1208 solution[iSequence] = lower;
1209 } else {
1210 thisStatus = AbcSimplex::atUpperBound;
1211 solution[iSequence] = upper;
1212 }
1213 } else {
1214 thisStatus = AbcSimplex::isFixed;
1215 solution[iSequence] = upper;
1216 }
1217 setInternalStatus(iSequence, thisStatus);
1218 }
1219 }
1220 }
1221 }
1222 int status = abcFactorization_->factorize(this, solveType, valuesPass);
1223 if (status) {
1224 handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
1225 << status
1226 << CoinMessageEol;
1227 return -1;
1228 } else if (!solveType) {
1229 // Initial basis - return number of singularities
1230 int numberSlacks = 0;
1231 for (int iRow = 0; iRow < numberRows_; iRow++) {
1232 if (getInternalStatus(iRow) == basic)
1233 numberSlacks++;
1234 }
1235 status = CoinMax(numberSlacks - numberRows_, 0);
1236 if (status)
1237 printf("%d singularities\n", status);
1238 // special case if all slack
1239 if (numberSlacks == numberRows_) {
1240 status = numberRows_ + 1;
1241 }
1242 }
1243 #ifdef ABC_USE_COIN_FACTORIZATION
1244 // sparse methods
1245 abcFactorization_->goSparse();
1246 #endif
1247 return status;
1248 }
1249 // Make sure no superbasic etc
cleanStatus(bool valuesPass)1250 void AbcSimplex::cleanStatus(bool valuesPass)
1251 {
1252 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1253 AbcSimplex::Status status = getInternalStatus(iSequence);
1254 if (status != basic && status != isFixed && abcUpper_[iSequence] == abcLower_[iSequence])
1255 setInternalStatus(iSequence, isFixed);
1256 }
1257 if (numberIterations_ == baseIteration_ && !valuesPass) {
1258 double largeValue = this->largeValue();
1259 double *COIN_RESTRICT solution = abcSolution_;
1260 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1261 AbcSimplex::Status status = getInternalStatus(iSequence);
1262 if (status == superBasic) {
1263 double lower = abcLower_[iSequence];
1264 double upper = abcUpper_[iSequence];
1265 double value = solution[iSequence];
1266 AbcSimplex::Status thisStatus = isFree;
1267 if (lower > -largeValue || upper < largeValue) {
1268 if (lower != upper) {
1269 if (fabs(value - lower) < fabs(value - upper)) {
1270 thisStatus = AbcSimplex::atLowerBound;
1271 solution[iSequence] = lower;
1272 } else {
1273 thisStatus = AbcSimplex::atUpperBound;
1274 solution[iSequence] = upper;
1275 }
1276 } else {
1277 thisStatus = AbcSimplex::isFixed;
1278 solution[iSequence] = upper;
1279 }
1280 setInternalStatus(iSequence, thisStatus);
1281 }
1282 }
1283 }
1284 }
1285 }
1286 // Sets objectiveValue_ from rawObjectiveValue_
setClpSimplexObjectiveValue()1287 void AbcSimplex::setClpSimplexObjectiveValue()
1288 {
1289 objectiveValue_ = rawObjectiveValue_ / (objectiveScale_ * rhsScale_);
1290 objectiveValue_ += objectiveOffset_;
1291 }
1292 /*
1293 This does basis housekeeping and does values for in/out variables.
1294 Can also decide to re-factorize
1295 */
housekeeping()1296 int AbcSimplex::housekeeping()
1297 {
1298 #define DELAYED_UPDATE
1299 #ifdef DELAYED_UPDATE
1300 if (algorithm_ < 0) {
1301 // modify dualout
1302 dualOut_ /= alpha_;
1303 dualOut_ *= -directionOut_;
1304 //double oldValue = valueIn_;
1305 if (directionIn_ == -1) {
1306 // as if from upper bound
1307 valueIn_ = upperIn_ + dualOut_;
1308 #if 0 //def ABC_DEBUG
1309 printf("In from upper dualout_ %g movement %g (old %g) valueIn_ %g using movement %g\n",
1310 dualOut_,movement_,movementOld,valueIn_,upperIn_+movement_);
1311 #endif
1312 } else {
1313 // as if from lower bound
1314 valueIn_ = lowerIn_ + dualOut_;
1315 #if 0 //def ABC_DEBUG
1316 printf("In from lower dualout_ %g movement %g (old %g) valueIn_ %g using movement %g\n",
1317 dualOut_,movement_,movementOld,valueIn_,lowerIn_+movement_);
1318 #endif
1319 }
1320 // outgoing
1321 if (directionOut_ > 0) {
1322 valueOut_ = lowerOut_;
1323 } else {
1324 valueOut_ = upperOut_;
1325 }
1326 #if ABC_NORMAL_DEBUG > 3
1327 if (rawObjectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16))
1328 printf("obj backwards %g %g\n", rawObjectiveValue_, oldobj);
1329 #endif
1330 }
1331 #endif
1332 #if 0 //ndef NDEBUG
1333 {
1334 int numberFlagged=0;
1335 for (int iRow = 0; iRow < numberRows_; iRow++) {
1336 int iPivot = abcPivotVariable_[iRow];
1337 if (flagged(iPivot))
1338 numberFlagged++;
1339 }
1340 assert (numberFlagged==numberFlagged_);
1341 }
1342 #endif
1343 // save value of incoming and outgoing
1344 numberIterations_++;
1345 changeMade_++; // something has happened
1346 // incoming variable
1347 if (handler_->logLevel() > 7) {
1348 //if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
1349 handler_->message(CLP_SIMPLEX_HOUSE1, messages_)
1350 << directionOut_
1351 << directionIn_ << theta_
1352 << dualOut_ << dualIn_ << alpha_
1353 << CoinMessageEol;
1354 if (getInternalStatus(sequenceIn_) == isFree) {
1355 handler_->message(CLP_SIMPLEX_FREEIN, messages_)
1356 << sequenceIn_
1357 << CoinMessageEol;
1358 }
1359 }
1360 // change of incoming
1361 char rowcol[] = { 'R', 'C' };
1362 if (abcUpper_[sequenceIn_] > 1.0e20 && abcLower_[sequenceIn_] < -1.0e20)
1363 progressFlag_ |= 2; // making real progress
1364 #ifndef DELAYED_UPDATE
1365 abcSolution_[sequenceIn_] = valueIn_;
1366 #endif
1367 if (abcUpper_[sequenceOut_] - abcLower_[sequenceOut_] < 1.0e-12)
1368 progressFlag_ |= 1; // making real progress
1369 if (sequenceIn_ != sequenceOut_) {
1370 if (alphaAccuracy_ > 0.0) {
1371 double value = fabs(alpha_);
1372 if (value > 1.0)
1373 alphaAccuracy_ *= value;
1374 else
1375 alphaAccuracy_ /= value;
1376 }
1377 setInternalStatus(sequenceIn_, basic);
1378 if (abcUpper_[sequenceOut_] - abcLower_[sequenceOut_] > 0) {
1379 // As Nonlinear costs may have moved bounds (to more feasible)
1380 // Redo using value
1381 if (fabs(valueOut_ - abcLower_[sequenceOut_]) < fabs(valueOut_ - abcUpper_[sequenceOut_])) {
1382 // going to lower
1383 setInternalStatus(sequenceOut_, atLowerBound);
1384 } else {
1385 // going to upper
1386 setInternalStatus(sequenceOut_, atUpperBound);
1387 }
1388 } else {
1389 // fixed
1390 setInternalStatus(sequenceOut_, isFixed);
1391 }
1392 #ifndef DELAYED_UPDATE
1393 abcSolution_[sequenceOut_] = valueOut_;
1394 #endif
1395 #if PARTITION_ROW_COPY == 1
1396 // move in row copy
1397 abcMatrix_->inOutUseful(sequenceIn_, sequenceOut_);
1398 #endif
1399 } else {
1400 //if (objective_->type()<2)
1401 //assert (fabs(theta_)>1.0e-13);
1402 // flip from bound to bound
1403 // As Nonlinear costs may have moved bounds (to more feasible)
1404 // Redo using value
1405 if (fabs(valueIn_ - abcLower_[sequenceIn_]) < fabs(valueIn_ - abcUpper_[sequenceIn_])) {
1406 // as if from upper bound
1407 setInternalStatus(sequenceIn_, atLowerBound);
1408 } else {
1409 // as if from lower bound
1410 setInternalStatus(sequenceIn_, atUpperBound);
1411 }
1412 }
1413 setClpSimplexObjectiveValue();
1414 if (handler_->logLevel() > 7) {
1415 //if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
1416 handler_->message(CLP_SIMPLEX_HOUSE2, messages_)
1417 << numberIterations_ << objectiveValue()
1418 << rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_)
1419 << rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_);
1420 handler_->printing(algorithm_ < 0) << dualOut_ << theta_;
1421 handler_->printing(algorithm_ > 0) << dualIn_ << theta_;
1422 handler_->message() << CoinMessageEol;
1423 }
1424 #if 0
1425 if (numberIterations_ > 10000)
1426 printf(" it %d %g %c%d %c%d\n"
1427 , numberIterations_, objectiveValue()
1428 , rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_)
1429 , rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_));
1430 #endif
1431 if (hitMaximumIterations())
1432 return 2;
1433 // check for small cycles
1434 int in = sequenceIn_;
1435 int out = sequenceOut_;
1436 int cycle = abcProgress_.cycle(in, out,
1437 directionIn_, directionOut_);
1438 if (cycle > 0) {
1439 #if ABC_NORMAL_DEBUG > 0
1440 if (handler_->logLevel() >= 63)
1441 printf("Cycle of %d\n", cycle);
1442 #endif
1443 // reset
1444 abcProgress_.startCheck();
1445 double random = randomNumberGenerator_.randomDouble();
1446 int extra = static_cast< int >(9.999 * random);
1447 int off[] = { 1, 1, 1, 1, 2, 2, 2, 3, 3, 4 };
1448 if (abcFactorization_->pivots() > cycle) {
1449 forceFactorization_ = CoinMax(1, cycle - off[extra]);
1450 } else {
1451 // need to reject something
1452 int iSequence;
1453 if (algorithm_ < 0) {
1454 iSequence = sequenceIn_;
1455 } else {
1456 /* should be better if don't reject incoming
1457 as it is in basis */
1458 iSequence = sequenceOut_;
1459 // so can't happen immediately again
1460 sequenceOut_ = -1;
1461 }
1462 char x = isColumn(iSequence) ? 'C' : 'R';
1463 if (handler_->logLevel() >= 63)
1464 handler_->message(CLP_SIMPLEX_FLAG, messages_)
1465 << x << sequenceWithin(iSequence)
1466 << CoinMessageEol;
1467 setFlagged(iSequence);
1468 //printf("flagging %d\n",iSequence);
1469 }
1470 return 1;
1471 }
1472 int invertNow = 0;
1473 // only time to re-factorize if one before real time
1474 // this is so user won't be surprised that maximumPivots has exact meaning
1475 int numberPivots = abcFactorization_->pivots();
1476 if (algorithm_ < 0)
1477 numberPivots++; // allow for update not done
1478 int maximumPivots = abcFactorization_->maximumPivots();
1479 int numberDense = 0; //abcFactorization_->numberDense();
1480 bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 > 2 * maximumIterations());
1481 if (numberPivots == maximumPivots || maximumPivots < 2) {
1482 // If dense then increase
1483 if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots && false) {
1484 abcFactorization_->maximumPivots(numberDense);
1485 }
1486 //printf("ZZ maxpivots %d its %d\n",numberPivots,maximumPivots);
1487 #if CLP_FACTORIZATION_NEW_TIMING > 1
1488 abcFactorization_->statsRefactor('M');
1489 #endif
1490 return 1;
1491 } else if ((abcFactorization_->timeToRefactorize() && !dontInvert)
1492 || invertNow) {
1493 //printf("ZZ time invertNow %s its %d\n",invertNow ? "yes":"no",numberPivots);
1494 #if CLP_FACTORIZATION_NEW_TIMING > 1
1495 abcFactorization_->statsRefactor('T');
1496 #endif
1497 return 1;
1498 } else if (forceFactorization_ > 0 &&
1499 #ifndef DELAYED_UPDATE
1500 abcFactorization_->pivots()
1501 #else
1502 abcFactorization_->pivots() + 1
1503 #endif
1504 >= forceFactorization_) {
1505 // relax
1506 forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
1507 if (forceFactorization_ > abcFactorization_->maximumPivots())
1508 forceFactorization_ = -1; //off
1509 //printf("ZZ forceFactor %d its %d\n",forceFactorization_,numberPivots);
1510 #if CLP_FACTORIZATION_NEW_TIMING > 1
1511 abcFactorization_->statsRefactor('F');
1512 #endif
1513 return 1;
1514 } else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
1515 // A bit worried problem may be cycling - lets factorize at random intervals for a short period
1516 int numberTooManyIterations = numberIterations_ - 10 * (numberRows_ + (numberColumns_ >> 2));
1517 double random = randomNumberGenerator_.randomDouble();
1518 int window = numberTooManyIterations % 5000;
1519 if (window < 2 * maximumPivots)
1520 random = 0.2 * random + 0.8; // randomly re-factorize but not too soon
1521 else
1522 random = 1.0; // switch off if not in window of opportunity
1523 int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
1524 if (abcFactorization_->pivots() >= random * maxNumber) {
1525 //printf("ZZ cycling a %d\n",numberPivots);
1526 return 1;
1527 } else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && numberIterations_ < 1000010 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
1528 //printf("ZZ cycling b %d\n",numberPivots);
1529 return 1;
1530 } else {
1531 // carry on iterating
1532 return 0;
1533 }
1534 } else {
1535 // carry on iterating
1536 return 0;
1537 }
1538 }
1539 // Swaps primal stuff
swapPrimalStuff()1540 void AbcSimplex::swapPrimalStuff()
1541 {
1542 if (sequenceOut_ < 0)
1543 return;
1544 assert(sequenceIn_ >= 0);
1545 abcSolution_[sequenceOut_] = valueOut_;
1546 abcSolution_[sequenceIn_] = valueIn_;
1547 solutionBasic_[pivotRow_] = valueIn_;
1548 if (algorithm_ < 0) {
1549 // and set bounds correctly
1550 static_cast< AbcSimplexDual * >(this)->originalBound(sequenceIn_);
1551 static_cast< AbcSimplexDual * >(this)->changeBound(sequenceOut_);
1552 } else {
1553 #if ABC_NORMAL_DEBUG > 2
1554 if (handler_->logLevel() == 63)
1555 printf("Code swapPrimalStuff for primal\n");
1556 #endif
1557 }
1558 if (pivotRow_ >= 0) { // may be flip in primal
1559 lowerBasic_[pivotRow_] = abcLower_[sequenceIn_];
1560 upperBasic_[pivotRow_] = abcUpper_[sequenceIn_];
1561 costBasic_[pivotRow_] = abcCost_[sequenceIn_];
1562 abcPivotVariable_[pivotRow_] = sequenceIn_;
1563 }
1564 }
1565 // Swaps dual stuff
swapDualStuff(int lastSequenceOut,int lastDirectionOut)1566 void AbcSimplex::swapDualStuff(int lastSequenceOut, int lastDirectionOut)
1567 {
1568 // outgoing
1569 // set dj to zero unless values pass
1570 if (lastSequenceOut >= 0) {
1571 if ((stateOfProblem_ & VALUES_PASS) == 0) {
1572 if (lastDirectionOut > 0) {
1573 abcDj_[lastSequenceOut] = theta_;
1574 } else {
1575 abcDj_[lastSequenceOut] = -theta_;
1576 }
1577 } else {
1578 if (lastDirectionOut > 0) {
1579 abcDj_[lastSequenceOut] -= theta_;
1580 } else {
1581 abcDj_[lastSequenceOut] += theta_;
1582 }
1583 }
1584 }
1585 if (sequenceIn_ >= 0) {
1586 abcDj_[sequenceIn_] = 0.0;
1587 //costBasic_[pivotRow_]=abcCost_[sequenceIn_];
1588 }
1589 }
solveMany(int number,ClpSimplex ** simplex)1590 static void solveMany(int number, ClpSimplex **simplex)
1591 {
1592 for (int i = 0; i < number - 1; i++)
1593 cilk_spawn simplex[i]->dual(0);
1594 simplex[number - 1]->dual(0);
1595 cilk_sync;
1596 }
crash(int type)1597 void AbcSimplex::crash(int type)
1598 {
1599 int i;
1600 for (i = 0; i < numberRows_; i++) {
1601 if (getInternalStatus(i) != basic)
1602 break;
1603 }
1604 if (i < numberRows_)
1605 return;
1606 // all slack
1607 if (type == 3) {
1608 // decomposition crash
1609 if (!abcMatrix_->gotRowCopy())
1610 abcMatrix_->createRowCopy();
1611 //const double * element = abcMatrix_->getPackedMatrix()->getElements();
1612 const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts();
1613 const int *row = abcMatrix_->getPackedMatrix()->getIndices();
1614 //const double * elementByRow = abcMatrix_->rowElements();
1615 const CoinBigIndex *rowStart = abcMatrix_->rowStart();
1616 const CoinBigIndex *rowEnd = abcMatrix_->rowEnd();
1617 const int *column = abcMatrix_->rowColumns();
1618 int *blockStart = usefulArray_[0].getIndices();
1619 int *columnBlock = blockStart + numberRows_;
1620 int *nextColumn = usefulArray_[1].getIndices();
1621 int *blockCount = nextColumn + numberColumns_;
1622 int direction[2] = { -1, 1 };
1623 int bestBreak = -1;
1624 double bestValue = 0.0;
1625 int iPass = 0;
1626 int halfway = (numberRows_ + 1) / 2;
1627 int firstMaster = -1;
1628 int lastMaster = -2;
1629 int numberBlocks = 0;
1630 int largestRows = 0;
1631 int largestColumns = 0;
1632 int numberEmpty = 0;
1633 int numberMaster = 0;
1634 int numberEmptyColumns = 0;
1635 int numberMasterColumns = 0;
1636 while (iPass < 2) {
1637 int increment = direction[iPass];
1638 int start = increment > 0 ? 0 : numberRows_ - 1;
1639 int stop = increment > 0 ? numberRows_ : -1;
1640 numberBlocks = 0;
1641 int thisBestBreak = -1;
1642 double thisBestValue = COIN_DBL_MAX;
1643 int numberRowsDone = 0;
1644 int numberMarkedColumns = 0;
1645 int maximumBlockSize = 0;
1646 for (int i = 0; i < numberRows_; i++) {
1647 blockStart[i] = -1;
1648 blockCount[i] = 0;
1649 }
1650 for (int i = 0; i < numberColumns_; i++) {
1651 columnBlock[i] = -1;
1652 nextColumn[i] = -1;
1653 }
1654 for (int iRow = start; iRow != stop; iRow += increment) {
1655 int iBlock = -1;
1656 for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
1657 int iColumn = column[j] - maximumAbcNumberRows_;
1658 int whichColumnBlock = columnBlock[iColumn];
1659 if (whichColumnBlock >= 0) {
1660 // column marked
1661 if (iBlock < 0) {
1662 // put row in that block
1663 iBlock = whichColumnBlock;
1664 } else if (iBlock != whichColumnBlock) {
1665 // merge
1666 blockCount[iBlock] += blockCount[whichColumnBlock];
1667 blockCount[whichColumnBlock] = 0;
1668 int jColumn = blockStart[whichColumnBlock];
1669 #ifndef NDEBUG
1670 int iiColumn = iColumn;
1671 #endif
1672 while (jColumn >= 0) {
1673 assert(columnBlock[jColumn] == whichColumnBlock);
1674 columnBlock[jColumn] = iBlock;
1675 #ifndef NDEBUG
1676 if (jColumn == iiColumn)
1677 iiColumn = -1;
1678 #endif
1679 iColumn = jColumn;
1680 jColumn = nextColumn[jColumn];
1681 }
1682 #ifndef NDEBUG
1683 assert(iiColumn == -1);
1684 #endif
1685 nextColumn[iColumn] = blockStart[iBlock];
1686 blockStart[iBlock] = blockStart[whichColumnBlock];
1687 blockStart[whichColumnBlock] = -1;
1688 }
1689 }
1690 }
1691 int n = numberMarkedColumns;
1692 if (iBlock < 0) {
1693 //new block
1694 if (rowEnd[iRow] > rowStart[iRow]) {
1695 numberBlocks++;
1696 iBlock = numberBlocks;
1697 int jColumn = column[rowStart[iRow]] - maximumAbcNumberRows_;
1698 columnBlock[jColumn] = iBlock;
1699 blockStart[iBlock] = jColumn;
1700 numberMarkedColumns++;
1701 for (CoinBigIndex j = rowStart[iRow] + 1; j < rowEnd[iRow]; j++) {
1702 int iColumn = column[j] - maximumAbcNumberRows_;
1703 columnBlock[iColumn] = iBlock;
1704 numberMarkedColumns++;
1705 nextColumn[jColumn] = iColumn;
1706 jColumn = iColumn;
1707 }
1708 blockCount[iBlock] = numberMarkedColumns - n;
1709 } else {
1710 // empty
1711 }
1712 } else {
1713 // put in existing block
1714 int jColumn = blockStart[iBlock];
1715 for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
1716 int iColumn = column[j] - maximumAbcNumberRows_;
1717 assert(columnBlock[iColumn] < 0 || columnBlock[iColumn] == iBlock);
1718 if (columnBlock[iColumn] < 0) {
1719 columnBlock[iColumn] = iBlock;
1720 numberMarkedColumns++;
1721 nextColumn[iColumn] = jColumn;
1722 jColumn = iColumn;
1723 }
1724 }
1725 blockStart[iBlock] = jColumn;
1726 blockCount[iBlock] += numberMarkedColumns - n;
1727 }
1728 if (iBlock >= 0)
1729 maximumBlockSize = CoinMax(maximumBlockSize, blockCount[iBlock]);
1730 numberRowsDone++;
1731 if (thisBestValue * numberRowsDone > maximumBlockSize && numberRowsDone > halfway) {
1732 thisBestBreak = iRow;
1733 thisBestValue = static_cast< double >(maximumBlockSize) / static_cast< double >(numberRowsDone);
1734 }
1735 }
1736 if (thisBestBreak == stop)
1737 thisBestValue = COIN_DBL_MAX;
1738 iPass++;
1739 if (iPass == 1) {
1740 bestBreak = thisBestBreak;
1741 bestValue = thisBestValue;
1742 } else {
1743 if (bestValue < thisBestValue) {
1744 firstMaster = 0;
1745 lastMaster = bestBreak;
1746 } else {
1747 firstMaster = thisBestBreak; // ? +1
1748 lastMaster = numberRows_;
1749 }
1750 }
1751 }
1752 bool useful = false;
1753 if (firstMaster < lastMaster) {
1754 for (int i = 0; i < numberRows_; i++)
1755 blockStart[i] = -1;
1756 for (int i = firstMaster; i < lastMaster; i++)
1757 blockStart[i] = -2;
1758 for (int i = 0; i < numberColumns_; i++)
1759 columnBlock[i] = -1;
1760 int firstRow = 0;
1761 numberBlocks = -1;
1762 while (true) {
1763 for (; firstRow < numberRows_; firstRow++) {
1764 if (blockStart[firstRow] == -1)
1765 break;
1766 }
1767 if (firstRow == numberRows_)
1768 break;
1769 int nRows = 0;
1770 numberBlocks++;
1771 int numberStack = 1;
1772 blockCount[0] = firstRow;
1773 while (numberStack) {
1774 int iRow = blockCount[--numberStack];
1775 for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
1776 int iColumn = column[j] - maximumAbcNumberRows_;
1777 int iBlock = columnBlock[iColumn];
1778 if (iBlock < 0) {
1779 columnBlock[iColumn] = numberBlocks;
1780 for (CoinBigIndex k = columnStart[iColumn];
1781 k < columnStart[iColumn + 1]; k++) {
1782 int jRow = row[k];
1783 int rowBlock = blockStart[jRow];
1784 if (rowBlock == -1) {
1785 nRows++;
1786 blockStart[jRow] = numberBlocks;
1787 blockCount[numberStack++] = jRow;
1788 }
1789 }
1790 }
1791 }
1792 }
1793 if (!nRows) {
1794 // empty!!
1795 numberBlocks--;
1796 }
1797 firstRow++;
1798 }
1799 // adjust
1800 numberBlocks++;
1801 for (int i = 0; i < numberBlocks; i++) {
1802 blockCount[i] = 0;
1803 nextColumn[i] = 0;
1804 }
1805 for (int iRow = 0; iRow < numberRows_; iRow++) {
1806 int iBlock = blockStart[iRow];
1807 if (iBlock >= 0) {
1808 blockCount[iBlock]++;
1809 } else {
1810 if (iBlock == -2)
1811 numberMaster++;
1812 else
1813 numberEmpty++;
1814 }
1815 }
1816 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1817 int iBlock = columnBlock[iColumn];
1818 if (iBlock >= 0) {
1819 nextColumn[iBlock]++;
1820 } else {
1821 if (columnStart[iColumn + 1] > columnStart[iColumn])
1822 numberMasterColumns++;
1823 else
1824 numberEmptyColumns++;
1825 }
1826 }
1827 for (int i = 0; i < numberBlocks; i++) {
1828 if (blockCount[i] + nextColumn[i] > largestRows + largestColumns) {
1829 largestRows = blockCount[i];
1830 largestColumns = nextColumn[i];
1831 }
1832 }
1833 useful = true;
1834 if (numberMaster > halfway || largestRows * 3 > numberRows_)
1835 useful = false;
1836 }
1837 if (useful) {
1838 #if ABC_NORMAL_DEBUG > 0
1839 if (logLevel() >= 2)
1840 printf("%d master rows %d <= < %d\n", lastMaster - firstMaster,
1841 firstMaster, lastMaster);
1842 printf("%s %d blocks (largest %d,%d), %d master rows (%d empty) out of %d, %d master columns (%d empty) out of %d\n",
1843 useful ? "**Useful" : "NoGood",
1844 numberBlocks, largestRows, largestColumns, numberMaster, numberEmpty, numberRows_,
1845 numberMasterColumns, numberEmptyColumns, numberColumns_);
1846 if (logLevel() >= 2) {
1847 for (int i = 0; i < numberBlocks; i++)
1848 printf("Block %d has %d rows and %d columns\n",
1849 i, blockCount[i], nextColumn[i]);
1850 }
1851 #endif
1852 #define NUMBER_DW_BLOCKS 20
1853 int minSize1 = (numberRows_ - numberMaster + NUMBER_DW_BLOCKS - 1) / NUMBER_DW_BLOCKS;
1854 int minSize2 = (numberRows_ - numberMaster + 2 * NUMBER_DW_BLOCKS - 1) / (2 * NUMBER_DW_BLOCKS);
1855 int *backRow = usefulArray_[1].getIndices();
1856 // first sort
1857 for (int i = 0; i < numberBlocks; i++) {
1858 backRow[i] = -(3 * blockCount[i] + 0 * nextColumn[i]);
1859 nextColumn[i] = i;
1860 }
1861 CoinSort_2(backRow, backRow + numberBlocks, nextColumn);
1862 // keep if >minSize2 or sum >minSize1;
1863 int n = 0;
1864 for (int i = 0; i < numberBlocks; i++) {
1865 int originalBlock = nextColumn[i];
1866 if (blockCount[originalBlock] < minSize2)
1867 break;
1868 n++;
1869 }
1870 int size = minSize1;
1871 for (int i = n; i < numberBlocks; i++) {
1872 int originalBlock = nextColumn[i];
1873 size -= blockCount[originalBlock];
1874 if (size > 0 && i < numberBlocks - 1) {
1875 blockCount[originalBlock] = -1;
1876 } else {
1877 size = minSize1;
1878 n++;
1879 }
1880 }
1881 int n2 = numberBlocks;
1882 numberBlocks = n;
1883 for (int i = n2 - 1; i >= 0; i--) {
1884 int originalBlock = nextColumn[i];
1885 if (blockCount[originalBlock] > 0)
1886 n--;
1887 blockCount[originalBlock] = n;
1888 }
1889 assert(!n);
1890 for (int iRow = 0; iRow < numberRows_; iRow++) {
1891 int iBlock = blockStart[iRow];
1892 if (iBlock >= 0)
1893 blockStart[iRow] = blockCount[iBlock];
1894 }
1895 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1896 int iBlock = columnBlock[iColumn];
1897 if (iBlock >= 0)
1898 columnBlock[iColumn] = blockCount[iBlock];
1899 }
1900 // stick to Clp for now
1901 ClpSimplex **simplex = new ClpSimplex *[numberBlocks];
1902 for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1903 int nRow = 0;
1904 int nColumn = 0;
1905 for (int iRow = 0; iRow < numberRows_; iRow++) {
1906 if (blockStart[iRow] == iBlock)
1907 blockCount[nRow++] = iRow;
1908 }
1909 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1910 if (columnBlock[iColumn] == iBlock)
1911 nextColumn[nColumn++] = iColumn;
1912 }
1913 simplex[iBlock] = new ClpSimplex(this, nRow, blockCount, nColumn, nextColumn);
1914 simplex[iBlock]->setSpecialOptions(simplex[iBlock]->specialOptions() & (~65536));
1915 if (logLevel() < 2)
1916 simplex[iBlock]->setLogLevel(0);
1917 }
1918 solveMany(numberBlocks, simplex);
1919 int numberBasic = numberMaster;
1920 int numberStructurals = 0;
1921 for (int i = 0; i < numberMaster; i++)
1922 abcPivotVariable_[i] = i + firstMaster;
1923 for (int iBlock = 0; iBlock < numberBlocks; iBlock++) {
1924 int nRow = 0;
1925 int nColumn = 0;
1926 // from Clp enum to Abc enum (and bound flip)
1927 unsigned char lookupToAbcSlack[6] = { 4, 6, 0 /*1*/, 1 /*0*/, 5, 7 };
1928 unsigned char *COIN_RESTRICT getStatus = simplex[iBlock]->statusArray() + simplex[iBlock]->numberColumns();
1929 double *COIN_RESTRICT solutionSaved = solutionSaved_;
1930 double *COIN_RESTRICT lowerSaved = lowerSaved_;
1931 double *COIN_RESTRICT upperSaved = upperSaved_;
1932 for (int iRow = 0; iRow < numberRows_; iRow++) {
1933 if (blockStart[iRow] == iBlock) {
1934 unsigned char status = getStatus[nRow++] & 7;
1935 AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbcSlack[status]);
1936 if (status != ClpSimplex::basic) {
1937 double lowerValue = lowerSaved[iRow];
1938 double upperValue = upperSaved[iRow];
1939 if (lowerValue == -COIN_DBL_MAX) {
1940 if (upperValue == COIN_DBL_MAX) {
1941 // free
1942 abcStatus = isFree;
1943 } else {
1944 abcStatus = atUpperBound;
1945 }
1946 } else if (upperValue == COIN_DBL_MAX) {
1947 abcStatus = atLowerBound;
1948 } else if (lowerValue == upperValue) {
1949 abcStatus = isFixed;
1950 }
1951 switch (abcStatus) {
1952 case isFixed:
1953 case atLowerBound:
1954 solutionSaved[iRow] = lowerValue;
1955 break;
1956 case atUpperBound:
1957 solutionSaved[iRow] = upperValue;
1958 break;
1959 default:
1960 break;
1961 }
1962 } else {
1963 // basic
1964 abcPivotVariable_[numberBasic++] = iRow;
1965 }
1966 internalStatus_[iRow] = abcStatus;
1967 }
1968 }
1969 // from Clp enum to Abc enum
1970 unsigned char lookupToAbc[6] = { 4, 6, 1, 0, 5, 7 };
1971 unsigned char *COIN_RESTRICT putStatus = internalStatus_ + maximumAbcNumberRows_;
1972 getStatus = simplex[iBlock]->statusArray();
1973 solutionSaved += maximumAbcNumberRows_;
1974 lowerSaved += maximumAbcNumberRows_;
1975 upperSaved += maximumAbcNumberRows_;
1976 int numberSaved = numberBasic;
1977 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1978 if (columnBlock[iColumn] == iBlock) {
1979 unsigned char status = getStatus[nColumn++] & 7;
1980 AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbc[status]);
1981 if (status != ClpSimplex::basic) {
1982 double lowerValue = lowerSaved[iColumn];
1983 double upperValue = upperSaved[iColumn];
1984 if (lowerValue == -COIN_DBL_MAX) {
1985 if (upperValue == COIN_DBL_MAX) {
1986 // free
1987 abcStatus = isFree;
1988 } else {
1989 abcStatus = atUpperBound;
1990 }
1991 } else if (upperValue == COIN_DBL_MAX) {
1992 abcStatus = atLowerBound;
1993 } else if (lowerValue == upperValue) {
1994 abcStatus = isFixed;
1995 } else if (abcStatus == isFree) {
1996 abcStatus = superBasic;
1997 }
1998 switch (abcStatus) {
1999 case isFixed:
2000 case atLowerBound:
2001 solutionSaved[iColumn] = lowerValue;
2002 break;
2003 case atUpperBound:
2004 solutionSaved[iColumn] = upperValue;
2005 break;
2006 default:
2007 break;
2008 }
2009 } else {
2010 // basic
2011 if (numberBasic < numberRows_)
2012 abcPivotVariable_[numberBasic++] = iColumn + maximumAbcNumberRows_;
2013 else
2014 abcStatus = superBasic;
2015 }
2016 putStatus[iColumn] = abcStatus;
2017 }
2018 }
2019 numberStructurals += numberBasic - numberSaved;
2020 delete simplex[iBlock];
2021 }
2022 #if ABC_NORMAL_DEBUG > 0
2023 printf("%d structurals put into basis\n", numberStructurals);
2024 #endif
2025 if (numberBasic < numberRows_) {
2026 for (int iRow = 0; iRow < numberRows_; iRow++) {
2027 AbcSimplex::Status status = getInternalStatus(iRow);
2028 if (status != AbcSimplex::basic) {
2029 abcPivotVariable_[numberBasic++] = iRow;
2030 setInternalStatus(iRow, basic);
2031 if (numberBasic == numberRows_)
2032 break;
2033 }
2034 }
2035 }
2036 assert(numberBasic == numberRows_);
2037 #if 0
2038 int n2=0;
2039 for (int i=0;i<numberTotal_;i++) {
2040 if (getInternalStatus(i)==basic)
2041 n2++;
2042 }
2043 assert (n2==numberRows_);
2044 std::sort(abcPivotVariable_,abcPivotVariable_+numberRows_);
2045 n2=-1;
2046 for (int i=0;i<numberRows_;i++) {
2047 assert (abcPivotVariable_[i]>n2);
2048 n2=abcPivotVariable_[i];
2049 }
2050 #endif
2051 delete[] simplex;
2052 return;
2053 } else {
2054 // try another
2055 type = 2;
2056 }
2057 }
2058 if (type == 1) {
2059 const double *linearObjective = abcCost_ + maximumAbcNumberRows_;
2060 double gamma = 0.0;
2061 int numberTotal = numberRows_ + numberColumns_;
2062 double *q = new double[numberTotal];
2063 double *v = q + numberColumns_;
2064 int *which = new int[numberTotal + 3 * numberRows_];
2065 int *ii = which + numberColumns_;
2066 int *r = ii + numberRows_;
2067 int *pivoted = r + numberRows_;
2068 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2069 gamma = CoinMax(gamma, linearObjective[iColumn]);
2070 }
2071 for (int iRow = 0; iRow < numberRows_; iRow++) {
2072 double lowerBound = abcLower_[iRow];
2073 double upperBound = abcUpper_[iRow];
2074 pivoted[iRow] = -1;
2075 ii[iRow] = 0;
2076 r[iRow] = 0;
2077 v[iRow] = COIN_DBL_MAX;
2078 if (lowerBound == upperBound)
2079 continue;
2080 if (lowerBound <= 0.0 && upperBound >= 0.0) {
2081 pivoted[iRow] = iRow;
2082 ii[iRow] = 1;
2083 r[iRow] = 1;
2084 }
2085 }
2086 int nPossible = 0;
2087 int lastPossible = 0;
2088 double cMaxDiv;
2089 if (gamma)
2090 cMaxDiv = 1.0 / (1000.0 * gamma);
2091 else
2092 cMaxDiv = 1.0;
2093 const double *columnLower = abcLower_ + maximumAbcNumberRows_;
2094 const double *columnUpper = abcUpper_ + maximumAbcNumberRows_;
2095 for (int iPass = 0; iPass < 3; iPass++) {
2096 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2097 double lowerBound = columnLower[iColumn];
2098 double upperBound = columnUpper[iColumn];
2099 if (lowerBound == upperBound)
2100 continue;
2101 double qValue;
2102 if (lowerBound > -1.0e20) {
2103 if (upperBound < 1.0e20) {
2104 // both
2105 qValue = lowerBound - upperBound;
2106 if (iPass != 2)
2107 qValue = COIN_DBL_MAX;
2108 } else {
2109 // just lower
2110 qValue = lowerBound;
2111 if (iPass != 1)
2112 qValue = COIN_DBL_MAX;
2113 }
2114 } else {
2115 if (upperBound < 1.0e20) {
2116 // just upper
2117 qValue = -upperBound;
2118 if (iPass != 1)
2119 qValue = COIN_DBL_MAX;
2120 } else {
2121 // free
2122 qValue = 0.0;
2123 if (iPass != 0)
2124 qValue = COIN_DBL_MAX;
2125 }
2126 }
2127 if (qValue != COIN_DBL_MAX) {
2128 which[nPossible] = iColumn;
2129 q[nPossible++] = qValue + linearObjective[iColumn] * cMaxDiv;
2130 ;
2131 }
2132 }
2133 CoinSort_2(q + lastPossible, q + nPossible, which + lastPossible);
2134 lastPossible = nPossible;
2135 }
2136 const double *element = abcMatrix_->getPackedMatrix()->getElements();
2137 const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts();
2138 //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths();
2139 const int *row = abcMatrix_->getPackedMatrix()->getIndices();
2140 int nPut = 0;
2141 for (int i = 0; i < nPossible; i++) {
2142 int iColumn = which[i];
2143 double maxAlpha = 0.0;
2144 int kRow = -1;
2145 double alternativeAlpha = 0.0;
2146 int jRow = -1;
2147 bool canTake = true;
2148 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
2149 int iRow = row[j];
2150 double alpha = fabs(element[j]);
2151 if (alpha > 0.01 * v[iRow]) {
2152 canTake = false;
2153 } else if (!ii[iRow] && alpha > alternativeAlpha) {
2154 alternativeAlpha = alpha;
2155 jRow = iRow;
2156 }
2157 if (!r[iRow] && alpha > maxAlpha) {
2158 maxAlpha = alpha;
2159 kRow = iRow;
2160 }
2161 }
2162 // only really works if scaled
2163 if (maxAlpha > 0.99) {
2164 pivoted[kRow] = iColumn + maximumAbcNumberRows_;
2165 v[kRow] = maxAlpha;
2166 ii[kRow] = 1;
2167 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
2168 int iRow = row[j];
2169 r[iRow]++;
2170 }
2171 nPut++;
2172 } else if (canTake && jRow >= 0) {
2173 pivoted[jRow] = iColumn + maximumAbcNumberRows_;
2174 v[jRow] = maxAlpha;
2175 ii[jRow] = 1;
2176 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
2177 int iRow = row[j];
2178 r[iRow]++;
2179 }
2180 nPut++;
2181 }
2182 }
2183 for (int iRow = 0; iRow < numberRows_; iRow++) {
2184 int iSequence = pivoted[iRow];
2185 if (iSequence >= 0 && iSequence < numberColumns_) {
2186 abcPivotVariable_[iRow] = iSequence;
2187 if (fabs(abcLower_[iRow]) < fabs(abcUpper_[iRow])) {
2188 setInternalStatus(iRow, atLowerBound);
2189 abcSolution_[iRow] = abcLower_[iRow];
2190 solutionSaved_[iRow] = abcLower_[iRow];
2191 } else {
2192 setInternalStatus(iRow, atUpperBound);
2193 abcSolution_[iRow] = abcUpper_[iRow];
2194 solutionSaved_[iRow] = abcUpper_[iRow];
2195 }
2196 setInternalStatus(iSequence, basic);
2197 }
2198 }
2199 #if ABC_NORMAL_DEBUG > 0
2200 printf("%d put into basis\n", nPut);
2201 #endif
2202 delete[] q;
2203 delete[] which;
2204 } else if (type == 2) {
2205 //return;
2206 int numberBad = 0;
2207 CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_);
2208 // Work on savedSolution and current
2209 int iVector = getAvailableArray();
2210 #define ALLOW_BAD_DJS
2211 #ifdef ALLOW_BAD_DJS
2212 double *modDj = usefulArray_[iVector].denseVector();
2213 #endif
2214 for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
2215 double dj = abcDj_[iSequence];
2216 modDj[iSequence] = 0.0;
2217 if (getInternalStatus(iSequence) == atLowerBound) {
2218 if (dj < -dualTolerance_) {
2219 double costThru = -dj * (abcUpper_[iSequence] - abcLower_[iSequence]);
2220 if (costThru < dualBound_) {
2221 // flip
2222 setInternalStatus(iSequence, atUpperBound);
2223 solutionSaved_[iSequence] = abcUpper_[iSequence];
2224 abcSolution_[iSequence] = abcUpper_[iSequence];
2225 } else {
2226 numberBad++;
2227 #ifdef ALLOW_BAD_DJS
2228 modDj[iSequence] = dj;
2229 dj = 0.0;
2230 #else
2231 break;
2232 #endif
2233 }
2234 } else {
2235 dj = CoinMax(dj, 0.0);
2236 }
2237 } else if (getInternalStatus(iSequence) == atLowerBound) {
2238 if (dj > dualTolerance_) {
2239 double costThru = dj * (abcUpper_[iSequence] - abcLower_[iSequence]);
2240 if (costThru < dualBound_) {
2241 assert(abcUpper_[iSequence] - abcLower_[iSequence] < 1.0e10);
2242 // flip
2243 setInternalStatus(iSequence, atLowerBound);
2244 solutionSaved_[iSequence] = abcLower_[iSequence];
2245 abcSolution_[iSequence] = abcLower_[iSequence];
2246 } else {
2247 numberBad++;
2248 #ifdef ALLOW_BAD_DJS
2249 modDj[iSequence] = dj;
2250 dj = 0.0;
2251 #else
2252 break;
2253 #endif
2254 }
2255 } else {
2256 dj = CoinMin(dj, 0.0);
2257 }
2258 } else {
2259 if (fabs(dj) < dualTolerance_) {
2260 dj = 0.0;
2261 } else {
2262 numberBad++;
2263 #ifdef ALLOW_BAD_DJS
2264 modDj[iSequence] = dj;
2265 dj = 0.0;
2266 #else
2267 break;
2268 #endif
2269 }
2270 }
2271 abcDj_[iSequence] = dj;
2272 }
2273 #ifndef ALLOW_BAD_DJS
2274 if (numberBad) {
2275 //CoinAbcMemset0(modDj+maximumAbcNumberRows_,numberColumns_);
2276 return;
2277 }
2278 #endif
2279 if (!abcMatrix_->gotRowCopy())
2280 abcMatrix_->createRowCopy();
2281 //const double * element = abcMatrix_->getPackedMatrix()->getElements();
2282 const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts() - maximumAbcNumberRows_;
2283 const int *row = abcMatrix_->getPackedMatrix()->getIndices();
2284 const double *elementByRow = abcMatrix_->rowElements();
2285 const CoinBigIndex *rowStart = abcMatrix_->rowStart();
2286 const CoinBigIndex *rowEnd = abcMatrix_->rowEnd();
2287 const int *column = abcMatrix_->rowColumns();
2288 CoinAbcMemset0(solutionBasic_, numberRows_);
2289 CoinAbcMemcpy(lowerBasic_, abcLower_, numberRows_);
2290 CoinAbcMemcpy(upperBasic_, abcUpper_, numberRows_);
2291 abcMatrix_->timesIncludingSlacks(-1.0, solutionSaved_, solutionBasic_);
2292 const double multiplier[] = { 1.0, -1.0 };
2293 int nBasic = 0;
2294 int *index = usefulArray_[iVector].getIndices();
2295 int iVector2 = getAvailableArray();
2296 int *index2 = usefulArray_[iVector2].getIndices();
2297 double *sort = usefulArray_[iVector2].denseVector();
2298 int average = CoinMax(5, rowEnd[numberRows_ - 1] / (8 * numberRows_));
2299 int nPossible = 0;
2300 if (numberRows_ > 10000) {
2301 // allow more
2302 for (int iRow = 0; iRow < numberRows_; iRow++) {
2303 double solutionValue = solutionBasic_[iRow];
2304 double infeasibility = 0.0;
2305 double lowerValue = lowerBasic_[iRow];
2306 double upperValue = upperBasic_[iRow];
2307 if (solutionValue < lowerValue - primalTolerance_) {
2308 infeasibility = -(lowerValue - solutionValue);
2309 } else if (solutionValue > upperValue + primalTolerance_) {
2310 infeasibility = upperValue - solutionValue;
2311 }
2312 int length = rowEnd[iRow] - rowStart[iRow];
2313 if (infeasibility)
2314 index2[nPossible++] = length;
2315 }
2316 std::sort(index2, index2 + nPossible);
2317 // see how much we need to get numberRows/10 or nPossible/3
2318 average = CoinMax(average, index2[CoinMin(numberRows_ / 10, nPossible / 3)]);
2319 nPossible = 0;
2320 }
2321 for (int iRow = 0; iRow < numberRows_; iRow++) {
2322 double solutionValue = solutionBasic_[iRow];
2323 double infeasibility = 0.0;
2324 double lowerValue = lowerBasic_[iRow];
2325 double upperValue = upperBasic_[iRow];
2326 if (solutionValue < lowerValue - primalTolerance_) {
2327 infeasibility = -(lowerValue - solutionValue);
2328 } else if (solutionValue > upperValue + primalTolerance_) {
2329 infeasibility = upperValue - solutionValue;
2330 }
2331 int length = rowEnd[iRow] - rowStart[iRow];
2332 if (infeasibility && length < average) {
2333 index2[nPossible] = iRow;
2334 sort[nPossible++] = 1.0e5 * infeasibility - iRow;
2335 //sort[nPossible++]=1.0e9*length-iRow;//infeasibility;
2336 }
2337 }
2338 CoinSort_2(sort, sort + nPossible, index2);
2339 for (int iWhich = 0; iWhich < nPossible; iWhich++) {
2340 int iRow = index2[iWhich];
2341 sort[iWhich] = 0.0;
2342 if (abcDj_[iRow])
2343 continue; // marked as invalid
2344 double solutionValue = solutionBasic_[iRow];
2345 double infeasibility = 0.0;
2346 double lowerValue = lowerBasic_[iRow];
2347 double upperValue = upperBasic_[iRow];
2348 if (solutionValue < lowerValue - primalTolerance_) {
2349 infeasibility = lowerValue - solutionValue;
2350 } else if (solutionValue > upperValue + primalTolerance_) {
2351 infeasibility = upperValue - solutionValue;
2352 }
2353 assert(infeasibility);
2354 double direction = infeasibility > 0 ? 1.0 : -1.0;
2355 infeasibility *= direction;
2356 int whichColumn = -1;
2357 double upperTheta = 1.0e30;
2358 int whichColumn2 = -1;
2359 double upperTheta2 = 1.0e30;
2360 double costThru = 0.0;
2361 int nThru = 0;
2362 for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2363 int iSequence = column[j];
2364 assert(iSequence >= maximumAbcNumberRows_);
2365 double dj = abcDj_[iSequence];
2366 double tableauValue = -elementByRow[j] * direction;
2367 unsigned char iStatus = internalStatus_[iSequence] & 7;
2368 if ((iStatus & 4) == 0) {
2369 double mult = multiplier[iStatus];
2370 double alpha = tableauValue * mult;
2371 double oldValue = dj * mult;
2372 assert(oldValue > -1.0e-2);
2373 double value = oldValue - upperTheta * alpha;
2374 if (value < 0.0) {
2375 upperTheta2 = upperTheta;
2376 whichColumn2 = whichColumn;
2377 costThru = alpha * (abcUpper_[iSequence] - abcLower_[iSequence]);
2378 nThru = 0;
2379 upperTheta = oldValue / alpha;
2380 whichColumn = iSequence;
2381 } else if (oldValue - upperTheta2 * alpha < 0.0) {
2382 costThru += alpha * (abcUpper_[iSequence] - abcLower_[iSequence]);
2383 index[nThru++] = iSequence;
2384 }
2385 } else if (iStatus < 6) {
2386 upperTheta = -1.0;
2387 upperTheta2 = elementByRow[j];
2388 whichColumn = iSequence;
2389 break;
2390 }
2391 }
2392 if (whichColumn < 0)
2393 continue;
2394 if (upperTheta != -1.0) {
2395 assert(upperTheta >= 0.0);
2396 if (costThru < infeasibility && whichColumn2 >= 0) {
2397 index[nThru++] = whichColumn;
2398 for (int i = 0; i < nThru; i++) {
2399 int iSequence = index[i];
2400 assert(abcUpper_[iSequence] - abcLower_[iSequence] < 1.0e10);
2401 // flip and use previous
2402 if (getInternalStatus(iSequence) == atLowerBound) {
2403 setInternalStatus(iSequence, atUpperBound);
2404 solutionSaved_[iSequence] = abcUpper_[iSequence];
2405 abcSolution_[iSequence] = abcUpper_[iSequence];
2406 } else {
2407 setInternalStatus(iSequence, atLowerBound);
2408 solutionSaved_[iSequence] = abcLower_[iSequence];
2409 abcSolution_[iSequence] = abcLower_[iSequence];
2410 }
2411 }
2412 whichColumn = whichColumn2;
2413 upperTheta = upperTheta2;
2414 }
2415 } else {
2416 // free coming in
2417 upperTheta = (abcDj_[whichColumn] * direction) / upperTheta2;
2418 }
2419 // update djs
2420 upperTheta *= -direction;
2421 for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2422 int iSequence = column[j];
2423 double dj = abcDj_[iSequence];
2424 double tableauValue = elementByRow[j];
2425 dj -= upperTheta * tableauValue;
2426 unsigned char iStatus = internalStatus_[iSequence] & 7;
2427 if ((iStatus & 4) == 0) {
2428 if (!iStatus) {
2429 assert(dj > -1.0e-2);
2430 dj = CoinMax(dj, 0.0);
2431 #ifdef ALLOW_BAD_DJS
2432 if (numberBad && modDj[iSequence]) {
2433 double bad = modDj[iSequence];
2434 if (dj + bad >= 0.0) {
2435 numberBad--;
2436 modDj[iSequence] = 0.0;
2437 dj += bad;
2438 } else {
2439 modDj[iSequence] += dj;
2440 dj = 0.0;
2441 }
2442 }
2443 #endif
2444 } else {
2445 assert(dj < 1.0e-2);
2446 dj = CoinMin(dj, 0.0);
2447 #ifdef ALLOW_BAD_DJS
2448 if (numberBad && modDj[iSequence]) {
2449 double bad = modDj[iSequence];
2450 if (dj + bad <= 0.0) {
2451 numberBad--;
2452 modDj[iSequence] = 0.0;
2453 dj += bad;
2454 } else {
2455 modDj[iSequence] += dj;
2456 dj = 0.0;
2457 }
2458 }
2459 #endif
2460 }
2461 } else if (iStatus < 6) {
2462 assert(fabs(dj) < 1.0e-4);
2463 dj = 0.0;
2464 }
2465 abcDj_[iSequence] = dj;
2466 }
2467 // do basis
2468 if (direction > 0.0) {
2469 if (upperBasic_[iRow] > lowerBasic_[iRow])
2470 setInternalStatus(iRow, atLowerBound);
2471 else
2472 setInternalStatus(iRow, isFixed);
2473 solutionSaved_[iRow] = abcLower_[iRow];
2474 abcSolution_[iRow] = abcLower_[iRow];
2475 } else {
2476 if (upperBasic_[iRow] > lowerBasic_[iRow])
2477 setInternalStatus(iRow, atUpperBound);
2478 else
2479 setInternalStatus(iRow, isFixed);
2480 solutionSaved_[iRow] = abcUpper_[iRow];
2481 abcSolution_[iRow] = abcUpper_[iRow];
2482 }
2483 setInternalStatus(whichColumn, basic);
2484 abcPivotVariable_[iRow] = whichColumn;
2485 nBasic++;
2486 // mark rows
2487 for (CoinBigIndex j = columnStart[whichColumn]; j < columnStart[whichColumn + 1]; j++) {
2488 int jRow = row[j];
2489 abcDj_[jRow] = 1.0;
2490 }
2491 }
2492 #ifdef ALLOW_BAD_DJS
2493 CoinAbcMemset0(modDj + maximumAbcNumberRows_, numberColumns_);
2494 #endif
2495 setAvailableArray(iVector);
2496 setAvailableArray(iVector2);
2497 #if ABC_NORMAL_DEBUG > 0
2498 printf("dual crash put %d in basis\n", nBasic);
2499 #endif
2500 } else {
2501 assert((stateOfProblem_ & VALUES_PASS2) != 0);
2502 // The idea is to put as many likely variables into basis as possible
2503 int n = 0;
2504 int iVector = getAvailableArray();
2505 int *index = usefulArray_[iVector].getIndices();
2506 double *array = usefulArray_[iVector].denseVector();
2507 int iVector2 = getAvailableArray();
2508 int *index2 = usefulArray_[iVector].getIndices();
2509 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2510 double dj = djSaved_[iSequence];
2511 double value = solution_[iSequence];
2512 double lower = abcLower_[iSequence];
2513 double upper = abcUpper_[iSequence];
2514 double gapUp = CoinMin(1.0e3, upper - value);
2515 assert(gapUp >= -1.0e-3);
2516 gapUp = CoinMax(gapUp, 0.0);
2517 double gapDown = CoinMin(1.0e3, value - lower);
2518 assert(gapDown >= -1.0e-3);
2519 gapDown = CoinMax(gapDown, 0.0);
2520 double measure = (CoinMin(gapUp, gapDown) + 1.0e-6) / (fabs(dj) + 1.0e-6);
2521 if (gapUp < primalTolerance_ * 10.0 && dj < dualTolerance_) {
2522 // set to ub
2523 setInternalStatus(iSequence, atUpperBound);
2524 solutionSaved_[iSequence] = abcUpper_[iSequence];
2525 abcSolution_[iSequence] = abcUpper_[iSequence];
2526 } else if (gapDown < primalTolerance_ * 10.0 && dj > -dualTolerance_) {
2527 // set to lb
2528 setInternalStatus(iSequence, atLowerBound);
2529 solutionSaved_[iSequence] = abcLower_[iSequence];
2530 abcSolution_[iSequence] = abcLower_[iSequence];
2531 } else if (upper > lower) {
2532 // set to nearest
2533 if (gapUp < gapDown) {
2534 // set to ub
2535 setInternalStatus(iSequence, atUpperBound);
2536 solutionSaved_[iSequence] = abcUpper_[iSequence];
2537 abcSolution_[iSequence] = abcUpper_[iSequence];
2538 } else {
2539 // set to lb
2540 setInternalStatus(iSequence, atLowerBound);
2541 solutionSaved_[iSequence] = abcLower_[iSequence];
2542 abcSolution_[iSequence] = abcLower_[iSequence];
2543 }
2544 array[n] = -measure;
2545 index[n++] = iSequence;
2546 }
2547 }
2548 // set slacks basic
2549 memset(internalStatus_, 6, numberRows_);
2550 CoinSort_2(solution_, solution_ + n, index);
2551 CoinAbcMemset0(array, n);
2552 for (int i = 0; i < numberRows_; i++)
2553 index2[i] = 0;
2554 // first put all possible slacks in
2555 int n2 = 0;
2556 for (int i = 0; i < n; i++) {
2557 int iSequence = index[i];
2558 if (iSequence < numberRows_) {
2559 index2[iSequence] = numberRows_;
2560 } else {
2561 index[n2++] = iSequence;
2562 }
2563 }
2564 n = n2;
2565 int numberIn = 0;
2566 const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts() - maximumAbcNumberRows_;
2567 //const int * columnLength = abcMatrix_->getPackedMatrix()->getVectorLengths();
2568 const int *row = abcMatrix_->getPackedMatrix()->getIndices();
2569 if (!abcMatrix_->gotRowCopy())
2570 abcMatrix_->createRowCopy();
2571 //const CoinBigIndex * rowStart = abcMatrix_->rowStart();
2572 //const CoinBigIndex * rowEnd = abcMatrix_->rowEnd();
2573 for (int i = 0; i < n; i++) {
2574 int iSequence = index[i];
2575 int bestRow = -1;
2576 int bestCount = numberRows_ + 1;
2577 for (CoinBigIndex j = columnStart[iSequence]; j < columnStart[iSequence + 1]; j++) {
2578 int iRow = row[j];
2579 if (!abcMatrix_->gotRowCopy())
2580 abcMatrix_->createRowCopy();
2581 const CoinBigIndex *rowStart = abcMatrix_->rowStart();
2582 const CoinBigIndex *rowEnd = abcMatrix_->rowEnd();
2583 if (!index2[iRow]) {
2584 int length = rowEnd[iRow] - rowStart[iRow];
2585 if (length < bestCount) {
2586 bestCount = length;
2587 bestRow = iRow;
2588 }
2589 }
2590 }
2591 if (bestRow >= 0) {
2592 numberIn++;
2593 for (CoinBigIndex j = columnStart[iSequence]; j < columnStart[iSequence + 1]; j++) {
2594 int iRow = row[j];
2595 index2[iRow]++;
2596 }
2597 setInternalStatus(iSequence, basic);
2598 abcPivotVariable_[bestRow] = iSequence;
2599 double dj = djSaved_[bestRow];
2600 double value = solution_[bestRow];
2601 double lower = abcLower_[bestRow];
2602 double upper = abcUpper_[bestRow];
2603 double gapUp = CoinMax(CoinMin(1.0e3, upper - value), 0.0);
2604 double gapDown = CoinMax(CoinMin(1.0e3, value - lower), 0.0);
2605 //double measure = (CoinMin(gapUp,gapDown)+1.0e-6)/(fabs(dj)+1.0e-6);
2606 if (gapUp < primalTolerance_ * 10.0 && dj < dualTolerance_) {
2607 // set to ub
2608 setInternalStatus(bestRow, atUpperBound);
2609 solutionSaved_[bestRow] = abcUpper_[bestRow];
2610 abcSolution_[bestRow] = abcUpper_[bestRow];
2611 } else if (gapDown < primalTolerance_ * 10.0 && dj > -dualTolerance_) {
2612 // set to lb
2613 setInternalStatus(bestRow, atLowerBound);
2614 solutionSaved_[bestRow] = abcLower_[bestRow];
2615 abcSolution_[bestRow] = abcLower_[bestRow];
2616 } else if (upper > lower) {
2617 // set to nearest
2618 if (gapUp < gapDown) {
2619 // set to ub
2620 setInternalStatus(bestRow, atUpperBound);
2621 solutionSaved_[bestRow] = abcUpper_[bestRow];
2622 abcSolution_[bestRow] = abcUpper_[bestRow];
2623 } else {
2624 // set to lb
2625 setInternalStatus(bestRow, atLowerBound);
2626 solutionSaved_[bestRow] = abcLower_[bestRow];
2627 abcSolution_[bestRow] = abcLower_[bestRow];
2628 }
2629 }
2630 }
2631 }
2632 #if ABC_NORMAL_DEBUG > 0
2633 printf("Idiot crash put %d in basis\n", numberIn);
2634 #endif
2635 setAvailableArray(iVector);
2636 setAvailableArray(iVector2);
2637 delete[] solution_;
2638 solution_ = NULL;
2639 }
2640 }
2641 /* Puts more stuff in basis
2642 1 bit set - do even if basis exists
2643 2 bit set - don't bother staying triangular
2644 */
putStuffInBasis(int type)2645 void AbcSimplex::putStuffInBasis(int type)
2646 {
2647 int i;
2648 for (i = 0; i < numberRows_; i++) {
2649 if (getInternalStatus(i) != basic)
2650 break;
2651 }
2652 if (i < numberRows_ && (type & 1) == 0)
2653 return;
2654 int iVector = getAvailableArray();
2655 // Column part is mustColumnIn
2656 int *COIN_RESTRICT mustRowOut = usefulArray_[iVector].getIndices();
2657 if (!abcMatrix_->gotRowCopy())
2658 abcMatrix_->createRowCopy();
2659 const double *elementByRow = abcMatrix_->rowElements();
2660 const CoinBigIndex *rowStart = abcMatrix_->rowStart();
2661 const CoinBigIndex *rowEnd = abcMatrix_->rowEnd();
2662 const int *column = abcMatrix_->rowColumns();
2663 for (int i = 0; i < numberTotal_; i++)
2664 mustRowOut[i] = -1;
2665 int nPossible = 0;
2666 // find equality rows with effective nonzero rhs
2667 for (int iRow = 0; iRow < numberRows_; iRow++) {
2668 if (abcUpper_[iRow] > abcLower_[iRow] || getInternalStatus(iRow) != basic) {
2669 mustRowOut[iRow] = -2;
2670 continue;
2671 }
2672 int chooseColumn[2] = { -1, -1 };
2673 for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2674 int iColumn = column[j];
2675 if (elementByRow[j] > 0.0) {
2676 if (chooseColumn[0] == -1)
2677 chooseColumn[0] = iColumn;
2678 else
2679 chooseColumn[0] = -2;
2680 } else {
2681 if (chooseColumn[1] == -1)
2682 chooseColumn[1] = iColumn;
2683 else
2684 chooseColumn[1] = -2;
2685 }
2686 }
2687 for (int iTry = 0; iTry < 2; iTry++) {
2688 int jColumn = chooseColumn[iTry];
2689 if (jColumn >= 0 && getInternalStatus(jColumn) != basic) {
2690 // see if has to be basic
2691 double lowerValue = -abcUpper_[iRow]; // check swap
2692 double upperValue = -abcLower_[iRow];
2693 int lowerInf = 0;
2694 int upperInf = 0;
2695 double alpha = 0.0;
2696 for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2697 int iColumn = column[j];
2698 if (iColumn != jColumn) {
2699 if (abcLower_[iColumn] > -1.0e20)
2700 lowerValue -= abcLower_[iColumn] * elementByRow[j];
2701 else
2702 lowerInf++;
2703 if (abcUpper_[iColumn] < 1.0e20)
2704 upperValue -= abcUpper_[iColumn] * elementByRow[j];
2705 else
2706 upperInf++;
2707 } else {
2708 alpha = elementByRow[j];
2709 }
2710 }
2711 // find values column must lie between (signs again)
2712 if (upperInf)
2713 upperValue = COIN_DBL_MAX;
2714 else
2715 upperValue /= alpha;
2716 if (lowerInf)
2717 lowerValue = -COIN_DBL_MAX;
2718 else
2719 lowerValue /= alpha;
2720 if (iTry) {
2721 // swap
2722 double temp = lowerValue;
2723 lowerValue = upperValue;
2724 upperValue = temp;
2725 }
2726 if (lowerValue > abcLower_[jColumn] + 10.0 * primalTolerance_ && upperValue < abcUpper_[jColumn] - 10.0 * primalTolerance_) {
2727 nPossible++;
2728 if (mustRowOut[jColumn] >= 0) {
2729 // choose one ???
2730 //printf("Column %d already selected on row %d now on %d\n",
2731 // jColumn,mustRowOut[jColumn],iRow);
2732 continue;
2733 }
2734 mustRowOut[jColumn] = iRow;
2735 mustRowOut[iRow] = jColumn;
2736 }
2737 }
2738 }
2739 }
2740 if (nPossible) {
2741 #if ABC_NORMAL_DEBUG > 0
2742 printf("%d possible candidates\n", nPossible);
2743 #endif
2744 if ((type & 2) == 0) {
2745 // triangular
2746 int iVector2 = getAvailableArray();
2747 int *COIN_RESTRICT counts = usefulArray_[iVector2].getIndices();
2748 CoinAbcMemset0(counts, numberRows_);
2749 for (int iRow = 0; iRow < numberRows_; iRow++) {
2750 int n = 0;
2751 for (CoinBigIndex j = rowStart[iRow]; j < rowEnd[iRow]; j++) {
2752 int iColumn = column[j];
2753 if (getInternalStatus(iColumn) == basic)
2754 n++;
2755 }
2756 counts[iRow] = n;
2757 }
2758 const CoinBigIndex *columnStart = abcMatrix_->getPackedMatrix()->getVectorStarts()
2759 - maximumAbcNumberRows_;
2760 const int *row = abcMatrix_->getPackedMatrix()->getIndices();
2761 for (int iRow = 0; iRow < numberRows_; iRow++) {
2762 if (!counts[iRow]) {
2763 int iColumn = mustRowOut[iRow];
2764 if (iColumn >= 0) {
2765 setInternalStatus(iRow, isFixed);
2766 solutionSaved_[iRow] = abcLower_[iRow];
2767 setInternalStatus(iColumn, basic);
2768 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
2769 int jRow = row[j];
2770 counts[jRow]++;
2771 }
2772 }
2773 }
2774 }
2775 setAvailableArray(iVector2);
2776 } else {
2777 // all
2778 for (int iRow = 0; iRow < numberRows_; iRow++) {
2779 int iColumn = mustRowOut[iRow];
2780 if (iColumn >= 0) {
2781 setInternalStatus(iRow, isFixed);
2782 solutionSaved_[iRow] = abcLower_[iRow];
2783 setInternalStatus(iColumn, basic);
2784 }
2785 }
2786 }
2787 // redo pivot array
2788 int numberBasic = 0;
2789 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2790 if (getInternalStatus(iSequence) == basic)
2791 abcPivotVariable_[numberBasic++] = iSequence;
2792 }
2793 assert(numberBasic == numberRows_);
2794 }
2795 setAvailableArray(iVector);
2796 }
2797 // Computes nonbasic cost and total cost
computeObjective()2798 void AbcSimplex::computeObjective()
2799 {
2800 sumNonBasicCosts_ = 0.0;
2801 rawObjectiveValue_ = 0.0;
2802 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2803 double value = abcSolution_[iSequence] * abcCost_[iSequence];
2804 rawObjectiveValue_ += value;
2805 if (getInternalStatus(iSequence) != basic)
2806 sumNonBasicCosts_ += value;
2807 }
2808 setClpSimplexObjectiveValue();
2809 }
2810 // This sets largest infeasibility and most infeasible
checkPrimalSolution(bool justBasic)2811 void AbcSimplex::checkPrimalSolution(bool justBasic)
2812 {
2813 rawObjectiveValue_ = CoinAbcInnerProduct(costBasic_, numberRows_, solutionBasic_);
2814 //rawObjectiveValue_ += sumNonBasicCosts_;
2815 setClpSimplexObjectiveValue();
2816 // now look at primal solution
2817 sumPrimalInfeasibilities_ = 0.0;
2818 numberPrimalInfeasibilities_ = 0;
2819 double primalTolerance = primalTolerance_;
2820 double relaxedTolerance = primalTolerance_;
2821 // we can't really trust infeasibilities if there is primal error
2822 double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_, 5.0 * primalTolerance_));
2823 // allow tolerance at least slightly bigger than standard
2824 relaxedTolerance = relaxedTolerance + error;
2825 sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2826 const double *COIN_RESTRICT lowerBasic = lowerBasic_;
2827 const double *COIN_RESTRICT upperBasic = upperBasic_;
2828 const double *COIN_RESTRICT solutionBasic = solutionBasic_;
2829 if (justBasic) {
2830 for (int iRow = 0; iRow < numberRows_; iRow++) {
2831 double infeasibility = 0.0;
2832 if (solutionBasic[iRow] > upperBasic[iRow]) {
2833 infeasibility = solutionBasic[iRow] - upperBasic[iRow];
2834 } else if (solutionBasic[iRow] < lowerBasic[iRow]) {
2835 infeasibility = lowerBasic[iRow] - solutionBasic[iRow];
2836 }
2837 if (infeasibility > primalTolerance) {
2838 sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2839 if (infeasibility > relaxedTolerance)
2840 sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2841 numberPrimalInfeasibilities_++;
2842 }
2843 }
2844 } else {
2845 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2846 double infeasibility = 0.0;
2847 if (abcSolution_[iSequence] > abcUpper_[iSequence]) {
2848 infeasibility = abcSolution_[iSequence] - abcUpper_[iSequence];
2849 } else if (abcSolution_[iSequence] < abcLower_[iSequence]) {
2850 infeasibility = abcLower_[iSequence] - abcSolution_[iSequence];
2851 }
2852 if (infeasibility > primalTolerance) {
2853 //assert (getInternalStatus(iSequence)==basic);
2854 sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2855 if (infeasibility > relaxedTolerance)
2856 sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2857 numberPrimalInfeasibilities_++;
2858 }
2859 }
2860 }
2861 }
checkDualSolution()2862 void AbcSimplex::checkDualSolution()
2863 {
2864
2865 sumDualInfeasibilities_ = 0.0;
2866 numberDualInfeasibilities_ = 0;
2867 int firstFreePrimal = -1;
2868 int firstFreeDual = -1;
2869 int numberSuperBasicWithDj = 0;
2870 bestPossibleImprovement_ = 0.0;
2871 // we can't really trust infeasibilities if there is dual error
2872 double error = CoinMin(1.0e-2, CoinMax(largestDualError_, 5.0 * dualTolerance_));
2873 // allow tolerance at least slightly bigger than standard
2874 double relaxedTolerance = currentDualTolerance_ + error;
2875 // allow bigger tolerance for possible improvement
2876 double possTolerance = 5.0 * relaxedTolerance;
2877 sumOfRelaxedDualInfeasibilities_ = 0.0;
2878 int numberNeedToMove = 0;
2879 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
2880 if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) {
2881 // not basic
2882 double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence];
2883 double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence];
2884 double value = abcDj_[iSequence];
2885 if (distanceUp > primalTolerance_) {
2886 // Check if "free"
2887 if (distanceDown <= primalTolerance_) {
2888 // should not be negative
2889 if (value < 0.0) {
2890 value = -value;
2891 if (value > currentDualTolerance_) {
2892 sumDualInfeasibilities_ += value - currentDualTolerance_;
2893 if (value > possTolerance)
2894 bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2895 if (value > relaxedTolerance)
2896 sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2897 numberDualInfeasibilities_++;
2898 }
2899 }
2900 } else {
2901 // free
2902 value = fabs(value);
2903 if (value > 1.0 * relaxedTolerance) {
2904 numberSuperBasicWithDj++;
2905 if (firstFreeDual < 0)
2906 firstFreeDual = iSequence;
2907 }
2908 if (value > currentDualTolerance_ || getInternalStatus(iSequence) != AbcSimplex::isFree) {
2909 numberNeedToMove++;
2910 if (firstFreePrimal < 0)
2911 firstFreePrimal = iSequence;
2912 sumDualInfeasibilities_ += value - currentDualTolerance_;
2913 if (value > possTolerance)
2914 bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2915 if (value > relaxedTolerance)
2916 sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2917 numberDualInfeasibilities_++;
2918 }
2919 }
2920 } else if (distanceDown > primalTolerance_) {
2921 // should not be positive
2922 if (value > 0.0) {
2923 if (value > currentDualTolerance_) {
2924 sumDualInfeasibilities_ += value - currentDualTolerance_;
2925 if (value > possTolerance)
2926 bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
2927 if (value > relaxedTolerance)
2928 sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2929 numberDualInfeasibilities_++;
2930 }
2931 }
2932 }
2933 }
2934 }
2935 numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ - numberSuperBasicWithDj;
2936 if (algorithm_ < 0 && firstFreeDual >= 0) {
2937 // dual
2938 firstFree_ = firstFreeDual;
2939 } else if (numberSuperBasicWithDj || numberNeedToMove) {
2940 //(abcProgress_.lastIterationNumber(0) <= 0)) {
2941 firstFree_ = firstFreePrimal;
2942 if (!sumDualInfeasibilities_)
2943 sumDualInfeasibilities_ = 1.0e-8;
2944 }
2945 }
2946 /* This sets sum and number of infeasibilities (Dual and Primal) */
checkBothSolutions()2947 void AbcSimplex::checkBothSolutions()
2948 {
2949 // old way
2950 checkPrimalSolution(true);
2951 checkDualSolution();
2952 }
2953 /*
2954 Unpacks one column of the matrix into indexed array
2955 */
unpack(CoinIndexedVector & rowArray,int sequence) const2956 void AbcSimplex::unpack(CoinIndexedVector &rowArray, int sequence) const
2957 {
2958 abcMatrix_->unpack(rowArray, sequence);
2959 }
2960 // Sets row pivot choice algorithm in dual
setDualRowPivotAlgorithm(AbcDualRowPivot & choice)2961 void AbcSimplex::setDualRowPivotAlgorithm(AbcDualRowPivot &choice)
2962 {
2963 delete abcDualRowPivot_;
2964 abcDualRowPivot_ = choice.clone(true);
2965 abcDualRowPivot_->setModel(this);
2966 }
2967 // Sets column pivot choice algorithm in primal
setPrimalColumnPivotAlgorithm(AbcPrimalColumnPivot & choice)2968 void AbcSimplex::setPrimalColumnPivotAlgorithm(AbcPrimalColumnPivot &choice)
2969 {
2970 delete abcPrimalColumnPivot_;
2971 abcPrimalColumnPivot_ = choice.clone(true);
2972 abcPrimalColumnPivot_->setModel(this);
2973 }
setFactorization(AbcSimplexFactorization & factorization)2974 void AbcSimplex::setFactorization(AbcSimplexFactorization &factorization)
2975 {
2976 if (abcFactorization_)
2977 abcFactorization_->setFactorization(factorization);
2978 else
2979 abcFactorization_ = new AbcSimplexFactorization(factorization,
2980 numberRows_);
2981 }
2982
2983 // Swaps factorization
2984 AbcSimplexFactorization *
swapFactorization(AbcSimplexFactorization * factorization)2985 AbcSimplex::swapFactorization(AbcSimplexFactorization *factorization)
2986 {
2987 AbcSimplexFactorization *swap = abcFactorization_;
2988 abcFactorization_ = factorization;
2989 return swap;
2990 }
2991
2992 /* Tightens primal bounds to make dual faster. Unless
2993 fixed, bounds are slightly looser than they could be.
2994 This is to make dual go faster and is probably not needed
2995 with a presolve. Returns non-zero if problem infeasible
2996 */
tightenPrimalBounds()2997 int AbcSimplex::tightenPrimalBounds()
2998 {
2999 int tightenType = 1;
3000 if (maximumIterations() >= 1000000 && maximumIterations() < 1000010)
3001 tightenType = maximumIterations() - 1000000;
3002 if (!tightenType)
3003 return 0;
3004 if (integerType_) {
3005 #if ABC_NORMAL_DEBUG > 0
3006 printf("Redo tighten to relax by 1 for integers (and don't be shocked by infeasibility)\n");
3007 #endif
3008 return 0;
3009 }
3010 // This needs to work on scaled matrix - then replace
3011 // Get a row copy in standard format
3012 CoinPackedMatrix copy;
3013 copy.setExtraGap(0.0);
3014 copy.setExtraMajor(0.0);
3015 copy.reverseOrderedCopyOf(*abcMatrix_->matrix());
3016 // get matrix data pointers
3017 const int *COIN_RESTRICT column = copy.getIndices();
3018 const CoinBigIndex *COIN_RESTRICT rowStart = copy.getVectorStarts();
3019 const int *COIN_RESTRICT rowLength = copy.getVectorLengths();
3020 const double *COIN_RESTRICT element = copy.getElements();
3021 int numberChanged = 1, iPass = 0;
3022 double large = largeValue(); // treat bounds > this as infinite
3023 #ifndef NDEBUG
3024 double large2 = 1.0e10 * large;
3025 #endif
3026 int numberInfeasible = 0;
3027 int totalTightened = 0;
3028
3029 double tolerance = primalTolerance();
3030
3031 // A copy of bounds is up at top
3032 translate(0); // move down
3033
3034 #define MAXPASS 10
3035
3036 // Loop round seeing if we can tighten bounds
3037 // Would be faster to have a stack of possible rows
3038 // and we put altered rows back on stack
3039 int numberCheck = -1;
3040 // temp swap signs so can use existing code
3041 double *COIN_RESTRICT rowLower = abcLower_;
3042 double *COIN_RESTRICT rowUpper = abcUpper_;
3043 for (int iRow = 0; iRow < numberRows_; iRow++) {
3044 double value = -rowLower[iRow];
3045 rowLower[iRow] = -rowUpper[iRow];
3046 rowUpper[iRow] = value;
3047 }
3048 // Keep lower and upper for rows
3049 int whichArray[5];
3050 for (int i = 0; i < 5; i++)
3051 whichArray[i] = getAvailableArray();
3052 double *COIN_RESTRICT minRhs = usefulArray_[whichArray[0]].denseVector();
3053 double *COIN_RESTRICT maxRhs = usefulArray_[whichArray[1]].denseVector();
3054 int *COIN_RESTRICT changedList = usefulArray_[whichArray[0]].getIndices();
3055 int *COIN_RESTRICT changedListNext = usefulArray_[whichArray[1]].getIndices();
3056 char *COIN_RESTRICT changed = reinterpret_cast< char * >(usefulArray_[whichArray[2]].getIndices());
3057 // -1 no infinite, -2 more than one infinite >=0 column
3058 int *COIN_RESTRICT whichInfiniteDown = usefulArray_[whichArray[3]].getIndices();
3059 int *COIN_RESTRICT whichInfiniteUp = usefulArray_[whichArray[3]].getIndices();
3060 int numberToDoNext = 0;
3061 for (int iRow = 0; iRow < numberRows_; iRow++) {
3062 if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
3063 changedListNext[numberToDoNext++] = iRow;
3064 whichInfiniteDown[iRow] = -1;
3065 whichInfiniteUp[iRow] = -1;
3066 } else {
3067 minRhs[iRow] = -COIN_DBL_MAX;
3068 maxRhs[iRow] = COIN_DBL_MAX;
3069 whichInfiniteDown[iRow] = -2;
3070 whichInfiniteUp[iRow] = -2;
3071 }
3072 }
3073 const int *COIN_RESTRICT row = abcMatrix_->matrix()->getIndices();
3074 const CoinBigIndex *COIN_RESTRICT columnStart = abcMatrix_->matrix()->getVectorStarts();
3075 const double *COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements();
3076
3077 double *COIN_RESTRICT columnLower = abcLower_ + maximumAbcNumberRows_;
3078 double *COIN_RESTRICT columnUpper = abcUpper_ + maximumAbcNumberRows_;
3079 while (numberChanged > numberCheck) {
3080 int numberToDo = numberToDoNext;
3081 numberToDoNext = 0;
3082 int *COIN_RESTRICT temp = changedList;
3083 changedList = changedListNext;
3084 changedListNext = temp;
3085 CoinAbcMemset0(changed, numberRows_);
3086
3087 numberChanged = 0; // Bounds tightened this pass
3088
3089 if (iPass == MAXPASS)
3090 break;
3091 iPass++;
3092 for (int k = 0; k < numberToDo; k++) {
3093 int iRow = changedList[k];
3094 // possible row
3095 int infiniteUpper = 0;
3096 int infiniteLower = 0;
3097 int firstInfiniteUpper = -1;
3098 int firstInfiniteLower = -1;
3099 double maximumUp = 0.0;
3100 double maximumDown = 0.0;
3101 double newBound;
3102 CoinBigIndex rStart = rowStart[iRow];
3103 CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
3104 CoinBigIndex j;
3105 // Compute possible lower and upper ranges
3106
3107 for (j = rStart; j < rEnd; ++j) {
3108 double value = element[j];
3109 int iColumn = column[j];
3110 if (value > 0.0) {
3111 if (columnUpper[iColumn] >= large) {
3112 firstInfiniteUpper = iColumn;
3113 ++infiniteUpper;
3114 } else {
3115 maximumUp += columnUpper[iColumn] * value;
3116 }
3117 if (columnLower[iColumn] <= -large) {
3118 firstInfiniteLower = iColumn;
3119 ++infiniteLower;
3120 } else {
3121 maximumDown += columnLower[iColumn] * value;
3122 }
3123 } else if (value < 0.0) {
3124 if (columnUpper[iColumn] >= large) {
3125 firstInfiniteLower = iColumn;
3126 ++infiniteLower;
3127 } else {
3128 maximumDown += columnUpper[iColumn] * value;
3129 }
3130 if (columnLower[iColumn] <= -large) {
3131 firstInfiniteUpper = iColumn;
3132 ++infiniteUpper;
3133 } else {
3134 maximumUp += columnLower[iColumn] * value;
3135 }
3136 }
3137 }
3138 // Build in a margin of error
3139 maximumUp += 1.0e-8 * fabs(maximumUp);
3140 maximumDown -= 1.0e-8 * fabs(maximumDown);
3141 double maxUp = maximumUp + infiniteUpper * 1.0e31;
3142 double maxDown = maximumDown - infiniteLower * 1.0e31;
3143 minRhs[iRow] = infiniteLower ? -COIN_DBL_MAX : maximumDown;
3144 maxRhs[iRow] = infiniteUpper ? COIN_DBL_MAX : maximumUp;
3145 if (infiniteLower == 0)
3146 whichInfiniteDown[iRow] = -1;
3147 else if (infiniteLower == 1)
3148 whichInfiniteDown[iRow] = firstInfiniteLower;
3149 else
3150 whichInfiniteDown[iRow] = -2;
3151 if (infiniteUpper == 0)
3152 whichInfiniteUp[iRow] = -1;
3153 else if (infiniteUpper == 1)
3154 whichInfiniteUp[iRow] = firstInfiniteUpper;
3155 else
3156 whichInfiniteUp[iRow] = -2;
3157 if (maxUp <= rowUpper[iRow] + tolerance && maxDown >= rowLower[iRow] - tolerance) {
3158
3159 // Row is redundant - make totally free
3160 // NO - can't do this for postsolve
3161 // rowLower[iRow]=-COIN_DBL_MAX;
3162 // rowUpper[iRow]=COIN_DBL_MAX;
3163 //printf("Redundant row in presolveX %d\n",iRow);
3164
3165 } else {
3166 if (maxUp < rowLower[iRow] - 100.0 * tolerance || maxDown > rowUpper[iRow] + 100.0 * tolerance) {
3167 // problem is infeasible - exit at once
3168 numberInfeasible++;
3169 break;
3170 }
3171 double lower = rowLower[iRow];
3172 double upper = rowUpper[iRow];
3173 for (j = rStart; j < rEnd; ++j) {
3174 double value = element[j];
3175 int iColumn = column[j];
3176 double nowLower = columnLower[iColumn];
3177 double nowUpper = columnUpper[iColumn];
3178 if (value > 0.0) {
3179 // positive value
3180 if (lower > -large) {
3181 if (!infiniteUpper) {
3182 assert(nowUpper < large2);
3183 newBound = nowUpper + (lower - maximumUp) / value;
3184 // relax if original was large
3185 if (fabs(maximumUp) > 1.0e8)
3186 newBound -= 1.0e-12 * fabs(maximumUp);
3187 } else if (infiniteUpper == 1 && nowUpper > large) {
3188 newBound = (lower - maximumUp) / value;
3189 // relax if original was large
3190 if (fabs(maximumUp) > 1.0e8)
3191 newBound -= 1.0e-12 * fabs(maximumUp);
3192 } else {
3193 newBound = -COIN_DBL_MAX;
3194 }
3195 if (newBound > nowLower + 1.0e-12 && newBound > -large) {
3196 // Tighten the lower bound
3197 numberChanged++;
3198 // check infeasible (relaxed)
3199 if (nowUpper < newBound) {
3200 if (nowUpper - newBound < -100.0 * tolerance)
3201 numberInfeasible++;
3202 else
3203 newBound = nowUpper;
3204 }
3205 columnLower[iColumn] = newBound;
3206 for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) {
3207 int jRow = row[j1];
3208 if (!changed[jRow]) {
3209 changed[jRow] = 1;
3210 changedListNext[numberToDoNext++] = jRow;
3211 }
3212 }
3213 // adjust
3214 double now;
3215 if (nowLower < -large) {
3216 now = 0.0;
3217 infiniteLower--;
3218 } else {
3219 now = nowLower;
3220 }
3221 maximumDown += (newBound - now) * value;
3222 nowLower = newBound;
3223 }
3224 }
3225 if (upper < large) {
3226 if (!infiniteLower) {
3227 assert(nowLower > -large2);
3228 newBound = nowLower + (upper - maximumDown) / value;
3229 // relax if original was large
3230 if (fabs(maximumDown) > 1.0e8)
3231 newBound += 1.0e-12 * fabs(maximumDown);
3232 } else if (infiniteLower == 1 && nowLower < -large) {
3233 newBound = (upper - maximumDown) / value;
3234 // relax if original was large
3235 if (fabs(maximumDown) > 1.0e8)
3236 newBound += 1.0e-12 * fabs(maximumDown);
3237 } else {
3238 newBound = COIN_DBL_MAX;
3239 }
3240 if (newBound < nowUpper - 1.0e-12 && newBound < large) {
3241 // Tighten the upper bound
3242 numberChanged++;
3243 // check infeasible (relaxed)
3244 if (nowLower > newBound) {
3245 if (newBound - nowLower < -100.0 * tolerance)
3246 numberInfeasible++;
3247 else
3248 newBound = nowLower;
3249 }
3250 columnUpper[iColumn] = newBound;
3251 for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) {
3252 int jRow = row[j1];
3253 if (!changed[jRow]) {
3254 changed[jRow] = 1;
3255 changedListNext[numberToDoNext++] = jRow;
3256 }
3257 }
3258 // adjust
3259 double now;
3260 if (nowUpper > large) {
3261 now = 0.0;
3262 infiniteUpper--;
3263 } else {
3264 now = nowUpper;
3265 }
3266 maximumUp += (newBound - now) * value;
3267 nowUpper = newBound;
3268 }
3269 }
3270 } else {
3271 // negative value
3272 if (lower > -large) {
3273 if (!infiniteUpper) {
3274 assert(nowLower < large2);
3275 newBound = nowLower + (lower - maximumUp) / value;
3276 // relax if original was large
3277 if (fabs(maximumUp) > 1.0e8)
3278 newBound += 1.0e-12 * fabs(maximumUp);
3279 } else if (infiniteUpper == 1 && nowLower < -large) {
3280 newBound = (lower - maximumUp) / value;
3281 // relax if original was large
3282 if (fabs(maximumUp) > 1.0e8)
3283 newBound += 1.0e-12 * fabs(maximumUp);
3284 } else {
3285 newBound = COIN_DBL_MAX;
3286 }
3287 if (newBound < nowUpper - 1.0e-12 && newBound < large) {
3288 // Tighten the upper bound
3289 numberChanged++;
3290 // check infeasible (relaxed)
3291 if (nowLower > newBound) {
3292 if (newBound - nowLower < -100.0 * tolerance)
3293 numberInfeasible++;
3294 else
3295 newBound = nowLower;
3296 }
3297 columnUpper[iColumn] = newBound;
3298 for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) {
3299 int jRow = row[j1];
3300 if (!changed[jRow]) {
3301 changed[jRow] = 1;
3302 changedListNext[numberToDoNext++] = jRow;
3303 }
3304 }
3305 // adjust
3306 double now;
3307 if (nowUpper > large) {
3308 now = 0.0;
3309 infiniteLower--;
3310 } else {
3311 now = nowUpper;
3312 }
3313 maximumDown += (newBound - now) * value;
3314 nowUpper = newBound;
3315 }
3316 }
3317 if (upper < large) {
3318 if (!infiniteLower) {
3319 assert(nowUpper < large2);
3320 newBound = nowUpper + (upper - maximumDown) / value;
3321 // relax if original was large
3322 if (fabs(maximumDown) > 1.0e8)
3323 newBound -= 1.0e-12 * fabs(maximumDown);
3324 } else if (infiniteLower == 1 && nowUpper > large) {
3325 newBound = (upper - maximumDown) / value;
3326 // relax if original was large
3327 if (fabs(maximumDown) > 1.0e8)
3328 newBound -= 1.0e-12 * fabs(maximumDown);
3329 } else {
3330 newBound = -COIN_DBL_MAX;
3331 }
3332 if (newBound > nowLower + 1.0e-12 && newBound > -large) {
3333 // Tighten the lower bound
3334 numberChanged++;
3335 // check infeasible (relaxed)
3336 if (nowUpper < newBound) {
3337 if (nowUpper - newBound < -100.0 * tolerance)
3338 numberInfeasible++;
3339 else
3340 newBound = nowUpper;
3341 }
3342 columnLower[iColumn] = newBound;
3343 for (CoinBigIndex j1 = columnStart[iColumn]; j1 < columnStart[iColumn + 1]; j1++) {
3344 int jRow = row[j1];
3345 if (!changed[jRow]) {
3346 changed[jRow] = 1;
3347 changedListNext[numberToDoNext++] = jRow;
3348 }
3349 }
3350 // adjust
3351 double now;
3352 if (nowLower < -large) {
3353 now = 0.0;
3354 infiniteUpper--;
3355 } else {
3356 now = nowLower;
3357 }
3358 maximumUp += (newBound - now) * value;
3359 nowLower = newBound;
3360 }
3361 }
3362 }
3363 }
3364 }
3365 }
3366 #if 0
3367 // see if we can tighten doubletons with infinite bounds
3368 if (iPass==1) {
3369 const double * COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements();
3370 for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
3371 CoinBigIndex start = columnStart[iColumn];
3372 if (start+2==columnStart[iColumn+1]&&
3373 (columnLower[iColumn]<-1.0e30||columnUpper[iColumn])) {
3374 int iRow0=row[start];
3375 int iRow1=row[start+1];
3376 printf("col %d row0 %d el0 %g whichDown0 %d whichUp0 %d row1 %d el1 %g whichDown1 %d whichUp1 %d\n",iColumn,
3377 iRow0,elementByColumn[start],whichInfiniteDown[iRow0],whichInfiniteUp[iRow0],
3378 iRow1,elementByColumn[start+1],whichInfiniteDown[iRow1],whichInfiniteUp[iRow1]);
3379 }
3380 }
3381 }
3382 #endif
3383 totalTightened += numberChanged;
3384 if (iPass == 1)
3385 numberCheck = numberChanged >> 4;
3386 if (numberInfeasible)
3387 break;
3388 }
3389 const double *COIN_RESTRICT saveLower = abcLower_ + maximumNumberTotal_ + maximumAbcNumberRows_;
3390 const double *COIN_RESTRICT saveUpper = abcUpper_ + maximumNumberTotal_ + maximumAbcNumberRows_;
3391 if (!numberInfeasible) {
3392 handler_->message(CLP_SIMPLEX_BOUNDTIGHTEN, messages_)
3393 << totalTightened
3394 << CoinMessageEol;
3395 int numberLowerChanged = 0;
3396 int numberUpperChanged = 0;
3397 if (!integerType_) {
3398 // Set bounds slightly loose
3399 // tighten rows as well
3400 //#define PRINT_TIGHTEN_PROGRESS 0
3401 for (int iRow = 0; iRow < numberRows_; iRow++) {
3402 assert(maxRhs[iRow] >= rowLower[iRow]);
3403 assert(minRhs[iRow] <= rowUpper[iRow]);
3404 rowLower[iRow] = CoinMax(rowLower[iRow], minRhs[iRow]);
3405 rowUpper[iRow] = CoinMin(rowUpper[iRow], maxRhs[iRow]);
3406 #if PRINT_TIGHTEN_PROGRESS > 3
3407 if (handler_->logLevel() > 5)
3408 printf("Row %d min %g max %g lower %g upper %g\n",
3409 iRow, minRhs[iRow], maxRhs[iRow], rowLower[iRow], rowUpper[iRow]);
3410 #endif
3411 }
3412 #if 0
3413 double useTolerance = 1.0e-4;
3414 double multiplier = 100.0;
3415 // To do this we need to compute min and max JUST for no costs?
3416 const double * COIN_RESTRICT cost = abcCost_+maximumAbcNumberRows_;
3417 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3418 if (saveUpper[iColumn] > saveLower[iColumn] + useTolerance) {
3419 double awayFromLower=0.0;
3420 double awayFromUpper=0.0;
3421 //double gap=columnUpper[iColumn]-columnLower[iColumn];
3422 for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
3423 int iRow=row[j];
3424 double value=elementByColumn[j];
3425 #if PRINT_TIGHTEN_PROGRESS > 3
3426 if (handler_->logLevel()>5)
3427 printf("row %d element %g\n",iRow,value);
3428 #endif
3429 if (value>0.0) {
3430 if (minRhs[iRow]+value*awayFromLower<rowLower[iRow]) {
3431 if (minRhs[iRow]>-1.0e10&&awayFromLower<1.0e10)
3432 awayFromLower = CoinMax(awayFromLower,(rowLower[iRow]-minRhs[iRow])/value);
3433 else
3434 awayFromLower=COIN_DBL_MAX;
3435 }
3436 if (maxRhs[iRow]-value*awayFromUpper>rowUpper[iRow]) {
3437 if (maxRhs[iRow]<1.0e10&&awayFromUpper<1.0e10)
3438 awayFromUpper = CoinMax(awayFromUpper,(maxRhs[iRow]-rowUpper[iRow])/value);
3439 else
3440 awayFromUpper=COIN_DBL_MAX;
3441 }
3442 } else {
3443 if (maxRhs[iRow]+value*awayFromLower>rowUpper[iRow]) {
3444 if (maxRhs[iRow]<1.0e10&&awayFromLower<1.0e10)
3445 awayFromLower = CoinMax(awayFromLower,(rowUpper[iRow]-maxRhs[iRow])/value);
3446 else
3447 awayFromLower=COIN_DBL_MAX;
3448 }
3449 if (minRhs[iRow]-value*awayFromUpper<rowLower[iRow]) {
3450 if (minRhs[iRow]>-1.0e10&&awayFromUpper<1.0e10)
3451 awayFromUpper = CoinMin(awayFromUpper,(minRhs[iRow]-rowLower[iRow])/value);
3452 else
3453 awayFromUpper=COIN_DBL_MAX;
3454 }
3455 }
3456 }
3457 // Might have to go as high as
3458 double upper = CoinMin(columnLower[iColumn]+awayFromLower,columnUpper[iColumn]);
3459 // and as low as
3460 double lower = CoinMax(columnUpper[iColumn]-awayFromUpper,columnLower[iColumn]);
3461 // but be sensible
3462 double gap=0.999*(CoinMin(columnUpper[iColumn]-columnLower[iColumn],1.0e10));
3463 if (awayFromLower>gap||awayFromUpper>gap) {
3464 if (fabs(columnUpper[iColumn]-upper)>1.0e-5) {
3465 #if PRINT_TIGHTEN_PROGRESS > 1
3466 printf("YYchange upper bound on %d from %g (original %g) to %g (awayFromLower %g) - cost %g\n",iColumn,
3467 columnUpper[iColumn],saveUpper[iColumn],upper,awayFromLower,cost[iColumn]);
3468 #endif
3469 upper=columnUpper[iColumn];
3470 }
3471 if (fabs(columnLower[iColumn]-lower)>1.0e-5) {
3472 #if PRINT_TIGHTEN_PROGRESS > 1
3473 printf("YYchange lower bound on %d from %g (original %g) to %g (awayFromUpper %g) - cost %g\n",iColumn,
3474 columnLower[iColumn],saveLower[iColumn],lower,awayFromUpper,cost[iColumn]);
3475 #endif
3476 lower=columnLower[iColumn];
3477 }
3478 }
3479 if (lower>upper) {
3480 #if PRINT_TIGHTEN_PROGRESS
3481 printf("XXchange upper bound on %d from %g (original %g) to %g (awayFromLower %g) - cost %g\n",iColumn,
3482 columnUpper[iColumn],saveUpper[iColumn],upper,awayFromLower,cost[iColumn]);
3483 printf("XXchange lower bound on %d from %g (original %g) to %g (awayFromUpper %g) - cost %g\n",iColumn,
3484 columnLower[iColumn],saveLower[iColumn],lower,awayFromUpper,cost[iColumn]);
3485 #endif
3486 lower=columnLower[iColumn];
3487 upper=columnUpper[iColumn];
3488 }
3489 //upper=CoinMax(upper,0.0);
3490 //upper=CoinMax(upper,0.0);
3491 if (cost[iColumn]>=0.0&&awayFromLower<1.0e10&&columnLower[iColumn]>-1.0e10) {
3492 // doesn't want to be too positive
3493 if (fabs(columnUpper[iColumn]-upper)>1.0e-5) {
3494 #if PRINT_TIGHTEN_PROGRESS > 2
3495 if (handler_->logLevel()>1)
3496 printf("change upper bound on %d from %g (original %g) to %g (awayFromLower %g) - cost %g\n",iColumn,
3497 columnUpper[iColumn],saveUpper[iColumn],upper,awayFromLower,cost[iColumn]);
3498 #endif
3499 double difference=columnUpper[iColumn]-upper;
3500 for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
3501 int iRow=row[j];
3502 double value=elementByColumn[j];
3503 if (value>0.0) {
3504 assert (difference*value>=0.0);
3505 maxRhs[iRow]-=difference*value;
3506 } else {
3507 assert (difference*value<=0.0);
3508 minRhs[iRow]-=difference*value;
3509 }
3510 }
3511 columnUpper[iColumn]=upper;
3512 numberUpperChanged++;
3513 }
3514 }
3515 if (cost[iColumn]<=0.0&&awayFromUpper<1.0e10&&columnUpper[iColumn]<1.0e10) {
3516 // doesn't want to be too negative
3517 if (fabs(columnLower[iColumn]-lower)>1.0e-5) {
3518 #if PRINT_TIGHTEN_PROGRESS > 2
3519 if (handler_->logLevel()>1)
3520 printf("change lower bound on %d from %g (original %g) to %g (awayFromUpper %g) - cost %g\n",iColumn,
3521 columnLower[iColumn],saveLower[iColumn],lower,awayFromUpper,cost[iColumn]);
3522 #endif
3523 double difference=columnLower[iColumn]-lower;
3524 for (CoinBigIndex j=columnStart[iColumn];j<columnStart[iColumn+1];j++) {
3525 int iRow=row[j];
3526 double value=elementByColumn[j];
3527 if (value>0.0) {
3528 assert (difference*value<=0.0);
3529 minRhs[iRow]-=difference*value;
3530 } else {
3531 assert (difference*value>=0.0);
3532 maxRhs[iRow]-=difference*value;
3533 }
3534 }
3535 columnLower[iColumn]=lower;
3536 numberLowerChanged++;
3537 }
3538 }
3539 // Make large bounds stay infinite
3540 if (saveUpper[iColumn] > 1.0e30 && columnUpper[iColumn] > 1.0e10) {
3541 columnUpper[iColumn] = COIN_DBL_MAX;
3542 }
3543 if (saveLower[iColumn] < -1.0e30 && columnLower[iColumn] < -1.0e10) {
3544 columnLower[iColumn] = -COIN_DBL_MAX;
3545 }
3546 if (columnUpper[iColumn] - columnLower[iColumn] < useTolerance + 1.0e-8) {
3547 // relax enough so will have correct dj
3548 columnLower[iColumn] = CoinMax(saveLower[iColumn],
3549 columnLower[iColumn] - multiplier * useTolerance);
3550 columnUpper[iColumn] = CoinMin(saveUpper[iColumn],
3551 columnUpper[iColumn] + multiplier * useTolerance);
3552 } else {
3553 if (columnUpper[iColumn] < saveUpper[iColumn]) {
3554 // relax a bit
3555 columnUpper[iColumn] = CoinMin(columnUpper[iColumn] + multiplier * useTolerance,
3556 saveUpper[iColumn]);
3557 }
3558 if (columnLower[iColumn] > saveLower[iColumn]) {
3559 // relax a bit
3560 columnLower[iColumn] = CoinMax(columnLower[iColumn] - multiplier * useTolerance,
3561 saveLower[iColumn]);
3562 }
3563 }
3564 if (0) {
3565 // debug
3566 // tighten rows as well
3567 for (int iRow = 0; iRow < numberRows_; iRow++) {
3568 if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
3569 // possible row
3570 int infiniteUpper = 0;
3571 int infiniteLower = 0;
3572 double maximumUp = 0.0;
3573 double maximumDown = 0.0;
3574 CoinBigIndex rStart = rowStart[iRow];
3575 CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
3576 CoinBigIndex j;
3577 // Compute possible lower and upper ranges
3578 for (j = rStart; j < rEnd; ++j) {
3579 double value = element[j];
3580 int iColumn = column[j];
3581 if (value > 0.0) {
3582 if (columnUpper[iColumn] >= large) {
3583 ++infiniteUpper;
3584 } else {
3585 maximumUp += columnUpper[iColumn] * value;
3586 }
3587 if (columnLower[iColumn] <= -large) {
3588 ++infiniteLower;
3589 } else {
3590 maximumDown += columnLower[iColumn] * value;
3591 }
3592 } else if (value < 0.0) {
3593 if (columnUpper[iColumn] >= large) {
3594 ++infiniteLower;
3595 } else {
3596 maximumDown += columnUpper[iColumn] * value;
3597 }
3598 if (columnLower[iColumn] <= -large) {
3599 ++infiniteUpper;
3600 } else {
3601 maximumUp += columnLower[iColumn] * value;
3602 }
3603 }
3604 }
3605 // Build in a margin of error
3606 double maxUp = maximumUp+1.0e-8 * fabs(maximumUp);
3607 double maxDown = maximumDown-1.0e-8 * fabs(maximumDown);
3608 double rLower=rowLower[iRow];
3609 double rUpper=rowUpper[iRow];
3610 if (!infiniteUpper&&maxUp < rUpper - tolerance) {
3611 if (rUpper!=COIN_DBL_MAX||maximumUp<1.0e5)
3612 rUpper=maximumUp;
3613 }
3614 if (!infiniteLower&&maxDown > rLower + tolerance) {
3615 if (rLower!=-COIN_DBL_MAX||maximumDown>-1.0e5)
3616 rLower=maximumDown;
3617 }
3618 if (infiniteUpper)
3619 maxUp=COIN_DBL_MAX;
3620 if (fabs(maxUp-maxRhs[iRow])>1.0e-3*(1.0+fabs(maxUp)))
3621 printf("bad Maxup row %d maxUp %g maxRhs %g\n",iRow,maxUp,maxRhs[iRow]);
3622 maxRhs[iRow]=maxUp;
3623 if (infiniteLower)
3624 maxDown=-COIN_DBL_MAX;
3625 if (fabs(maxDown-minRhs[iRow])>1.0e-3*(1.0+fabs(maxDown)))
3626 printf("bad Maxdown row %d maxDown %g minRhs %g\n",iRow,maxDown,minRhs[iRow]);
3627 minRhs[iRow]=maxDown;
3628 assert(rLower<=rUpper);
3629 }
3630 }
3631 // end debug
3632 }
3633 }
3634 }
3635 if (tightenType>1) {
3636 // relax
3637 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3638 // relax enough so will have correct dj
3639 double lower=saveLower[iColumn];
3640 double upper=saveUpper[iColumn];
3641 columnLower[iColumn] = CoinMax(lower,columnLower[iColumn] - multiplier);
3642 columnUpper[iColumn] = CoinMin(upper,columnUpper[iColumn] + multiplier);
3643 }
3644 }
3645 #endif
3646 } else {
3647 #if ABC_NORMAL_DEBUG > 0
3648 printf("Redo tighten to relax by 1 for integers\n");
3649 #endif
3650 }
3651 #if ABC_NORMAL_DEBUG > 0
3652 printf("Tighten - phase1 %d bounds changed, phase2 %d lower, %d upper\n",
3653 totalTightened, numberLowerChanged, numberUpperChanged);
3654 #endif
3655 if (integerType_) {
3656 const double *columnScale = scaleToExternal_ + maximumAbcNumberRows_;
3657 const double *inverseColumnScale = scaleFromExternal_ + maximumAbcNumberRows_;
3658 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3659 if (integerType_[iColumn]) {
3660 double value;
3661 double lower = columnLower[iColumn] * inverseColumnScale[iColumn];
3662 double upper = columnUpper[iColumn] * inverseColumnScale[iColumn];
3663 value = floor(lower + 0.5);
3664 if (fabs(value - lower) > primalTolerance_)
3665 value = ceil(lower);
3666 columnLower_[iColumn] = value;
3667 columnLower[iColumn] = value * columnScale[iColumn];
3668 value = floor(upper + 0.5);
3669 if (fabs(value - upper) > primalTolerance_)
3670 value = floor(upper);
3671 columnUpper_[iColumn] = value;
3672 columnUpper[iColumn] = value * columnScale[iColumn];
3673 if (columnLower_[iColumn] > columnUpper_[iColumn])
3674 numberInfeasible++;
3675 }
3676 }
3677 if (numberInfeasible) {
3678 handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
3679 << numberInfeasible
3680 << CoinMessageEol;
3681 // restore column bounds
3682 CoinMemcpyN(saveLower, numberColumns_, columnLower_);
3683 CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
3684 }
3685 }
3686 } else {
3687 handler_->message(CLP_SIMPLEX_INFEASIBILITIES, messages_)
3688 << numberInfeasible
3689 << CoinMessageEol;
3690 // restore column bounds
3691 CoinMemcpyN(saveLower, numberColumns_, columnLower_);
3692 CoinMemcpyN(saveUpper, numberColumns_, columnUpper_);
3693 }
3694 if (!numberInfeasible) {
3695 // tighten rows as well
3696 for (int iRow = 0; iRow < numberRows_; iRow++) {
3697 if (rowLower[iRow] > -large || rowUpper[iRow] < large) {
3698 // possible row
3699 int infiniteUpper = 0;
3700 int infiniteLower = 0;
3701 double maximumUp = 0.0;
3702 double maximumDown = 0.0;
3703 CoinBigIndex rStart = rowStart[iRow];
3704 CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
3705 CoinBigIndex j;
3706 // Compute possible lower and upper ranges
3707 for (j = rStart; j < rEnd; ++j) {
3708 double value = element[j];
3709 int iColumn = column[j];
3710 if (value > 0.0) {
3711 if (columnUpper[iColumn] >= large) {
3712 ++infiniteUpper;
3713 } else {
3714 maximumUp += columnUpper[iColumn] * value;
3715 }
3716 if (columnLower[iColumn] <= -large) {
3717 ++infiniteLower;
3718 } else {
3719 maximumDown += columnLower[iColumn] * value;
3720 }
3721 } else if (value < 0.0) {
3722 if (columnUpper[iColumn] >= large) {
3723 ++infiniteLower;
3724 } else {
3725 maximumDown += columnUpper[iColumn] * value;
3726 }
3727 if (columnLower[iColumn] <= -large) {
3728 ++infiniteUpper;
3729 } else {
3730 maximumUp += columnLower[iColumn] * value;
3731 }
3732 }
3733 }
3734 // Build in a margin of error
3735 double maxUp = maximumUp + 1.0e-8 * fabs(maximumUp);
3736 double maxDown = maximumDown - 1.0e-8 * fabs(maximumDown);
3737 if (!infiniteUpper && maxUp < rowUpper[iRow] - tolerance) {
3738 if (rowUpper[iRow] != COIN_DBL_MAX || maximumUp < 1.0e5)
3739 rowUpper[iRow] = maximumUp;
3740 }
3741 if (!infiniteLower && maxDown > rowLower[iRow] + tolerance) {
3742 if (rowLower[iRow] != -COIN_DBL_MAX || maximumDown > -1.0e5)
3743 rowLower[iRow] = maximumDown;
3744 }
3745 assert(rowLower[iRow] <= rowUpper[iRow]);
3746 if (rowUpper[iRow] - rowLower[iRow] < 1.0e-12) {
3747 if (fabs(rowLower[iRow]) > fabs(rowUpper[iRow]))
3748 rowLower[iRow] = rowUpper[iRow];
3749 else
3750 rowUpper[iRow] = rowLower[iRow];
3751 }
3752 }
3753 }
3754 // flip back
3755 for (int iRow = 0; iRow < numberRows_; iRow++) {
3756 double value = -rowLower[iRow];
3757 rowLower[iRow] = -rowUpper[iRow];
3758 rowUpper[iRow] = value;
3759 }
3760 // move back - relaxed unless integer
3761 if (!integerType_) {
3762 double *COIN_RESTRICT lowerSaved = abcLower_ + maximumNumberTotal_;
3763 double *COIN_RESTRICT upperSaved = abcUpper_ + maximumNumberTotal_;
3764 double tolerance = 2.0 * primalTolerance_;
3765 for (int i = 0; i < numberTotal_; i++) {
3766 double lower = abcLower_[i];
3767 double upper = abcUpper_[i];
3768 if (lower != upper) {
3769 lower = CoinMax(lower - tolerance, lowerSaved[i]);
3770 upper = CoinMin(upper + tolerance, upperSaved[i]);
3771 }
3772 abcLower_[i] = lower;
3773 lowerSaved[i] = lower;
3774 abcUpper_[i] = upper;
3775 upperSaved[i] = upper;
3776 }
3777 } else {
3778 CoinAbcMemcpy(abcLower_ + maximumNumberTotal_, abcLower_, numberTotal_);
3779 CoinAbcMemcpy(abcUpper_ + maximumNumberTotal_, abcUpper_, numberTotal_);
3780 }
3781 #if 0
3782 for (int i=0;i<numberTotal_;i++) {
3783 if (abcSolution_[i+maximumNumberTotal_]>abcUpper_[i+maximumNumberTotal_]+1.0e-5)
3784 printf ("above %d %g %g\n",i,
3785 abcSolution_[i+maximumNumberTotal_],abcUpper_[i+maximumNumberTotal_]);
3786 if (abcSolution_[i+maximumNumberTotal_]<abcLower_[i+maximumNumberTotal_]-1.0e-5)
3787 printf ("below %d %g %g\n",i,
3788 abcSolution_[i+maximumNumberTotal_],abcLower_[i+maximumNumberTotal_]);
3789 }
3790 #endif
3791 permuteBasis();
3792 // move solution
3793 CoinAbcMemcpy(abcSolution_, solutionSaved_, numberTotal_);
3794 }
3795 CoinAbcMemset0(minRhs, numberRows_);
3796 CoinAbcMemset0(maxRhs, numberRows_);
3797 for (int i = 0; i < 5; i++)
3798 setAvailableArray(whichArray[i]);
3799 return (numberInfeasible);
3800 }
3801 // dual
3802 #include "AbcSimplexDual.hpp"
3803 #ifdef ABC_INHERIT
3804 // Returns 0 if dual can be skipped
doAbcDual()3805 int ClpSimplex::doAbcDual()
3806 {
3807 if ((abcState_ & CLP_ABC_WANTED) == 0) {
3808 return 1; // not wanted
3809 } else {
3810 assert(abcSimplex_);
3811 int crashState = (abcState_ >> 8) & 7;
3812 if (crashState && crashState < 4) {
3813 switch (crashState) {
3814 case 1:
3815 crash(1000.0, 1);
3816 break;
3817 case 2:
3818 crash(abcSimplex_->dualBound(), 0);
3819 break;
3820 case 3:
3821 crash(0.0, 3);
3822 break;
3823 }
3824 }
3825 // this and abcSimplex_ are approximately same pointer
3826 if ((abcState_ & CLP_ABC_FULL_DONE) == 0) {
3827 abcSimplex_->gutsOfResize(numberRows_, numberColumns_);
3828 abcSimplex_->translate(DO_SCALE_AND_MATRIX | DO_BASIS_AND_ORDER | DO_STATUS | DO_SOLUTION);
3829 //abcSimplex_->permuteIn();
3830 int maximumPivotsAbc = CoinMin(abcSimplex_->factorization()->maximumPivots(), numberColumns_ + 200);
3831 abcSimplex_->factorization()->maximumPivots(maximumPivotsAbc);
3832 if (numberColumns_)
3833 abcSimplex_->tightenPrimalBounds();
3834 }
3835 if ((abcSimplex_->stateOfProblem() & VALUES_PASS2) != 0) {
3836 if (solution_)
3837 crashState = 6;
3838 }
3839 if (crashState && crashState > 3) {
3840 abcSimplex_->crash(crashState - 2);
3841 }
3842 if (crashState && crashState < 4) {
3843 abcSimplex_->putStuffInBasis(1 + 2);
3844 }
3845 int returnCode = abcSimplex_->doAbcDual();
3846 //set to force crossover test returnCode=1;
3847 abcSimplex_->permuteOut(ROW_PRIMAL_OK | ROW_DUAL_OK | COLUMN_PRIMAL_OK | COLUMN_DUAL_OK | ALL_STATUS_OK);
3848 if (returnCode) {
3849 // fix as this model is all messed up e.g. scaling
3850 scalingFlag_ = abs(scalingFlag_);
3851 //delete [] rowScale_;
3852 //delete [] columnScale_;
3853 rowScale_ = NULL;
3854 columnScale_ = NULL;
3855 #if 0
3856 delete [] savedRowScale_;
3857 delete [] savedColumnScale_;
3858 savedRowScale_ = NULL;
3859 savedColumnScale_ = NULL;
3860 inverseRowScale_ = NULL;
3861 inverseColumnScale_ = NULL;
3862 #endif
3863 if (perturbation_ == 50)
3864 perturbation_ = 51;
3865 //perturbation_=100;
3866 #if ABC_NORMAL_DEBUG > 0
3867 printf("Bad exit with status of %d\n", problemStatus_);
3868 #endif
3869 //problemStatus_=10;
3870 int status = problemStatus_;
3871 //problemStatus_=-1;
3872 if (status != 10) {
3873 dual(); // Do ClpSimplexDual
3874 } else {
3875 moreSpecialOptions_ |= 256;
3876 primal(1);
3877 }
3878 } else {
3879 // double check objective
3880 //printf("%g %g\n",objectiveValue(),abcSimplex_->objectiveValue());
3881 perturbation_ = 100;
3882 moreSpecialOptions_ |= 256;
3883 //primal(1);
3884 }
3885 return returnCode;
3886 }
3887 }
3888 #endif
3889 #include "AbcSimplexPrimal.hpp"
3890 // Do dual (return 1 if cleanup needed)
doAbcDual()3891 int AbcSimplex::doAbcDual()
3892 {
3893 if (objective_) {
3894 objective_->setActivated(0);
3895 } else {
3896 // create dummy stuff
3897 assert(!numberColumns_);
3898 if (!numberRows_)
3899 problemStatus_ = 0; // say optimal
3900 return 0;
3901 }
3902 /* Note use of "down casting". The only class the user sees is AbcSimplex.
3903 Classes AbcSimplexDual, AbcSimplexPrimal, (AbcSimplexNonlinear)
3904 and AbcSimplexOther all exist and inherit from AbcSimplex but have no
3905 additional data and have no destructor or (non-default) constructor.
3906
3907 This is to stop classes becoming too unwieldy and so I (JJHF) can use e.g. "perturb"
3908 in primal and dual.
3909
3910 As far as I can see this is perfectly safe.
3911 */
3912 algorithm_ = -1;
3913 baseIteration_ = numberIterations_;
3914 if (!abcMatrix_->gotRowCopy())
3915 abcMatrix_->createRowCopy();
3916 #ifdef EARLY_FACTORIZE
3917 if (maximumIterations() > 1000000 && maximumIterations() < 1000999) {
3918 numberEarly_ = maximumIterations() - 1000000;
3919 #if ABC_NORMAL_DEBUG > 0
3920 printf("Setting numberEarly_ to %d\n", numberEarly_);
3921 #endif
3922 numberEarly_ += numberEarly_ << 16;
3923 }
3924 #endif
3925 minimumThetaMovement_ = 1.0e-10;
3926 if (fabs(infeasibilityCost_ - 1.0e10) < 1.1e9
3927 && fabs(infeasibilityCost_ - 1.0e10) > 1.0e5 && false) {
3928 minimumThetaMovement_ = 1.0e-3 / fabs(infeasibilityCost_ - 1.0e10);
3929 printf("trying minimum theta movement of %g\n", minimumThetaMovement_);
3930 infeasibilityCost_ = 1.0e10;
3931 }
3932 int savePerturbation = perturbation_;
3933 static_cast< AbcSimplexDual * >(this)->dual();
3934 minimumThetaMovement_ = 1.0e-10;
3935 if ((specialOptions_ & 2048) != 0 && problemStatus_ == 10 && !numberPrimalInfeasibilities_
3936 && sumDualInfeasibilities_ < 1000.0 * dualTolerance_)
3937 problemStatus_ = 0; // ignore
3938 onStopped(); // set secondary status if stopped
3939 if (problemStatus_ == 1 && numberFlagged_) {
3940 problemStatus_ = 10;
3941 static_cast< AbcSimplexPrimal * >(this)->unflag();
3942 }
3943 if (perturbation_ < 101) {
3944 //printf("Perturbed djs?\n");
3945 double *save = abcCost_;
3946 abcCost_ = costSaved_;
3947 computeInternalObjectiveValue();
3948 abcCost_ = save;
3949 }
3950 int totalNumberIterations = numberIterations_;
3951 if (problemStatus_ == 10 && (moreSpecialOptions_ & 32768) != 0 && sumDualInfeasibilities_ < 0.1) {
3952 problemStatus_ = 0;
3953 }
3954 #ifndef TRY_ABC_GUS
3955 if (problemStatus_ == 10) {
3956 int savePerturbation = perturbation_;
3957 perturbation_ = 100;
3958 copyFromSaved(2);
3959 minimumThetaMovement_ = 1.0e-12;
3960 baseIteration_ = numberIterations_;
3961 static_cast< AbcSimplexPrimal * >(this)->primal(0);
3962 totalNumberIterations += numberIterations_ - baseIteration_;
3963 perturbation_ = savePerturbation;
3964 }
3965 #else
3966 int iPass = 0;
3967 while (problemStatus_ == 10 && minimumThetaMovement_ > 0.99999e-12) {
3968 iPass++;
3969 if (initialSumInfeasibilities_ >= 0.0) {
3970 if (savePerturbation >= 100) {
3971 perturbation_ = 100;
3972 } else {
3973 if (savePerturbation == 50)
3974 perturbation_ = CoinMax(51, HEAVY_PERTURBATION - 4 - iPass); //smaller
3975 else
3976 perturbation_ = CoinMax(51, savePerturbation - 1 - iPass); //smaller
3977 }
3978 } else {
3979 perturbation_ = savePerturbation;
3980 abcFactorization_->setPivots(100000); // force factorization
3981 initialSumInfeasibilities_ = 1.23456789e-5;
3982 // clean pivots
3983 int numberBasic = 0;
3984 int *pivot = pivotVariable();
3985 for (int i = 0; i < numberTotal_; i++) {
3986 if (getInternalStatus(i) == basic)
3987 pivot[numberBasic++] = i;
3988 }
3989 assert(numberBasic == numberRows_);
3990 }
3991 copyFromSaved(14);
3992 // Say second call
3993 if (iPass > 3)
3994 moreSpecialOptions_ |= 256;
3995 if (iPass > 2)
3996 perturbation_ = 101;
3997 baseIteration_ = numberIterations_;
3998 static_cast< AbcSimplexPrimal * >(this)->primal(0);
3999 totalNumberIterations += numberIterations_ - baseIteration_;
4000 if (problemStatus_ == 10) {
4001 // reduce
4002 minimumThetaMovement_ *= 1.0e-1;
4003 copyFromSaved(14);
4004 baseIteration_ = numberIterations_;
4005 static_cast< AbcSimplexDual * >(this)->dual();
4006 totalNumberIterations += numberIterations_ - baseIteration_;
4007 }
4008 perturbation_ = savePerturbation;
4009 }
4010 #endif
4011 numberIterations_ = totalNumberIterations;
4012 return (problemStatus_ < 0 || problemStatus_ == 4 || problemStatus_ == 10) ? 1 : 0;
4013 }
dual()4014 int AbcSimplex::dual()
4015 {
4016 if (!abcMatrix_->gotRowCopy())
4017 abcMatrix_->createRowCopy();
4018 assert(!numberFlagged_);
4019 baseIteration_ = numberIterations_;
4020 //double savedPivotTolerance = abcFactorization_->pivotTolerance();
4021 if (objective_) {
4022 objective_->setActivated(0);
4023 } else {
4024 // create dummy stuff
4025 assert(!numberColumns_);
4026 if (!numberRows_)
4027 problemStatus_ = 0; // say optimal
4028 return 0;
4029 }
4030 /* Note use of "down casting". The only class the user sees is AbcSimplex.
4031 Classes AbcSimplexDual, AbcSimplexPrimal, (AbcSimplexNonlinear)
4032 and AbcSimplexOther all exist and inherit from AbcSimplex but have no
4033 additional data and have no destructor or (non-default) constructor.
4034
4035 This is to stop classes becoming too unwieldy and so I (JJF) can use e.g. "perturb"
4036 in primal and dual.
4037
4038 As far as I can see this is perfectly safe.
4039 */
4040 minimumThetaMovement_ = 1.0e-12;
4041 int returnCode = static_cast< AbcSimplexDual * >(this)->dual();
4042 if ((specialOptions_ & 2048) != 0 && problemStatus_ == 10 && !numberPrimalInfeasibilities_
4043 && sumDualInfeasibilities_ < 1000.0 * dualTolerance_ && perturbation_ >= 100)
4044 problemStatus_ = 0; // ignore
4045 if (problemStatus_ == 10) {
4046 #if 1 //ABC_NORMAL_DEBUG>0
4047 printf("return and use ClpSimplexPrimal\n");
4048 abort();
4049 #else
4050 //printf("Cleaning up with primal\n");
4051 //lastAlgorithm=1;
4052 int savePerturbation = perturbation_;
4053 int saveLog = handler_->logLevel();
4054 //handler_->setLogLevel(63);
4055 perturbation_ = 101;
4056 bool denseFactorization = initialDenseFactorization();
4057 // It will be safe to allow dense
4058 setInitialDenseFactorization(true);
4059 // Allow for catastrophe
4060 int saveMax = intParam_[ClpMaxNumIteration];
4061 if (numberIterations_) {
4062 // normal
4063 if (intParam_[ClpMaxNumIteration] > 100000 + numberIterations_)
4064 intParam_[ClpMaxNumIteration]
4065 = numberIterations_ + 1000 + 2 * numberRows_ + numberColumns_;
4066 } else {
4067 // Not normal allow more
4068 baseIteration_ += 2 * (numberRows_ + numberColumns_);
4069 }
4070 // check which algorithms allowed
4071 int dummy;
4072 if (problemStatus_ == 10 && saveObjective == objective_)
4073 startFinishOptions |= 2;
4074 baseIteration_ = numberIterations_;
4075 // Say second call
4076 moreSpecialOptions_ |= 256;
4077 minimumThetaMovement_ = 1.0e-12;
4078 returnCode = static_cast< AbcSimplexPrimal * >(this)->primal(1, startFinishOptions);
4079 // Say not second call
4080 moreSpecialOptions_ &= ~256;
4081 baseIteration_ = 0;
4082 if (saveObjective != objective_) {
4083 // We changed objective to see if infeasible
4084 delete objective_;
4085 objective_ = saveObjective;
4086 if (!problemStatus_) {
4087 // carry on
4088 returnCode = static_cast< AbcSimplexPrimal * >(this)->primal(1, startFinishOptions);
4089 }
4090 }
4091 if (problemStatus_ == 3 && numberIterations_ < saveMax) {
4092 // flatten solution and try again
4093 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
4094 if (getInternalStatus(iSequence) != basic) {
4095 setInternalStatus(iSequence, superBasic);
4096 // but put to bound if close
4097 if (fabs(rowActivity_[iSequence] - rowLower_[iSequence])
4098 <= primalTolerance_) {
4099 abcSolution_[iSequence] = abcLower_[iSequence];
4100 setInternalStatus(iSequence, atLowerBound);
4101 } else if (fabs(rowActivity_[iSequence] - rowUpper_[iSequence])
4102 <= primalTolerance_) {
4103 abcSolution_[iSequence] = abcUpper_[iSequence];
4104 setInternalStatus(iSequence, atUpperBound);
4105 }
4106 }
4107 }
4108 problemStatus_ = -1;
4109 intParam_[ClpMaxNumIteration] = CoinMin(numberIterations_ + 1000 + 2 * numberRows_ + numberColumns_, saveMax);
4110 perturbation_ = savePerturbation;
4111 baseIteration_ = numberIterations_;
4112 // Say second call
4113 moreSpecialOptions_ |= 256;
4114 returnCode = static_cast< AbcSimplexPrimal * >(this)->primal(0, startFinishOptions);
4115 // Say not second call
4116 moreSpecialOptions_ &= ~256;
4117 baseIteration_ = 0;
4118 computeObjectiveValue();
4119 // can't rely on djs either
4120 memset(reducedCost_, 0, numberColumns_ * sizeof(double));
4121 }
4122 intParam_[ClpMaxNumIteration] = saveMax;
4123
4124 setInitialDenseFactorization(denseFactorization);
4125 perturbation_ = savePerturbation;
4126 if (problemStatus_ == 10) {
4127 if (!numberPrimalInfeasibilities_)
4128 problemStatus_ = 0;
4129 else
4130 problemStatus_ = 4;
4131 }
4132 handler_->setLogLevel(saveLog);
4133 #endif
4134 }
4135 onStopped(); // set secondary status if stopped
4136 return returnCode;
4137 }
4138 #ifdef ABC_INHERIT
4139 // primal
4140 // Returns 0 if primal can be skipped
doAbcPrimal(int ifValuesPass)4141 int ClpSimplex::doAbcPrimal(int ifValuesPass)
4142 {
4143 if ((abcState_ & CLP_ABC_WANTED) == 0) {
4144 return 1; // not wanted
4145 } else {
4146 assert(abcSimplex_);
4147 if (ifValuesPass)
4148 abcSimplex_->setStateOfProblem(abcSimplex_->stateOfProblem() | VALUES_PASS);
4149 int crashState = (abcState_ >> 8) & 7;
4150 if (crashState && crashState < 4) {
4151 switch (crashState) {
4152 case 1:
4153 crash(1000.0, 1);
4154 break;
4155 case 2:
4156 crash(abcSimplex_->dualBound(), 0);
4157 break;
4158 case 3:
4159 crash(0.0, 3);
4160 break;
4161 }
4162 }
4163 // this and abcSimplex_ are approximately same pointer
4164 if ((abcState_ & CLP_ABC_FULL_DONE) == 0) {
4165 abcSimplex_->gutsOfResize(numberRows_, numberColumns_);
4166 abcSimplex_->translate(DO_SCALE_AND_MATRIX | DO_BASIS_AND_ORDER | DO_STATUS | DO_SOLUTION);
4167 //abcSimplex_->permuteIn();
4168 int maximumPivotsAbc = CoinMin(abcSimplex_->factorization()->maximumPivots(), numberColumns_ + 200);
4169 abcSimplex_->factorization()->maximumPivots(maximumPivotsAbc);
4170 abcSimplex_->copyFromSaved(1);
4171 }
4172 if ((abcSimplex_->stateOfProblem() & VALUES_PASS2) != 0) {
4173 if (solution_)
4174 crashState = 6;
4175 }
4176 if (crashState && crashState > 3) {
4177 abcSimplex_->crash(crashState - 2);
4178 }
4179 if (crashState && crashState < 4) {
4180 abcSimplex_->putStuffInBasis(1 + 2);
4181 }
4182 int returnCode = abcSimplex_->doAbcPrimal(ifValuesPass);
4183 //set to force crossover test returnCode=1;
4184 abcSimplex_->permuteOut(ROW_PRIMAL_OK | ROW_DUAL_OK | COLUMN_PRIMAL_OK | COLUMN_DUAL_OK | ALL_STATUS_OK);
4185 if (returnCode && problemStatus_ < 3) {
4186 // fix as this model is all messed up e.g. scaling
4187 scalingFlag_ = abs(scalingFlag_);
4188 //delete [] rowScale_;
4189 //delete [] columnScale_;
4190 rowScale_ = NULL;
4191 columnScale_ = NULL;
4192 #if 0
4193 delete [] savedRowScale_;
4194 delete [] savedColumnScale_;
4195 savedRowScale_ = NULL;
4196 savedColumnScale_ = NULL;
4197 inverseRowScale_ = NULL;
4198 inverseColumnScale_ = NULL;
4199 #endif
4200 if (perturbation_ == 50)
4201 perturbation_ = 51;
4202 //perturbation_=100;
4203 #if ABC_NORMAL_DEBUG > 0
4204 if (problemStatus_)
4205 printf("Bad exit with status of %d\n", problemStatus_);
4206 #endif
4207 //problemStatus_=10;
4208 int status = problemStatus_;
4209 // should not be using
4210 for (int i = 0; i < 6; i++) {
4211 if (rowArray_[i])
4212 rowArray_[i]->clear();
4213 if (columnArray_[i])
4214 columnArray_[i]->clear();
4215 }
4216 //problemStatus_=-1;
4217 if (status < 3 && status != 0) {
4218 if (status != 10) {
4219 primal(); // Do ClpSimplexPrimal
4220 } else {
4221 moreSpecialOptions_ |= 256;
4222 dual();
4223 }
4224 }
4225 } else {
4226 // double check objective
4227 //printf("%g %g\n",objectiveValue(),abcSimplex_->objectiveValue());
4228 }
4229 numberIterations_ = abcSimplex_->numberIterations();
4230 return returnCode;
4231 }
4232 }
4233 #endif
4234 // Do primal (return 1 if cleanup needed)
doAbcPrimal(int ifValuesPass)4235 int AbcSimplex::doAbcPrimal(int ifValuesPass)
4236 {
4237 if (objective_) {
4238 objective_->setActivated(0);
4239 } else {
4240 // create dummy stuff
4241 assert(!numberColumns_);
4242 if (!numberRows_)
4243 problemStatus_ = 0; // say optimal
4244 return 0;
4245 }
4246 /* Note use of "down casting". The only class the user sees is AbcSimplex.
4247 Classes AbcSimplexDual, AbcSimplexPrimal, (AbcSimplexNonlinear)
4248 and AbcSimplexOther all exist and inherit from AbcSimplex but have no
4249 additional data and have no destructor or (non-default) constructor.
4250
4251 This is to stop classes becoming too unwieldy and so I (JJHF) can use e.g. "perturb"
4252 in primal and dual.
4253
4254 As far as I can see this is perfectly safe.
4255 */
4256 algorithm_ = -1;
4257 int savePerturbation = perturbation_;
4258 baseIteration_ = numberIterations_;
4259 if (!abcMatrix_->gotRowCopy())
4260 abcMatrix_->createRowCopy();
4261 #ifdef EARLY_FACTORIZE
4262 if (maximumIterations() > 1000000 && maximumIterations() < 1000999) {
4263 numberEarly_ = maximumIterations() - 1000000;
4264 #if ABC_NORMAL_DEBUG > 0
4265 printf("Setting numberEarly_ to %d\n", numberEarly_);
4266 #endif
4267 numberEarly_ += numberEarly_ << 16;
4268 }
4269 #endif
4270 // add values pass options
4271 minimumThetaMovement_ = 1.0e-12;
4272 if ((specialOptions_ & 8192) != 0)
4273 minimumThetaMovement_ = 1.0e-15;
4274 int returnCode = static_cast< AbcSimplexPrimal * >(this)->primal(ifValuesPass);
4275 stateOfProblem_ &= ~VALUES_PASS;
4276 #ifndef TRY_ABC_GUS
4277 if (returnCode == 10) {
4278 copyFromSaved(14);
4279 int savePerturbation = perturbation_;
4280 perturbation_ = 51;
4281 minimumThetaMovement_ = 1.0e-12;
4282 returnCode = static_cast< AbcSimplexDual * >(this)->dual();
4283 perturbation_ = savePerturbation;
4284 }
4285 #else
4286 minimumThetaMovement_ = 1.0e-12;
4287 int iPass = 0;
4288 while (problemStatus_ == 10 && minimumThetaMovement_ > 1.0e-15) {
4289 iPass++;
4290 if (minimumThetaMovement_ == 1.0e-12)
4291 perturbation_ = CoinMin(savePerturbation, 55);
4292 else
4293 perturbation_ = 100;
4294 copyFromSaved(14);
4295 // reduce ?
4296 // Say second call
4297 // Say second call
4298 if (iPass > 3)
4299 moreSpecialOptions_ |= 256;
4300 if (iPass > 2)
4301 perturbation_ = 101;
4302 baseIteration_ = numberIterations_;
4303 // make sure no superbasic
4304 cleanStatus();
4305 if ((specialOptions_ & 8192) == 0)
4306 static_cast< AbcSimplexDual * >(this)->dual();
4307 baseIteration_ = numberIterations_;
4308 if (problemStatus_ == 10) {
4309 if (initialSumInfeasibilities_ >= 0.0) {
4310 if (savePerturbation >= 100) {
4311 perturbation_ = 100;
4312 } else {
4313 if (savePerturbation == 50)
4314 perturbation_ = CoinMax(51, HEAVY_PERTURBATION - 4 - iPass); //smaller
4315 else
4316 perturbation_ = CoinMax(51, savePerturbation - 1 - iPass); //smaller
4317 }
4318 } else {
4319 specialOptions_ |= 8192; // stop going to dual
4320 perturbation_ = savePerturbation;
4321 abcFactorization_->setPivots(100000); // force factorization
4322 initialSumInfeasibilities_ = 1.23456789e-5;
4323 // clean pivots
4324 int numberBasic = 0;
4325 int *pivot = pivotVariable();
4326 for (int i = 0; i < numberTotal_; i++) {
4327 if (getInternalStatus(i) == basic)
4328 pivot[numberBasic++] = i;
4329 }
4330 assert(numberBasic == numberRows_);
4331 }
4332 if (iPass > 2)
4333 perturbation_ = 100;
4334 copyFromSaved(14);
4335 baseIteration_ = numberIterations_;
4336 static_cast< AbcSimplexPrimal * >(this)->primal(0);
4337 }
4338 minimumThetaMovement_ *= 0.5;
4339 perturbation_ = savePerturbation;
4340 }
4341 #endif
4342 return returnCode;
4343 }
4344 /* Sets up all slack basis and resets solution to
4345 as it was after initial load or readMps */
allSlackBasis()4346 void AbcSimplex::allSlackBasis()
4347 {
4348 // set column status to one nearest zero
4349 CoinFillN(internalStatus_, numberRows_, static_cast< unsigned char >(basic));
4350 const double *COIN_RESTRICT lower = abcLower_;
4351 const double *COIN_RESTRICT upper = abcUpper_;
4352 double *COIN_RESTRICT solution = abcSolution_;
4353 // But set value to zero if lb <0.0 and ub>0.0
4354 for (int i = maximumAbcNumberRows_; i < numberTotal_; i++) {
4355 if (lower[i] >= 0.0) {
4356 solution[i] = lower[i];
4357 setInternalStatus(i, atLowerBound);
4358 } else if (upper[i] <= 0.0) {
4359 solution[i] = upper[i];
4360 setInternalStatus(i, atUpperBound);
4361 } else if (lower[i] < -1.0e20 && upper[i] > 1.0e20) {
4362 // free
4363 solution[i] = 0.0;
4364 setInternalStatus(i, isFree);
4365 } else if (fabs(lower[i]) < fabs(upper[i])) {
4366 solution[i] = 0.0;
4367 setInternalStatus(i, atLowerBound);
4368 } else {
4369 solution[i] = 0.0;
4370 setInternalStatus(i, atUpperBound);
4371 }
4372 }
4373 }
4374 #if 0 //ndef SLIM_NOIO
4375 // Read an mps file from the given filename
4376 int
4377 AbcSimplex::readMps(const char *filename,
4378 bool keepNames,
4379 bool ignoreErrors)
4380 {
4381 int status = ClpSimplex::readMps(filename, keepNames, ignoreErrors);
4382 ClpSimplex * thisModel=this;
4383 gutsOfResize(thisModel->numberRows(),thisModel->numberColumns());
4384 translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
4385 return status;
4386 }
4387 // Read GMPL files from the given filenames
4388 int
4389 AbcSimplex::readGMPL(const char *filename, const char * dataName,
4390 bool keepNames)
4391 {
4392 int status = ClpSimplex::readGMPL(filename, dataName, keepNames);
4393 translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
4394 return status;
4395 }
4396 // Read file in LP format from file with name filename.
4397 int
4398 AbcSimplex::readLp(const char *filename, const double epsilon )
4399 {
4400 FILE *fp = fopen(filename, "r");
4401
4402 if(!fp) {
4403 printf("### ERROR: AbcSimplex::readLp(): Unable to open file %s for reading\n",
4404 filename);
4405 return(1);
4406 }
4407 CoinLpIO m;
4408 m.readLp(fp, epsilon);
4409 fclose(fp);
4410
4411 // set problem name
4412 setStrParam(ClpProbName, m.getProblemName());
4413 // no errors
4414 loadProblem(*m.getMatrixByRow(), m.getColLower(), m.getColUpper(),
4415 m.getObjCoefficients(), m.getRowLower(), m.getRowUpper());
4416
4417 if (m.integerColumns()) {
4418 integerType_ = new char[numberColumns_];
4419 CoinMemcpyN(m.integerColumns(), numberColumns_, integerType_);
4420 } else {
4421 integerType_ = NULL;
4422 }
4423 translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
4424 unsigned int maxLength = 0;
4425 int iRow;
4426 rowNames_ = std::vector<std::string> ();
4427 columnNames_ = std::vector<std::string> ();
4428 rowNames_.reserve(numberRows_);
4429 for (iRow = 0; iRow < numberRows_; iRow++) {
4430 const char * name = m.rowName(iRow);
4431 if (name) {
4432 maxLength = CoinMax(maxLength, static_cast<unsigned int> (strlen(name)));
4433 rowNames_.push_back(name);
4434 } else {
4435 rowNames_.push_back("");
4436 }
4437 }
4438
4439 int iColumn;
4440 columnNames_.reserve(numberColumns_);
4441 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4442 const char * name = m.columnName(iColumn);
4443 if (name) {
4444 maxLength = CoinMax(maxLength, static_cast<unsigned int> (strlen(name)));
4445 columnNames_.push_back(name);
4446 } else {
4447 columnNames_.push_back("");
4448 }
4449 }
4450 lengthNames_ = static_cast<int> (maxLength);
4451
4452 return 0;
4453 }
4454 #endif
4455 // Factorization frequency
factorizationFrequency() const4456 int AbcSimplex::factorizationFrequency() const
4457 {
4458 if (abcFactorization_)
4459 return abcFactorization_->maximumPivots();
4460 else
4461 return -1;
4462 }
setFactorizationFrequency(int value)4463 void AbcSimplex::setFactorizationFrequency(int value)
4464 {
4465 if (abcFactorization_)
4466 abcFactorization_->maximumPivots(value);
4467 }
4468
4469 /* Get a clean factorization - i.e. throw out singularities
4470 may do more later */
cleanFactorization(int ifValuesPass)4471 int AbcSimplex::cleanFactorization(int ifValuesPass)
4472 {
4473 int status = internalFactorize(ifValuesPass ? 10 : 0);
4474 if (status < 0) {
4475 return 1; // some error
4476 } else {
4477 firstFree_ = 0;
4478 return 0;
4479 }
4480 }
4481 // Save data
4482 ClpDataSave
saveData()4483 AbcSimplex::saveData()
4484 {
4485 ClpDataSave saved;
4486 saved.dualBound_ = dualBound_;
4487 saved.infeasibilityCost_ = infeasibilityCost_;
4488 saved.pivotTolerance_ = abcFactorization_->pivotTolerance();
4489 saved.zeroFactorizationTolerance_ = abcFactorization_->zeroTolerance();
4490 saved.zeroSimplexTolerance_ = zeroTolerance_;
4491 saved.perturbation_ = perturbation_;
4492 saved.forceFactorization_ = forceFactorization_;
4493 saved.acceptablePivot_ = acceptablePivot_;
4494 saved.objectiveScale_ = objectiveScale_;
4495 // Progress indicator
4496 abcProgress_.fillFromModel(this);
4497 return saved;
4498 }
4499 // Restore data
restoreData(ClpDataSave saved)4500 void AbcSimplex::restoreData(ClpDataSave saved)
4501 {
4502 //abcFactorization_->sparseThreshold(saved.sparseThreshold_);
4503 abcFactorization_->pivotTolerance(saved.pivotTolerance_);
4504 abcFactorization_->zeroTolerance(saved.zeroFactorizationTolerance_);
4505 zeroTolerance_ = saved.zeroSimplexTolerance_;
4506 perturbation_ = saved.perturbation_;
4507 infeasibilityCost_ = saved.infeasibilityCost_;
4508 dualBound_ = saved.dualBound_;
4509 forceFactorization_ = saved.forceFactorization_;
4510 objectiveScale_ = saved.objectiveScale_;
4511 acceptablePivot_ = saved.acceptablePivot_;
4512 }
4513 // To flag a variable
setFlagged(int sequence)4514 void AbcSimplex::setFlagged(int sequence)
4515 {
4516 assert(sequence >= 0 && sequence < numberTotal_);
4517 if (sequence >= 0) {
4518 internalStatus_[sequence] |= 64;
4519 // use for real disaster lastFlaggedIteration_ = numberIterations_;
4520 numberFlagged_++;
4521 }
4522 }
4523 #ifndef NDEBUG
4524 // For errors to make sure print to screen
4525 // only called in debug mode
indexError(int index,std::string methodName)4526 static void indexError(int index,
4527 std::string methodName)
4528 {
4529 std::cerr << "Illegal index " << index << " in AbcSimplex::" << methodName << std::endl;
4530 throw CoinError("Illegal index", methodName, "AbcSimplex");
4531 }
4532 #endif
4533 // After modifying first copy refreshes second copy and marks as updated
refreshCosts()4534 void AbcSimplex::refreshCosts()
4535 {
4536 whatsChanged_ &= ~OBJECTIVE_SAME;
4537 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
4538 }
refreshLower(unsigned int type)4539 void AbcSimplex::refreshLower(unsigned int type)
4540 {
4541 if (type)
4542 whatsChanged_ &= type;
4543 CoinAbcMemcpy(abcLower_, lowerSaved_, numberTotal_);
4544 }
refreshUpper(unsigned int type)4545 void AbcSimplex::refreshUpper(unsigned int type)
4546 {
4547 if (type)
4548 whatsChanged_ &= type;
4549 CoinAbcMemcpy(abcUpper_, upperSaved_, numberTotal_);
4550 }
4551 /* Set an objective function coefficient */
setObjectiveCoefficient(int elementIndex,double elementValue)4552 void AbcSimplex::setObjectiveCoefficient(int elementIndex, double elementValue)
4553 {
4554 #ifndef NDEBUG
4555 if (elementIndex < 0 || elementIndex >= numberColumns_) {
4556 indexError(elementIndex, "setObjectiveCoefficient");
4557 }
4558 #endif
4559 if (objective()[elementIndex] != elementValue) {
4560 objective()[elementIndex] = elementValue;
4561 // work arrays exist - update as well
4562 whatsChanged_ &= ~OBJECTIVE_SAME;
4563 double direction = optimizationDirection_ * objectiveScale_;
4564 costSaved_[elementIndex + maximumAbcNumberRows_] = direction * elementValue
4565 * columnUseScale_[elementIndex + maximumAbcNumberRows_];
4566 }
4567 }
4568 /* Set a single row lower bound<br>
4569 Use -DBL_MAX for -infinity. */
setRowLower(int elementIndex,double elementValue)4570 void AbcSimplex::setRowLower(int elementIndex, double elementValue)
4571 {
4572 #ifndef NDEBUG
4573 int n = numberRows_;
4574 if (elementIndex < 0 || elementIndex >= n) {
4575 indexError(elementIndex, "setRowLower");
4576 }
4577 #endif
4578 if (elementValue < -1.0e27)
4579 elementValue = -COIN_DBL_MAX;
4580 if (rowLower_[elementIndex] != elementValue) {
4581 rowLower_[elementIndex] = elementValue;
4582 // work arrays exist - update as well
4583 whatsChanged_ &= ~ROW_LOWER_SAME;
4584 if (rowLower_[elementIndex] == -COIN_DBL_MAX) {
4585 lowerSaved_[elementIndex] = -COIN_DBL_MAX;
4586 } else {
4587 lowerSaved_[elementIndex] = elementValue * rhsScale_
4588 * rowUseScale_[elementIndex];
4589 }
4590 }
4591 }
4592
4593 /* Set a single row upper bound<br>
4594 Use DBL_MAX for infinity. */
setRowUpper(int elementIndex,double elementValue)4595 void AbcSimplex::setRowUpper(int elementIndex, double elementValue)
4596 {
4597 #ifndef NDEBUG
4598 int n = numberRows_;
4599 if (elementIndex < 0 || elementIndex >= n) {
4600 indexError(elementIndex, "setRowUpper");
4601 }
4602 #endif
4603 if (elementValue > 1.0e27)
4604 elementValue = COIN_DBL_MAX;
4605 if (rowUpper_[elementIndex] != elementValue) {
4606 rowUpper_[elementIndex] = elementValue;
4607 // work arrays exist - update as well
4608 whatsChanged_ &= ~ROW_UPPER_SAME;
4609 if (rowUpper_[elementIndex] == COIN_DBL_MAX) {
4610 upperSaved_[elementIndex] = COIN_DBL_MAX;
4611 } else {
4612 upperSaved_[elementIndex] = elementValue * rhsScale_
4613 * rowUseScale_[elementIndex];
4614 }
4615 }
4616 }
4617
4618 /* Set a single row lower and upper bound */
setRowBounds(int elementIndex,double lowerValue,double upperValue)4619 void AbcSimplex::setRowBounds(int elementIndex,
4620 double lowerValue, double upperValue)
4621 {
4622 #ifndef NDEBUG
4623 int n = numberRows_;
4624 if (elementIndex < 0 || elementIndex >= n) {
4625 indexError(elementIndex, "setRowBounds");
4626 }
4627 #endif
4628 if (lowerValue < -1.0e27)
4629 lowerValue = -COIN_DBL_MAX;
4630 if (upperValue > 1.0e27)
4631 upperValue = COIN_DBL_MAX;
4632 //CoinAssert (upperValue>=lowerValue);
4633 if (rowLower_[elementIndex] != lowerValue) {
4634 rowLower_[elementIndex] = lowerValue;
4635 // work arrays exist - update as well
4636 whatsChanged_ &= ~ROW_LOWER_SAME;
4637 if (rowLower_[elementIndex] == -COIN_DBL_MAX) {
4638 lowerSaved_[elementIndex] = -COIN_DBL_MAX;
4639 } else {
4640 lowerSaved_[elementIndex] = lowerValue * rhsScale_
4641 * rowUseScale_[elementIndex];
4642 }
4643 }
4644 if (rowUpper_[elementIndex] != upperValue) {
4645 rowUpper_[elementIndex] = upperValue;
4646 // work arrays exist - update as well
4647 whatsChanged_ &= ~ROW_UPPER_SAME;
4648 if (rowUpper_[elementIndex] == COIN_DBL_MAX) {
4649 upperSaved_[elementIndex] = COIN_DBL_MAX;
4650 } else {
4651 upperSaved_[elementIndex] = upperValue * rhsScale_
4652 * rowUseScale_[elementIndex];
4653 }
4654 }
4655 }
setRowSetBounds(const int * indexFirst,const int * indexLast,const double * boundList)4656 void AbcSimplex::setRowSetBounds(const int *indexFirst,
4657 const int *indexLast,
4658 const double *boundList)
4659 {
4660 #ifndef NDEBUG
4661 int n = numberRows_;
4662 #endif
4663 int numberChanged = 0;
4664 const int *saveFirst = indexFirst;
4665 while (indexFirst != indexLast) {
4666 const int iRow = *indexFirst++;
4667 #ifndef NDEBUG
4668 if (iRow < 0 || iRow >= n) {
4669 indexError(iRow, "setRowSetBounds");
4670 }
4671 #endif
4672 double lowerValue = *boundList++;
4673 double upperValue = *boundList++;
4674 if (lowerValue < -1.0e27)
4675 lowerValue = -COIN_DBL_MAX;
4676 if (upperValue > 1.0e27)
4677 upperValue = COIN_DBL_MAX;
4678 //CoinAssert (upperValue>=lowerValue);
4679 if (rowLower_[iRow] != lowerValue) {
4680 rowLower_[iRow] = lowerValue;
4681 whatsChanged_ &= ~ROW_LOWER_SAME;
4682 numberChanged++;
4683 }
4684 if (rowUpper_[iRow] != upperValue) {
4685 rowUpper_[iRow] = upperValue;
4686 whatsChanged_ &= ~ROW_UPPER_SAME;
4687 numberChanged++;
4688 }
4689 }
4690 if (numberChanged) {
4691 indexFirst = saveFirst;
4692 while (indexFirst != indexLast) {
4693 const int iRow = *indexFirst++;
4694 if (rowLower_[iRow] == -COIN_DBL_MAX) {
4695 lowerSaved_[iRow] = -COIN_DBL_MAX;
4696 } else {
4697 lowerSaved_[iRow] = rowLower_[iRow] * rhsScale_
4698 * rowUseScale_[iRow];
4699 }
4700 if (rowUpper_[iRow] == COIN_DBL_MAX) {
4701 upperSaved_[iRow] = COIN_DBL_MAX;
4702 } else {
4703 upperSaved_[iRow] = rowUpper_[iRow] * rhsScale_
4704 * rowUseScale_[iRow];
4705 }
4706 }
4707 }
4708 }
4709 //-----------------------------------------------------------------------------
4710 /* Set a single column lower bound<br>
4711 Use -DBL_MAX for -infinity. */
setColumnLower(int elementIndex,double elementValue)4712 void AbcSimplex::setColumnLower(int elementIndex, double elementValue)
4713 {
4714 #ifndef NDEBUG
4715 int n = numberColumns_;
4716 if (elementIndex < 0 || elementIndex >= n) {
4717 indexError(elementIndex, "setColumnLower");
4718 }
4719 #endif
4720 if (elementValue < -1.0e27)
4721 elementValue = -COIN_DBL_MAX;
4722 if (columnLower_[elementIndex] != elementValue) {
4723 columnLower_[elementIndex] = elementValue;
4724 // work arrays exist - update as well
4725 whatsChanged_ &= ~COLUMN_LOWER_SAME;
4726 double value;
4727 if (columnLower_[elementIndex] == -COIN_DBL_MAX) {
4728 value = -COIN_DBL_MAX;
4729 } else {
4730 value = elementValue * rhsScale_
4731 * inverseColumnUseScale_[elementIndex];
4732 }
4733 lowerSaved_[elementIndex + maximumAbcNumberRows_] = value;
4734 }
4735 }
4736
4737 /* Set a single column upper bound<br>
4738 Use DBL_MAX for infinity. */
setColumnUpper(int elementIndex,double elementValue)4739 void AbcSimplex::setColumnUpper(int elementIndex, double elementValue)
4740 {
4741 #ifndef NDEBUG
4742 int n = numberColumns_;
4743 if (elementIndex < 0 || elementIndex >= n) {
4744 indexError(elementIndex, "setColumnUpper");
4745 }
4746 #endif
4747 if (elementValue > 1.0e27)
4748 elementValue = COIN_DBL_MAX;
4749 if (columnUpper_[elementIndex] != elementValue) {
4750 columnUpper_[elementIndex] = elementValue;
4751 // work arrays exist - update as well
4752 whatsChanged_ &= ~COLUMN_UPPER_SAME;
4753 double value;
4754 if (columnUpper_[elementIndex] == COIN_DBL_MAX) {
4755 value = COIN_DBL_MAX;
4756 } else {
4757 value = elementValue * rhsScale_
4758 * inverseColumnUseScale_[elementIndex];
4759 }
4760 upperSaved_[elementIndex + maximumAbcNumberRows_] = value;
4761 }
4762 }
4763
4764 /* Set a single column lower and upper bound */
setColumnBounds(int elementIndex,double lowerValue,double upperValue)4765 void AbcSimplex::setColumnBounds(int elementIndex,
4766 double lowerValue, double upperValue)
4767 {
4768 #ifndef NDEBUG
4769 int n = numberColumns_;
4770 if (elementIndex < 0 || elementIndex >= n) {
4771 indexError(elementIndex, "setColumnBounds");
4772 }
4773 #endif
4774 if (lowerValue < -1.0e27)
4775 lowerValue = -COIN_DBL_MAX;
4776 if (columnLower_[elementIndex] != lowerValue) {
4777 columnLower_[elementIndex] = lowerValue;
4778 // work arrays exist - update as well
4779 whatsChanged_ &= ~COLUMN_LOWER_SAME;
4780 if (columnLower_[elementIndex] == -COIN_DBL_MAX) {
4781 lowerSaved_[elementIndex + maximumAbcNumberRows_] = -COIN_DBL_MAX;
4782 } else {
4783 lowerSaved_[elementIndex + maximumAbcNumberRows_] = lowerValue * rhsScale_
4784 * inverseColumnUseScale_[elementIndex];
4785 }
4786 }
4787 if (upperValue > 1.0e27)
4788 upperValue = COIN_DBL_MAX;
4789 if (columnUpper_[elementIndex] != upperValue) {
4790 columnUpper_[elementIndex] = upperValue;
4791 // work arrays exist - update as well
4792 whatsChanged_ &= ~COLUMN_UPPER_SAME;
4793 if (columnUpper_[elementIndex] == COIN_DBL_MAX) {
4794 upperSaved_[elementIndex + maximumAbcNumberRows_] = COIN_DBL_MAX;
4795 } else {
4796 upperSaved_[elementIndex + maximumAbcNumberRows_] = upperValue * rhsScale_
4797 * inverseColumnUseScale_[elementIndex];
4798 }
4799 }
4800 }
setColumnSetBounds(const int * indexFirst,const int * indexLast,const double * boundList)4801 void AbcSimplex::setColumnSetBounds(const int *indexFirst,
4802 const int *indexLast,
4803 const double *boundList)
4804 {
4805 #ifndef NDEBUG
4806 int n = numberColumns_;
4807 #endif
4808 int numberChanged = 0;
4809 const int *saveFirst = indexFirst;
4810 while (indexFirst != indexLast) {
4811 const int iColumn = *indexFirst++;
4812 #ifndef NDEBUG
4813 if (iColumn < 0 || iColumn >= n) {
4814 indexError(iColumn, "setColumnSetBounds");
4815 }
4816 #endif
4817 double lowerValue = *boundList++;
4818 double upperValue = *boundList++;
4819 if (lowerValue < -1.0e27)
4820 lowerValue = -COIN_DBL_MAX;
4821 if (upperValue > 1.0e27)
4822 upperValue = COIN_DBL_MAX;
4823 //CoinAssert (upperValue>=lowerValue);
4824 if (columnLower_[iColumn] != lowerValue) {
4825 columnLower_[iColumn] = lowerValue;
4826 whatsChanged_ &= ~COLUMN_LOWER_SAME;
4827 numberChanged++;
4828 }
4829 if (columnUpper_[iColumn] != upperValue) {
4830 columnUpper_[iColumn] = upperValue;
4831 whatsChanged_ &= ~COLUMN_UPPER_SAME;
4832 numberChanged++;
4833 }
4834 }
4835 if (numberChanged) {
4836 indexFirst = saveFirst;
4837 while (indexFirst != indexLast) {
4838 const int iColumn = *indexFirst++;
4839 if (columnLower_[iColumn] == -COIN_DBL_MAX) {
4840 lowerSaved_[iColumn + maximumAbcNumberRows_] = -COIN_DBL_MAX;
4841 } else {
4842 lowerSaved_[iColumn + maximumAbcNumberRows_] = columnLower_[iColumn] * rhsScale_
4843 * inverseColumnUseScale_[iColumn];
4844 }
4845 if (columnUpper_[iColumn] == COIN_DBL_MAX) {
4846 upperSaved_[iColumn + maximumAbcNumberRows_] = COIN_DBL_MAX;
4847 } else {
4848 upperSaved_[iColumn + maximumAbcNumberRows_] = columnUpper_[iColumn] * rhsScale_
4849 * inverseColumnUseScale_[iColumn];
4850 }
4851 }
4852 }
4853 }
4854 #ifndef SLIM_CLP
4855 #include "CoinWarmStartBasis.hpp"
4856
4857 // Returns a basis (to be deleted by user)
4858 CoinWarmStartBasis *
getBasis() const4859 AbcSimplex::getBasis() const
4860 {
4861 CoinWarmStartBasis *basis = new CoinWarmStartBasis();
4862 basis->setSize(numberColumns_, numberRows_);
4863
4864 unsigned char lookup[8] = { 2, 3, 255, 255, 0, 0, 1, 3 };
4865 for (int iRow = 0; iRow < numberRows_; iRow++) {
4866 int iStatus = getInternalStatus(iRow);
4867 iStatus = lookup[iStatus];
4868 basis->setArtifStatus(iRow, static_cast< CoinWarmStartBasis::Status >(iStatus));
4869 }
4870 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
4871 int iStatus = getInternalStatus(iColumn + maximumAbcNumberRows_);
4872 iStatus = lookup[iStatus];
4873 basis->setStructStatus(iColumn, static_cast< CoinWarmStartBasis::Status >(iStatus));
4874 }
4875 return basis;
4876 }
4877 #endif
4878 // Compute objective value from solution
computeObjectiveValue(bool useInternalArrays)4879 void AbcSimplex::computeObjectiveValue(bool useInternalArrays)
4880 {
4881 const double *obj = objective();
4882 if (!useInternalArrays) {
4883 objectiveValue_ = 0.0;
4884 for (int iSequence = 0; iSequence < numberColumns_; iSequence++) {
4885 double value = columnActivity_[iSequence];
4886 objectiveValue_ += value * obj[iSequence];
4887 }
4888 // But remember direction as we are using external objective
4889 objectiveValue_ *= optimizationDirection_;
4890 } else {
4891 rawObjectiveValue_ = 0.0;
4892 for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4893 double value = abcSolution_[iSequence];
4894 rawObjectiveValue_ += value * abcCost_[iSequence];
4895 }
4896 setClpSimplexObjectiveValue();
4897 }
4898 }
4899 // Compute minimization objective value from internal solution
4900 double
computeInternalObjectiveValue()4901 AbcSimplex::computeInternalObjectiveValue()
4902 {
4903 rawObjectiveValue_ = 0.0;
4904 for (int iSequence = maximumAbcNumberRows_; iSequence < numberTotal_; iSequence++) {
4905 double value = abcSolution_[iSequence];
4906 rawObjectiveValue_ += value * abcCost_[iSequence];
4907 }
4908 setClpSimplexObjectiveValue();
4909 return objectiveValue_ - dblParam_[ClpObjOffset];
4910 }
4911 // If user left factorization frequency then compute
defaultFactorizationFrequency()4912 void AbcSimplex::defaultFactorizationFrequency()
4913 {
4914 if (factorizationFrequency() == 200) {
4915 // User did not touch preset
4916 const int cutoff1 = 10000;
4917 const int cutoff2 = 100000;
4918 const int base = 75;
4919 const int freq0 = 50;
4920 const int freq1 = 150;
4921 const int maximum = 10000;
4922 int frequency;
4923 if (numberRows_ < cutoff1)
4924 frequency = base + numberRows_ / freq0;
4925 else
4926 frequency = base + cutoff1 / freq0 + (numberRows_ - cutoff1) / freq1;
4927 setFactorizationFrequency(CoinMin(maximum, frequency));
4928 }
4929 }
4930 // Gets clean and emptyish factorization
4931 AbcSimplexFactorization *
getEmptyFactorization()4932 AbcSimplex::getEmptyFactorization()
4933 {
4934 if ((specialOptions_ & 65536) == 0) {
4935 assert(!abcFactorization_);
4936 abcFactorization_ = new AbcSimplexFactorization();
4937 } else if (!abcFactorization_) {
4938 abcFactorization_ = new AbcSimplexFactorization();
4939 }
4940 return abcFactorization_;
4941 }
4942 // Resizes rim part of model
resize(int newNumberRows,int newNumberColumns)4943 void AbcSimplex::resize(int newNumberRows, int newNumberColumns)
4944 {
4945 assert(maximumAbcNumberRows_ >= 0);
4946 // save
4947 int numberRows = numberRows_;
4948 int numberColumns = numberColumns_;
4949 ClpSimplex::resize(newNumberRows, newNumberColumns);
4950 // back out
4951 numberRows_ = numberRows;
4952 numberColumns_ = numberColumns;
4953 gutsOfResize(newNumberRows, newNumberColumns);
4954 }
4955 // Return true if the objective limit test can be relied upon
isObjectiveLimitTestValid() const4956 bool AbcSimplex::isObjectiveLimitTestValid() const
4957 {
4958 if (problemStatus_ == 0) {
4959 return true;
4960 } else if (problemStatus_ == 1) {
4961 // ok if dual
4962 return (algorithm_ < 0);
4963 } else if (problemStatus_ == 2) {
4964 // ok if primal
4965 return (algorithm_ > 0);
4966 } else {
4967 return false;
4968 }
4969 }
4970 /*
4971 Permutes in from ClpModel data - assumes scale factors done
4972 and AbcMatrix exists but is in original order (including slacks)
4973 For now just add basicArray at end
4974 ==
4975 But could partition into
4976 normal (i.e. reasonable lower/upper)
4977 abnormal - free, odd bounds
4978 fixed
4979 ==
4980 sets a valid pivotVariable
4981 Slacks always shifted by offset
4982 Fixed variables always shifted by offset
4983 Recode to allow row objective so can use pi from idiot etc
4984 */
permuteIn()4985 void AbcSimplex::permuteIn()
4986 {
4987 // for now only if maximumAbcNumberRows_==numberRows_
4988 //assert(maximumAbcNumberRows_==numberRows_);
4989 numberTotal_ = maximumAbcNumberRows_ + numberColumns_;
4990 double direction = optimizationDirection_;
4991 // all this could use cilk
4992 // move primal stuff to end
4993 const double *COIN_RESTRICT rowScale = scaleFromExternal_;
4994 const double *COIN_RESTRICT inverseRowScale = scaleToExternal_;
4995 const double *COIN_RESTRICT columnScale = scaleToExternal_ + maximumAbcNumberRows_;
4996 const double *COIN_RESTRICT inverseColumnScale = scaleFromExternal_ + maximumAbcNumberRows_;
4997 double *COIN_RESTRICT offsetRhs = offsetRhs_;
4998 double *COIN_RESTRICT saveLower = lowerSaved_ + maximumAbcNumberRows_;
4999 double *COIN_RESTRICT saveUpper = upperSaved_ + maximumAbcNumberRows_;
5000 double *COIN_RESTRICT saveCost = costSaved_ + maximumAbcNumberRows_;
5001 double *COIN_RESTRICT saveSolution = solutionSaved_ + maximumAbcNumberRows_;
5002 CoinAbcMemset0(offsetRhs, maximumAbcNumberRows_);
5003 const double *COIN_RESTRICT objective = this->objective();
5004 objectiveOffset_ = 0.0;
5005 double *COIN_RESTRICT offset = offset_ + maximumAbcNumberRows_;
5006 const int *COIN_RESTRICT row = abcMatrix_->matrix()->getIndices();
5007 const CoinBigIndex *COIN_RESTRICT columnStart = abcMatrix_->matrix()->getVectorStarts();
5008 const double *COIN_RESTRICT elementByColumn = abcMatrix_->matrix()->getElements();
5009 largestGap_ = 1.0e-12;
5010 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
5011 double scale = inverseColumnScale[iColumn];
5012 double lowerValue = columnLower_[iColumn];
5013 double upperValue = columnUpper_[iColumn];
5014 double thisOffset = 0.0;
5015 #if 1 //ndef CLP_USER_DRIVEN
5016 double lowerValue2 = fabs(lowerValue);
5017 double upperValue2 = fabs(upperValue);
5018 if (CoinMin(lowerValue2, upperValue2) < 1000.0) {
5019 // move to zero
5020 if (lowerValue2 > upperValue2)
5021 thisOffset = upperValue;
5022 else
5023 thisOffset = lowerValue;
5024 }
5025 #endif
5026 offset[iColumn] = thisOffset;
5027 if (thisOffset) {
5028 objectiveOffset_ += thisOffset * objective[iColumn] * optimizationDirection_;
5029 double scaledOffset = thisOffset * scale;
5030 for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn + 1]; j++) {
5031 int iRow = row[j];
5032 offsetRhs[iRow] += scaledOffset * elementByColumn[j];
5033 }
5034 }
5035 lowerValue -= thisOffset;
5036 if (lowerValue > -1.0e30)
5037 lowerValue *= scale;
5038 saveLower[iColumn] = lowerValue;
5039 upperValue -= thisOffset;
5040 if (upperValue < 1.0e30)
5041 upperValue *= scale;
5042 saveUpper[iColumn] = upperValue;
5043 largestGap_ = CoinMax(largestGap_, upperValue - lowerValue);
5044 saveSolution[iColumn] = scale * (columnActivity_[iColumn] - thisOffset);
5045 saveCost[iColumn] = objective[iColumn] * direction * columnScale[iColumn];
5046 }
5047 CoinAbcMemset0(offset_, maximumAbcNumberRows_);
5048 saveLower -= maximumAbcNumberRows_;
5049 saveUpper -= maximumAbcNumberRows_;
5050 saveCost -= maximumAbcNumberRows_;
5051 saveSolution -= maximumAbcNumberRows_;
5052 for (int iRow = 0; iRow < numberRows_; iRow++) {
5053 // note switch of lower to upper
5054 double upperValue = -rowLower_[iRow];
5055 double lowerValue = -rowUpper_[iRow];
5056 double thisOffset = offsetRhs_[iRow];
5057 double scale = rowScale[iRow];
5058 if (lowerValue > -1.0e30)
5059 lowerValue = lowerValue * scale + thisOffset;
5060 saveLower[iRow] = lowerValue;
5061 if (upperValue < 1.0e30)
5062 upperValue = upperValue * scale + thisOffset;
5063 saveUpper[iRow] = upperValue;
5064 largestGap_ = CoinMax(largestGap_, upperValue - lowerValue);
5065 saveCost[iRow] = 0.0;
5066 dual_[iRow] *= direction * inverseRowScale[iRow];
5067 saveSolution[iRow] = 0.0; // not necessary
5068 }
5069 dualBound_ = CoinMin(dualBound_, largestGap_);
5070 // Compute rhsScale_ and objectiveScale_
5071 double minValue = COIN_DBL_MAX;
5072 double maxValue = 0.0;
5073 CoinAbcMinMaxAbsNormalValues(costSaved_ + maximumAbcNumberRows_, numberTotal_ - maximumAbcNumberRows_, minValue, maxValue);
5074 // scale to 1000.0 ?
5075 if (minValue && false) {
5076 objectiveScale_ = 1000.0 / sqrt(minValue * maxValue);
5077 objectiveScale_ = CoinMin(1.0, 1000.0 / maxValue);
5078 #ifndef NDEBUG
5079 double smallestNormal = COIN_DBL_MAX;
5080 double smallestAny = COIN_DBL_MAX;
5081 double largestAny = 0.0;
5082 for (int i = 0; i < numberTotal_; i++) {
5083 double value = fabs(costSaved_[i]);
5084 if (value) {
5085 if (value > 1.0e-8)
5086 smallestNormal = CoinMin(smallestNormal, value);
5087 smallestAny = CoinMin(smallestAny, value);
5088 largestAny = CoinMax(largestAny, value);
5089 }
5090 }
5091 printf("objectiveScale_ %g min_used %g (min_reasonable %g, min_any %g) max_used %g (max_any %g)\n",
5092 objectiveScale_, minValue, smallestNormal, smallestAny, maxValue, largestAny);
5093 #endif
5094 } else {
5095 //objectiveScale_=1.0;
5096 }
5097 CoinAbcScale(costSaved_, objectiveScale_, numberTotal_);
5098 minValue = COIN_DBL_MAX;
5099 maxValue = 0.0;
5100 CoinAbcMinMaxAbsNormalValues(lowerSaved_, numberTotal_, minValue, maxValue);
5101 CoinAbcMinMaxAbsNormalValues(upperSaved_, numberTotal_, minValue, maxValue);
5102 // scale to 100.0 ?
5103 if (minValue && false) {
5104 rhsScale_ = 100.0 / sqrt(minValue * maxValue);
5105 #ifndef NDEBUG
5106 double smallestNormal = COIN_DBL_MAX;
5107 double smallestAny = COIN_DBL_MAX;
5108 double largestAny = 0.0;
5109 for (int i = 0; i < numberTotal_; i++) {
5110 double value = fabs(lowerSaved_[i]);
5111 if (value && value != COIN_DBL_MAX) {
5112 if (value > 1.0e-8)
5113 smallestNormal = CoinMin(smallestNormal, value);
5114 smallestAny = CoinMin(smallestAny, value);
5115 largestAny = CoinMax(largestAny, value);
5116 }
5117 }
5118 for (int i = 0; i < numberTotal_; i++) {
5119 double value = fabs(upperSaved_[i]);
5120 if (value && value != COIN_DBL_MAX) {
5121 if (value > 1.0e-8)
5122 smallestNormal = CoinMin(smallestNormal, value);
5123 smallestAny = CoinMin(smallestAny, value);
5124 largestAny = CoinMax(largestAny, value);
5125 }
5126 }
5127 printf("rhsScale_ %g min_used %g (min_reasonable %g, min_any %g) max_used %g (max_any %g)\n",
5128 rhsScale_, minValue, smallestNormal, smallestAny, maxValue, largestAny);
5129 #endif
5130 } else {
5131 rhsScale_ = 1.0;
5132 }
5133 CoinAbcScaleNormalValues(lowerSaved_, rhsScale_, 1.0e-13, numberTotal_);
5134 CoinAbcScaleNormalValues(upperSaved_, rhsScale_, 1.0e-13, numberTotal_);
5135 // copy
5136 CoinAbcMemcpy(abcLower_, abcLower_ + maximumNumberTotal_, numberTotal_);
5137 CoinAbcMemcpy(abcUpper_, abcUpper_ + maximumNumberTotal_, numberTotal_);
5138 CoinAbcMemcpy(abcSolution_, abcSolution_ + maximumNumberTotal_, numberTotal_);
5139 CoinAbcMemcpy(abcCost_, abcCost_ + maximumNumberTotal_, numberTotal_);
5140 }
permuteBasis()5141 void AbcSimplex::permuteBasis()
5142 {
5143 assert(abcPivotVariable_ || (!numberRows_ && !numberColumns_));
5144 int numberBasic = 0;
5145 // from Clp enum to Abc enum (and bound flip)
5146 unsigned char lookupToAbcSlack[6] = { 4, 6, 0 /*1*/, 1 /*0*/, 5, 7 };
5147 unsigned char *COIN_RESTRICT getStatus = status_ + numberColumns_;
5148 double *COIN_RESTRICT solutionSaved = solutionSaved_;
5149 double *COIN_RESTRICT lowerSaved = lowerSaved_;
5150 double *COIN_RESTRICT upperSaved = upperSaved_;
5151 bool ordinaryVariables = true;
5152 bool valuesPass = (stateOfProblem_ & VALUES_PASS) != 0;
5153 if (valuesPass) {
5154 // get solution
5155 CoinAbcMemset0(abcSolution_, numberRows_);
5156 abcMatrix_->timesIncludingSlacks(-1.0, abcSolution_, abcSolution_);
5157 //double * temp = new double[numberRows_];
5158 //memset(temp,0,numberRows_*sizeof(double));
5159 //matrix_->times(1.0,columnActivity_,temp);
5160 CoinAbcMemcpy(solutionSaved_, abcSolution_, numberRows_);
5161 int n = 0;
5162 for (int i = 0; i < numberTotal_; i++) {
5163 if (solutionSaved_[i] > upperSaved_[i] + 1.0e-5)
5164 n++;
5165 else if (solutionSaved_[i] < lowerSaved_[i] - 1.0e-5)
5166 n++;
5167 }
5168 #if ABC_NORMAL_DEBUG > 0
5169 if (n)
5170 printf("%d infeasibilities\n", n);
5171 #endif
5172 }
5173 // dual at present does not like superbasic
5174 for (int iRow = 0; iRow < numberRows_; iRow++) {
5175 unsigned char status = getStatus[iRow] & 7;
5176 AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbcSlack[status]);
5177 if (status != ClpSimplex::basic) {
5178 double lowerValue = lowerSaved[iRow];
5179 double upperValue = upperSaved[iRow];
5180 if (lowerValue == -COIN_DBL_MAX) {
5181 if (upperValue == COIN_DBL_MAX) {
5182 // free
5183 abcStatus = isFree;
5184 ordinaryVariables = false;
5185 } else {
5186 abcStatus = atUpperBound;
5187 }
5188 } else if (upperValue == COIN_DBL_MAX) {
5189 abcStatus = atLowerBound;
5190 } else if (lowerValue == upperValue) {
5191 abcStatus = isFixed;
5192 } else if (abcStatus == isFixed) {
5193 double value = solutionSaved[iRow];
5194 if (value - lowerValue < upperValue - value)
5195 abcStatus = atLowerBound;
5196 else
5197 abcStatus = atUpperBound;
5198 }
5199 if (valuesPass) {
5200 double value = solutionSaved[iRow];
5201 if (value > lowerValue + primalTolerance_ && value < upperValue - primalTolerance_ && (abcStatus == atLowerBound || abcStatus == atUpperBound))
5202 abcStatus = superBasic;
5203 }
5204 switch (abcStatus) {
5205 case isFixed:
5206 case atLowerBound:
5207 solutionSaved[iRow] = lowerValue;
5208 break;
5209 case atUpperBound:
5210 solutionSaved[iRow] = upperValue;
5211 break;
5212 default:
5213 ordinaryVariables = false;
5214 break;
5215 }
5216 } else {
5217 // basic
5218 abcPivotVariable_[numberBasic++] = iRow;
5219 }
5220 internalStatus_[iRow] = abcStatus;
5221 }
5222 // from Clp enum to Abc enum
5223 unsigned char lookupToAbc[6] = { 4, 6, 1, 0, 5, 7 };
5224 unsigned char *COIN_RESTRICT putStatus = internalStatus_ + maximumAbcNumberRows_;
5225 getStatus = status_;
5226 solutionSaved += maximumAbcNumberRows_;
5227 lowerSaved += maximumAbcNumberRows_;
5228 upperSaved += maximumAbcNumberRows_;
5229 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
5230 unsigned char status = getStatus[iColumn] & 7;
5231 AbcSimplex::Status abcStatus = static_cast< AbcSimplex::Status >(lookupToAbc[status]);
5232 if (status != ClpSimplex::basic) {
5233 double lowerValue = lowerSaved[iColumn];
5234 double upperValue = upperSaved[iColumn];
5235 if (lowerValue == -COIN_DBL_MAX) {
5236 if (upperValue == COIN_DBL_MAX) {
5237 // free
5238 abcStatus = isFree;
5239 ordinaryVariables = false;
5240 } else {
5241 abcStatus = atUpperBound;
5242 }
5243 } else if (upperValue == COIN_DBL_MAX) {
5244 abcStatus = atLowerBound;
5245 } else if (lowerValue == upperValue) {
5246 abcStatus = isFixed;
5247 } else if (abcStatus == isFixed) {
5248 double value = solutionSaved[iColumn];
5249 if (value - lowerValue < upperValue - value)
5250 abcStatus = atLowerBound;
5251 else
5252 abcStatus = atUpperBound;
5253 } else if (abcStatus == isFree) {
5254 abcStatus = superBasic;
5255 }
5256 if (valuesPass && (abcStatus == atLowerBound || abcStatus == atUpperBound)) {
5257 double value = solutionSaved[iColumn];
5258 if (value > lowerValue + primalTolerance_) {
5259 if (value < upperValue - primalTolerance_) {
5260 abcStatus = superBasic;
5261 } else {
5262 abcStatus = atUpperBound;
5263 }
5264 } else {
5265 abcStatus = atLowerBound;
5266 }
5267 }
5268 switch (abcStatus) {
5269 case isFixed:
5270 case atLowerBound:
5271 solutionSaved[iColumn] = lowerValue;
5272 break;
5273 case atUpperBound:
5274 solutionSaved[iColumn] = upperValue;
5275 break;
5276 default:
5277 ordinaryVariables = false;
5278 break;
5279 }
5280 } else {
5281 // basic
5282 if (numberBasic < numberRows_)
5283 abcPivotVariable_[numberBasic++] = iColumn + maximumAbcNumberRows_;
5284 else
5285 abcStatus = superBasic;
5286 }
5287 putStatus[iColumn] = abcStatus;
5288 }
5289 ordinaryVariables_ = ordinaryVariables ? 1 : 0;
5290 if (numberBasic < numberRows_) {
5291 for (int iRow = 0; iRow < numberRows_; iRow++) {
5292 AbcSimplex::Status status = getInternalStatus(iRow);
5293 if (status != AbcSimplex::basic) {
5294 abcPivotVariable_[numberBasic++] = iRow;
5295 setInternalStatus(iRow, basic);
5296 if (numberBasic == numberRows_)
5297 break;
5298 }
5299 }
5300 }
5301 // copy
5302 CoinAbcMemcpy(internalStatus_ + maximumNumberTotal_, internalStatus_, numberTotal_);
5303 }
5304 // Permutes out - bit settings same as stateOfProblem
permuteOut(int whatsWanted)5305 void AbcSimplex::permuteOut(int whatsWanted)
5306 {
5307 assert(whatsWanted);
5308 if ((whatsWanted & ALL_STATUS_OK) != 0 && (stateOfProblem_ & ALL_STATUS_OK) == 0) {
5309 stateOfProblem_ |= ALL_STATUS_OK;
5310 // from Abc enum to Clp enum (and bound flip)
5311 unsigned char lookupToClpSlack[8] = { 2, 3, 255, 255, 0, 0, 1, 5 };
5312 unsigned char *COIN_RESTRICT putStatus = status_ + numberColumns_;
5313 const unsigned char *COIN_RESTRICT getStatus = internalStatus_;
5314 for (int iRow = 0; iRow < numberRows_; iRow++) {
5315 putStatus[iRow] = lookupToClpSlack[getStatus[iRow] & 7];
5316 }
5317 // from Abc enum to Clp enum
5318 unsigned char lookupToClp[8] = { 3, 2, 255, 255, 0, 0, 1, 5 };
5319 putStatus = status_;
5320 getStatus = internalStatus_ + maximumAbcNumberRows_;
5321 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
5322 putStatus[iColumn] = lookupToClp[getStatus[iColumn] & 7];
5323 }
5324 }
5325 const double *COIN_RESTRICT rowScale = scaleFromExternal_;
5326 const double *COIN_RESTRICT inverseRowScale = scaleToExternal_;
5327 const double *COIN_RESTRICT columnScale = scaleToExternal_; // allow for offset in loop +maximumAbcNumberRows_;
5328 const double *COIN_RESTRICT inverseColumnScale = scaleFromExternal_; // allow for offset in loop +maximumAbcNumberRows_;
5329 double scaleC = 1.0 / objectiveScale_;
5330 double scaleR = 1.0 / rhsScale_;
5331 int numberPrimalScaled = 0;
5332 int numberPrimalUnscaled = 0;
5333 int numberDualScaled = 0;
5334 int numberDualUnscaled = 0;
5335 bool filledInSolution = false;
5336 if ((whatsWanted & ROW_PRIMAL_OK) != 0 && (stateOfProblem_ & ROW_PRIMAL_OK) == 0) {
5337 stateOfProblem_ |= ROW_PRIMAL_OK;
5338 if (!filledInSolution) {
5339 filledInSolution = true;
5340 CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_);
5341 CoinAbcScatterZeroTo(abcDj_, abcPivotVariable_, numberRows_);
5342 }
5343 // Collect infeasibilities
5344 double *COIN_RESTRICT putSolution = rowActivity_;
5345 const double *COIN_RESTRICT rowLower = rowLower_;
5346 const double *COIN_RESTRICT rowUpper = rowUpper_;
5347 const double *COIN_RESTRICT abcLower = abcLower_;
5348 const double *COIN_RESTRICT abcUpper = abcUpper_;
5349 const double *COIN_RESTRICT abcSolution = abcSolution_;
5350 const double *COIN_RESTRICT offsetRhs = offsetRhs_;
5351 for (int i = 0; i < numberRows_; i++) {
5352 double scaleFactor = inverseRowScale[i];
5353 double valueScaled = abcSolution[i];
5354 double lowerScaled = abcLower[i];
5355 double upperScaled = abcUpper[i];
5356 if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
5357 if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
5358 numberPrimalScaled++;
5359 else
5360 upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
5361 }
5362 double value = (offsetRhs[i] - valueScaled * scaleR) * scaleFactor;
5363 putSolution[i] = value;
5364 if (value < rowLower[i] - primalTolerance_)
5365 numberPrimalUnscaled++;
5366 else if (value > rowUpper[i] + primalTolerance_)
5367 numberPrimalUnscaled++;
5368 }
5369 }
5370 if ((whatsWanted & COLUMN_PRIMAL_OK) != 0 && (stateOfProblem_ & COLUMN_PRIMAL_OK) == 0) {
5371 stateOfProblem_ |= COLUMN_PRIMAL_OK;
5372 // Collect infeasibilities
5373 if (!filledInSolution) {
5374 filledInSolution = true;
5375 CoinAbcScatterTo(solutionBasic_, abcSolution_, abcPivotVariable_, numberRows_);
5376 CoinAbcScatterZeroTo(abcDj_, abcPivotVariable_, numberRows_);
5377 }
5378 // Collect infeasibilities
5379 double *COIN_RESTRICT putSolution = columnActivity_ - maximumAbcNumberRows_;
5380 const double *COIN_RESTRICT columnLower = columnLower_ - maximumAbcNumberRows_;
5381 const double *COIN_RESTRICT columnUpper = columnUpper_ - maximumAbcNumberRows_;
5382 for (int i = maximumAbcNumberRows_; i < numberTotal_; i++) {
5383 double scaleFactor = columnScale[i];
5384 double valueScaled = abcSolution_[i];
5385 double lowerScaled = abcLower_[i];
5386 double upperScaled = abcUpper_[i];
5387 if (lowerScaled > -1.0e20 || upperScaled < 1.0e20) {
5388 if (valueScaled < lowerScaled - primalTolerance_ || valueScaled > upperScaled + primalTolerance_)
5389 numberPrimalScaled++;
5390 else
5391 upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
5392 }
5393 double value = (valueScaled * scaleR) * scaleFactor + offset_[i];
5394 putSolution[i] = value;
5395 if (value < columnLower[i] - primalTolerance_)
5396 numberPrimalUnscaled++;
5397 else if (value > columnUpper[i] + primalTolerance_)
5398 numberPrimalUnscaled++;
5399 }
5400 }
5401 if ((whatsWanted & COLUMN_DUAL_OK) != 0 && (stateOfProblem_ & COLUMN_DUAL_OK) == 0) {
5402 // Get fixed djs
5403 CoinAbcMemcpy(abcDj_, abcCost_, numberTotal_);
5404 abcMatrix_->transposeTimesAll(dual_, abcDj_);
5405 stateOfProblem_ |= COLUMN_DUAL_OK;
5406 // Collect infeasibilities
5407 double *COIN_RESTRICT putDj = reducedCost_ - maximumAbcNumberRows_;
5408 for (int i = maximumAbcNumberRows_; i < numberTotal_; i++) {
5409 double scaleFactor = inverseColumnScale[i];
5410 double valueDual = abcDj_[i];
5411 double value = abcSolution_[i];
5412 bool aboveLower = value > abcLower_[i] + primalTolerance_;
5413 bool belowUpper = value < abcUpper_[i] - primalTolerance_;
5414 if (aboveLower && valueDual > dualTolerance_)
5415 numberDualScaled++;
5416 if (belowUpper && valueDual < -dualTolerance_)
5417 numberDualScaled++;
5418 valueDual *= scaleFactor * scaleC;
5419 putDj[i] = valueDual;
5420 if (aboveLower && valueDual > dualTolerance_)
5421 numberDualUnscaled++;
5422 if (belowUpper && valueDual < -dualTolerance_)
5423 numberDualUnscaled++;
5424 }
5425 }
5426 if ((whatsWanted & ROW_DUAL_OK) != 0 && (stateOfProblem_ & ROW_DUAL_OK) == 0) {
5427 stateOfProblem_ |= ROW_DUAL_OK;
5428 // Collect infeasibilities
5429 for (int i = 0; i < numberRows_; i++) {
5430 double scaleFactor = rowScale[i];
5431 double valueDual = abcDj_[i]; // ? +- and direction
5432 double value = abcSolution_[i];
5433 bool aboveLower = value > abcLower_[i] + primalTolerance_;
5434 bool belowUpper = value < abcUpper_[i] - primalTolerance_;
5435 if (aboveLower && valueDual > dualTolerance_)
5436 numberDualScaled++;
5437 if (belowUpper && valueDual < -dualTolerance_)
5438 numberDualScaled++;
5439 valueDual *= scaleFactor * scaleC;
5440 dual_[i] = -(valueDual - abcCost_[i]); // sign
5441 if (aboveLower && valueDual > dualTolerance_)
5442 numberDualUnscaled++;
5443 if (belowUpper && valueDual < -dualTolerance_)
5444 numberDualUnscaled++;
5445 }
5446 }
5447 // do after djs
5448 if (!problemStatus_ && (!secondaryStatus_ || secondaryStatus_ == 2 || secondaryStatus_ == 3)) {
5449 // See if we need to set secondary status
5450 if (numberPrimalUnscaled) {
5451 if (numberDualUnscaled || secondaryStatus_ == 3)
5452 secondaryStatus_ = 4;
5453 else
5454 secondaryStatus_ = 2;
5455 } else if (numberDualUnscaled) {
5456 if (secondaryStatus_ == 0)
5457 secondaryStatus_ = 3;
5458 else
5459 secondaryStatus_ = 4;
5460 }
5461 }
5462 if (scalingFlag_) {
5463 if (problemStatus_ == 2) {
5464 for (int i = 0; i < numberColumns_; i++) {
5465 ray_[i] *= columnScale[i];
5466 }
5467 } else if (problemStatus_ == 1 && ray_) {
5468 for (int i = 0; i < numberRows_; i++) {
5469 ray_[i] *= rowScale[i];
5470 }
5471 }
5472 }
5473 }
5474 #if ABC_DEBUG > 1
5475 // For debug - moves solution back to external and computes stuff
checkMoveBack(bool checkDuals)5476 void AbcSimplex::checkMoveBack(bool checkDuals)
5477 {
5478 stateOfProblem_ &= ~(ROW_PRIMAL_OK | ROW_DUAL_OK | COLUMN_PRIMAL_OK | COLUMN_DUAL_OK | ALL_STATUS_OK);
5479 permuteOut(ROW_PRIMAL_OK | ROW_DUAL_OK | COLUMN_PRIMAL_OK | COLUMN_DUAL_OK | ALL_STATUS_OK);
5480 ClpSimplex::computeObjectiveValue(false);
5481 #if ABC_NORMAL_DEBUG > 0
5482 printf("Check objective %g\n", objectiveValue() - objectiveOffset_);
5483 #endif
5484 double *region = new double[numberRows_ + numberColumns_];
5485 CoinAbcMemset0(region, numberRows_);
5486 ClpModel::times(1.0, columnActivity_, region);
5487 int numberInf;
5488 double sumInf;
5489 numberInf = 0;
5490 sumInf = 0.0;
5491 for (int i = 0; i < numberColumns_; i++) {
5492 if (columnActivity_[i] > columnUpper_[i] + 1.0e-7) {
5493 numberInf++;
5494 sumInf += columnActivity_[i] - columnUpper_[i];
5495 } else if (columnActivity_[i] < columnLower_[i] - 1.0e-7) {
5496 numberInf++;
5497 sumInf -= columnActivity_[i] - columnLower_[i];
5498 }
5499 }
5500 #if ABC_NORMAL_DEBUG > 0
5501 if (numberInf)
5502 printf("Check column infeasibilities %d sum %g\n", numberInf, sumInf);
5503 #endif
5504 numberInf = 0;
5505 sumInf = 0.0;
5506 for (int i = 0; i < numberRows_; i++) {
5507 if (rowActivity_[i] > rowUpper_[i] + 1.0e-7) {
5508 numberInf++;
5509 sumInf += rowActivity_[i] - rowUpper_[i];
5510 } else if (rowActivity_[i] < rowLower_[i] - 1.0e-7) {
5511 numberInf++;
5512 sumInf -= rowActivity_[i] - rowLower_[i];
5513 }
5514 }
5515 #if ABC_NORMAL_DEBUG > 0
5516 if (numberInf)
5517 printf("Check row infeasibilities %d sum %g\n", numberInf, sumInf);
5518 #endif
5519 CoinAbcMemcpy(region, objective(), numberColumns_);
5520 ClpModel::transposeTimes(-1.0, dual_, region);
5521 numberInf = 0;
5522 sumInf = 0.0;
5523 for (int i = 0; i < numberColumns_; i++) {
5524 if (columnUpper_[i] > columnLower_[i]) {
5525 if (getColumnStatus(i) == ClpSimplex::atLowerBound) {
5526 if (region[i] < -1.0e-7) {
5527 numberInf++;
5528 sumInf -= region[i];
5529 }
5530 } else if (getColumnStatus(i) == ClpSimplex::atUpperBound) {
5531 if (region[i] > 1.0e-7) {
5532 numberInf++;
5533 sumInf += region[i];
5534 }
5535 }
5536 }
5537 }
5538 #if ABC_NORMAL_DEBUG > 0
5539 if (numberInf)
5540 printf("Check column dj infeasibilities %d sum %g\n", numberInf, sumInf);
5541 #endif
5542 if (checkDuals) {
5543 numberInf = 0;
5544 sumInf = 0.0;
5545 for (int i = 0; i < numberRows_; i++) {
5546 if (rowUpper_[i] > rowLower_[i]) {
5547 if (getRowStatus(i) == ClpSimplex::atLowerBound) {
5548 if (dual_[i] < -1.0e-7) {
5549 numberInf++;
5550 sumInf -= region[i];
5551 }
5552 } else if (getRowStatus(i) == ClpSimplex::atUpperBound) {
5553 if (dual_[i] > 1.0e-7) {
5554 numberInf++;
5555 sumInf += region[i];
5556 }
5557 }
5558 }
5559 }
5560 #if ABC_NORMAL_DEBUG > 0
5561 if (numberInf)
5562 printf("Check row dual infeasibilities %d sum %g\n", numberInf, sumInf);
5563 #endif
5564 }
5565 delete[] region;
5566 }
5567 #if 0
5568 void xxxxxx(const char * where) {
5569 printf("xxxxx %s\n",where);
5570 }
5571 #endif
5572 // For debug - checks solutionBasic
checkSolutionBasic() const5573 void AbcSimplex::checkSolutionBasic() const
5574 {
5575 //work space
5576 int whichArray[2];
5577 for (int i = 0; i < 2; i++)
5578 whichArray[i] = getAvailableArray();
5579 CoinIndexedVector *arrayVector = &usefulArray_[whichArray[0]];
5580 double *solution = usefulArray_[whichArray[1]].denseVector();
5581 CoinAbcMemcpy(solution, abcSolution_, numberTotal_);
5582 // accumulate non basic stuff
5583 double *array = arrayVector->denseVector();
5584 CoinAbcScatterZeroTo(solution, abcPivotVariable_, numberRows_);
5585 abcMatrix_->timesIncludingSlacks(-1.0, solution, array);
5586 //arrayVector->scan(0,numberRows_,zeroTolerance_);
5587 // Ftran adjusted RHS
5588 //if (arrayVector->getNumElements())
5589 abcFactorization_->updateFullColumn(*arrayVector);
5590 CoinAbcScatterTo(array, solution, abcPivotVariable_, numberRows_);
5591 double largestDifference = 0.0;
5592 int whichDifference = -1;
5593 for (int i = 0; i < numberRows_; i++) {
5594 double difference = fabs(solutionBasic_[i] - array[i]);
5595 if (difference > 1.0e-5 && numberRows_ < 100)
5596 printf("solutionBasic difference is %g on row %d solutionBasic_ %g computed %g\n",
5597 difference, i, solutionBasic_[i], array[i]);
5598 if (difference > largestDifference) {
5599 largestDifference = difference;
5600 whichDifference = i;
5601 }
5602 }
5603 if (largestDifference > 1.0e-9)
5604 printf("Debug largest solutionBasic difference is %g on row %d solutionBasic_ %g computed %g\n",
5605 largestDifference, whichDifference, solutionBasic_[whichDifference], array[whichDifference]);
5606 arrayVector->clear();
5607 CoinAbcMemset0(solution, numberTotal_);
5608 for (int i = 0; i < 2; i++)
5609 setAvailableArray(whichArray[i]);
5610 }
5611
5612 // For debug - summarizes dj situation
checkDjs(int type) const5613 void AbcSimplex::checkDjs(int type) const
5614 {
5615 if (type) {
5616 //work space
5617 int whichArrays[2];
5618 for (int i = 0; i < 2; i++)
5619 whichArrays[i] = getAvailableArray();
5620
5621 CoinIndexedVector *arrayVector = &usefulArray_[whichArrays[0]];
5622 double *array = arrayVector->denseVector();
5623 int *index = arrayVector->getIndices();
5624 int number = 0;
5625 for (int iRow = 0; iRow < numberRows_; iRow++) {
5626 double value = costBasic_[iRow];
5627 if (value) {
5628 array[iRow] = value;
5629 index[number++] = iRow;
5630 }
5631 }
5632 arrayVector->setNumElements(number);
5633 // Btran basic costs
5634 abcFactorization_->updateFullColumnTranspose(*arrayVector);
5635 double largestDifference = 0.0;
5636 int whichDifference = -1;
5637 if (type == 2) {
5638 for (int i = 0; i < numberRows_; i++) {
5639 double difference = fabs(dual_[i] - array[i]);
5640 if (difference > largestDifference) {
5641 largestDifference = difference;
5642 whichDifference = i;
5643 }
5644 }
5645 if (largestDifference > 1.0e-9)
5646 printf("Debug largest dual difference is %g on row %d dual_ %g computed %g\n",
5647 largestDifference, whichDifference, dual_[whichDifference], array[whichDifference]);
5648 }
5649 // now look at djs
5650 double *djs = usefulArray_[whichArrays[1]].denseVector();
5651 CoinAbcMemcpy(djs, abcCost_, numberTotal_);
5652 abcMatrix_->transposeTimesNonBasic(-1.0, array, djs);
5653 largestDifference = 0.0;
5654 whichDifference = -1;
5655 for (int i = 0; i < numberTotal_; i++) {
5656 Status status = getInternalStatus(i);
5657 double difference = fabs(abcDj_[i] - djs[i]);
5658 switch (status) {
5659 case basic:
5660 // probably not an error
5661 break;
5662 case AbcSimplex::isFixed:
5663 break;
5664 case isFree:
5665 case superBasic:
5666 case atUpperBound:
5667 case atLowerBound:
5668 if (difference > 1.0e-5 && numberTotal_ < 200)
5669 printf("Debug dj difference is %g on column %d abcDj_ %g computed %g\n",
5670 difference, i, abcDj_[i], djs[i]);
5671 if (difference > largestDifference) {
5672 largestDifference = difference;
5673 whichDifference = i;
5674 }
5675 break;
5676 }
5677 }
5678 if (largestDifference > 1.0e-9)
5679 printf("Debug largest dj difference is %g on column %d abcDj_ %g computed %g\n",
5680 largestDifference, whichDifference, abcDj_[whichDifference], djs[whichDifference]);
5681 CoinAbcMemset0(djs, numberTotal_);
5682 arrayVector->clear();
5683 for (int i = 0; i < 2; i++)
5684 setAvailableArray(whichArrays[i]);
5685 }
5686 int state[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
5687 int badT[8] = { 0, 0, 0, 0, 0, 0, 0, 0 };
5688 //int badP[8]={0,0,0,0,0,0,0,0};
5689 int numberInfeasibilities = 0;
5690 for (int i = 0; i < numberTotal_; i++) {
5691 int iStatus = internalStatus_[i] & 7;
5692 state[iStatus]++;
5693 Status status = getInternalStatus(i);
5694 double djValue = abcDj_[i];
5695 double value = abcSolution_[i];
5696 double lowerValue = abcLower_[i];
5697 double upperValue = abcUpper_[i];
5698 //double perturbationValue=abcPerturbation_[i];
5699 switch (status) {
5700 case basic:
5701 // probably not an error
5702 if (fabs(djValue) > currentDualTolerance_)
5703 badT[iStatus]++;
5704 //if(fabs(djValue)>perturbationValue)
5705 //badP[iStatus]++;
5706 break;
5707 case AbcSimplex::isFixed:
5708 break;
5709 case isFree:
5710 case superBasic:
5711 if (fabs(djValue) > currentDualTolerance_)
5712 badT[iStatus]++;
5713 //if(fabs(djValue)>perturbationValue)
5714 //badP[iStatus]++;
5715 break;
5716 case atUpperBound:
5717 if (fabs(value - upperValue) > primalTolerance_)
5718 numberInfeasibilities++;
5719 if (djValue > currentDualTolerance_)
5720 badT[iStatus]++;
5721 //if(djValue>perturbationValue)
5722 //badP[iStatus]++;
5723 break;
5724 case atLowerBound:
5725 if (fabs(value - lowerValue) > primalTolerance_)
5726 numberInfeasibilities++;
5727 if (-djValue > currentDualTolerance_)
5728 badT[iStatus]++;
5729 //if(-djValue>perturbationValue)
5730 //badP[iStatus]++;
5731 break;
5732 }
5733 }
5734 if (numberInfeasibilities)
5735 printf("Debug %d variables away from bound when should be\n", numberInfeasibilities);
5736 int numberBad = 0;
5737 for (int i = 0; i < 8; i++) {
5738 numberBad += badT[i] /*+badP[i]*/;
5739 }
5740 if (numberBad) {
5741 const char *type[] = { "atLowerBound",
5742 "atUpperBound",
5743 "??",
5744 "??",
5745 "isFree",
5746 "superBasic",
5747 "basic",
5748 "isFixed" };
5749 for (int i = 0; i < 8; i++) {
5750 if (state[i])
5751 printf("Debug - %s %d total %d bad with tolerance %d bad with perturbation\n",
5752 type[i], state[i], badT[i], 0 /*badP[i]*/);
5753 }
5754 }
5755 }
5756 #else
5757 // For debug - moves solution back to external and computes stuff
checkMoveBack(bool)5758 void AbcSimplex::checkMoveBack(bool)
5759 {
5760 }
5761 // For debug - checks solutionBasic
checkSolutionBasic() const5762 void AbcSimplex::checkSolutionBasic() const
5763 {
5764 }
5765
5766 // For debug - summarizes dj situation
checkDjs(int) const5767 void AbcSimplex::checkDjs(int) const
5768 {
5769 }
5770 #endif
5771 #ifndef EARLY_FACTORIZE
5772 #define ABC_NUMBER_USEFUL_NORMAL ABC_NUMBER_USEFUL
5773 #else
5774 #define ABC_NUMBER_USEFUL_NORMAL ABC_NUMBER_USEFUL - 1
5775 #endif
5776 static double *elAddress[ABC_NUMBER_USEFUL_NORMAL];
5777 // For debug - prints summary of arrays which are out of kilter
checkArrays(int ignoreEmpty) const5778 void AbcSimplex::checkArrays(int ignoreEmpty) const
5779 {
5780 if (!numberIterations_ || !elAddress[0]) {
5781 for (int i = 0; i < ABC_NUMBER_USEFUL_NORMAL; i++)
5782 elAddress[i] = usefulArray_[i].denseVector();
5783 } else {
5784 if (elAddress[0] != usefulArray_[0].denseVector()) {
5785 printf("elAddress not zero and does not match??\n");
5786 for (int i = 0; i < ABC_NUMBER_USEFUL_NORMAL; i++)
5787 elAddress[i] = usefulArray_[i].denseVector();
5788 }
5789 for (int i = 0; i < ABC_NUMBER_USEFUL_NORMAL; i++)
5790 assert(elAddress[i] == usefulArray_[i].denseVector());
5791 }
5792 for (int i = 0; i < ABC_NUMBER_USEFUL_NORMAL; i++) {
5793 int check = 1 << i;
5794 if (usefulArray_[i].getNumElements()) {
5795 if ((stateOfProblem_ & check) == 0)
5796 printf("Array %d has %d elements but says it is available!\n",
5797 i, usefulArray_[i].getNumElements());
5798 if (usefulArray_[i].getNumPartitions())
5799 usefulArray_[i].checkClean();
5800 else
5801 usefulArray_[i].CoinIndexedVector::checkClean();
5802 } else {
5803 if ((stateOfProblem_ & check) != 0 && (ignoreEmpty & check) == 0)
5804 printf("Array %d has no elements but says it is not available\n", i);
5805 if (usefulArray_[i].getNumPartitions())
5806 usefulArray_[i].checkClear();
5807 else
5808 usefulArray_[i].CoinIndexedVector::checkClear();
5809 }
5810 }
5811 }
5812 #include "CoinAbcCommon.hpp"
5813 #if 0
5814 #include "ClpFactorization.hpp"
5815 // Loads tolerances etc
5816 void
5817 ClpSimplex::loadTolerancesEtc(const AbcTolerancesEtc & data)
5818 {
5819 zeroTolerance_ = data.zeroTolerance_;
5820 primalToleranceToGetOptimal_ = data.primalToleranceToGetOptimal_;
5821 largeValue_ = data.largeValue_;
5822 alphaAccuracy_ = data.alphaAccuracy_;
5823 dualBound_ = data.dualBound_;
5824 dualTolerance_ = data.dualTolerance_;
5825 primalTolerance_ = data.primalTolerance_;
5826 infeasibilityCost_ = data.infeasibilityCost_;
5827 incomingInfeasibility_ = data.incomingInfeasibility_;
5828 allowedInfeasibility_ = data.allowedInfeasibility_;
5829 baseIteration_ = data.baseIteration_;
5830 numberRefinements_ = data.numberRefinements_;
5831 forceFactorization_ = data.forceFactorization_;
5832 perturbation_ = data.perturbation_;
5833 dontFactorizePivots_ = data.dontFactorizePivots_;
5834 /// For factorization
5835 abcFactorization_->maximumPivots(data.maximumPivots_);
5836 }
5837 // Loads tolerances etc
5838 void
5839 ClpSimplex::unloadTolerancesEtc(AbcTolerancesEtc & data)
5840 {
5841 data.zeroTolerance_ = zeroTolerance_;
5842 data.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
5843 data.largeValue_ = largeValue_;
5844 data.alphaAccuracy_ = alphaAccuracy_;
5845 data.dualBound_ = dualBound_;
5846 data.dualTolerance_ = dualTolerance_;
5847 data.primalTolerance_ = primalTolerance_;
5848 data.infeasibilityCost_ = infeasibilityCost_;
5849 data.incomingInfeasibility_ = incomingInfeasibility_;
5850 data.allowedInfeasibility_ = allowedInfeasibility_;
5851 data.baseIteration_ = baseIteration_;
5852 data.numberRefinements_ = numberRefinements_;
5853 data.forceFactorization_ = forceFactorization_;
5854 data.perturbation_ = perturbation_;
5855 data.dontFactorizePivots_ = dontFactorizePivots_;
5856 /// For factorization
5857 data.maximumPivots_ = abcFactorization_->maximumPivots();
5858 }
5859 // Loads tolerances etc
5860 void
5861 AbcSimplex::loadTolerancesEtc(const AbcTolerancesEtc & data)
5862 {
5863 zeroTolerance_ = data.zeroTolerance_;
5864 primalToleranceToGetOptimal_ = data.primalToleranceToGetOptimal_;
5865 largeValue_ = data.largeValue_;
5866 alphaAccuracy_ = data.alphaAccuracy_;
5867 dualBound_ = data.dualBound_;
5868 dualTolerance_ = data.dualTolerance_;
5869 primalTolerance_ = data.primalTolerance_;
5870 infeasibilityCost_ = data.infeasibilityCost_;
5871 incomingInfeasibility_ = data.incomingInfeasibility_;
5872 allowedInfeasibility_ = data.allowedInfeasibility_;
5873 baseIteration_ = data.baseIteration_;
5874 numberRefinements_ = data.numberRefinements_;
5875 forceFactorization_ = data.forceFactorization_;
5876 perturbation_ = data.perturbation_;
5877 dontFactorizePivots_ = data.dontFactorizePivots_;
5878 /// For factorization
5879 abcFactorization_->maximumPivots(data.maximumPivots_);
5880 }
5881 // Loads tolerances etc
5882 void
5883 AbcSimplex::unloadTolerancesEtc(AbcTolerancesEtc & data)
5884 {
5885 data.zeroTolerance_ = zeroTolerance_;
5886 data.primalToleranceToGetOptimal_ = primalToleranceToGetOptimal_;
5887 data.largeValue_ = largeValue_;
5888 data.alphaAccuracy_ = alphaAccuracy_;
5889 data.dualBound_ = dualBound_;
5890 data.dualTolerance_ = dualTolerance_;
5891 data.primalTolerance_ = primalTolerance_;
5892 data.infeasibilityCost_ = infeasibilityCost_;
5893 data.incomingInfeasibility_ = incomingInfeasibility_;
5894 data.allowedInfeasibility_ = allowedInfeasibility_;
5895 data.baseIteration_ = baseIteration_;
5896 data.numberRefinements_ = numberRefinements_;
5897 data.forceFactorization_ = forceFactorization_;
5898 data.perturbation_ = perturbation_;
5899 data.dontFactorizePivots_ = dontFactorizePivots_;
5900 /// For factorization
5901 data.maximumPivots_ = abcFactorization_->maximumPivots();
5902 }
5903 #endif
5904 // Swaps two variables (for now just updates basic list) and sets status
swap(int pivotRow,int nonBasicPosition,Status newStatus)5905 void AbcSimplex::swap(int pivotRow, int nonBasicPosition, Status newStatus)
5906 {
5907 // set status
5908 setInternalStatus(nonBasicPosition, newStatus);
5909 solutionBasic_[pivotRow] = abcSolution_[nonBasicPosition];
5910 lowerBasic_[pivotRow] = abcLower_[nonBasicPosition];
5911 upperBasic_[pivotRow] = abcUpper_[nonBasicPosition];
5912 costBasic_[pivotRow] = abcCost_[nonBasicPosition];
5913 }
5914 // Swaps two variables (for now just updates basic list)
swap(int pivotRow,int nonBasicPosition)5915 void AbcSimplex::swap(int pivotRow, int nonBasicPosition)
5916 {
5917 solutionBasic_[pivotRow] = abcSolution_[nonBasicPosition];
5918 lowerBasic_[pivotRow] = lowerSaved_[nonBasicPosition];
5919 upperBasic_[pivotRow] = upperSaved_[nonBasicPosition];
5920 costBasic_[pivotRow] = abcCost_[nonBasicPosition];
5921 }
5922 // Move status and solution to ClpSimplex
moveStatusToClp(ClpSimplex * clpModel)5923 void AbcSimplex::moveStatusToClp(ClpSimplex *clpModel)
5924 {
5925 assert(clpModel);
5926 if (algorithm_ < 0)
5927 clpModel->setObjectiveValue(clpObjectiveValue());
5928 else
5929 clpModel->setObjectiveValue(objectiveValue());
5930 clpModel->setProblemStatus(problemStatus_);
5931 clpModel->setSecondaryStatus(secondaryStatus_);
5932 clpModel->setNumberIterations(numberIterations_);
5933 clpModel->setSumDualInfeasibilities(sumDualInfeasibilities_);
5934 clpModel->setSumOfRelaxedDualInfeasibilities(sumOfRelaxedDualInfeasibilities_);
5935 clpModel->setNumberDualInfeasibilities(numberDualInfeasibilities_);
5936 clpModel->setSumPrimalInfeasibilities(sumPrimalInfeasibilities_);
5937 clpModel->setSumOfRelaxedPrimalInfeasibilities(sumOfRelaxedPrimalInfeasibilities_);
5938 clpModel->setNumberPrimalInfeasibilities(numberPrimalInfeasibilities_);
5939 permuteOut(ROW_PRIMAL_OK | ROW_DUAL_OK | COLUMN_PRIMAL_OK | COLUMN_DUAL_OK | ALL_STATUS_OK);
5940 CoinAbcMemcpy(clpModel->primalColumnSolution(), primalColumnSolution(), numberColumns_);
5941 CoinAbcMemcpy(clpModel->dualColumnSolution(), dualColumnSolution(), numberColumns_);
5942 CoinAbcMemcpy(clpModel->primalRowSolution(), primalRowSolution(), numberRows_);
5943 CoinAbcMemcpy(clpModel->dualRowSolution(), dualRowSolution(), numberRows_);
5944 CoinAbcMemcpy(clpModel->statusArray(), statusArray(), numberTotal_);
5945 }
5946 // Move status and solution from ClpSimplex
moveStatusFromClp(ClpSimplex * clpModel)5947 void AbcSimplex::moveStatusFromClp(ClpSimplex *clpModel)
5948 {
5949 assert(clpModel);
5950 problemStatus_ = clpModel->problemStatus();
5951 secondaryStatus_ = clpModel->secondaryStatus();
5952 numberIterations_ = clpModel->numberIterations();
5953 sumDualInfeasibilities_ = clpModel->sumDualInfeasibilities();
5954 sumOfRelaxedDualInfeasibilities_ = clpModel->sumOfRelaxedDualInfeasibilities();
5955 numberDualInfeasibilities_ = clpModel->numberDualInfeasibilities();
5956 sumPrimalInfeasibilities_ = clpModel->sumPrimalInfeasibilities();
5957 sumOfRelaxedPrimalInfeasibilities_ = clpModel->sumOfRelaxedPrimalInfeasibilities();
5958 numberPrimalInfeasibilities_ = clpModel->numberPrimalInfeasibilities();
5959 CoinAbcMemcpy(primalColumnSolution(), clpModel->primalColumnSolution(), numberColumns_);
5960 CoinAbcMemcpy(dualColumnSolution(), clpModel->dualColumnSolution(), numberColumns_);
5961 CoinAbcMemcpy(primalRowSolution(), clpModel->primalRowSolution(), numberRows_);
5962 CoinAbcMemcpy(dualRowSolution(), clpModel->dualRowSolution(), numberRows_);
5963 CoinAbcMemcpy(statusArray(), clpModel->statusArray(), numberTotal_);
5964 translate(DO_SCALE_AND_MATRIX | DO_BASIS_AND_ORDER | DO_STATUS | DO_SOLUTION);
5965 }
5966 // Clears an array and says available (-1 does all)
clearArrays(int which)5967 void AbcSimplex::clearArrays(int which)
5968 {
5969 if (which >= 0) {
5970 if (usefulArray_[which].getNumElements())
5971 usefulArray_[which].clear();
5972 int check = 1 << which;
5973 stateOfProblem_ &= ~check;
5974 } else if (which == -1) {
5975 for (int i = 0; i < ABC_NUMBER_USEFUL_NORMAL; i++) {
5976 if (usefulArray_[i].getNumElements()) {
5977 usefulArray_[i].clear();
5978 usefulArray_[i].clearAndReset();
5979 }
5980 }
5981 stateOfProblem_ &= ~255;
5982 } else {
5983 for (int i = 0; i < ABC_NUMBER_USEFUL_NORMAL; i++) {
5984 usefulArray_[i].clearAndReset();
5985 }
5986 stateOfProblem_ &= ~255;
5987 }
5988 }
5989 // Clears an array and says available
clearArrays(CoinPartitionedVector * which)5990 void AbcSimplex::clearArrays(CoinPartitionedVector *which)
5991 {
5992 int check = 1;
5993 int i;
5994 for (i = 0; i < ABC_NUMBER_USEFUL_NORMAL; i++) {
5995 if (&usefulArray_[i] == which) {
5996 usefulArray_[i].clear();
5997 stateOfProblem_ &= ~check;
5998 break;
5999 }
6000 check = check << 1;
6001 }
6002 assert(i < ABC_NUMBER_USEFUL_NORMAL);
6003 }
6004 // set multiple sequence in
setMultipleSequenceIn(int sequenceIn[4])6005 void AbcSimplex::setMultipleSequenceIn(int sequenceIn[4])
6006 {
6007 memcpy(multipleSequenceIn_, sequenceIn, sizeof(multipleSequenceIn_));
6008 }
6009 // Returns first available empty array (and sets flag)
getAvailableArray() const6010 int AbcSimplex::getAvailableArray() const
6011 {
6012 int which;
6013 int status = 1;
6014 for (which = 0; which < ABC_NUMBER_USEFUL_NORMAL; which++) {
6015 if ((stateOfProblem_ & status) == 0)
6016 break;
6017 status = status << 1;
6018 }
6019 assert(which < ABC_NUMBER_USEFUL_NORMAL);
6020 assert(!usefulArray_[which].getNumElements());
6021 stateOfProblem_ |= status;
6022 return which;
6023 }
atFakeBound(int sequence) const6024 bool AbcSimplex::atFakeBound(int sequence) const
6025 {
6026 FakeBound status = getFakeBound(sequence);
6027 bool atBound = false;
6028 switch (status) {
6029 case noFake:
6030 break;
6031 case lowerFake:
6032 atBound = (getInternalStatus(sequence) == atLowerBound);
6033 break;
6034 case upperFake:
6035 atBound = (getInternalStatus(sequence) == atUpperBound);
6036 break;
6037 case bothFake:
6038 atBound = true;
6039 break;
6040 }
6041 return atBound;
6042 }
AbcSimplexProgress()6043 AbcSimplexProgress::AbcSimplexProgress()
6044 : ClpSimplexProgress()
6045 {
6046 }
6047
6048 //-----------------------------------------------------------------------------
6049
~AbcSimplexProgress()6050 AbcSimplexProgress::~AbcSimplexProgress()
6051 {
6052 }
6053 // Copy constructor from model
AbcSimplexProgress(ClpSimplex * rhs)6054 AbcSimplexProgress::AbcSimplexProgress(ClpSimplex *rhs)
6055 : ClpSimplexProgress(rhs)
6056 {
6057 }
6058 // Assignment operator. This copies the data
6059 AbcSimplexProgress &
operator =(const AbcSimplexProgress & rhs)6060 AbcSimplexProgress::operator=(const AbcSimplexProgress &rhs)
6061 {
6062 if (this != &rhs) {
6063 ClpSimplexProgress::operator=(rhs);
6064 }
6065 return *this;
6066 }
looping()6067 int AbcSimplexProgress::looping()
6068 {
6069 if (!model_)
6070 return -1;
6071 #ifdef ABC_INHERIT
6072 AbcSimplex *model = model_->abcSimplex();
6073 #else
6074 AbcSimplex *model;
6075 return -1;
6076 #endif
6077 double objective;
6078 if (model_->algorithm() < 0) {
6079 objective = model_->rawObjectiveValue();
6080 objective -= model_->bestPossibleImprovement();
6081 } else {
6082 objective = model->abcNonLinearCost()->feasibleReportCost();
6083 }
6084 double infeasibility;
6085 double realInfeasibility = 0.0;
6086 int numberInfeasibilities;
6087 int iterationNumber = model->numberIterations();
6088 //numberTimesFlagged_ = 0;
6089 if (model->algorithm() < 0) {
6090 // dual
6091 infeasibility = model->sumPrimalInfeasibilities();
6092 numberInfeasibilities = model->numberPrimalInfeasibilities();
6093 } else {
6094 //primal
6095 infeasibility = model->sumDualInfeasibilities();
6096 realInfeasibility = model->abcNonLinearCost()->sumInfeasibilities();
6097 numberInfeasibilities = model->numberDualInfeasibilities();
6098 }
6099 int i;
6100 int numberMatched = 0;
6101 int matched = 0;
6102 int nsame = 0;
6103 for (i = 0; i < CLP_PROGRESS; i++) {
6104 bool matchedOnObjective = objective == objective_[i];
6105 bool matchedOnInfeasibility = infeasibility == infeasibility_[i];
6106 bool matchedOnInfeasibilities = (numberInfeasibilities == numberInfeasibilities_[i]);
6107
6108 if (matchedOnObjective && matchedOnInfeasibility && matchedOnInfeasibilities) {
6109 matched |= (1 << i);
6110 // Check not same iteration
6111 if (iterationNumber != iterationNumber_[i]) {
6112 numberMatched++;
6113 #if ABC_NORMAL_DEBUG > 0
6114 // here mainly to get over compiler bug?
6115 if (model->messageHandler()->logLevel() > 10)
6116 printf("%d %d %d %d %d loop check\n", i, numberMatched,
6117 matchedOnObjective, matchedOnInfeasibility,
6118 matchedOnInfeasibilities);
6119 #endif
6120 } else {
6121 // stuck but code should notice
6122 nsame++;
6123 }
6124 }
6125 if (i) {
6126 objective_[i - 1] = objective_[i];
6127 infeasibility_[i - 1] = infeasibility_[i];
6128 realInfeasibility_[i - 1] = realInfeasibility_[i];
6129 numberInfeasibilities_[i - 1] = numberInfeasibilities_[i];
6130 iterationNumber_[i - 1] = iterationNumber_[i];
6131 }
6132 }
6133 objective_[CLP_PROGRESS - 1] = objective;
6134 infeasibility_[CLP_PROGRESS - 1] = infeasibility;
6135 realInfeasibility_[CLP_PROGRESS - 1] = realInfeasibility;
6136 numberInfeasibilities_[CLP_PROGRESS - 1] = numberInfeasibilities;
6137 iterationNumber_[CLP_PROGRESS - 1] = iterationNumber;
6138 if (nsame == CLP_PROGRESS)
6139 numberMatched = CLP_PROGRESS; // really stuck
6140 if (model->progressFlag())
6141 numberMatched = 0;
6142 numberTimes_++;
6143 if (numberTimes_ < 10)
6144 numberMatched = 0;
6145 // skip if just last time as may be checking something
6146 if (matched == (1 << (CLP_PROGRESS - 1)))
6147 numberMatched = 0;
6148 if (numberMatched) {
6149 model->messageHandler()->message(CLP_POSSIBLELOOP, model->messages())
6150 << numberMatched
6151 << matched
6152 << numberTimes_
6153 << CoinMessageEol;
6154 printf("loop detected %d times out of %d\n", numberBadTimes_, numberTimes_);
6155 numberBadTimes_++;
6156 if (numberBadTimes_ < 10) {
6157 // make factorize every iteration
6158 model->forceFactorization(1);
6159 if (numberBadTimes_ < 2) {
6160 startCheck(); // clear other loop check
6161 if (model->algorithm() < 0) {
6162 // dual - change tolerance
6163 model->setCurrentDualTolerance(model->currentDualTolerance() * 1.05);
6164 // if infeasible increase dual bound
6165 if (model->currentDualBound() < 1.0e17) {
6166 model->setDualBound(model->currentDualBound() * 1.1);
6167 static_cast< AbcSimplexDual * >(model)->bounceTolerances(100);
6168 }
6169 } else {
6170 // primal - change tolerance
6171 if (numberBadTimes_ > 3)
6172 model->setCurrentPrimalTolerance(model->currentPrimalTolerance() * 1.05);
6173 // if infeasible increase infeasibility cost
6174 //if (model->nonLinearCost()->numberInfeasibilities() &&
6175 // model->infeasibilityCost() < 1.0e17) {
6176 // model->setInfeasibilityCost(model->infeasibilityCost() * 1.1);
6177 //}
6178 }
6179 } else {
6180 // flag
6181 int iSequence;
6182 if (model->algorithm() < 0) {
6183 // dual
6184 if (model->currentDualBound() > 1.0e14)
6185 model->setDualBound(1.0e14);
6186 iSequence = in_[CLP_CYCLE - 1];
6187 } else {
6188 // primal
6189 if (model->infeasibilityCost() > 1.0e14)
6190 model->setInfeasibilityCost(1.0e14);
6191 iSequence = out_[CLP_CYCLE - 1];
6192 }
6193 if (iSequence >= 0) {
6194 char x = model->isColumn(iSequence) ? 'C' : 'R';
6195 if (model->messageHandler()->logLevel() >= 63)
6196 model->messageHandler()->message(CLP_SIMPLEX_FLAG, model->messages())
6197 << x << model->sequenceWithin(iSequence)
6198 << CoinMessageEol;
6199 // if Gub then needs to be sequenceIn_
6200 int save = model->sequenceIn();
6201 model->setSequenceIn(iSequence);
6202 model->setFlagged(iSequence);
6203 model->setLastBadIteration(model->numberIterations());
6204 model->setSequenceIn(save);
6205 //printf("flagging %d from loop\n",iSequence);
6206 startCheck();
6207 } else {
6208 // Give up
6209 #if ABC_NORMAL_DEBUG > 0
6210 if (model->messageHandler()->logLevel() >= 63)
6211 printf("***** All flagged?\n");
6212 #endif
6213 return 4;
6214 }
6215 // reset
6216 numberBadTimes_ = 2;
6217 }
6218 return -2;
6219 } else {
6220 // look at solution and maybe declare victory
6221 if (infeasibility < 1.0e-4) {
6222 return 0;
6223 } else {
6224 model->messageHandler()->message(CLP_LOOP, model->messages())
6225 << CoinMessageEol;
6226 #ifndef NDEBUG
6227 printf("debug loop AbcSimplex A\n");
6228 abort();
6229 #endif
6230 return 3;
6231 }
6232 }
6233 }
6234 return -1;
6235 }
6236 #if 0
6237 void
6238 AbcSimplex::loadProblem ( const CoinPackedMatrix& matrix,
6239 const double* collb, const double* colub,
6240 const double* obj,
6241 const double* rowlb, const double* rowub,
6242 const double * rowObjective)
6243 {
6244 ClpSimplex::loadProblem(matrix, collb, colub, obj, rowlb, rowub,
6245 rowObjective);
6246 translate(DO_SCALE_AND_MATRIX|DO_BASIS_AND_ORDER|DO_STATUS|DO_SOLUTION);
6247 }
6248 #endif
6249 // For debug - check pivotVariable consistent
checkConsistentPivots() const6250 void AbcSimplex::checkConsistentPivots() const
6251 {
6252 unsigned char *copyStatus = CoinCopyOfArray(internalStatus_, numberTotal_);
6253 int nBad = 0;
6254 for (int i = 0; i < numberRows_; i++) {
6255 int k = abcPivotVariable_[i];
6256 if (k >= 0 && k < numberTotal_) {
6257 if ((copyStatus[k] & 7) != 6) {
6258 nBad++;
6259 printf("%d pivoting on row %d is not basic - status %d\n", k, i, copyStatus[k] & 7);
6260 }
6261 copyStatus[k] = 16;
6262 } else {
6263 nBad++;
6264 printf("%d pivoting on row %d is bad\n", k, i);
6265 }
6266 }
6267 for (int i = 0; i < numberTotal_; i++) {
6268 if ((copyStatus[i] & 7) == 6) {
6269 nBad++;
6270 printf("%d without pivot is basic\n", i);
6271 }
6272 }
6273 assert(!nBad);
6274 delete[] copyStatus;
6275 }
6276 // Print stuff
printStuff() const6277 void AbcSimplex::printStuff() const
6278 {
6279 if (numberIterations_ != 2821)
6280 return;
6281 printf("XXStart\n");
6282 for (int i = 0; i < this->numberTotal(); i++) {
6283 if (this->getInternalStatus(i) != AbcSimplex::basic)
6284 printf("%d status %d primal %g dual %g lb %g ub %g\n",
6285 i, this->getInternalStatus(i), this->solutionRegion()[i],
6286 this->djRegion()[i], this->lowerRegion()[i], this->upperRegion()[i]);
6287 }
6288 for (int i = 0; i < this->numberRows(); i++) {
6289 printf("%d %g <= %g <= %g -pivot %d cost %g\n",
6290 i, this->lowerBasic()[i], this->solutionBasic()[i],
6291 this->upperBasic()[i], this->pivotVariable()[i], this->costBasic()[i]);
6292 }
6293 printf("XXend\n");
6294 }
6295 #if ABC_PARALLEL == 1
6296 // so thread can find out which one it is
whichThread() const6297 int AbcSimplex::whichThread() const
6298 {
6299 pthread_t thisThread = pthread_self();
6300 int whichThread;
6301 for (whichThread = 0; whichThread < NUMBER_THREADS; whichThread++) {
6302 if (pthread_equal(thisThread, abcThread_[whichThread]))
6303 break;
6304 }
6305 assert(whichThread < NUMBER_THREADS);
6306 return whichThread;
6307 }
6308 #endif
6309
6310 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
6311 */
6312