1 /* $Id: AbcSimplexPrimal.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
2 // Copyright (C) 2002, 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 /* Notes on implementation of primal simplex algorithm.
7
8 When primal feasible(A):
9
10 If dual feasible, we are optimal. Otherwise choose an infeasible
11 basic variable to enter basis from a bound (B). We now need to find an
12 outgoing variable which will leave problem primal feasible so we get
13 the column of the tableau corresponding to the incoming variable
14 (with the correct sign depending if variable will go up or down).
15
16 We now perform a ratio test to determine which outgoing variable will
17 preserve primal feasibility (C). If no variable found then problem
18 is unbounded (in primal sense). If there is a variable, we then
19 perform pivot and repeat. Trivial?
20
21 -------------------------------------------
22
23 A) How do we get primal feasible? All variables have fake costs
24 outside their feasible region so it is trivial to declare problem
25 feasible. OSL did not have a phase 1/phase 2 approach but
26 instead effectively put an extra cost on infeasible basic variables
27 I am taking the same approach here, although it is generalized
28 to allow for non-linear costs and dual information.
29
30 In OSL, this weight was changed heuristically, here at present
31 it is only increased if problem looks finished. If problem is
32 feasible I check for unboundedness. If not unbounded we
33 could play with going into dual. As long as weights increase
34 any algorithm would be finite.
35
36 B) Which incoming variable to choose is a virtual base class.
37 For difficult problems steepest edge is preferred while for
38 very easy (large) problems we will need partial scan.
39
40 C) Sounds easy, but this is hardest part of algorithm.
41 1) Instead of stopping at first choice, we may be able
42 to allow that variable to go through bound and if objective
43 still improving choose again. These mini iterations can
44 increase speed by orders of magnitude but we may need to
45 go to more of a bucket choice of variable rather than looking
46 at them one by one (for speed).
47 2) Accuracy. Basic infeasibilities may be less than
48 tolerance. Pivoting on these makes objective go backwards.
49 OSL modified cost so a zero move was made, Gill et al
50 modified so a strictly positive move was made.
51 The two problems are that re-factorizations can
52 change rinfeasibilities above and below tolerances and that when
53 finished we need to reset costs and try again.
54 3) Degeneracy. Gill et al helps but may not be enough. We
55 may need more. Also it can improve speed a lot if we perturb
56 the rhs and bounds significantly.
57
58 References:
59 Forrest and Goldfarb, Steepest-edge simplex algorithms for
60 linear programming - Mathematical Programming 1992
61 Forrest and Tomlin, Implementing the simplex method for
62 the Optimization Subroutine Library - IBM Systems Journal 1992
63 Gill, Murray, Saunders, Wright A Practical Anti-Cycling
64 Procedure for Linear and Nonlinear Programming SOL report 1988
65
66
67 TODO:
68
69 a) Better recovery procedures. At present I never check on forward
70 progress. There is checkpoint/restart with reducing
71 re-factorization frequency, but this is only on singular
72 factorizations.
73 b) Fast methods for large easy problems (and also the option for
74 the code to automatically choose which method).
75 c) We need to be able to stop in various ways for OSI - this
76 is fairly easy.
77
78 */
79
80 #include "CoinPragma.hpp"
81
82 #include <math.h>
83
84 #include "CoinHelperFunctions.hpp"
85 #include "CoinAbcHelperFunctions.hpp"
86 #include "AbcSimplexPrimal.hpp"
87 #include "AbcSimplexFactorization.hpp"
88 #include "AbcNonLinearCost.hpp"
89 #include "CoinPackedMatrix.hpp"
90 #include "CoinIndexedVector.hpp"
91 #include "AbcPrimalColumnPivot.hpp"
92 #include "ClpMessage.hpp"
93 #include "ClpEventHandler.hpp"
94 #include <cfloat>
95 #include <cassert>
96 #include <string>
97 #include <stdio.h>
98 #include <iostream>
99 #ifdef CLP_USER_DRIVEN1
100 /* Returns true if variable sequenceOut can leave basis when
101 model->sequenceIn() enters.
102 This function may be entered several times for each sequenceOut.
103 The first time realAlpha will be positive if going to lower bound
104 and negative if going to upper bound (scaled bounds in lower,upper) - then will be zero.
105 currentValue is distance to bound.
106 currentTheta is current theta.
107 alpha is fabs(pivot element).
108 Variable will change theta if currentValue - currentTheta*alpha < 0.0
109 */
110 bool userChoiceValid1(const AbcSimplex *model,
111 int sequenceOut,
112 double currentValue,
113 double currentTheta,
114 double alpha,
115 double realAlpha);
116 /* This returns true if chosen in/out pair valid.
117 The main thing to check would be variable flipping bounds may be
118 OK. This would be signaled by reasonable theta_ and valueOut_.
119 If you return false sequenceIn_ will be flagged as ineligible.
120 */
121 bool userChoiceValid2(const AbcSimplex *model);
122 /* If a good pivot then you may wish to unflag some variables.
123 */
124 void userChoiceWasGood(AbcSimplex *model);
125 #endif
126 #ifdef TRY_SPLIT_VALUES_PASS
127 static double valuesChunk = -10.0;
128 static double valuesRatio = 0.1;
129 static int valuesStop = -1;
130 static int keepValuesPass = -1;
131 static int doOrdinaryVariables = -1;
132 #endif
133 // primal
primal(int ifValuesPass,int)134 int AbcSimplexPrimal::primal(int ifValuesPass, int /*startFinishOptions*/)
135 {
136
137 /*
138 Method
139
140 It tries to be a single phase approach with a weight of 1.0 being
141 given to getting optimal and a weight of infeasibilityCost_ being
142 given to getting primal feasible. In this version I have tried to
143 be clever in a stupid way. The idea of fake bounds in dual
144 seems to work so the primal analogue would be that of getting
145 bounds on reduced costs (by a presolve approach) and using
146 these for being above or below feasible region. I decided to waste
147 memory and keep these explicitly. This allows for non-linear
148 costs!
149
150 The code is designed to take advantage of sparsity so arrays are
151 seldom zeroed out from scratch or gone over in their entirety.
152 The only exception is a full scan to find incoming variable for
153 Dantzig row choice. For steepest edge we keep an updated list
154 of dual infeasibilities (actually squares).
155 On easy problems we don't need full scan - just
156 pick first reasonable.
157
158 One problem is how to tackle degeneracy and accuracy. At present
159 I am using the modification of costs which I put in OSL and which was
160 extended by Gill et al. I am still not sure of the exact details.
161
162 The flow of primal is three while loops as follows:
163
164 while (not finished) {
165
166 while (not clean solution) {
167
168 Factorize and/or clean up solution by changing bounds so
169 primal feasible. If looks finished check fake primal bounds.
170 Repeat until status is iterating (-1) or finished (0,1,2)
171
172 }
173
174 while (status==-1) {
175
176 Iterate until no pivot in or out or time to re-factorize.
177
178 Flow is:
179
180 choose pivot column (incoming variable). if none then
181 we are primal feasible so looks as if done but we need to
182 break and check bounds etc.
183
184 Get pivot column in tableau
185
186 Choose outgoing row. If we don't find one then we look
187 primal unbounded so break and check bounds etc. (Also the
188 pivot tolerance is larger after any iterations so that may be
189 reason)
190
191 If we do find outgoing row, we may have to adjust costs to
192 keep going forwards (anti-degeneracy). Check pivot will be stable
193 and if unstable throw away iteration and break to re-factorize.
194 If minor error re-factorize after iteration.
195
196 Update everything (this may involve changing bounds on
197 variables to stay primal feasible.
198
199 }
200
201 }
202
203 At present we never check we are going forwards. I overdid that in
204 OSL so will try and make a last resort.
205
206 Needs partial scan pivot in option.
207
208 May need other anti-degeneracy measures, especially if we try and use
209 loose tolerances as a way to solve in fewer iterations.
210
211 I like idea of dynamic scaling. This gives opportunity to decouple
212 different implications of scaling for accuracy, iteration count and
213 feasibility tolerance.
214
215 */
216
217 algorithm_ = +1;
218 moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
219 #if ABC_PARALLEL > 0
220 if (parallelMode()) {
221 // extra regions
222 // for moment allow for ordered factorization
223 int length = 2 * numberRows_ + abcFactorization_->maximumPivots();
224 for (int i = 0; i < 6; i++) {
225 delete rowArray_[i];
226 rowArray_[i] = new CoinIndexedVector();
227 rowArray_[i]->reserve(length);
228 delete columnArray_[i];
229 columnArray_[i] = new CoinIndexedVector();
230 columnArray_[i]->reserve(length);
231 }
232 }
233 #endif
234 // save data
235 ClpDataSave data = saveData();
236 dualTolerance_ = dblParam_[ClpDualTolerance];
237 primalTolerance_ = dblParam_[ClpPrimalTolerance];
238 currentDualTolerance_ = dualTolerance_;
239 //nextCleanNonBasicIteration_=baseIteration_+numberRows_;
240 // Save so can see if doing after dual
241 int initialStatus = problemStatus_;
242 int initialIterations = numberIterations_;
243 int initialNegDjs = -1;
244 // copy bounds to perturbation
245 CoinAbcMemcpy(abcPerturbation_, abcLower_, numberTotal_);
246 CoinAbcMemcpy(abcPerturbation_ + numberTotal_, abcUpper_, numberTotal_);
247 #if 0
248 if (numberRows_>80000&&numberRows_<90000) {
249 FILE * fp = fopen("save.stuff", "rb");
250 if (fp) {
251 fread(internalStatus_,sizeof(char),numberTotal_,fp);
252 fread(abcSolution_,sizeof(double),numberTotal_,fp);
253 } else {
254 printf("can't open save.stuff\n");
255 abort();
256 }
257 fclose(fp);
258 }
259 #endif
260 // initialize - maybe values pass and algorithm_ is +1
261 /// Initial coding here
262 if (nonLinearCost_ == NULL && algorithm_ > 0) {
263 // get a valid nonlinear cost function
264 abcNonLinearCost_ = new AbcNonLinearCost(this);
265 }
266 statusOfProblemInPrimal(0);
267 /*if (!startup(ifValuesPass))*/ {
268
269 // Set average theta
270 abcNonLinearCost_->setAverageTheta(1.0e3);
271
272 // Say no pivot has occurred (for steepest edge and updates)
273 pivotRow_ = -2;
274
275 // This says whether to restore things etc
276 int factorType = 0;
277 if (problemStatus_ < 0 && perturbation_ < 100 && !ifValuesPass) {
278 perturb(0);
279 if (perturbation_ != 100) {
280 // Can't get here if values pass
281 assert(!ifValuesPass);
282 gutsOfPrimalSolution(3);
283 if (handler_->logLevel() > 2) {
284 handler_->message(CLP_SIMPLEX_STATUS, messages_)
285 << numberIterations_ << objectiveValue();
286 handler_->printing(sumPrimalInfeasibilities_ > 0.0)
287 << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
288 handler_->printing(sumDualInfeasibilities_ > 0.0)
289 << sumDualInfeasibilities_ << numberDualInfeasibilities_;
290 handler_->printing(numberDualInfeasibilitiesWithoutFree_
291 < numberDualInfeasibilities_)
292 << numberDualInfeasibilitiesWithoutFree_;
293 handler_->message() << CoinMessageEol;
294 }
295 }
296 }
297 // Start check for cycles
298 abcProgress_.fillFromModel(this);
299 abcProgress_.startCheck();
300 bool needNewWeights = false;
301 double pivotTolerance = abcFactorization_->pivotTolerance();
302 /*
303 Status of problem:
304 0 - optimal
305 1 - infeasible
306 2 - unbounded
307 -1 - iterating
308 -2 - factorization wanted
309 -3 - redo checking without factorization
310 -4 - looks infeasible
311 -5 - looks unbounded
312 */
313 while (problemStatus_ < 0) {
314 // clear all arrays
315 clearArrays(-1);
316 // If getting nowhere - why not give it a kick
317 #if 1
318 if (perturbation_ < 101 && numberIterations_ == 6078 /*> 2 * (numberRows_ + numberColumns_)*/ && (specialOptions_ & 4) == 0
319 && initialStatus != 10) {
320 perturb(1);
321 }
322 #endif
323 // If we have done no iterations - special
324 if (lastGoodIteration_ == numberIterations_ && factorType)
325 factorType = 3;
326 if (pivotTolerance < abcFactorization_->pivotTolerance()) {
327 // force factorization
328 pivotTolerance = abcFactorization_->pivotTolerance();
329 factorType = 5;
330 }
331
332 // may factorize, checks if problem finished
333 if (factorType)
334 statusOfProblemInPrimal(factorType + (needNewWeights ? 10 : 0));
335 needNewWeights = false;
336 if (problemStatus_ == 10 && (moreSpecialOptions_ & 2048) != 0) {
337 problemStatus_ = 0;
338 }
339 if (initialStatus == 10) {
340 initialStatus = -1;
341 // cleanup phase
342 if (initialIterations != numberIterations_) {
343 if (numberDualInfeasibilities_ > 10000 && numberDualInfeasibilities_ > 10 * initialNegDjs) {
344 // getting worse - try perturbing
345 if (perturbation_ < 101 && (specialOptions_ & 4) == 0) {
346 perturb(1);
347 statusOfProblemInPrimal(factorType);
348 }
349 }
350 } else {
351 // save number of negative djs
352 if (!numberPrimalInfeasibilities_)
353 initialNegDjs = numberDualInfeasibilities_;
354 // make sure weight won't be changed
355 if (infeasibilityCost_ == 1.0e10)
356 infeasibilityCost_ = 1.000001e10;
357 }
358 }
359
360 // Say good factorization
361 factorType = 1;
362
363 // Say no pivot has occurred (for steepest edge and updates)
364 pivotRow_ = -2;
365
366 // exit if victory declared
367 if (problemStatus_ >= 0) {
368 #ifdef CLP_USER_DRIVEN
369 int status = eventHandler_->event(ClpEventHandler::endInPrimal);
370 if (status >= 0 && status < 10) {
371 // carry on
372 problemStatus_ = -1;
373 if (status == 0)
374 continue; // re-factorize
375 } else if (status >= 10) {
376 problemStatus_ = status - 10;
377 break;
378 } else {
379 break;
380 }
381 #else
382 break;
383 #endif
384 }
385
386 // test for maximum iterations
387 if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) {
388 problemStatus_ = 3;
389 break;
390 }
391
392 if (firstFree_ < 0) {
393 if (ifValuesPass) {
394 // end of values pass
395 ifValuesPass = 0;
396 needNewWeights = true;
397 int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
398 if (status >= 0) {
399 problemStatus_ = 5;
400 secondaryStatus_ = ClpEventHandler::endOfValuesPass;
401 break;
402 }
403 //#define FEB_TRY
404 #if 1 //def FEB_TRY
405 if (perturbation_ < 100
406 #ifdef TRY_SPLIT_VALUES_PASS
407 && valuesStop < 0
408 #endif
409 )
410 perturb(0);
411 #endif
412 }
413 }
414 if (abcNonLinearCost_->numberInfeasibilities() > 0 && !abcProgress_.initialWeight_ && !ifValuesPass && infeasibilityCost_ == 1.0e10
415 #ifdef TRY_SPLIT_VALUES_PASS
416 && valuesStop < 0
417 #endif
418 ) {
419 // first time infeasible - start up weight computation
420 // compute with original costs
421 CoinAbcMemcpy(djSaved_, abcDj_, numberTotal_);
422 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
423 CoinAbcGatherFrom(abcCost_, costBasic_, abcPivotVariable_, numberRows_);
424 gutsOfPrimalSolution(1);
425 int numberSame = 0;
426 int numberDifferent = 0;
427 int numberZero = 0;
428 int numberFreeSame = 0;
429 int numberFreeDifferent = 0;
430 int numberFreeZero = 0;
431 int n = 0;
432 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
433 if (getInternalStatus(iSequence) != basic && !flagged(iSequence)) {
434 // not basic
435 double distanceUp = abcUpper_[iSequence] - abcSolution_[iSequence];
436 double distanceDown = abcSolution_[iSequence] - abcLower_[iSequence];
437 double feasibleDj = abcDj_[iSequence];
438 double infeasibleDj = djSaved_[iSequence] - feasibleDj;
439 double value = feasibleDj * infeasibleDj;
440 if (distanceUp > primalTolerance_) {
441 // Check if "free"
442 if (distanceDown > primalTolerance_) {
443 // free
444 if (value > dualTolerance_) {
445 numberFreeSame++;
446 } else if (value < -dualTolerance_) {
447 numberFreeDifferent++;
448 abcDj_[n++] = feasibleDj / infeasibleDj;
449 } else {
450 numberFreeZero++;
451 }
452 } else {
453 // should not be negative
454 if (value > dualTolerance_) {
455 numberSame++;
456 } else if (value < -dualTolerance_) {
457 numberDifferent++;
458 abcDj_[n++] = feasibleDj / infeasibleDj;
459 } else {
460 numberZero++;
461 }
462 }
463 } else if (distanceDown > primalTolerance_) {
464 // should not be positive
465 if (value > dualTolerance_) {
466 numberSame++;
467 } else if (value < -dualTolerance_) {
468 numberDifferent++;
469 abcDj_[n++] = feasibleDj / infeasibleDj;
470 } else {
471 numberZero++;
472 }
473 }
474 }
475 abcProgress_.initialWeight_ = -1.0;
476 }
477 #ifdef CLP_USEFUL_PRINTOUT
478 printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
479 numberSame, numberDifferent, numberZero,
480 numberFreeSame, numberFreeDifferent, numberFreeZero);
481 #endif
482 // we want most to be same
483 if (n) {
484 std::sort(abcDj_, abcDj_ + n);
485 #ifdef CLP_USEFUL_PRINTOUT
486 double most = 0.95;
487 int which = static_cast< int >((1.0 - most) * static_cast< double >(n));
488 double take2 = -abcDj_[which] * infeasibilityCost_;
489 printf("XXXX inf cost %g take %g (range %g %g)\n", infeasibilityCost_, take2, -abcDj_[0] * infeasibilityCost_, -abcDj_[n - 1] * infeasibilityCost_);
490 #endif
491 double take = -abcDj_[0] * infeasibilityCost_;
492 // was infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10);
493 infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e3), 1.0000001e10);
494 #ifdef CLP_USEFUL_PRINTOUT
495 printf("XXXX changing weight to %g\n", infeasibilityCost_);
496 #endif
497 // may need to force redo of weights
498 needNewWeights = true;
499 }
500 abcNonLinearCost_->checkInfeasibilities(0.0);
501 gutsOfPrimalSolution(3);
502 }
503 // Check event
504 {
505 int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
506 if (status >= 0) {
507 problemStatus_ = 5;
508 secondaryStatus_ = ClpEventHandler::endOfFactorization;
509 break;
510 }
511 }
512 // Iterate
513 whileIterating(ifValuesPass ? 1 : 0);
514 }
515 }
516 // if infeasible get real values
517 //printf("XXXXY final cost %g\n",infeasibilityCost_);
518 abcProgress_.initialWeight_ = 0.0;
519 if (problemStatus_ == 1 && secondaryStatus_ != 6) {
520 infeasibilityCost_ = 0.0;
521 copyFromSaved();
522 delete abcNonLinearCost_;
523 abcNonLinearCost_ = new AbcNonLinearCost(this);
524 abcNonLinearCost_->checkInfeasibilities(0.0);
525 sumPrimalInfeasibilities_ = abcNonLinearCost_->sumInfeasibilities();
526 numberPrimalInfeasibilities_ = abcNonLinearCost_->numberInfeasibilities();
527 // and get good feasible duals
528 gutsOfPrimalSolution(1);
529 }
530 // clean up
531 unflag();
532 abcProgress_.clearTimesFlagged();
533 #if ABC_PARALLEL > 0
534 if (parallelMode()) {
535 for (int i = 0; i < 6; i++) {
536 delete rowArray_[i];
537 rowArray_[i] = NULL;
538 delete columnArray_[i];
539 columnArray_[i] = NULL;
540 }
541 }
542 #endif
543 //finish(startFinishOptions);
544 restoreData(data);
545 setObjectiveValue(abcNonLinearCost_->feasibleReportCost() + objectiveOffset_);
546 // clean up
547 if (problemStatus_ == 10)
548 abcNonLinearCost_->checkInfeasibilities(0.0);
549 delete abcNonLinearCost_;
550 abcNonLinearCost_ = NULL;
551 #if 0
552 if (numberRows_>80000&&numberRows_<90000) {
553 FILE * fp = fopen("save.stuff", "wb");
554 if (fp) {
555 fwrite(internalStatus_,sizeof(char),numberTotal_,fp);
556 fwrite(abcSolution_,sizeof(double),numberTotal_,fp);
557 } else {
558 printf("can't open save.stuff\n");
559 abort();
560 }
561 fclose(fp);
562 }
563 #endif
564 return problemStatus_;
565 }
566 /*
567 Reasons to come out:
568 -1 iterations etc
569 -2 inaccuracy
570 -3 slight inaccuracy (and done iterations)
571 -4 end of values pass and done iterations
572 +0 looks optimal (might be infeasible - but we will investigate)
573 +2 looks unbounded
574 +3 max iterations
575 */
whileIterating(int valuesOption)576 int AbcSimplexPrimal::whileIterating(int valuesOption)
577 {
578 #if 1
579 #define DELAYED_UPDATE
580 arrayForBtran_ = 0;
581 //setUsedArray(arrayForBtran_);
582 arrayForFtran_ = 1;
583 setUsedArray(arrayForFtran_);
584 arrayForFlipBounds_ = 2;
585 setUsedArray(arrayForFlipBounds_);
586 arrayForTableauRow_ = 3;
587 setUsedArray(arrayForTableauRow_);
588 //arrayForDualColumn_=4;
589 //setUsedArray(arrayForDualColumn_);
590 #if ABC_PARALLEL < 2
591 arrayForReplaceColumn_ = 4; //4;
592 #else
593 arrayForReplaceColumn_ = 6; //4;
594 setUsedArray(arrayForReplaceColumn_);
595 #endif
596 //arrayForFlipRhs_=5;
597 //setUsedArray(arrayForFlipRhs_);
598 #endif
599 // Say if values pass
600 int ifValuesPass = (firstFree_ >= 0) ? 1 : 0;
601 int returnCode = -1;
602 int superBasicType = 1;
603 if (valuesOption > 1)
604 superBasicType = 3;
605 int numberStartingInfeasibilities = abcNonLinearCost_->numberInfeasibilities();
606 // status stays at -1 while iterating, >=0 finished, -2 to invert
607 // status -3 to go to top without an invert
608 int numberFlaggedStart = abcProgress_.timesFlagged();
609 while (problemStatus_ == -1) {
610 if (!ifValuesPass) {
611 // choose column to come in
612 // can use pivotRow_ to update weights
613 // pass in list of cost changes so can do row updates (rowArray_[1])
614 // NOTE rowArray_[0] is used by computeDuals which is a
615 // slow way of getting duals but might be used
616 int saveSequence = sequenceIn_;
617 primalColumn(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForTableauRow_],
618 &usefulArray_[arrayForFlipBounds_]);
619 if (saveSequence >= 0 && saveSequence != sequenceOut_) {
620 if (getInternalStatus(saveSequence) == basic)
621 abcDj_[saveSequence] = 0.0;
622 }
623 #if ABC_NORMAL_DEBUG > 0
624 if (handler_->logLevel() == 63) {
625 for (int i = 0; i < numberTotal_; i++) {
626 if (fabs(abcDj_[i]) > 1.0e2)
627 assert(getInternalStatus(i) != basic);
628 }
629 double largestCost = 0.0;
630 double largestDj = 0.0;
631 double largestGoodDj = 0.0;
632 int iLargest = -1;
633 int numberInfeasibilities = 0;
634 double sum = 0.0;
635 for (int i = 0; i < numberTotal_; i++) {
636 if (getInternalStatus(i) == isFixed)
637 continue;
638 largestCost = CoinMax(largestCost, fabs(abcCost_[i]));
639 double dj = abcDj_[i];
640 if (getInternalStatus(i) == atLowerBound)
641 dj = -dj;
642 largestDj = CoinMax(largestDj, fabs(dj));
643 if (largestGoodDj < dj) {
644 largestGoodDj = dj;
645 iLargest = i;
646 }
647 if (getInternalStatus(i) == basic) {
648 if (abcSolution_[i] > abcUpper_[i] + primalTolerance_) {
649 numberInfeasibilities++;
650 sum += abcSolution_[i] - abcUpper_[i] - primalTolerance_;
651 } else if (abcSolution_[i] < abcLower_[i] - primalTolerance_) {
652 numberInfeasibilities++;
653 sum -= abcSolution_[i] - abcLower_[i] + primalTolerance_;
654 }
655 }
656 }
657 char xx;
658 int kLargest;
659 if (iLargest >= numberRows_) {
660 xx = 'C';
661 kLargest = iLargest - numberRows_;
662 } else {
663 xx = 'R';
664 kLargest = iLargest;
665 }
666 printf("largest cost %g, largest dj %g best %g (%d==%c%d) - %d infeasible (sum %g) nonlininf %d\n",
667 largestCost, largestDj, largestGoodDj, iLargest, xx, kLargest, numberInfeasibilities, sum,
668 abcNonLinearCost_->numberInfeasibilities());
669 assert(getInternalStatus(iLargest) != basic);
670 }
671 #endif
672 } else {
673 // in values pass
674 for (int i = 0; i < 4; i++)
675 multipleSequenceIn_[i] = -1;
676 int sequenceIn = nextSuperBasic(superBasicType, &usefulArray_[arrayForFlipBounds_]);
677 if (valuesOption > 1)
678 superBasicType = 2;
679 if (sequenceIn < 0) {
680 // end of values pass - initialize weights etc
681 sequenceIn_ = -1;
682 handler_->message(CLP_END_VALUES_PASS, messages_)
683 << numberIterations_;
684 stateOfProblem_ &= ~VALUES_PASS;
685 #ifdef TRY_SPLIT_VALUES_PASS
686 valuesStop = numberIterations_ + doOrdinaryVariables;
687 #endif
688 abcPrimalColumnPivot_->saveWeights(this, 5);
689 problemStatus_ = -2; // factorize now
690 pivotRow_ = -1; // say no weights update
691 returnCode = -4;
692 // Clean up
693 for (int i = 0; i < numberTotal_; i++) {
694 if (getInternalStatus(i) == atLowerBound || getInternalStatus(i) == isFixed)
695 abcSolution_[i] = abcLower_[i];
696 else if (getInternalStatus(i) == atUpperBound)
697 abcSolution_[i] = abcUpper_[i];
698 }
699 break;
700 } else {
701 // normal
702 sequenceIn_ = sequenceIn;
703 valueIn_ = abcSolution_[sequenceIn_];
704 lowerIn_ = abcLower_[sequenceIn_];
705 upperIn_ = abcUpper_[sequenceIn_];
706 dualIn_ = abcDj_[sequenceIn_];
707 // see if any more
708 if (maximumIterations() == 100000) {
709 multipleSequenceIn_[0] = sequenceIn;
710 for (int i = 1; i < 4; i++) {
711 int sequenceIn = nextSuperBasic(superBasicType, &usefulArray_[arrayForFlipBounds_]);
712 if (sequenceIn >= 0)
713 multipleSequenceIn_[i] = sequenceIn;
714 else
715 break;
716 }
717 }
718 }
719 }
720 pivotRow_ = -1;
721 sequenceOut_ = -1;
722 usefulArray_[arrayForFtran_].clear();
723 if (sequenceIn_ >= 0) {
724 // we found a pivot column
725 assert(!flagged(sequenceIn_));
726 //#define MULTIPLE_PRICE
727 // do second half of iteration
728 if (multipleSequenceIn_[1] == -1 || maximumIterations() != 100000) {
729 returnCode = pivotResult(ifValuesPass);
730 } else {
731 if (multipleSequenceIn_[0] < 0)
732 multipleSequenceIn_[0] = sequenceIn_;
733 returnCode = pivotResult4(ifValuesPass);
734 #ifdef MULTIPLE_PRICE
735 if (sequenceIn_ >= 0)
736 returnCode = pivotResult(ifValuesPass);
737 #endif
738 }
739 if (numberStartingInfeasibilities && !abcNonLinearCost_->numberInfeasibilities()) {
740 //if (abcFactorization_->pivots()>200&&numberIterations_>2*(numberRows_+numberColumns_))
741 if (abcFactorization_->pivots() > 2 && numberIterations_ > (numberRows_ + numberColumns_) && (stateOfProblem_ & PESSIMISTIC) != 0)
742 returnCode = -2; // refactorize - maybe just after n iterations
743 }
744 if (returnCode < -1 && returnCode > -5) {
745 problemStatus_ = -2; //
746 } else if (returnCode == -5) {
747 if (abcProgress_.timesFlagged() > 10 + numberFlaggedStart)
748 problemStatus_ = -2;
749 if ((moreSpecialOptions_ & 16) == 0 && abcFactorization_->pivots()) {
750 moreSpecialOptions_ |= 16;
751 problemStatus_ = -2;
752 }
753 // otherwise something flagged - continue;
754 } else if (returnCode == 2) {
755 problemStatus_ = -5; // looks unbounded
756 } else if (returnCode == 4) {
757 problemStatus_ = -2; // looks unbounded but has iterated
758 } else if (returnCode != -1) {
759 assert(returnCode == 3);
760 if (problemStatus_ != 5)
761 problemStatus_ = 3;
762 }
763 } else {
764 // no pivot column
765 #if ABC_NORMAL_DEBUG > 3
766 if (handler_->logLevel() & 32)
767 printf("** no column pivot\n");
768 #endif
769 if (abcNonLinearCost_->numberInfeasibilities())
770 problemStatus_ = -4; // might be infeasible
771 // Force to re-factorize early next time
772 int numberPivots = abcFactorization_->pivots();
773 returnCode = 0;
774 #ifdef CLP_USER_DRIVEN
775 // If large number of pivots trap later?
776 if (problemStatus_ == -1 && numberPivots < 13000) {
777 int status = eventHandler_->event(ClpEventHandler::noCandidateInPrimal);
778 if (status >= 0 && status < 10) {
779 // carry on
780 problemStatus_ = -1;
781 if (status == 0)
782 break;
783 } else if (status >= 10) {
784 problemStatus_ = status - 10;
785 break;
786 } else {
787 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
788 break;
789 }
790 }
791 #else
792 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
793 break;
794 #endif
795 }
796 }
797 if (valuesOption > 1)
798 usefulArray_[arrayForFlipBounds_].setNumElements(0);
799 return returnCode;
800 }
801 /* Checks if finished. Updates status */
statusOfProblemInPrimal(int type)802 void AbcSimplexPrimal::statusOfProblemInPrimal(int type)
803 {
804 int saveFirstFree = firstFree_;
805 // number of pivots done
806 int numberPivots = abcFactorization_->pivots();
807 #if 0
808 printf("statusOf %d pivots type %d pstatus %d forceFac %d dont %d\n",
809 numberPivots,type,problemStatus_,forceFactorization_,dontFactorizePivots_);
810 #endif
811 //int saveType=type;
812 if (type > 3) {
813 // force factorization
814 numberPivots = 9999999;
815 if (type > 10)
816 type -= 10;
817 else
818 type &= 3;
819 }
820 #ifndef TRY_ABC_GUS
821 bool doFactorization = (type != 3 && (numberPivots > dontFactorizePivots_ || numberIterations_ == baseIteration_ || problemStatus_ == 10));
822 #else
823 bool doFactorization = (type != 3 && (numberPivots > dontFactorizePivots_ || numberIterations_ == baseIteration_));
824 #endif
825 bool doWeights = doFactorization || problemStatus_ == 10;
826 if (type == 2) {
827 // trouble - restore solution
828 restoreGoodStatus(1);
829 forceFactorization_ = 1; // a bit drastic but ..
830 pivotRow_ = -1; // say no weights update
831 changeMade_++; // say change made
832 }
833 int tentativeStatus = problemStatus_;
834 double lastSumInfeasibility = COIN_DBL_MAX;
835 int lastNumberInfeasibility = 1;
836 #ifndef CLP_CAUTION
837 #define CLP_CAUTION 1
838 #endif
839 #if CLP_CAUTION
840 double lastAverageInfeasibility = sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_ + 1);
841 #endif
842 if (numberIterations_ && type) {
843 lastSumInfeasibility = abcNonLinearCost_->sumInfeasibilities();
844 lastNumberInfeasibility = abcNonLinearCost_->numberInfeasibilities();
845 } else {
846 lastAverageInfeasibility = 1.0e10;
847 }
848 bool ifValuesPass = (stateOfProblem_ & VALUES_PASS) != 0;
849 bool takenAction = false;
850 double sumInfeasibility = 0.0;
851 if (problemStatus_ > -3 || problemStatus_ == -4) {
852 // factorize
853 // later on we will need to recover from singularities
854 // also we could skip if first time
855 // do weights
856 // This may save pivotRow_ for use
857 if (doFactorization)
858 abcPrimalColumnPivot_->saveWeights(this, 1);
859
860 if (!type) {
861 // be optimistic
862 stateOfProblem_ &= ~PESSIMISTIC;
863 // but use 0.1
864 double newTolerance = CoinMax(0.1, saveData_.pivotTolerance_);
865 abcFactorization_->pivotTolerance(newTolerance);
866 }
867 if (doFactorization) {
868 // is factorization okay?
869 int solveType = ifValuesPass ? 11 : 1;
870 if (!type)
871 solveType++;
872 int factorStatus = internalFactorize(solveType);
873 if (factorStatus) {
874 if (type != 1 || largestPrimalError_ > 1.0e3
875 || largestDualError_ > 1.0e3) {
876 #if ABC_NORMAL_DEBUG > 0
877 printf("Bad initial basis\n");
878 #endif
879 internalFactorize(2);
880 } else {
881 // no - restore previous basis
882 // Keep any flagged variables
883 for (int i = 0; i < numberTotal_; i++) {
884 if (flagged(i))
885 internalStatusSaved_[i] |= 64; //say flagged
886 }
887 restoreGoodStatus(1);
888 if (numberPivots <= 1) {
889 // throw out something
890 if (sequenceIn_ >= 0 && getInternalStatus(sequenceIn_) != basic) {
891 setFlagged(sequenceIn_);
892 } else if (sequenceOut_ >= 0 && getInternalStatus(sequenceOut_) != basic) {
893 setFlagged(sequenceOut_);
894 }
895 abcProgress_.incrementTimesFlagged();
896 double newTolerance = CoinMax(0.5 + 0.499 * randomNumberGenerator_.randomDouble(), abcFactorization_->pivotTolerance());
897 abcFactorization_->pivotTolerance(newTolerance);
898 } else {
899 // Go to safe
900 abcFactorization_->pivotTolerance(0.99);
901 }
902 forceFactorization_ = 1; // a bit drastic but ..
903 type = 2;
904 abcNonLinearCost_->checkInfeasibilities();
905 if (internalFactorize(2) != 0) {
906 largestPrimalError_ = 1.0e4; // force other type
907 }
908 }
909 changeMade_++; // say change made
910 }
911 }
912 if (problemStatus_ != -4)
913 problemStatus_ = -3;
914 if (abcProgress_.realInfeasibility_[0] < 1.0e-1 && primalTolerance_ == 1.0e-7 && abcProgress_.iterationNumber_[0] > 0 && abcProgress_.iterationNumber_[CLP_PROGRESS - 1] - abcProgress_.iterationNumber_[0] > 25) {
915 // default - so user did not set
916 int iP;
917 double minAverage = COIN_DBL_MAX;
918 double maxAverage = 0.0;
919 for (iP = 0; iP < CLP_PROGRESS; iP++) {
920 int n = abcProgress_.numberInfeasibilities_[iP];
921 if (!n) {
922 break;
923 } else {
924 double average = abcProgress_.realInfeasibility_[iP];
925 if (average > 0.1)
926 break;
927 average /= static_cast< double >(n);
928 minAverage = CoinMin(minAverage, average);
929 maxAverage = CoinMax(maxAverage, average);
930 }
931 }
932 if (iP == CLP_PROGRESS && minAverage < 1.0e-5 && maxAverage < 1.0e-3) {
933 // change tolerance
934 #if CBC_USEFUL_PRINTING > 0
935 printf("CCchanging tolerance\n");
936 #endif
937 primalTolerance_ = 1.0e-6;
938 dblParam_[ClpPrimalTolerance] = 1.0e-6;
939 moreSpecialOptions_ |= 4194304;
940 }
941 }
942 // at this stage status is -3 or -5 if looks unbounded
943 // get primal and dual solutions
944 // put back original costs and then check
945 // createRim(4); // costs do not change
946 if (ifValuesPass && numberIterations_ == baseIteration_) {
947 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
948 lastSumInfeasibility = abcNonLinearCost_->largestInfeasibility();
949 }
950 #ifdef CLP_USER_DRIVEN
951 int status = eventHandler_->event(ClpEventHandler::goodFactorization);
952 if (status >= 0) {
953 lastSumInfeasibility = COIN_DBL_MAX;
954 }
955 #endif
956 if (ifValuesPass && numberIterations_ == baseIteration_) {
957 double *save = new double[numberRows_];
958 int numberOut = 1;
959 while (numberOut) {
960 for (int iRow = 0; iRow < numberRows_; iRow++) {
961 int iPivot = abcPivotVariable_[iRow];
962 save[iRow] = abcSolution_[iPivot];
963 }
964 gutsOfPrimalSolution(2);
965 double badInfeasibility = abcNonLinearCost_->largestInfeasibility();
966 numberOut = 0;
967 // But may be very large rhs etc
968 double useError = CoinMin(largestPrimalError_,
969 1.0e5 / CoinAbcMaximumAbsElement(abcSolution_, numberTotal_));
970 if ((lastSumInfeasibility < incomingInfeasibility_ || badInfeasibility > (CoinMax(10.0, 100.0 * lastSumInfeasibility)))
971 && (badInfeasibility > 10.0 || useError > 1.0e-3)) {
972 //printf("Original largest infeas %g, now %g, primalError %g\n",
973 // lastSumInfeasibility,abcNonLinearCost_->largestInfeasibility(),
974 // largestPrimalError_);
975 // throw out up to 1000 structurals
976 int *sort = new int[numberRows_];
977 // first put back solution and store difference
978 for (int iRow = 0; iRow < numberRows_; iRow++) {
979 int iPivot = abcPivotVariable_[iRow];
980 double difference = fabs(abcSolution_[iPivot] - save[iRow]);
981 abcSolution_[iPivot] = save[iRow];
982 save[iRow] = difference;
983 }
984 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
985 if (handler_->logLevel() > 1)
986 printf("Largest infeasibility %g\n", abcNonLinearCost_->largestInfeasibility());
987 int numberBasic = 0;
988 for (int iRow = 0; iRow < numberRows_; iRow++) {
989 int iPivot = abcPivotVariable_[iRow];
990 if (iPivot >= numberRows_) {
991 // column
992 double difference = save[iRow];
993 if (difference > 1.0e-4) {
994 sort[numberOut] = iRow;
995 save[numberOut++] = -difference;
996 if (getInternalStatus(iPivot) == basic)
997 numberBasic++;
998 }
999 }
1000 }
1001 if (!numberBasic) {
1002 //printf("no errors on basic - going to all slack - numberOut %d\n",numberOut);
1003 // allow
1004 numberOut = 0;
1005 }
1006 CoinSort_2(save, save + numberOut, sort);
1007 numberOut = CoinMin(1000, numberOut);
1008 // for now bring in any slack
1009 int jRow = 0;
1010 for (int i = 0; i < numberOut; i++) {
1011 int iRow = sort[i];
1012 int iColumn = abcPivotVariable_[iRow];
1013 assert(getInternalStatus(iColumn) == basic);
1014 setInternalStatus(iColumn, superBasic);
1015 while (getInternalStatus(jRow) == basic)
1016 jRow++;
1017 setInternalStatus(jRow, basic);
1018 abcPivotVariable_[iRow] = jRow;
1019 if (fabs(abcSolution_[iColumn]) > 1.0e10) {
1020 if (abcUpper_[iColumn] < 0.0) {
1021 abcSolution_[iColumn] = abcUpper_[iColumn];
1022 } else if (abcLower_[iColumn] > 0.0) {
1023 abcSolution_[iColumn] = abcLower_[iColumn];
1024 } else {
1025 abcSolution_[iColumn] = 0.0;
1026 }
1027 }
1028 }
1029 delete[] sort;
1030 }
1031 if (numberOut) {
1032 double savePivotTolerance = abcFactorization_->pivotTolerance();
1033 abcFactorization_->pivotTolerance(0.99);
1034 int factorStatus = internalFactorize(12);
1035 abcFactorization_->pivotTolerance(savePivotTolerance);
1036 assert(!factorStatus);
1037 }
1038 }
1039 delete[] save;
1040 }
1041 gutsOfPrimalSolution(3);
1042 sumInfeasibility = abcNonLinearCost_->sumInfeasibilities();
1043 int reason2 = 0;
1044 #if CLP_CAUTION
1045 #if CLP_CAUTION == 2
1046 double test2 = 1.0e5;
1047 #else
1048 double test2 = 1.0e-1;
1049 #endif
1050 if (!lastSumInfeasibility && sumInfeasibility > 1.0e3 * primalTolerance_ && lastAverageInfeasibility < test2 && numberPivots > 10 && !ifValuesPass)
1051 reason2 = 3;
1052 if (lastSumInfeasibility < 1.0e-6 && sumInfeasibility > 1.0e-3 && numberPivots > 10 && !ifValuesPass)
1053 reason2 = 4;
1054 #endif
1055 if ((sumInfeasibility > 1.0e7 && sumInfeasibility > 100.0 * lastSumInfeasibility
1056 && abcFactorization_->pivotTolerance() < 0.11)
1057 || (largestPrimalError_ > 1.0e10 && largestDualError_ > 1.0e10))
1058 reason2 = 2;
1059 if (reason2) {
1060 takenAction = true;
1061 problemStatus_ = tentativeStatus;
1062 doFactorization = true;
1063 if (numberPivots) {
1064 // go back
1065 // trouble - restore solution
1066 restoreGoodStatus(1);
1067 sequenceIn_ = -1;
1068 sequenceOut_ = -1;
1069 abcNonLinearCost_->checkInfeasibilities();
1070 if (reason2 < 3) {
1071 // Go to safe
1072 abcFactorization_->pivotTolerance(CoinMin(0.99, 1.01 * abcFactorization_->pivotTolerance()));
1073 forceFactorization_ = 1; // a bit drastic but ..
1074 } else if (forceFactorization_ < 0) {
1075 forceFactorization_ = CoinMin(numberPivots / 4, 100);
1076 } else {
1077 forceFactorization_ = CoinMin(forceFactorization_,
1078 CoinMax(3, numberPivots / 4));
1079 }
1080 pivotRow_ = -1; // say no weights update
1081 changeMade_++; // say change made
1082 if (numberPivots == 1) {
1083 // throw out something
1084 if (sequenceIn_ >= 0 && getInternalStatus(sequenceIn_) != basic) {
1085 setFlagged(sequenceIn_);
1086 } else if (sequenceOut_ >= 0 && getInternalStatus(sequenceOut_) != basic) {
1087 setFlagged(sequenceOut_);
1088 }
1089 abcProgress_.incrementTimesFlagged();
1090 }
1091 type = 2; // so will restore weights
1092 if (internalFactorize(2) != 0) {
1093 largestPrimalError_ = 1.0e4; // force other type
1094 }
1095 numberPivots = 0;
1096 gutsOfPrimalSolution(3);
1097 sumInfeasibility = abcNonLinearCost_->sumInfeasibilities();
1098 }
1099 }
1100 }
1101 // Double check reduced costs if no action
1102 if (abcProgress_.lastIterationNumber(0) == numberIterations_) {
1103 if (abcPrimalColumnPivot_->looksOptimal()) {
1104 numberDualInfeasibilities_ = 0;
1105 sumDualInfeasibilities_ = 0.0;
1106 }
1107 }
1108 // If in primal and small dj give up
1109 if ((specialOptions_ & 1024) != 0 && !numberPrimalInfeasibilities_ && numberDualInfeasibilities_) {
1110 double average = sumDualInfeasibilities_ / (static_cast< double >(numberDualInfeasibilities_));
1111 if (numberIterations_ > 300 && average < 1.0e-4) {
1112 numberDualInfeasibilities_ = 0;
1113 sumDualInfeasibilities_ = 0.0;
1114 }
1115 }
1116 // Check if looping
1117 int loop;
1118 ifValuesPass = firstFree_ >= 0;
1119 if (type != 2 && !ifValuesPass)
1120 loop = abcProgress_.looping();
1121 else
1122 loop = -1;
1123 if ((moreSpecialOptions_ & 2048) != 0 && !numberPrimalInfeasibilities_ && numberDualInfeasibilities_) {
1124 double average = sumDualInfeasibilities_ / (static_cast< double >(numberDualInfeasibilities_));
1125 if (abcProgress_.lastIterationNumber(2) == numberIterations_ && ((abcProgress_.timesFlagged() > 2 && average < 1.0e-1) || abcProgress_.timesFlagged() > 20)) {
1126 numberDualInfeasibilities_ = 0;
1127 sumDualInfeasibilities_ = 0.0;
1128 problemStatus_ = 3;
1129 loop = 0;
1130 }
1131 }
1132 if (loop >= 0) {
1133 if (!problemStatus_) {
1134 // declaring victory
1135 numberPrimalInfeasibilities_ = 0;
1136 sumPrimalInfeasibilities_ = 0.0;
1137 // say been feasible
1138 abcState_ |= CLP_ABC_BEEN_FEASIBLE;
1139 } else {
1140 problemStatus_ = loop; //exit if in loop
1141 problemStatus_ = 10; // instead - try other algorithm
1142 numberPrimalInfeasibilities_ = abcNonLinearCost_->numberInfeasibilities();
1143 }
1144 problemStatus_ = 10; // instead - try other algorithm
1145 return;
1146 } else if (loop < -1) {
1147 // Is it time for drastic measures
1148 if (abcNonLinearCost_->numberInfeasibilities() && abcProgress_.badTimes() > 5 && abcProgress_.oddState() < 10 && abcProgress_.oddState() >= 0) {
1149 abcProgress_.newOddState();
1150 abcNonLinearCost_->zapCosts();
1151 }
1152 // something may have changed
1153 gutsOfPrimalSolution(3);
1154 }
1155 // If progress then reset costs
1156 if (loop == -1 && !abcNonLinearCost_->numberInfeasibilities() && abcProgress_.oddState() < 0) {
1157 copyFromSaved(); // costs back
1158 delete abcNonLinearCost_;
1159 abcNonLinearCost_ = new AbcNonLinearCost(this);
1160 abcProgress_.endOddState();
1161 gutsOfPrimalSolution(3);
1162 }
1163 abcProgress_.modifyObjective(abcNonLinearCost_->feasibleReportCost());
1164 if (!lastNumberInfeasibility && sumInfeasibility && numberPivots > 1 && !ifValuesPass && !takenAction && abcProgress_.lastObjective(0) >= abcProgress_.lastObjective(CLP_PROGRESS - 1)) {
1165 // first increase minimumThetaMovement_;
1166 // be more careful
1167 //abcFactorization_->saferTolerances (-0.99, -1.002);
1168 #if ABC_NORMAL_DEBUG > 0
1169 if (handler_->logLevel() == 63)
1170 printf("thought feasible but sum is %g force %d minimum theta %g\n", sumInfeasibility,
1171 forceFactorization_, minimumThetaMovement_);
1172 #endif
1173 stateOfProblem_ |= PESSIMISTIC;
1174 if (minimumThetaMovement_ < 1.0e-10) {
1175 minimumThetaMovement_ *= 1.1;
1176 } else if (0) {
1177 int maxFactor = abcFactorization_->maximumPivots();
1178 if (maxFactor > 10) {
1179 if (forceFactorization_ < 0)
1180 forceFactorization_ = maxFactor;
1181 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 2));
1182 if (handler_->logLevel() == 63)
1183 printf("Reducing factorization frequency to %d\n", forceFactorization_);
1184 }
1185 }
1186 }
1187 // Flag to say whether to go to dual to clean up
1188 bool goToDual = false;
1189 // really for free variables in
1190 //if((progressFlag_&2)!=0)
1191 //problemStatus_=-1;;
1192 progressFlag_ = 0; //reset progress flag
1193
1194 handler_->message(CLP_SIMPLEX_STATUS, messages_)
1195 << numberIterations_ << abcNonLinearCost_->feasibleReportCost() + objectiveOffset_;
1196 handler_->printing(abcNonLinearCost_->numberInfeasibilities() > 0)
1197 << abcNonLinearCost_->sumInfeasibilities() << abcNonLinearCost_->numberInfeasibilities();
1198 handler_->printing(sumDualInfeasibilities_ > 0.0)
1199 << sumDualInfeasibilities_ << numberDualInfeasibilities_;
1200 handler_->printing(numberDualInfeasibilitiesWithoutFree_
1201 < numberDualInfeasibilities_)
1202 << numberDualInfeasibilitiesWithoutFree_;
1203 handler_->message() << CoinMessageEol;
1204 if (!primalFeasible()) {
1205 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
1206 gutsOfPrimalSolution(2);
1207 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
1208 }
1209 double trueInfeasibility = abcNonLinearCost_->sumInfeasibilities();
1210 if (!abcNonLinearCost_->numberInfeasibilities() && infeasibilityCost_ == 1.0e10 && !ifValuesPass && true) {
1211 // say been feasible
1212 abcState_ |= CLP_ABC_BEEN_FEASIBLE;
1213 // relax if default
1214 infeasibilityCost_ = CoinMin(CoinMax(100.0 * sumDualInfeasibilities_, 1.0e8), 1.00000001e10);
1215 // reset looping criterion
1216 abcProgress_.reset();
1217 trueInfeasibility = 1.123456e10;
1218 }
1219 if (trueInfeasibility > 1.0) {
1220 // If infeasibility going up may change weights
1221 double testValue = trueInfeasibility - 1.0e-4 * (10.0 + trueInfeasibility);
1222 double lastInf = abcProgress_.lastInfeasibility(1);
1223 double lastInf3 = abcProgress_.lastInfeasibility(3);
1224 double thisObj = abcProgress_.lastObjective(0);
1225 double thisInf = abcProgress_.lastInfeasibility(0);
1226 thisObj += infeasibilityCost_ * 2.0 * thisInf;
1227 double lastObj = abcProgress_.lastObjective(1);
1228 lastObj += infeasibilityCost_ * 2.0 * lastInf;
1229 double lastObj3 = abcProgress_.lastObjective(3);
1230 lastObj3 += infeasibilityCost_ * 2.0 * lastInf3;
1231 if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-7
1232 && firstFree_ < 0 && thisInf >= lastInf) {
1233 #if ABC_NORMAL_DEBUG > 0
1234 if (handler_->logLevel() == 63)
1235 printf("lastobj %g this %g force %d\n", lastObj, thisObj, forceFactorization_);
1236 #endif
1237 int maxFactor = abcFactorization_->maximumPivots();
1238 if (maxFactor > 10) {
1239 if (forceFactorization_ < 0)
1240 forceFactorization_ = maxFactor;
1241 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 2));
1242 #if ABC_NORMAL_DEBUG > 0
1243 if (handler_->logLevel() == 63)
1244 printf("Reducing factorization frequency to %d\n", forceFactorization_);
1245 #endif
1246 }
1247 } else if (lastObj3 < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj3)) - 1.0e-7
1248 && firstFree_ < 0 && thisInf >= lastInf) {
1249 #if ABC_NORMAL_DEBUG > 0
1250 if (handler_->logLevel() == 63)
1251 printf("lastobj3 %g this3 %g force %d\n", lastObj3, thisObj, forceFactorization_);
1252 #endif
1253 int maxFactor = abcFactorization_->maximumPivots();
1254 if (maxFactor > 10) {
1255 if (forceFactorization_ < 0)
1256 forceFactorization_ = maxFactor;
1257 forceFactorization_ = CoinMax(1, (forceFactorization_ * 2) / 3);
1258 #if ABC_NORMAL_DEBUG > 0
1259 if (handler_->logLevel() == 63)
1260 printf("Reducing factorization frequency to %d\n", forceFactorization_);
1261 #endif
1262 }
1263 } else if (lastInf < testValue || trueInfeasibility == 1.123456e10) {
1264 if (infeasibilityCost_ < 1.0e14) {
1265 infeasibilityCost_ *= 1.5;
1266 // reset looping criterion
1267 abcProgress_.reset();
1268 #if ABC_NORMAL_DEBUG > 0
1269 if (handler_->logLevel() == 63)
1270 printf("increasing weight to %g\n", infeasibilityCost_);
1271 #endif
1272 gutsOfPrimalSolution(3);
1273 }
1274 }
1275 }
1276 // we may wish to say it is optimal even if infeasible
1277 bool alwaysOptimal = (specialOptions_ & 1) != 0;
1278 #if CLP_CAUTION
1279 // If twice nearly there ...
1280 if (lastAverageInfeasibility < 2.0 * dualTolerance_) {
1281 double averageInfeasibility = sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_ + 1);
1282 printf("last av %g now %g\n", lastAverageInfeasibility,
1283 averageInfeasibility);
1284 if (averageInfeasibility < 2.0 * dualTolerance_)
1285 sumOfRelaxedDualInfeasibilities_ = 0.0;
1286 }
1287 #endif
1288 // give code benefit of doubt
1289 if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1290 // say been feasible
1291 abcState_ |= CLP_ABC_BEEN_FEASIBLE;
1292 // say optimal (with these bounds etc)
1293 numberDualInfeasibilities_ = 0;
1294 sumDualInfeasibilities_ = 0.0;
1295 numberPrimalInfeasibilities_ = 0;
1296 sumPrimalInfeasibilities_ = 0.0;
1297 // But if real primal infeasibilities nonzero carry on
1298 if (abcNonLinearCost_->numberInfeasibilities()) {
1299 // most likely to happen if infeasible
1300 double relaxedToleranceP = primalTolerance_;
1301 // we can't really trust infeasibilities if there is primal error
1302 double error = CoinMin(1.0e-2, largestPrimalError_);
1303 // allow tolerance at least slightly bigger than standard
1304 relaxedToleranceP = relaxedToleranceP + error;
1305 int ninfeas = abcNonLinearCost_->numberInfeasibilities();
1306 double sum = abcNonLinearCost_->sumInfeasibilities();
1307 double average = sum / static_cast< double >(ninfeas);
1308 #if ABC_NORMAL_DEBUG > 3
1309 if (handler_->logLevel() > 0)
1310 printf("nonLinearCost says infeasible %d summing to %g\n",
1311 ninfeas, sum);
1312 #endif
1313 if (average > relaxedToleranceP) {
1314 sumOfRelaxedPrimalInfeasibilities_ = sum;
1315 numberPrimalInfeasibilities_ = ninfeas;
1316 sumPrimalInfeasibilities_ = sum;
1317 #if ABC_NORMAL_DEBUG > 3
1318 bool unflagged =
1319 #endif
1320 unflag();
1321 abcProgress_.clearTimesFlagged();
1322 #if ABC_NORMAL_DEBUG > 3
1323 if (unflagged && handler_->logLevel() > 0)
1324 printf(" - but flagged variables\n");
1325 #endif
1326 }
1327 }
1328 }
1329 //double saveSum=sumOfRelaxedDualInfeasibilities_;
1330 // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
1331 if ((dualFeasible() || problemStatus_ == -4) && !ifValuesPass) {
1332 // see if extra helps
1333 if (abcNonLinearCost_->numberInfeasibilities() && (abcNonLinearCost_->sumInfeasibilities() > 1.0e-3 || sumOfRelaxedPrimalInfeasibilities_)
1334 && !alwaysOptimal) {
1335 //may need infeasiblity cost changed
1336 // we can see if we can construct a ray
1337 // do twice to make sure Primal solution has settled
1338 // put non-basics to bounds in case tolerance moved
1339 // put back original costs
1340 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
1341 ;
1342 abcNonLinearCost_->checkInfeasibilities(0.0);
1343 gutsOfPrimalSolution(3);
1344 // put back original costs
1345 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
1346 ;
1347 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
1348 gutsOfPrimalSolution(3);
1349 // may have fixed infeasibilities - double check
1350 if (abcNonLinearCost_->numberInfeasibilities() == 0) {
1351 // carry on
1352 problemStatus_ = -1;
1353 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
1354 } else {
1355 if ((infeasibilityCost_ >= 1.0e18 || numberDualInfeasibilities_ == 0) && perturbation_ == 101
1356 && (specialOptions_ & 8192) == 0) {
1357 goToDual = unPerturb(); // stop any further perturbation
1358 #ifndef TRY_ABC_GUS
1359 if (abcNonLinearCost_->sumInfeasibilities() > 1.0e-1)
1360 goToDual = false;
1361 #endif
1362 numberDualInfeasibilities_ = 1; // carry on
1363 problemStatus_ = -1;
1364 } else if (numberDualInfeasibilities_ == 0 && largestDualError_ > 1.0e-2
1365 #ifndef TRY_ABC_GUS
1366 && ((moreSpecialOptions_ & 256) == 0 && (specialOptions_ & 8192) == 0)
1367 #else
1368 && (specialOptions_ & 8192) == 0
1369 #endif
1370 ) {
1371 goToDual = true;
1372 abcFactorization_->pivotTolerance(CoinMax(0.5, abcFactorization_->pivotTolerance()));
1373 }
1374 if (!goToDual) {
1375 if (infeasibilityCost_ >= 1.0e20 || numberDualInfeasibilities_ == 0) {
1376 // we are infeasible - use as ray
1377 delete[] ray_;
1378 ray_ = new double[numberRows_];
1379 CoinMemcpyN(dual_, numberRows_, ray_);
1380 // and get feasible duals
1381 infeasibilityCost_ = 0.0;
1382 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
1383 ;
1384 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
1385 gutsOfPrimalSolution(3);
1386 // so will exit
1387 infeasibilityCost_ = 1.0e30;
1388 // reset infeasibilities
1389 sumPrimalInfeasibilities_ = abcNonLinearCost_->sumInfeasibilities();
1390 ;
1391 numberPrimalInfeasibilities_ = abcNonLinearCost_->numberInfeasibilities();
1392 }
1393 if (infeasibilityCost_ < 1.0e20) {
1394 infeasibilityCost_ *= 5.0;
1395 // reset looping criterion
1396 abcProgress_.reset();
1397 changeMade_++; // say change made
1398 handler_->message(CLP_PRIMAL_WEIGHT, messages_)
1399 << infeasibilityCost_
1400 << CoinMessageEol;
1401 // put back original costs and then check
1402 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
1403 ;
1404 abcNonLinearCost_->checkInfeasibilities(0.0);
1405 gutsOfPrimalSolution(3);
1406 problemStatus_ = -1; //continue
1407 #ifndef TRY_ABC_GUS
1408 goToDual = false;
1409 #else
1410 if ((specialOptions_ & 8192) == 0 && !sumOfRelaxedDualInfeasibilities_)
1411 goToDual = true;
1412 #endif
1413 } else {
1414 // say infeasible
1415 problemStatus_ = 1;
1416 }
1417 }
1418 }
1419 } else {
1420 // may be optimal
1421 if (perturbation_ == 101) {
1422 goToDual = unPerturb(); // stop any further perturbation
1423 #ifndef TRY_ABC_GUS
1424 if ((numberRows_ > 20000 || numberDualInfeasibilities_) && !numberTimesOptimal_)
1425 #else
1426 if ((specialOptions_ & 8192) != 0)
1427 #endif
1428 goToDual = false; // Better to carry on a bit longer
1429 lastCleaned_ = -1; // carry on
1430 }
1431 bool unflagged = (unflag() != 0);
1432 abcProgress_.clearTimesFlagged();
1433 if (lastCleaned_ != numberIterations_ || unflagged) {
1434 handler_->message(CLP_PRIMAL_OPTIMAL, messages_)
1435 << primalTolerance_
1436 << CoinMessageEol;
1437 if (numberTimesOptimal_ < 4) {
1438 numberTimesOptimal_++;
1439 changeMade_++; // say change made
1440 if (numberTimesOptimal_ == 1) {
1441 // better to have small tolerance even if slower
1442 abcFactorization_->zeroTolerance(CoinMin(abcFactorization_->zeroTolerance(), 1.0e-15));
1443 }
1444 lastCleaned_ = numberIterations_;
1445 if (primalTolerance_ != dblParam_[ClpPrimalTolerance])
1446 handler_->message(CLP_PRIMAL_ORIGINAL, messages_)
1447 << CoinMessageEol;
1448 double oldTolerance = primalTolerance_;
1449 primalTolerance_ = dblParam_[ClpPrimalTolerance];
1450 // put back original costs and then check
1451 #if 0 //ndef NDEBUG
1452 double largestDifference=0.0;
1453 for (int i=0;i<numberTotal_;i++)
1454 largestDifference=CoinMax(largestDifference,fabs(abcCost_[i]-costSaved_[i]));
1455 if (handler_->logLevel()==63)
1456 printf("largest change in cost %g\n",largestDifference);
1457 #endif
1458 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
1459 ;
1460 abcNonLinearCost_->checkInfeasibilities(oldTolerance);
1461 gutsOfPrimalSolution(3);
1462 if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1463 // say optimal (with these bounds etc)
1464 numberDualInfeasibilities_ = 0;
1465 sumDualInfeasibilities_ = 0.0;
1466 numberPrimalInfeasibilities_ = 0;
1467 sumPrimalInfeasibilities_ = 0.0;
1468 goToDual = false;
1469 }
1470 if (dualFeasible() && !abcNonLinearCost_->numberInfeasibilities() && lastCleaned_ >= 0)
1471 problemStatus_ = 0;
1472 else
1473 problemStatus_ = -1;
1474 } else {
1475 problemStatus_ = 0; // optimal
1476 if (lastCleaned_ < numberIterations_) {
1477 handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
1478 << CoinMessageEol;
1479 }
1480 }
1481 } else {
1482 if (!alwaysOptimal || !sumOfRelaxedPrimalInfeasibilities_)
1483 problemStatus_ = 0; // optimal
1484 else
1485 problemStatus_ = 1; // infeasible
1486 }
1487 }
1488 } else {
1489 // see if looks unbounded
1490 if (problemStatus_ == -5) {
1491 if (abcNonLinearCost_->numberInfeasibilities()) {
1492 if (infeasibilityCost_ > 1.0e18 && perturbation_ == 101) {
1493 // back off weight
1494 infeasibilityCost_ = 1.0e13;
1495 // reset looping criterion
1496 abcProgress_.reset();
1497 unPerturb(); // stop any further perturbation
1498 }
1499 //we need infeasiblity cost changed
1500 if (infeasibilityCost_ < 1.0e20) {
1501 infeasibilityCost_ *= 5.0;
1502 // reset looping criterion
1503 abcProgress_.reset();
1504 changeMade_++; // say change made
1505 handler_->message(CLP_PRIMAL_WEIGHT, messages_)
1506 << infeasibilityCost_
1507 << CoinMessageEol;
1508 // put back original costs and then check
1509 CoinAbcMemcpy(abcCost_, costSaved_, numberTotal_);
1510 ;
1511 gutsOfPrimalSolution(3);
1512 problemStatus_ = -1; //continue
1513 } else {
1514 // say infeasible
1515 problemStatus_ = 1;
1516 // we are infeasible - use as ray
1517 delete[] ray_;
1518 ray_ = new double[numberRows_];
1519 CoinMemcpyN(dual_, numberRows_, ray_);
1520 }
1521 } else {
1522 // say unbounded
1523 problemStatus_ = 2;
1524 }
1525 } else {
1526 // carry on
1527 problemStatus_ = -1;
1528 if (type == 3 && !ifValuesPass) {
1529 //bool unflagged =
1530 unflag();
1531 abcProgress_.clearTimesFlagged();
1532 if (sumDualInfeasibilities_ < 1.0e-3 || (sumDualInfeasibilities_ / static_cast< double >(numberDualInfeasibilities_)) < 1.0e-5 || abcProgress_.lastIterationNumber(0) == numberIterations_) {
1533 if (!numberPrimalInfeasibilities_) {
1534 if (numberTimesOptimal_ < 4) {
1535 numberTimesOptimal_++;
1536 changeMade_++; // say change made
1537 } else {
1538 problemStatus_ = 0;
1539 secondaryStatus_ = 5;
1540 }
1541 }
1542 }
1543 }
1544 }
1545 }
1546 if (problemStatus_ == 0) {
1547 double objVal = (abcNonLinearCost_->feasibleCost()
1548 + objective_->nonlinearOffset());
1549 objVal /= (objectiveScale_ * rhsScale_);
1550 double tol = 1.0e-10 * CoinMax(fabs(objVal), fabs(objectiveValue_)) + 1.0e-8;
1551 if (fabs(objVal - objectiveValue_) > tol) {
1552 #if ABC_NORMAL_DEBUG > 3
1553 if (handler_->logLevel() > 0)
1554 printf("nonLinearCost has feasible obj of %g, objectiveValue_ is %g\n",
1555 objVal, objectiveValue_);
1556 #endif
1557 objectiveValue_ = objVal;
1558 }
1559 }
1560 if (type == 0 || type == 1) {
1561 saveGoodStatus();
1562 }
1563 // see if in Cbc etc
1564 bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
1565 bool disaster = false;
1566 if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
1567 disasterArea_->saveInfo();
1568 disaster = true;
1569 }
1570 if (disaster)
1571 problemStatus_ = 3;
1572 if (problemStatus_ < 0 && !changeMade_) {
1573 problemStatus_ = 4; // unknown
1574 }
1575 lastGoodIteration_ = numberIterations_;
1576 if (numberIterations_ > lastBadIteration_ + 100)
1577 moreSpecialOptions_ &= ~16; // clear check accuracy flag
1578 #ifndef TRY_ABC_GUS
1579 if ((moreSpecialOptions_ & 256) != 0 || saveSum || (specialOptions_ & 8192) != 0)
1580 goToDual = false;
1581 #else
1582 if ((specialOptions_ & 8192) != 0)
1583 goToDual = false;
1584 #endif
1585 if (goToDual || (numberIterations_ > 1000 + baseIteration_ && largestPrimalError_ > 1.0e6 && largestDualError_ > 1.0e6)) {
1586 problemStatus_ = 10; // try dual
1587 // See if second call
1588 if ((moreSpecialOptions_ & 256) != 0 || abcNonLinearCost_->sumInfeasibilities() > 1.0e20) {
1589 numberPrimalInfeasibilities_ = abcNonLinearCost_->numberInfeasibilities();
1590 sumPrimalInfeasibilities_ = abcNonLinearCost_->sumInfeasibilities();
1591 // say infeasible
1592 if (numberPrimalInfeasibilities_ && (abcState_ & CLP_ABC_BEEN_FEASIBLE) == 0)
1593 problemStatus_ = 1;
1594 }
1595 }
1596 // make sure first free monotonic
1597 if (firstFree_ >= 0 && saveFirstFree >= 0) {
1598 firstFree_ = (numberIterations_) ? saveFirstFree : -1;
1599 nextSuperBasic(1, NULL);
1600 }
1601 #ifdef TRY_SPLIT_VALUES_PASS
1602 if (valuesChunk > 0) {
1603 bool doFixing = firstFree_ < 0;
1604 if (numberIterations_ == baseIteration_ && numberDualInfeasibilitiesWithoutFree_ + 1000 < numberDualInfeasibilities_) {
1605 doFixing = true;
1606 int numberSuperBasic = 0;
1607 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1608 if (getInternalStatus(iSequence) == superBasic)
1609 numberSuperBasic++;
1610 }
1611 keepValuesPass = static_cast< int >((1.01 * numberSuperBasic) / valuesChunk);
1612 doOrdinaryVariables = static_cast< int >(valuesRatio * keepValuesPass);
1613 } else if (valuesStop > 0) {
1614 if (numberIterations_ >= valuesStop || problemStatus_ >= 0) {
1615 gutsOfSolution(3);
1616 abcNonLinearCost_->refreshFromPerturbed(primalTolerance_);
1617 gutsOfSolution(3);
1618 int numberSuperBasic = 0;
1619 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1620 if (getInternalStatus(iSequence) == isFixed) {
1621 if (abcUpper_[iSequence] > abcLower_[iSequence]) {
1622 setInternalStatus(iSequence, superBasic);
1623 numberSuperBasic++;
1624 }
1625 }
1626 }
1627 if (numberSuperBasic) {
1628 stateOfProblem_ |= VALUES_PASS;
1629 problemStatus_ = -1;
1630 gutsOfSolution(3);
1631 } else {
1632 doFixing = false;
1633 }
1634 valuesStop = -1;
1635 } else {
1636 doFixing = false;
1637 }
1638 }
1639 if (doFixing) {
1640 abcNonLinearCost_->refreshFromPerturbed(primalTolerance_);
1641 int numberSuperBasic = 0;
1642 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1643 if (getInternalStatus(iSequence) == isFixed) {
1644 if (abcUpper_[iSequence] > abcLower_[iSequence])
1645 setInternalStatus(iSequence, superBasic);
1646 }
1647 if (getInternalStatus(iSequence) == superBasic)
1648 numberSuperBasic++;
1649 }
1650 int numberFixed = 0;
1651 firstFree_ = -1;
1652 if (numberSuperBasic) {
1653 double threshold = static_cast< double >(keepValuesPass) / numberSuperBasic;
1654 if (1.1 * keepValuesPass > numberSuperBasic)
1655 threshold = 1.1;
1656 for (int iSequence = 0; iSequence < numberTotal_; iSequence++) {
1657 if (getInternalStatus(iSequence) == superBasic) {
1658 if (randomNumberGenerator_.randomDouble() < threshold) {
1659 if (firstFree_ < 0)
1660 firstFree_ = iSequence;
1661 } else {
1662 numberFixed++;
1663 abcLower_[iSequence] = abcSolution_[iSequence];
1664 abcUpper_[iSequence] = abcSolution_[iSequence];
1665 setInternalStatus(iSequence, isFixed);
1666 }
1667 }
1668 }
1669 }
1670 if (!numberFixed) {
1671 // end
1672 valuesChunk = -1;
1673 }
1674 }
1675 }
1676 #endif
1677 if (doFactorization || doWeights) {
1678 // restore weights (if saved) - also recompute infeasibility list
1679 if (tentativeStatus > -3)
1680 abcPrimalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4);
1681 else
1682 abcPrimalColumnPivot_->saveWeights(this, 3);
1683 }
1684 }
1685 // Computes solutions - 1 do duals, 2 do primals, 3 both
gutsOfPrimalSolution(int type)1686 int AbcSimplex::gutsOfPrimalSolution(int type)
1687 {
1688 //work space
1689 int whichArray[2];
1690 for (int i = 0; i < 2; i++)
1691 whichArray[i] = getAvailableArray();
1692 CoinIndexedVector *array1 = &usefulArray_[whichArray[0]];
1693 CoinIndexedVector *array2 = &usefulArray_[whichArray[1]];
1694 // do work and check
1695 int numberRefinements = 0;
1696 if ((type & 2) != 0) {
1697 numberRefinements = computePrimals(array1, array2);
1698 if (algorithm_ > 0 && abcNonLinearCost_ != NULL) {
1699 // primal algorithm
1700 // get correct bounds on all variables
1701 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
1702 if (abcNonLinearCost_->numberInfeasibilities())
1703 if (handler_->detail(CLP_SIMPLEX_NONLINEAR, messages_) < 100) {
1704 handler_->message(CLP_SIMPLEX_NONLINEAR, messages_)
1705 << abcNonLinearCost_->changeInCost()
1706 << abcNonLinearCost_->numberInfeasibilities()
1707 << CoinMessageEol;
1708 }
1709 }
1710 checkPrimalSolution(true);
1711 }
1712 if ((type & 1) != 0
1713 #if CAN_HAVE_ZERO_OBJ > 1
1714 && (specialOptions_ & 2097152) == 0
1715 #endif
1716 ) {
1717 numberRefinements += computeDuals(NULL, array1, array2);
1718 checkDualSolution();
1719 }
1720 for (int i = 0; i < 2; i++)
1721 setAvailableArray(whichArray[i]);
1722 rawObjectiveValue_ += sumNonBasicCosts_;
1723 //computeObjective(); // ? done in checkDualSolution??
1724 setClpSimplexObjectiveValue();
1725 if (handler_->logLevel() > 3 || (largestPrimalError_ > 1.0e-2 || largestDualError_ > 1.0e-2))
1726 handler_->message(CLP_SIMPLEX_ACCURACY, messages_)
1727 << largestPrimalError_
1728 << largestDualError_
1729 << CoinMessageEol;
1730 if ((largestPrimalError_ > 1.0e1 || largestDualError_ > 1.0e1)
1731 && numberRows_ > 100 && numberIterations_) {
1732 #if ABC_NORMAL_DEBUG > 0
1733 printf("Large errors - primal %g dual %g\n",
1734 largestPrimalError_, largestDualError_);
1735 #endif
1736 // Change factorization tolerance
1737 //if (abcFactorization_->zeroTolerance() > 1.0e-18)
1738 //abcFactorization_->zeroTolerance(1.0e-18);
1739 }
1740 return numberRefinements;
1741 }
1742 /*
1743 Row array has pivot column
1744 This chooses pivot row.
1745 For speed, we may need to go to a bucket approach when many
1746 variables go through bounds
1747 On exit rhsArray will have changes in costs of basic variables
1748 */
primalRow(CoinIndexedVector * rowArray,CoinIndexedVector * rhsArray,CoinIndexedVector * spareArray,int valuesPass)1749 void AbcSimplexPrimal::primalRow(CoinIndexedVector *rowArray,
1750 CoinIndexedVector *rhsArray,
1751 CoinIndexedVector *spareArray,
1752 int valuesPass)
1753 {
1754 #if 1
1755 for (int iRow = 0; iRow < numberRows_; iRow++) {
1756 int iPivot = abcPivotVariable_[iRow];
1757 assert(costBasic_[iRow] == abcCost_[iPivot]);
1758 assert(lowerBasic_[iRow] == abcLower_[iPivot]);
1759 assert(upperBasic_[iRow] == abcUpper_[iPivot]);
1760 assert(solutionBasic_[iRow] == abcSolution_[iPivot]);
1761 }
1762 #endif
1763 double saveDj = dualIn_;
1764 if (valuesPass) {
1765 dualIn_ = abcCost_[sequenceIn_];
1766
1767 double *work = rowArray->denseVector();
1768 int number = rowArray->getNumElements();
1769 int *which = rowArray->getIndices();
1770
1771 int iIndex;
1772 for (iIndex = 0; iIndex < number; iIndex++) {
1773
1774 int iRow = which[iIndex];
1775 double alpha = work[iRow];
1776 dualIn_ -= alpha * costBasic_[iRow];
1777 }
1778 // determine direction here
1779 if (dualIn_ < -dualTolerance_) {
1780 directionIn_ = 1;
1781 } else if (dualIn_ > dualTolerance_) {
1782 directionIn_ = -1;
1783 } else {
1784 // towards nearest bound
1785 if (valueIn_ - lowerIn_ < upperIn_ - valueIn_) {
1786 directionIn_ = -1;
1787 dualIn_ = dualTolerance_;
1788 } else {
1789 directionIn_ = 1;
1790 dualIn_ = -dualTolerance_;
1791 }
1792 }
1793 }
1794
1795 // sequence stays as row number until end
1796 pivotRow_ = -1;
1797 int numberRemaining = 0;
1798
1799 double totalThru = 0.0; // for when variables flip
1800 // Allow first few iterations to take tiny
1801 double acceptablePivot = 1.0e-1 * acceptablePivot_;
1802 if (numberIterations_ > 100)
1803 acceptablePivot = acceptablePivot_;
1804 if (abcFactorization_->pivots() > 10)
1805 acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
1806 else if (abcFactorization_->pivots() > 5)
1807 acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
1808 else if (abcFactorization_->pivots())
1809 acceptablePivot = acceptablePivot_; // relax
1810 double bestEverPivot = acceptablePivot;
1811 int lastPivotRow = -1;
1812 double lastPivot = 0.0;
1813 double lastTheta = 1.0e50;
1814
1815 // use spareArrays to put ones looked at in
1816 // First one is list of candidates
1817 // We could compress if we really know we won't need any more
1818 // Second array has current set of pivot candidates
1819 // with a backup list saved in double * part of indexed vector
1820
1821 // pivot elements
1822 double *spare;
1823 // indices
1824 int *index;
1825 spareArray->clear();
1826 spare = spareArray->denseVector();
1827 index = spareArray->getIndices();
1828
1829 // we also need somewhere for effective rhs
1830 double *rhs = rhsArray->denseVector();
1831 // and we can use indices to point to alpha
1832 // that way we can store fabs(alpha)
1833 int *indexPoint = rhsArray->getIndices();
1834 //int numberFlip=0; // Those which may change if flips
1835
1836 /*
1837 First we get a list of possible pivots. We can also see if the
1838 problem looks unbounded.
1839
1840 At first we increase theta and see what happens. We start
1841 theta at a reasonable guess. If in right area then we do bit by bit.
1842 We save possible pivot candidates
1843
1844 */
1845
1846 // do first pass to get possibles
1847 // We can also see if unbounded
1848
1849 double *work = rowArray->denseVector();
1850 int number = rowArray->getNumElements();
1851 int *which = rowArray->getIndices();
1852
1853 // we need to swap sign if coming in from ub
1854 double way = directionIn_;
1855 double maximumMovement;
1856 if (way > 0.0)
1857 maximumMovement = CoinMin(1.0e30, upperIn_ - valueIn_);
1858 else
1859 maximumMovement = CoinMin(1.0e30, valueIn_ - lowerIn_);
1860
1861 double averageTheta = abcNonLinearCost_->averageTheta();
1862 double tentativeTheta = CoinMin(10.0 * averageTheta, maximumMovement);
1863 double upperTheta = maximumMovement;
1864 if (tentativeTheta > 0.5 * maximumMovement)
1865 tentativeTheta = maximumMovement;
1866 bool thetaAtMaximum = tentativeTheta == maximumMovement;
1867 // In case tiny bounds increase
1868 if (maximumMovement < 1.0)
1869 tentativeTheta *= 1.1;
1870 double dualCheck = fabs(dualIn_);
1871 // but make a bit more pessimistic
1872 dualCheck = CoinMax(dualCheck - 100.0 * dualTolerance_, 0.99 * dualCheck);
1873
1874 int iIndex;
1875 //#define CLP_DEBUG
1876 #if ABC_NORMAL_DEBUG > 3
1877 if (numberIterations_ == 1369 || numberIterations_ == -3840) {
1878 double dj = abcCost_[sequenceIn_];
1879 printf("cost in on %d is %g, dual in %g\n", sequenceIn_, dj, dualIn_);
1880 for (iIndex = 0; iIndex < number; iIndex++) {
1881
1882 int iRow = which[iIndex];
1883 double alpha = work[iRow];
1884 int iPivot = abcPivotVariable_[iRow];
1885 dj -= alpha * costBasic_[iRow];
1886 printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n",
1887 iRow, iPivot, lowerBasic_[iRow], solutionBasic_[iRow], upperBasic_[iRow],
1888 alpha, solutionBasic_[iRow] - 1.0e9 * alpha, costBasic_[iRow], dj);
1889 }
1890 }
1891 #endif
1892 while (true) {
1893 totalThru = 0.0;
1894 // We also re-compute reduced cost
1895 numberRemaining = 0;
1896 dualIn_ = abcCost_[sequenceIn_];
1897 #ifndef NDEBUG
1898 //double tolerance = primalTolerance_ * 1.002;
1899 #endif
1900 for (iIndex = 0; iIndex < number; iIndex++) {
1901
1902 int iRow = which[iIndex];
1903 double alpha = work[iRow];
1904 dualIn_ -= alpha * costBasic_[iRow];
1905 alpha *= way;
1906 double oldValue = solutionBasic_[iRow];
1907 // get where in bound sequence
1908 // note that after this alpha is actually fabs(alpha)
1909 bool possible;
1910 // do computation same way as later on in primal
1911 if (alpha > 0.0) {
1912 // basic variable going towards lower bound
1913 double bound = lowerBasic_[iRow];
1914 // must be exactly same as when used
1915 double change = tentativeTheta * alpha;
1916 possible = (oldValue - change) <= bound + primalTolerance_;
1917 oldValue -= bound;
1918 } else {
1919 // basic variable going towards upper bound
1920 double bound = upperBasic_[iRow];
1921 // must be exactly same as when used
1922 double change = tentativeTheta * alpha;
1923 possible = (oldValue - change) >= bound - primalTolerance_;
1924 oldValue = bound - oldValue;
1925 alpha = -alpha;
1926 }
1927 double value;
1928 //assert (oldValue >= -10.0*tolerance);
1929 if (possible) {
1930 value = oldValue - upperTheta * alpha;
1931 #ifdef CLP_USER_DRIVEN1
1932 if (!userChoiceValid1(this, iPivot, oldValue,
1933 upperTheta, alpha, work[iRow] * way))
1934 value = 0.0; // say can't use
1935 #endif
1936 if (value < -primalTolerance_ && alpha >= acceptablePivot) {
1937 upperTheta = (oldValue + primalTolerance_) / alpha;
1938 }
1939 // add to list
1940 spare[numberRemaining] = alpha;
1941 rhs[numberRemaining] = oldValue;
1942 indexPoint[numberRemaining] = iIndex;
1943 index[numberRemaining++] = iRow;
1944 totalThru += alpha;
1945 setActive(iRow);
1946 //} else if (value<primalTolerance_*1.002) {
1947 // May change if is a flip
1948 //indexRhs[numberFlip++]=iRow;
1949 }
1950 }
1951 if (upperTheta < maximumMovement && totalThru * infeasibilityCost_ >= 1.0001 * dualCheck) {
1952 // Can pivot here
1953 break;
1954 } else if (!thetaAtMaximum) {
1955 //printf("Going round with average theta of %g\n",averageTheta);
1956 tentativeTheta = maximumMovement;
1957 thetaAtMaximum = true; // seems to be odd compiler error
1958 } else {
1959 break;
1960 }
1961 }
1962 totalThru = 0.0;
1963
1964 theta_ = maximumMovement;
1965
1966 bool goBackOne = false;
1967 if (objective_->type() > 1)
1968 dualIn_ = saveDj;
1969
1970 //printf("%d remain out of %d\n",numberRemaining,number);
1971 int iTry = 0;
1972 #define MAXTRY 1000
1973 if (numberRemaining && upperTheta < maximumMovement) {
1974 // First check if previously chosen one will work
1975
1976 // first get ratio with tolerance
1977 for (; iTry < MAXTRY; iTry++) {
1978
1979 upperTheta = maximumMovement;
1980 int iBest = -1;
1981 for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
1982
1983 double alpha = spare[iIndex];
1984 double oldValue = rhs[iIndex];
1985 double value = oldValue - upperTheta * alpha;
1986
1987 #ifdef CLP_USER_DRIVEN1
1988 int sequenceOut = abcPivotVariable_[index[iIndex]];
1989 if (!userChoiceValid1(this, sequenceOut, oldValue,
1990 upperTheta, alpha, 0.0))
1991 value = 0.0; // say can't use
1992 #endif
1993 if (value < -primalTolerance_ && alpha >= acceptablePivot) {
1994 upperTheta = (oldValue + primalTolerance_) / alpha;
1995 iBest = iIndex; // just in case weird numbers
1996 }
1997 }
1998
1999 // now look at best in this lot
2000 // But also see how infeasible small pivots will make
2001 double sumInfeasibilities = 0.0;
2002 double bestPivot = acceptablePivot;
2003 pivotRow_ = -1;
2004 for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
2005
2006 int iRow = index[iIndex];
2007 double alpha = spare[iIndex];
2008 double oldValue = rhs[iIndex];
2009 double value = oldValue - upperTheta * alpha;
2010
2011 if (value <= 0 || iBest == iIndex) {
2012 // how much would it cost to go thru and modify bound
2013 double trueAlpha = way * work[iRow];
2014 totalThru += abcNonLinearCost_->changeInCost(iRow, trueAlpha, rhs[iIndex]);
2015 setActive(iRow);
2016 if (alpha > bestPivot) {
2017 bestPivot = alpha;
2018 theta_ = oldValue / bestPivot;
2019 pivotRow_ = iRow;
2020 } else if (alpha < acceptablePivot
2021 #ifdef CLP_USER_DRIVEN1
2022 || !userChoiceValid1(this, abcPivotVariable_[index[iIndex]],
2023 oldValue, upperTheta, alpha, 0.0)
2024 #endif
2025 ) {
2026 if (value < -primalTolerance_)
2027 sumInfeasibilities += -value - primalTolerance_;
2028 }
2029 }
2030 }
2031 if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e-6 && bestPivot < 1.0e-3) {
2032 // back to previous one
2033 goBackOne = true;
2034 break;
2035 } else if (pivotRow_ == -1 && upperTheta > largeValue_) {
2036 if (lastPivot > acceptablePivot) {
2037 // back to previous one
2038 goBackOne = true;
2039 } else {
2040 // can only get here if all pivots so far too small
2041 }
2042 break;
2043 } else if (totalThru >= dualCheck) {
2044 if (sumInfeasibilities > primalTolerance_ && !abcNonLinearCost_->numberInfeasibilities()) {
2045 // Looks a bad choice
2046 if (lastPivot > acceptablePivot) {
2047 goBackOne = true;
2048 } else {
2049 // say no good
2050 dualIn_ = 0.0;
2051 }
2052 }
2053 break; // no point trying another loop
2054 } else {
2055 lastPivotRow = pivotRow_;
2056 lastTheta = theta_;
2057 if (bestPivot > bestEverPivot)
2058 bestEverPivot = bestPivot;
2059 }
2060 }
2061 // can get here without pivotRow_ set but with lastPivotRow
2062 if (goBackOne || (pivotRow_ < 0 && lastPivotRow >= 0)) {
2063 // back to previous one
2064 pivotRow_ = lastPivotRow;
2065 theta_ = lastTheta;
2066 }
2067 } else if (pivotRow_ < 0 && maximumMovement > 1.0e20) {
2068 // looks unbounded
2069 valueOut_ = COIN_DBL_MAX; // say odd
2070 if (abcNonLinearCost_->numberInfeasibilities()) {
2071 // but infeasible??
2072 // move variable but don't pivot
2073 tentativeTheta = 1.0e50;
2074 for (iIndex = 0; iIndex < number; iIndex++) {
2075 int iRow = which[iIndex];
2076 double alpha = work[iRow];
2077 alpha *= way;
2078 double oldValue = solutionBasic_[iRow];
2079 // get where in bound sequence
2080 // note that after this alpha is actually fabs(alpha)
2081 if (alpha > 0.0) {
2082 // basic variable going towards lower bound
2083 double bound = lowerBasic_[iRow];
2084 oldValue -= bound;
2085 } else {
2086 // basic variable going towards upper bound
2087 double bound = upperBasic_[iRow];
2088 oldValue = bound - oldValue;
2089 alpha = -alpha;
2090 }
2091 if (oldValue - tentativeTheta * alpha < 0.0) {
2092 tentativeTheta = oldValue / alpha;
2093 }
2094 }
2095 // If free in then see if we can get to 0.0
2096 if (lowerIn_ < -1.0e20 && upperIn_ > 1.0e20) {
2097 if (dualIn_ * valueIn_ > 0.0) {
2098 if (fabs(valueIn_) < 1.0e-2 && (tentativeTheta < fabs(valueIn_) || tentativeTheta > 1.0e20)) {
2099 tentativeTheta = fabs(valueIn_);
2100 }
2101 }
2102 }
2103 if (tentativeTheta < 1.0e10)
2104 valueOut_ = valueIn_ + way * tentativeTheta;
2105 }
2106 }
2107 //if (iTry>50)
2108 //printf("** %d tries\n",iTry);
2109 if (pivotRow_ >= 0) {
2110 alpha_ = work[pivotRow_];
2111 // translate to sequence
2112 sequenceOut_ = abcPivotVariable_[pivotRow_];
2113 valueOut_ = solution(sequenceOut_);
2114 lowerOut_ = abcLower_[sequenceOut_];
2115 upperOut_ = abcUpper_[sequenceOut_];
2116 //#define MINIMUMTHETA 1.0e-11
2117 // Movement should be minimum for anti-degeneracy - unless
2118 // fixed variable out
2119 double minimumTheta;
2120 if (upperOut_ > lowerOut_)
2121 minimumTheta = minimumThetaMovement_;
2122 else
2123 minimumTheta = 0.0;
2124 // But can't go infeasible
2125 double distance;
2126 if (alpha_ * way > 0.0)
2127 distance = valueOut_ - lowerOut_;
2128 else
2129 distance = upperOut_ - valueOut_;
2130 if (distance - minimumTheta * fabs(alpha_) < -primalTolerance_)
2131 minimumTheta = CoinMax(0.0, (distance + 0.5 * primalTolerance_) / fabs(alpha_));
2132 // will we need to increase tolerance
2133 //#define CLP_DEBUG
2134 double largestInfeasibility = primalTolerance_;
2135 if (theta_ < minimumTheta && (specialOptions_ & 4) == 0 && !valuesPass) {
2136 theta_ = minimumTheta;
2137 for (iIndex = 0; iIndex < numberRemaining - numberRemaining; iIndex++) {
2138 largestInfeasibility = CoinMax(largestInfeasibility,
2139 -(rhs[iIndex] - spare[iIndex] * theta_));
2140 }
2141 //#define CLP_DEBUG
2142 #if ABC_NORMAL_DEBUG > 3
2143 if (largestInfeasibility > primalTolerance_ && (handler_->logLevel() & 32) > -1)
2144 printf("Primal tolerance increased from %g to %g\n",
2145 primalTolerance_, largestInfeasibility);
2146 #endif
2147 //#undef CLP_DEBUG
2148 primalTolerance_ = CoinMax(primalTolerance_, largestInfeasibility);
2149 }
2150 // Need to look at all in some cases
2151 if (theta_ > tentativeTheta) {
2152 for (iIndex = 0; iIndex < number; iIndex++)
2153 setActive(which[iIndex]);
2154 }
2155 if (way < 0.0)
2156 theta_ = -theta_;
2157 double newValue = valueOut_ - theta_ * alpha_;
2158 // If 4 bit set - Force outgoing variables to exact bound (primal)
2159 if (alpha_ * way < 0.0) {
2160 directionOut_ = -1; // to upper bound
2161 if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
2162 upperOut_ = abcNonLinearCost_->nearest(pivotRow_, newValue);
2163 } else {
2164 upperOut_ = newValue;
2165 }
2166 } else {
2167 directionOut_ = 1; // to lower bound
2168 if (fabs(theta_) > 1.0e-6 || (specialOptions_ & 4) != 0) {
2169 lowerOut_ = abcNonLinearCost_->nearest(pivotRow_, newValue);
2170 } else {
2171 lowerOut_ = newValue;
2172 }
2173 }
2174 dualOut_ = reducedCost(sequenceOut_);
2175 } else if (maximumMovement < 1.0e20) {
2176 // flip
2177 pivotRow_ = -2; // so we can tell its a flip
2178 sequenceOut_ = sequenceIn_;
2179 valueOut_ = valueIn_;
2180 dualOut_ = dualIn_;
2181 lowerOut_ = lowerIn_;
2182 upperOut_ = upperIn_;
2183 alpha_ = 0.0;
2184 if (way < 0.0) {
2185 directionOut_ = 1; // to lower bound
2186 theta_ = lowerOut_ - valueOut_;
2187 } else {
2188 directionOut_ = -1; // to upper bound
2189 theta_ = upperOut_ - valueOut_;
2190 }
2191 }
2192
2193 double theta1 = CoinMax(theta_, 1.0e-12);
2194 double theta2 = numberIterations_ * abcNonLinearCost_->averageTheta();
2195 // Set average theta
2196 abcNonLinearCost_->setAverageTheta((theta1 + theta2) / (static_cast< double >(numberIterations_ + 1)));
2197 //if (numberIterations_%1000==0)
2198 //printf("average theta is %g\n",abcNonLinearCost_->averageTheta());
2199
2200 // clear arrays
2201 CoinZeroN(spare, numberRemaining);
2202 CoinZeroN(rhs, numberRemaining);
2203
2204 // put back original bounds etc
2205 CoinMemcpyN(index, numberRemaining, rhsArray->getIndices());
2206 rhsArray->setNumElements(numberRemaining);
2207 //rhsArray->setPacked();
2208 abcNonLinearCost_->goBackAll(rhsArray);
2209 rhsArray->clear();
2210 }
2211 /*
2212 Row array has pivot column
2213 This chooses pivot row.
2214 For speed, we may need to go to a bucket approach when many
2215 variables go through bounds
2216 On exit rhsArray will have changes in costs of basic variables
2217 */
primalRow(CoinIndexedVector * rowArray,CoinIndexedVector * rhsArray,CoinIndexedVector * spareArray,pivotStruct & stuff)2218 void AbcSimplexPrimal::primalRow(CoinIndexedVector *rowArray,
2219 CoinIndexedVector *rhsArray,
2220 CoinIndexedVector *spareArray,
2221 pivotStruct &stuff)
2222 {
2223 int valuesPass = stuff.valuesPass_;
2224 double dualIn = stuff.dualIn_;
2225 double lowerIn = stuff.lowerIn_;
2226 double upperIn = stuff.upperIn_;
2227 double valueIn = stuff.valueIn_;
2228 int sequenceIn = stuff.sequenceIn_;
2229 int directionIn = stuff.directionIn_;
2230 int pivotRow = -1;
2231 int sequenceOut = -1;
2232 double dualOut = 0.0;
2233 double lowerOut = 0.0;
2234 double upperOut = 0.0;
2235 double valueOut = 0.0;
2236 double theta;
2237 double alpha;
2238 int directionOut = 0;
2239 if (valuesPass) {
2240 dualIn = abcCost_[sequenceIn];
2241
2242 double *work = rowArray->denseVector();
2243 int number = rowArray->getNumElements();
2244 int *which = rowArray->getIndices();
2245
2246 int iIndex;
2247 for (iIndex = 0; iIndex < number; iIndex++) {
2248
2249 int iRow = which[iIndex];
2250 double alpha = work[iRow];
2251 dualIn -= alpha * costBasic_[iRow];
2252 }
2253 // determine direction here
2254 if (dualIn < -dualTolerance_) {
2255 directionIn = 1;
2256 } else if (dualIn > dualTolerance_) {
2257 directionIn = -1;
2258 } else {
2259 // towards nearest bound
2260 if (valueIn - lowerIn < upperIn - valueIn) {
2261 directionIn = -1;
2262 dualIn = dualTolerance_;
2263 } else {
2264 directionIn = 1;
2265 dualIn = -dualTolerance_;
2266 }
2267 }
2268 }
2269
2270 // sequence stays as row number until end
2271 pivotRow = -1;
2272 int numberRemaining = 0;
2273
2274 double totalThru = 0.0; // for when variables flip
2275 // Allow first few iterations to take tiny
2276 double acceptablePivot = 1.0e-1 * acceptablePivot_;
2277 if (numberIterations_ > 100)
2278 acceptablePivot = acceptablePivot_;
2279 if (abcFactorization_->pivots() > 10)
2280 acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
2281 else if (abcFactorization_->pivots() > 5)
2282 acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
2283 else if (abcFactorization_->pivots())
2284 acceptablePivot = acceptablePivot_; // relax
2285 double bestEverPivot = acceptablePivot;
2286 int lastPivotRow = -1;
2287 double lastPivot = 0.0;
2288 double lastTheta = 1.0e50;
2289
2290 // use spareArrays to put ones looked at in
2291 // First one is list of candidates
2292 // We could compress if we really know we won't need any more
2293 // Second array has current set of pivot candidates
2294 // with a backup list saved in double * part of indexed vector
2295
2296 // pivot elements
2297 double *spare;
2298 // indices
2299 int *index;
2300 spareArray->clear();
2301 spare = spareArray->denseVector();
2302 index = spareArray->getIndices();
2303
2304 // we also need somewhere for effective rhs
2305 double *rhs = rhsArray->denseVector();
2306 // and we can use indices to point to alpha
2307 // that way we can store fabs(alpha)
2308 int *indexPoint = rhsArray->getIndices();
2309 //int numberFlip=0; // Those which may change if flips
2310
2311 /*
2312 First we get a list of possible pivots. We can also see if the
2313 problem looks unbounded.
2314
2315 At first we increase theta and see what happens. We start
2316 theta at a reasonable guess. If in right area then we do bit by bit.
2317 We save possible pivot candidates
2318
2319 */
2320
2321 // do first pass to get possibles
2322 // We can also see if unbounded
2323
2324 double *work = rowArray->denseVector();
2325 int number = rowArray->getNumElements();
2326 int *which = rowArray->getIndices();
2327
2328 // we need to swap sign if coming in from ub
2329 double way = directionIn;
2330 double maximumMovement;
2331 if (way > 0.0)
2332 maximumMovement = CoinMin(1.0e30, upperIn - valueIn);
2333 else
2334 maximumMovement = CoinMin(1.0e30, valueIn - lowerIn);
2335
2336 double averageTheta = abcNonLinearCost_->averageTheta();
2337 double tentativeTheta = CoinMin(10.0 * averageTheta, maximumMovement);
2338 double upperTheta = maximumMovement;
2339 if (tentativeTheta > 0.5 * maximumMovement)
2340 tentativeTheta = maximumMovement;
2341 bool thetaAtMaximum = tentativeTheta == maximumMovement;
2342 // In case tiny bounds increase
2343 if (maximumMovement < 1.0)
2344 tentativeTheta *= 1.1;
2345 double dualCheck = fabs(dualIn);
2346 // but make a bit more pessimistic
2347 dualCheck = CoinMax(dualCheck - 100.0 * dualTolerance_, 0.99 * dualCheck);
2348
2349 int iIndex;
2350 while (true) {
2351 totalThru = 0.0;
2352 // We also re-compute reduced cost
2353 numberRemaining = 0;
2354 dualIn = abcCost_[sequenceIn];
2355 for (iIndex = 0; iIndex < number; iIndex++) {
2356 int iRow = which[iIndex];
2357 double alpha = work[iRow];
2358 dualIn -= alpha * costBasic_[iRow];
2359 alpha *= way;
2360 double oldValue = solutionBasic_[iRow];
2361 // get where in bound sequence
2362 // note that after this alpha is actually fabs(alpha)
2363 bool possible;
2364 // do computation same way as later on in primal
2365 if (alpha > 0.0) {
2366 // basic variable going towards lower bound
2367 double bound = lowerBasic_[iRow];
2368 // must be exactly same as when used
2369 double change = tentativeTheta * alpha;
2370 possible = (oldValue - change) <= bound + primalTolerance_;
2371 oldValue -= bound;
2372 } else {
2373 // basic variable going towards upper bound
2374 double bound = upperBasic_[iRow];
2375 // must be exactly same as when used
2376 double change = tentativeTheta * alpha;
2377 possible = (oldValue - change) >= bound - primalTolerance_;
2378 oldValue = bound - oldValue;
2379 alpha = -alpha;
2380 }
2381 double value;
2382 if (possible) {
2383 value = oldValue - upperTheta * alpha;
2384 #ifdef CLP_USER_DRIVEN1
2385 if (!userChoiceValid1(this, iPivot, oldValue,
2386 upperTheta, alpha, work[iRow] * way))
2387 value = 0.0; // say can't use
2388 #endif
2389 if (value < -primalTolerance_ && alpha >= acceptablePivot) {
2390 upperTheta = (oldValue + primalTolerance_) / alpha;
2391 }
2392 // add to list
2393 spare[numberRemaining] = alpha;
2394 rhs[numberRemaining] = oldValue;
2395 indexPoint[numberRemaining] = iIndex;
2396 index[numberRemaining++] = iRow;
2397 totalThru += alpha;
2398 setActive(iRow);
2399 }
2400 }
2401 if (upperTheta < maximumMovement && totalThru * infeasibilityCost_ >= 1.0001 * dualCheck) {
2402 // Can pivot here
2403 break;
2404 } else if (!thetaAtMaximum) {
2405 //printf("Going round with average theta of %g\n",averageTheta);
2406 tentativeTheta = maximumMovement;
2407 thetaAtMaximum = true; // seems to be odd compiler error
2408 } else {
2409 break;
2410 }
2411 }
2412 totalThru = 0.0;
2413
2414 theta = maximumMovement;
2415
2416 bool goBackOne = false;
2417
2418 //printf("%d remain out of %d\n",numberRemaining,number);
2419 int iTry = 0;
2420 if (numberRemaining && upperTheta < maximumMovement) {
2421 // First check if previously chosen one will work
2422
2423 // first get ratio with tolerance
2424 for (; iTry < MAXTRY; iTry++) {
2425
2426 upperTheta = maximumMovement;
2427 int iBest = -1;
2428 for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
2429
2430 double alpha = spare[iIndex];
2431 double oldValue = rhs[iIndex];
2432 double value = oldValue - upperTheta * alpha;
2433
2434 #ifdef CLP_USER_DRIVEN1
2435 int sequenceOut = abcPivotVariable_[index[iIndex]];
2436 if (!userChoiceValid1(this, sequenceOut, oldValue,
2437 upperTheta, alpha, 0.0))
2438 value = 0.0; // say can't use
2439 #endif
2440 if (value < -primalTolerance_ && alpha >= acceptablePivot) {
2441 upperTheta = (oldValue + primalTolerance_) / alpha;
2442 iBest = iIndex; // just in case weird numbers
2443 }
2444 }
2445
2446 // now look at best in this lot
2447 // But also see how infeasible small pivots will make
2448 double sumInfeasibilities = 0.0;
2449 double bestPivot = acceptablePivot;
2450 pivotRow = -1;
2451 for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
2452
2453 int iRow = index[iIndex];
2454 double alpha = spare[iIndex];
2455 double oldValue = rhs[iIndex];
2456 double value = oldValue - upperTheta * alpha;
2457
2458 if (value <= 0 || iBest == iIndex) {
2459 // how much would it cost to go thru and modify bound
2460 double trueAlpha = way * work[iRow];
2461 totalThru += abcNonLinearCost_->changeInCost(iRow, trueAlpha, rhs[iIndex]);
2462 setActive(iRow);
2463 if (alpha > bestPivot) {
2464 bestPivot = alpha;
2465 theta = oldValue / bestPivot;
2466 pivotRow = iRow;
2467 } else if (alpha < acceptablePivot
2468 #ifdef CLP_USER_DRIVEN1
2469 || !userChoiceValid1(this, abcPivotVariable_[index[iIndex]],
2470 oldValue, upperTheta, alpha, 0.0)
2471 #endif
2472 ) {
2473 if (value < -primalTolerance_)
2474 sumInfeasibilities += -value - primalTolerance_;
2475 }
2476 }
2477 }
2478 if (bestPivot < 0.1 * bestEverPivot && bestEverPivot > 1.0e-6 && bestPivot < 1.0e-3) {
2479 // back to previous one
2480 goBackOne = true;
2481 break;
2482 } else if (pivotRow == -1 && upperTheta > largeValue_) {
2483 if (lastPivot > acceptablePivot) {
2484 // back to previous one
2485 goBackOne = true;
2486 } else {
2487 // can only get here if all pivots so far too small
2488 }
2489 break;
2490 } else if (totalThru >= dualCheck) {
2491 if (sumInfeasibilities > primalTolerance_ && !abcNonLinearCost_->numberInfeasibilities()) {
2492 // Looks a bad choice
2493 if (lastPivot > acceptablePivot) {
2494 goBackOne = true;
2495 } else {
2496 // say no good
2497 dualIn = 0.0;
2498 }
2499 }
2500 break; // no point trying another loop
2501 } else {
2502 lastPivotRow = pivotRow;
2503 lastTheta = theta;
2504 if (bestPivot > bestEverPivot)
2505 bestEverPivot = bestPivot;
2506 }
2507 }
2508 // can get here without pivotRow set but with lastPivotRow
2509 if (goBackOne || (pivotRow < 0 && lastPivotRow >= 0)) {
2510 // back to previous one
2511 pivotRow = lastPivotRow;
2512 theta = lastTheta;
2513 }
2514 } else if (pivotRow < 0 && maximumMovement > 1.0e20) {
2515 // looks unbounded
2516 valueOut = COIN_DBL_MAX; // say odd
2517 if (abcNonLinearCost_->numberInfeasibilities()) {
2518 // but infeasible??
2519 // move variable but don't pivot
2520 tentativeTheta = 1.0e50;
2521 for (iIndex = 0; iIndex < number; iIndex++) {
2522 int iRow = which[iIndex];
2523 double alpha = work[iRow];
2524 alpha *= way;
2525 double oldValue = solutionBasic_[iRow];
2526 // get where in bound sequence
2527 // note that after this alpha is actually fabs(alpha)
2528 if (alpha > 0.0) {
2529 // basic variable going towards lower bound
2530 double bound = lowerBasic_[iRow];
2531 oldValue -= bound;
2532 } else {
2533 // basic variable going towards upper bound
2534 double bound = upperBasic_[iRow];
2535 oldValue = bound - oldValue;
2536 alpha = -alpha;
2537 }
2538 if (oldValue - tentativeTheta * alpha < 0.0) {
2539 tentativeTheta = oldValue / alpha;
2540 }
2541 }
2542 // If free in then see if we can get to 0.0
2543 if (lowerIn < -1.0e20 && upperIn > 1.0e20) {
2544 if (dualIn * valueIn > 0.0) {
2545 if (fabs(valueIn) < 1.0e-2 && (tentativeTheta < fabs(valueIn) || tentativeTheta > 1.0e20)) {
2546 tentativeTheta = fabs(valueIn);
2547 }
2548 }
2549 }
2550 if (tentativeTheta < 1.0e10)
2551 valueOut = valueIn + way * tentativeTheta;
2552 }
2553 }
2554 if (pivotRow >= 0) {
2555 alpha = work[pivotRow];
2556 // translate to sequence
2557 sequenceOut = abcPivotVariable_[pivotRow];
2558 valueOut = solutionBasic_[pivotRow];
2559 lowerOut = lowerBasic_[pivotRow];
2560 upperOut = upperBasic_[pivotRow];
2561 // Movement should be minimum for anti-degeneracy - unless
2562 // fixed variable out
2563 double minimumTheta;
2564 if (upperOut > lowerOut)
2565 minimumTheta = minimumThetaMovement_;
2566 else
2567 minimumTheta = 0.0;
2568 // But can't go infeasible
2569 double distance;
2570 if (alpha * way > 0.0)
2571 distance = valueOut - lowerOut;
2572 else
2573 distance = upperOut - valueOut;
2574 if (distance - minimumTheta * fabs(alpha) < -primalTolerance_)
2575 minimumTheta = CoinMax(0.0, (distance + 0.5 * primalTolerance_) / fabs(alpha));
2576 // will we need to increase tolerance
2577 //double largestInfeasibility = primalTolerance_;
2578 if (theta < minimumTheta && (specialOptions_ & 4) == 0 && !valuesPass) {
2579 theta = minimumTheta;
2580 //for (iIndex = 0; iIndex < numberRemaining - numberRemaining; iIndex++) {
2581 //largestInfeasibility = CoinMax(largestInfeasibility,
2582 // -(rhs[iIndex] - spare[iIndex] * theta));
2583 //}
2584 //primalTolerance_ = CoinMax(primalTolerance_, largestInfeasibility);
2585 }
2586 // Need to look at all in some cases
2587 if (theta > tentativeTheta) {
2588 for (iIndex = 0; iIndex < number; iIndex++)
2589 setActive(which[iIndex]);
2590 }
2591 if (way < 0.0)
2592 theta = -theta;
2593 double newValue = valueOut - theta * alpha;
2594 // If 4 bit set - Force outgoing variables to exact bound (primal)
2595 if (alpha * way < 0.0) {
2596 directionOut = -1; // to upper bound
2597 if (fabs(theta) > 1.0e-6 || (specialOptions_ & 4) != 0) {
2598 upperOut = abcNonLinearCost_->nearest(pivotRow, newValue);
2599 } else {
2600 upperOut = newValue;
2601 }
2602 } else {
2603 directionOut = 1; // to lower bound
2604 if (fabs(theta) > 1.0e-6 || (specialOptions_ & 4) != 0) {
2605 lowerOut = abcNonLinearCost_->nearest(pivotRow, newValue);
2606 } else {
2607 lowerOut = newValue;
2608 }
2609 }
2610 dualOut = reducedCost(sequenceOut);
2611 } else if (maximumMovement < 1.0e20) {
2612 // flip
2613 pivotRow = -2; // so we can tell its a flip
2614 sequenceOut = sequenceIn;
2615 valueOut = valueIn;
2616 dualOut = dualIn;
2617 lowerOut = lowerIn;
2618 upperOut = upperIn;
2619 alpha = 0.0;
2620 if (way < 0.0) {
2621 directionOut = 1; // to lower bound
2622 theta = lowerOut - valueOut;
2623 } else {
2624 directionOut = -1; // to upper bound
2625 theta = upperOut - valueOut;
2626 }
2627 }
2628
2629 // clear arrays
2630 CoinZeroN(spare, numberRemaining);
2631 CoinZeroN(rhs, numberRemaining);
2632
2633 // put back original bounds etc
2634 CoinMemcpyN(index, numberRemaining, rhsArray->getIndices());
2635 rhsArray->setNumElements(numberRemaining);
2636 abcNonLinearCost_->goBackAll(rhsArray);
2637 rhsArray->clear();
2638 stuff.theta_ = theta;
2639 stuff.alpha_ = alpha;
2640 stuff.dualIn_ = dualIn;
2641 stuff.dualOut_ = dualOut;
2642 stuff.lowerOut_ = lowerOut;
2643 stuff.upperOut_ = upperOut;
2644 stuff.valueOut_ = valueOut;
2645 stuff.sequenceOut_ = sequenceOut;
2646 stuff.directionOut_ = directionOut;
2647 stuff.pivotRow_ = pivotRow;
2648 }
2649 /*
2650 Chooses primal pivot column
2651 updateArray has cost updates (also use pivotRow_ from last iteration)
2652 Would be faster with separate region to scan
2653 and will have this (with square of infeasibility) when steepest
2654 For easy problems we can just choose one of the first columns we look at
2655 */
primalColumn(CoinPartitionedVector * updates,CoinPartitionedVector * spareRow2,CoinPartitionedVector * spareColumn1)2656 void AbcSimplexPrimal::primalColumn(CoinPartitionedVector *updates,
2657 CoinPartitionedVector *spareRow2,
2658 CoinPartitionedVector *spareColumn1)
2659 {
2660 for (int i = 0; i < 4; i++)
2661 multipleSequenceIn_[i] = -1;
2662 sequenceIn_ = abcPrimalColumnPivot_->pivotColumn(updates,
2663 spareRow2, spareColumn1);
2664 if (sequenceIn_ >= 0) {
2665 valueIn_ = abcSolution_[sequenceIn_];
2666 dualIn_ = abcDj_[sequenceIn_];
2667 lowerIn_ = abcLower_[sequenceIn_];
2668 upperIn_ = abcUpper_[sequenceIn_];
2669 if (dualIn_ > 0.0)
2670 directionIn_ = -1;
2671 else
2672 directionIn_ = 1;
2673 } else {
2674 sequenceIn_ = -1;
2675 }
2676 }
2677 /* The primals are updated by the given array.
2678 Returns number of infeasibilities.
2679 After rowArray will have list of cost changes
2680 */
updatePrimalsInPrimal(CoinIndexedVector * rowArray,double theta,double & objectiveChange,int valuesPass)2681 int AbcSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector *rowArray,
2682 double theta,
2683 double &objectiveChange,
2684 int valuesPass)
2685 {
2686 // Cost on pivot row may change - may need to change dualIn
2687 double oldCost = 0.0;
2688 if (pivotRow_ >= 0)
2689 oldCost = abcCost_[sequenceOut_];
2690 double *work = rowArray->denseVector();
2691 int number = rowArray->getNumElements();
2692 int *which = rowArray->getIndices();
2693
2694 int newNumber = 0;
2695 abcNonLinearCost_->setChangeInCost(0.0);
2696 // allow for case where bound+tolerance == bound
2697 //double tolerance = 0.999999*primalTolerance_;
2698 double relaxedTolerance = 1.001 * primalTolerance_;
2699 int iIndex;
2700 if (!valuesPass) {
2701 for (iIndex = 0; iIndex < number; iIndex++) {
2702
2703 int iRow = which[iIndex];
2704 double alpha = work[iRow];
2705 work[iRow] = 0.0;
2706 int iPivot = abcPivotVariable_[iRow];
2707 double change = theta * alpha;
2708 double value = abcSolution_[iPivot] - change;
2709 abcSolution_[iPivot] = value;
2710 value = solutionBasic_[iRow] - change;
2711 solutionBasic_[iRow] = value;
2712 #ifndef NDEBUG
2713 // check if not active then okay
2714 if (!active(iRow) && (specialOptions_ & 4) == 0 && pivotRow_ != -1) {
2715 // But make sure one going out is feasible
2716 if (change > 0.0) {
2717 // going down
2718 if (value <= abcLower_[iPivot] + primalTolerance_) {
2719 if (iPivot == sequenceOut_ && value > abcLower_[iPivot] - relaxedTolerance)
2720 value = abcLower_[iPivot];
2721 //double difference = abcNonLinearCost_->setOneBasic(iRow, value);
2722 //assert (!difference || fabs(change) > 1.0e9);
2723 }
2724 } else {
2725 // going up
2726 if (value >= abcUpper_[iPivot] - primalTolerance_) {
2727 if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance)
2728 value = abcUpper_[iPivot];
2729 //double difference = abcNonLinearCost_->setOneBasic(iRow, value);
2730 //assert (!difference || fabs(change) > 1.0e9);
2731 }
2732 }
2733 }
2734 #endif
2735 if (active(iRow) || theta_ < 0.0) {
2736 clearActive(iRow);
2737 // But make sure one going out is feasible
2738 if (change > 0.0) {
2739 // going down
2740 if (value <= abcLower_[iPivot] + primalTolerance_) {
2741 if (iPivot == sequenceOut_ && value >= abcLower_[iPivot] - relaxedTolerance)
2742 value = abcLower_[iPivot];
2743 double difference = abcNonLinearCost_->setOneBasic(iRow, value);
2744 if (difference) {
2745 work[iRow] = difference;
2746 //change reduced cost on this
2747 //abcDj_[iPivot] = -difference;
2748 which[newNumber++] = iRow;
2749 }
2750 }
2751 } else {
2752 // going up
2753 if (value >= abcUpper_[iPivot] - primalTolerance_) {
2754 if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance)
2755 value = abcUpper_[iPivot];
2756 double difference = abcNonLinearCost_->setOneBasic(iRow, value);
2757 if (difference) {
2758 work[iRow] = difference;
2759 //change reduced cost on this
2760 //abcDj_[iPivot] = -difference;
2761 which[newNumber++] = iRow;
2762 }
2763 }
2764 }
2765 }
2766 }
2767 } else {
2768 // values pass so look at all
2769 for (iIndex = 0; iIndex < number; iIndex++) {
2770
2771 int iRow = which[iIndex];
2772 double alpha = work[iRow];
2773 work[iRow] = 0.0;
2774 int iPivot = abcPivotVariable_[iRow];
2775 double change = theta * alpha;
2776 double value = abcSolution_[iPivot] - change;
2777 abcSolution_[iPivot] = value;
2778 value = solutionBasic_[iRow] - change;
2779 solutionBasic_[iRow] = value;
2780 clearActive(iRow);
2781 // But make sure one going out is feasible
2782 if (change > 0.0) {
2783 // going down
2784 if (value <= abcLower_[iPivot] + primalTolerance_) {
2785 if (iPivot == sequenceOut_ && value > abcLower_[iPivot] - relaxedTolerance)
2786 value = abcLower_[iPivot];
2787 double difference = abcNonLinearCost_->setOneBasic(iRow, value);
2788 if (difference) {
2789 work[iRow] = difference;
2790 //change reduced cost on this
2791 abcDj_[iPivot] = -difference;
2792 which[newNumber++] = iRow;
2793 }
2794 }
2795 } else {
2796 // going up
2797 if (value >= abcUpper_[iPivot] - primalTolerance_) {
2798 if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance)
2799 value = abcUpper_[iPivot];
2800 double difference = abcNonLinearCost_->setOneBasic(iRow, value);
2801 if (difference) {
2802 work[iRow] = difference;
2803 //change reduced cost on this
2804 abcDj_[iPivot] = -difference;
2805 which[newNumber++] = iRow;
2806 }
2807 }
2808 }
2809 }
2810 }
2811 objectiveChange += abcNonLinearCost_->changeInCost();
2812 //rowArray->setPacked();
2813 if (pivotRow_ >= 0) {
2814 double dualIn = dualIn_ + (oldCost - abcCost_[sequenceOut_]);
2815 // update change vector to include pivot
2816 if (!work[pivotRow_])
2817 which[newNumber++] = pivotRow_;
2818 work[pivotRow_] -= dualIn;
2819 // above seems marginally better for updates - below should also work
2820 //work[pivotRow_] = -dualIn_;
2821 }
2822 rowArray->setNumElements(newNumber);
2823 return 0;
2824 }
2825 // Perturbs problem
perturb(int)2826 void AbcSimplexPrimal::perturb(int /*type*/)
2827 {
2828 if (perturbation_ > 100)
2829 return; //perturbed already
2830 if (perturbation_ == 100)
2831 perturbation_ = 50; // treat as normal
2832 int savePerturbation = perturbation_;
2833 int i;
2834 copyFromSaved(14); // copy bounds and costs
2835 // look at element range
2836 double smallestNegative;
2837 double largestNegative;
2838 double smallestPositive;
2839 double largestPositive;
2840 matrix_->rangeOfElements(smallestNegative, largestNegative,
2841 smallestPositive, largestPositive);
2842 smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
2843 largestPositive = CoinMax(fabs(largestNegative), largestPositive);
2844 double elementRatio = largestPositive / smallestPositive;
2845 if ((!numberIterations_ || initialSumInfeasibilities_ == 1.23456789e-5) && perturbation_ == 50) {
2846 // See if we need to perturb
2847 int numberTotal = CoinMax(numberRows_, numberColumns_);
2848 double *sort = new double[numberTotal];
2849 int nFixed = 0;
2850 for (i = 0; i < numberRows_; i++) {
2851 double lo = fabs(rowLower_[i]);
2852 double up = fabs(rowUpper_[i]);
2853 double value = 0.0;
2854 if (lo && lo < 1.0e20) {
2855 if (up && up < 1.0e20) {
2856 value = 0.5 * (lo + up);
2857 if (lo == up)
2858 nFixed++;
2859 } else {
2860 value = lo;
2861 }
2862 } else {
2863 if (up && up < 1.0e20)
2864 value = up;
2865 }
2866 sort[i] = value;
2867 }
2868 std::sort(sort, sort + numberRows_);
2869 int number = 1;
2870 double last = sort[0];
2871 for (i = 1; i < numberRows_; i++) {
2872 if (last != sort[i])
2873 number++;
2874 last = sort[i];
2875 }
2876 #ifdef KEEP_GOING_IF_FIXED
2877 //printf("ratio number diff rhs %g (%d %d fixed), element ratio %g\n",((double)number)/((double) numberRows_),
2878 // numberRows_,nFixed,elementRatio);
2879 #endif
2880 if (number * 4 > numberRows_ || elementRatio > 1.0e12) {
2881 perturbation_ = 100;
2882 delete[] sort;
2883 // Make sure feasible bounds
2884 if (abcNonLinearCost_) {
2885 abcNonLinearCost_->refresh();
2886 abcNonLinearCost_->checkInfeasibilities();
2887 //abcNonLinearCost_->feasibleBounds();
2888 }
2889 moveToBasic();
2890 return; // good enough
2891 }
2892 number = 0;
2893 #ifdef KEEP_GOING_IF_FIXED
2894 if (!integerType_) {
2895 // look at columns
2896 nFixed = 0;
2897 for (i = 0; i < numberColumns_; i++) {
2898 double lo = fabs(columnLower_[i]);
2899 double up = fabs(columnUpper_[i]);
2900 double value = 0.0;
2901 if (lo && lo < 1.0e20) {
2902 if (up && up < 1.0e20) {
2903 value = 0.5 * (lo + up);
2904 if (lo == up)
2905 nFixed++;
2906 } else {
2907 value = lo;
2908 }
2909 } else {
2910 if (up && up < 1.0e20)
2911 value = up;
2912 }
2913 sort[i] = value;
2914 }
2915 std::sort(sort, sort + numberColumns_);
2916 number = 1;
2917 last = sort[0];
2918 for (i = 1; i < numberColumns_; i++) {
2919 if (last != sort[i])
2920 number++;
2921 last = sort[i];
2922 }
2923 //printf("cratio number diff bounds %g (%d %d fixed)\n",((double)number)/((double) numberColumns_),
2924 // numberColumns_,nFixed);
2925 }
2926 #endif
2927 delete[] sort;
2928 if (number * 4 > numberColumns_) {
2929 perturbation_ = 100;
2930 // Make sure feasible bounds
2931 if (abcNonLinearCost_) {
2932 abcNonLinearCost_->refresh();
2933 abcNonLinearCost_->checkInfeasibilities();
2934 //abcNonLinearCost_->feasibleBounds();
2935 }
2936 moveToBasic();
2937 return; // good enough
2938 }
2939 }
2940 // primal perturbation
2941 double perturbation = 1.0e-20;
2942 double bias = 1.0;
2943 int numberNonZero = 0;
2944 // maximum fraction of rhs/bounds to perturb
2945 double maximumFraction = 1.0e-5;
2946 double overallMultiplier = (perturbation_ == 50 || perturbation_ > 54) ? 2.0 : 0.2;
2947 #ifdef HEAVY_PERTURBATION
2948 if (perturbation_ == 50)
2949 perturbation_ = HEAVY_PERTURBATION;
2950 #endif
2951 if (perturbation_ >= 50) {
2952 perturbation = 1.0e-4;
2953 for (i = 0; i < numberColumns_ + numberRows_; i++) {
2954 if (abcUpper_[i] > abcLower_[i] + primalTolerance_) {
2955 double lowerValue, upperValue;
2956 if (abcLower_[i] > -1.0e20)
2957 lowerValue = fabs(abcLower_[i]);
2958 else
2959 lowerValue = 0.0;
2960 if (abcUpper_[i] < 1.0e20)
2961 upperValue = fabs(abcUpper_[i]);
2962 else
2963 upperValue = 0.0;
2964 double value = CoinMax(fabs(lowerValue), fabs(upperValue));
2965 value = CoinMin(value, abcUpper_[i] - abcLower_[i]);
2966 #if 1
2967 if (value) {
2968 perturbation += value;
2969 numberNonZero++;
2970 }
2971 #else
2972 perturbation = CoinMax(perturbation, value);
2973 #endif
2974 }
2975 }
2976 if (numberNonZero)
2977 perturbation /= static_cast< double >(numberNonZero);
2978 else
2979 perturbation = 1.0e-1;
2980 if (perturbation_ > 50 && perturbation_ < 55) {
2981 // reduce
2982 while (perturbation_ < 55) {
2983 perturbation_++;
2984 perturbation *= 0.25;
2985 bias *= 0.25;
2986 }
2987 perturbation_ = 50;
2988 } else if (perturbation_ >= 55 && perturbation_ < 60) {
2989 // increase
2990 while (perturbation_ > 55) {
2991 overallMultiplier *= 1.2;
2992 perturbation_--;
2993 perturbation *= 4.0;
2994 }
2995 perturbation_ = 50;
2996 }
2997 } else if (perturbation_ < 100) {
2998 perturbation = pow(10.0, perturbation_);
2999 // user is in charge
3000 maximumFraction = 1.0;
3001 }
3002 double largestZero = 0.0;
3003 double largest = 0.0;
3004 double largestPerCent = 0.0;
3005 bool printOut = (handler_->logLevel() == 63);
3006 printOut = false; //off
3007 // Check if all slack
3008 int number = 0;
3009 int iSequence;
3010 for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3011 if (getInternalStatus(iSequence) == basic)
3012 number++;
3013 }
3014 if (rhsScale_ > 100.0) {
3015 // tone down perturbation
3016 maximumFraction *= 0.1;
3017 }
3018 if (savePerturbation == 51) {
3019 perturbation = CoinMin(0.1, perturbation);
3020 maximumFraction *= 0.1;
3021 }
3022 //if (number != numberRows_)
3023 //type = 1;
3024 // modify bounds
3025 // Change so at least 1.0e-5 and no more than 0.1
3026 // For now just no more than 0.1
3027 // printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
3028 // seems much slower???
3029 //const double * COIN_RESTRICT perturbationArray = perturbationSaved_;
3030 // Make sure feasible bounds
3031 if (abcNonLinearCost_ && true) {
3032 abcNonLinearCost_->refresh();
3033 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
3034 //abcNonLinearCost_->feasibleBounds();
3035 }
3036 double tolerance = 100.0 * primalTolerance_;
3037 int numberChanged = 0;
3038 // Set bit if fixed
3039 for (int i = 0; i < numberRows_; i++) {
3040 if (rowLower_[i] != rowUpper_[i])
3041 internalStatus_[i] &= ~128;
3042 else
3043 internalStatus_[i] |= 128;
3044 }
3045 for (int i = 0; i < numberColumns_; i++) {
3046 if (columnLower_[i] != columnUpper_[i])
3047 internalStatus_[i + numberRows_] &= ~128;
3048 else
3049 internalStatus_[i + numberRows_] |= 128;
3050 }
3051 //double multiplier = perturbation*maximumFraction;
3052 for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
3053 if (getInternalStatus(iSequence) == basic) {
3054 if ((internalStatus_[i] & 128) != 0)
3055 continue;
3056 double lowerValue = abcLower_[iSequence];
3057 double upperValue = abcUpper_[iSequence];
3058 if (upperValue > lowerValue + tolerance) {
3059 double solutionValue = abcSolution_[iSequence];
3060 double difference = upperValue - lowerValue;
3061 difference = CoinMin(difference, perturbation);
3062 difference = CoinMin(difference, fabs(solutionValue) + 1.0);
3063 double value = maximumFraction * (difference + bias);
3064 value = CoinMin(value, 0.1);
3065 value = CoinMax(value, primalTolerance_);
3066 double perturbationValue = overallMultiplier * randomNumberGenerator_.randomDouble();
3067 value *= perturbationValue;
3068 if (solutionValue - lowerValue <= primalTolerance_) {
3069 abcLower_[iSequence] -= value;
3070 // change correct saved value
3071 if (abcNonLinearCost_->getCurrentStatus(iSequence) == CLP_FEASIBLE)
3072 abcPerturbation_[iSequence] = abcLower_[iSequence];
3073 else
3074 abcPerturbation_[iSequence + numberTotal_] = abcLower_[iSequence];
3075 } else if (upperValue - solutionValue <= primalTolerance_) {
3076 abcUpper_[iSequence] += value;
3077 // change correct saved value
3078 if (abcNonLinearCost_->getCurrentStatus(iSequence) == CLP_FEASIBLE)
3079 abcPerturbation_[iSequence + numberTotal_] = abcUpper_[iSequence];
3080 else
3081 abcPerturbation_[iSequence] = abcUpper_[iSequence];
3082 } else {
3083 #if 0
3084 if (iSequence >= numberColumns_) {
3085 // may not be at bound - but still perturb (unless free)
3086 if (upperValue > 1.0e30 && lowerValue < -1.0e30)
3087 value = 0.0;
3088 else
3089 value = - value; // as -1.0 in matrix
3090 } else {
3091 value = 0.0;
3092 }
3093 #else
3094 value = 0.0;
3095 #endif
3096 }
3097 if (value) {
3098 numberChanged++;
3099 if (printOut)
3100 printf("col %d lower from %g to %g, upper from %g to %g\n",
3101 iSequence, abcLower_[iSequence], lowerValue, abcUpper_[iSequence], upperValue);
3102 if (solutionValue) {
3103 largest = CoinMax(largest, value);
3104 if (value > (fabs(solutionValue) + 1.0) * largestPerCent)
3105 largestPerCent = value / (fabs(solutionValue) + 1.0);
3106 } else {
3107 largestZero = CoinMax(largestZero, value);
3108 }
3109 }
3110 }
3111 }
3112 }
3113 if (!numberChanged) {
3114 // do non basic columns?
3115 for (iSequence = 0; iSequence < maximumAbcNumberRows_ + numberColumns_; iSequence++) {
3116 if (getInternalStatus(iSequence) != basic) {
3117 if ((internalStatus_[i] & 128) != 0)
3118 continue;
3119 double lowerValue = abcLower_[iSequence];
3120 double upperValue = abcUpper_[iSequence];
3121 if (upperValue > lowerValue + tolerance) {
3122 double solutionValue = abcSolution_[iSequence];
3123 double difference = upperValue - lowerValue;
3124 difference = CoinMin(difference, perturbation);
3125 difference = CoinMin(difference, fabs(solutionValue) + 1.0);
3126 double value = maximumFraction * (difference + bias);
3127 value = CoinMin(value, 0.1);
3128 value = CoinMax(value, primalTolerance_);
3129 double perturbationValue = overallMultiplier * randomNumberGenerator_.randomDouble();
3130 value *= perturbationValue;
3131 if (solutionValue - lowerValue <= primalTolerance_) {
3132 abcLower_[iSequence] -= value;
3133 // change correct saved value
3134 if (abcNonLinearCost_->getCurrentStatus(iSequence) == CLP_FEASIBLE)
3135 abcPerturbation_[iSequence] = abcLower_[iSequence];
3136 else
3137 abcPerturbation_[iSequence + numberTotal_] = abcLower_[iSequence];
3138 } else if (upperValue - solutionValue <= primalTolerance_) {
3139 abcUpper_[iSequence] += value;
3140 // change correct saved value
3141 if (abcNonLinearCost_->getCurrentStatus(iSequence) == CLP_FEASIBLE)
3142 abcPerturbation_[iSequence + numberTotal_] = abcUpper_[iSequence];
3143 else
3144 abcPerturbation_[iSequence] = abcUpper_[iSequence];
3145 } else {
3146 value = 0.0;
3147 }
3148 if (value) {
3149 if (printOut)
3150 printf("col %d lower from %g to %g, upper from %g to %g\n",
3151 iSequence, abcLower_[iSequence], lowerValue, abcUpper_[iSequence], upperValue);
3152 if (solutionValue) {
3153 largest = CoinMax(largest, value);
3154 if (value > (fabs(solutionValue) + 1.0) * largestPerCent)
3155 largestPerCent = value / (fabs(solutionValue) + 1.0);
3156 } else {
3157 largestZero = CoinMax(largestZero, value);
3158 }
3159 }
3160 }
3161 }
3162 }
3163 }
3164 // Clean up
3165 for (i = 0; i < numberColumns_ + numberRows_; i++) {
3166 internalStatus_[i] &= ~128;
3167 switch (getInternalStatus(i)) {
3168
3169 case basic:
3170 break;
3171 case atUpperBound:
3172 abcSolution_[i] = abcUpper_[i];
3173 break;
3174 case isFixed:
3175 case atLowerBound:
3176 abcSolution_[i] = abcLower_[i];
3177 break;
3178 case isFree:
3179 break;
3180 case superBasic:
3181 break;
3182 }
3183 }
3184 handler_->message(CLP_SIMPLEX_PERTURB, messages_)
3185 << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
3186 << CoinMessageEol;
3187 // redo nonlinear costs
3188 //delete abcNonLinearCost_;abcNonLinearCost_=new AbcNonLinearCost(this);//abort();// something elseabcNonLinearCost_->refresh();
3189 moveToBasic();
3190 if (!numberChanged) {
3191 // we changed nonbasic
3192 gutsOfPrimalSolution(3);
3193 // Make sure feasible bounds
3194 if (abcNonLinearCost_) {
3195 //abcNonLinearCost_->refresh();
3196 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
3197 //abcNonLinearCost_->feasibleBounds();
3198 }
3199 abcPrimalColumnPivot_->saveWeights(this, 3);
3200 }
3201 // say perturbed
3202 perturbation_ = 101;
3203 }
3204 // un perturb
unPerturb()3205 bool AbcSimplexPrimal::unPerturb()
3206 {
3207 if (perturbation_ != 101)
3208 return false;
3209 // put back original bounds and costs
3210 copyFromSaved();
3211 // copy bounds to perturbation
3212 CoinAbcMemcpy(abcPerturbation_, abcLower_, numberTotal_);
3213 CoinAbcMemcpy(abcPerturbation_ + numberTotal_, abcUpper_, numberTotal_);
3214 //sanityCheck();
3215 // unflag
3216 unflag();
3217 abcProgress_.clearTimesFlagged();
3218 // get a valid nonlinear cost function
3219 delete abcNonLinearCost_;
3220 abcNonLinearCost_ = new AbcNonLinearCost(this);
3221 perturbation_ = 102; // stop any further perturbation
3222 // move non basic variables to new bounds
3223 abcNonLinearCost_->checkInfeasibilities(0.0);
3224 gutsOfSolution(3);
3225 abcNonLinearCost_->checkInfeasibilities(primalTolerance_);
3226 // Try using dual
3227 return abcNonLinearCost_->sumInfeasibilities() != 0.0;
3228 }
3229 // Unflag all variables and return number unflagged
unflag()3230 int AbcSimplexPrimal::unflag()
3231 {
3232 int i;
3233 int number = numberRows_ + numberColumns_;
3234 int numberFlagged = 0;
3235 // we can't really trust infeasibilities if there is dual error
3236 // allow tolerance bigger than standard to check on duals
3237 double relaxedToleranceD = dualTolerance_ + CoinMin(1.0e-2, 10.0 * largestDualError_);
3238 for (i = 0; i < number; i++) {
3239 if (flagged(i)) {
3240 clearFlagged(i);
3241 // only say if reasonable dj
3242 if (fabs(abcDj_[i]) > relaxedToleranceD)
3243 numberFlagged++;
3244 }
3245 }
3246 #if ABC_NORMAL_DEBUG > 0
3247 if (handler_->logLevel() > 2 && numberFlagged && objective_->type() > 1)
3248 printf("%d unflagged\n", numberFlagged);
3249 #endif
3250 return numberFlagged;
3251 }
3252 // Do not change infeasibility cost and always say optimal
alwaysOptimal(bool onOff)3253 void AbcSimplexPrimal::alwaysOptimal(bool onOff)
3254 {
3255 if (onOff)
3256 specialOptions_ |= 1;
3257 else
3258 specialOptions_ &= ~1;
3259 }
alwaysOptimal() const3260 bool AbcSimplexPrimal::alwaysOptimal() const
3261 {
3262 return (specialOptions_ & 1) != 0;
3263 }
3264 // Flatten outgoing variables i.e. - always to exact bound
exactOutgoing(bool onOff)3265 void AbcSimplexPrimal::exactOutgoing(bool onOff)
3266 {
3267 if (onOff)
3268 specialOptions_ |= 4;
3269 else
3270 specialOptions_ &= ~4;
3271 }
exactOutgoing() const3272 bool AbcSimplexPrimal::exactOutgoing() const
3273 {
3274 return (specialOptions_ & 4) != 0;
3275 }
3276 /*
3277 Reasons to come out (normal mode/user mode):
3278 -1 normal
3279 -2 factorize now - good iteration/ NA
3280 -3 slight inaccuracy - refactorize - iteration done/ same but factor done
3281 -4 inaccuracy - refactorize - no iteration/ NA
3282 -5 something flagged - go round again/ pivot not possible
3283 +2 looks unbounded
3284 +3 max iterations (iteration done)
3285 */
pivotResult(int ifValuesPass)3286 int AbcSimplexPrimal::pivotResult(int ifValuesPass)
3287 {
3288
3289 bool roundAgain = true;
3290 int returnCode = -1;
3291
3292 // loop round if user setting and doing refactorization
3293 while (roundAgain) {
3294 roundAgain = false;
3295 returnCode = -1;
3296 pivotRow_ = -1;
3297 sequenceOut_ = -1;
3298 usefulArray_[arrayForFtran_].clear();
3299
3300 // we found a pivot column
3301 // update the incoming column
3302 unpack(usefulArray_[arrayForFtran_]);
3303 // save reduced cost
3304 double saveDj = dualIn_;
3305 abcFactorization_->updateColumnFT(usefulArray_[arrayForFtran_]);
3306 // do ratio test and re-compute dj
3307 #ifdef CLP_USER_DRIVEN
3308 if ((moreSpecialOptions_ & 512) == 0) {
3309 #endif
3310 primalRow(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForFlipBounds_],
3311 &usefulArray_[arrayForTableauRow_],
3312 ifValuesPass);
3313 #ifdef CLP_USER_DRIVEN
3314 // user can tell which use it is
3315 int status = eventHandler_->event(ClpEventHandler::pivotRow);
3316 if (status >= 0) {
3317 problemStatus_ = 5;
3318 secondaryStatus_ = ClpEventHandler::pivotRow;
3319 break;
3320 }
3321 } else {
3322 int status = eventHandler_->event(ClpEventHandler::pivotRow);
3323 if (status >= 0) {
3324 problemStatus_ = 5;
3325 secondaryStatus_ = ClpEventHandler::pivotRow;
3326 break;
3327 }
3328 }
3329 #endif
3330 if (ifValuesPass) {
3331 saveDj = dualIn_;
3332 //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
3333 if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3334 if (fabs(dualIn_) < 1.0e2 * dualTolerance_ && objective_->type() < 2) {
3335 // try other way
3336 directionIn_ = -directionIn_;
3337 primalRow(&usefulArray_[arrayForFtran_], &usefulArray_[arrayForFlipBounds_],
3338 &usefulArray_[arrayForTableauRow_],
3339 0);
3340 }
3341 if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3342 // reject it
3343 char x = isColumn(sequenceIn_) ? 'C' : 'R';
3344 handler_->message(CLP_SIMPLEX_FLAG, messages_)
3345 << x << sequenceWithin(sequenceIn_)
3346 << CoinMessageEol;
3347 setFlagged(sequenceIn_);
3348 abcProgress_.incrementTimesFlagged();
3349 abcProgress_.clearBadTimes();
3350 lastBadIteration_ = numberIterations_; // say be more cautious
3351 clearAll();
3352 pivotRow_ = -1;
3353 returnCode = -5;
3354 break;
3355 }
3356 }
3357 }
3358 double checkValue = 1.0e-2;
3359 if (largestDualError_ > 1.0e-5)
3360 checkValue = 1.0e-1;
3361 double test2 = dualTolerance_;
3362 double test1 = 1.0e-20;
3363 #if 0 //def FEB_TRY
3364 if (abcFactorization_->pivots() < 1) {
3365 test1 = -1.0e-4;
3366 if ((saveDj < 0.0 && dualIn_ < -1.0e-5 * dualTolerance_) ||
3367 (saveDj > 0.0 && dualIn_ > 1.0e-5 * dualTolerance_))
3368 test2 = 0.0; // allow through
3369 }
3370 #endif
3371 if (!ifValuesPass && (saveDj * dualIn_ < test1 || fabs(saveDj - dualIn_) > checkValue * (1.0 + fabs(saveDj)) || fabs(dualIn_) < test2)) {
3372 if (!(saveDj * dualIn_ > 0.0 && CoinMin(fabs(saveDj), fabs(dualIn_)) > 1.0e5)) {
3373 char x = isColumn(sequenceIn_) ? 'C' : 'R';
3374 handler_->message(CLP_PRIMAL_DJ, messages_)
3375 << x << sequenceWithin(sequenceIn_)
3376 << saveDj << dualIn_
3377 << CoinMessageEol;
3378 if (lastGoodIteration_ != numberIterations_) {
3379 clearAll();
3380 pivotRow_ = -1; // say no weights update
3381 returnCode = -4;
3382 if (lastGoodIteration_ + 1 == numberIterations_) {
3383 // not looking wonderful - try cleaning bounds
3384 // put non-basics to bounds in case tolerance moved
3385 abcNonLinearCost_->checkInfeasibilities(0.0);
3386 }
3387 sequenceOut_ = -1;
3388 sequenceIn_ = -1;
3389 if (abcFactorization_->pivots() < 10 && abcFactorization_->pivotTolerance() < 0.25)
3390 abcFactorization_->saferTolerances(1.0, -1.03);
3391 break;
3392 } else {
3393 // take on more relaxed criterion
3394 if (saveDj * dualIn_ < test1 || fabs(saveDj - dualIn_) > 2.0e-1 * (1.0 + fabs(dualIn_)) || fabs(dualIn_) < test2) {
3395 // need to reject something
3396 char x = isColumn(sequenceIn_) ? 'C' : 'R';
3397 handler_->message(CLP_SIMPLEX_FLAG, messages_)
3398 << x << sequenceWithin(sequenceIn_)
3399 << CoinMessageEol;
3400 setFlagged(sequenceIn_);
3401 abcProgress_.incrementTimesFlagged();
3402 #if 1 //def FEB_TRY \
3403 // Make safer?
3404 double tolerance = abcFactorization_->pivotTolerance();
3405 abcFactorization_->saferTolerances(1.0, -1.03);
3406 #endif
3407 abcProgress_.clearBadTimes();
3408 lastBadIteration_ = numberIterations_; // say be more cautious
3409 clearAll();
3410 pivotRow_ = -1;
3411 returnCode = -5;
3412 if (tolerance < abcFactorization_->pivotTolerance())
3413 returnCode = -4;
3414 sequenceOut_ = -1;
3415 sequenceIn_ = -1;
3416 break;
3417 }
3418 }
3419 } else {
3420 //printf("%d %g %g\n",numberIterations_,saveDj,dualIn_);
3421 }
3422 }
3423 if (pivotRow_ >= 0) {
3424 #ifdef CLP_USER_DRIVEN1
3425 // Got good pivot - may need to unflag stuff
3426 userChoiceWasGood(this);
3427 #endif
3428 // if stable replace in basis
3429 // check update
3430 //abcFactorization_->checkReplacePart1a(&usefulArray_[arrayForReplaceColumn_],pivotRow_);
3431 //usefulArray_[arrayForReplaceColumn_].print();
3432 ftAlpha_ = abcFactorization_->checkReplacePart1(&usefulArray_[arrayForReplaceColumn_], pivotRow_);
3433 int updateStatus = abcFactorization_->checkReplacePart2(pivotRow_, btranAlpha_, alpha_, ftAlpha_);
3434 abcFactorization_->replaceColumnPart3(this,
3435 &usefulArray_[arrayForReplaceColumn_],
3436 &usefulArray_[arrayForFtran_],
3437 pivotRow_,
3438 ftAlpha_ ? ftAlpha_ : alpha_);
3439
3440 // if no pivots, bad update but reasonable alpha - take and invert
3441 if (updateStatus == 2 && lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
3442 updateStatus = 4;
3443 if (updateStatus == 1 || updateStatus == 4) {
3444 // slight error
3445 if (abcFactorization_->pivots() > 5 || updateStatus == 4) {
3446 returnCode = -3;
3447 }
3448 } else if (updateStatus == 2) {
3449 // major error
3450 // better to have small tolerance even if slower
3451 abcFactorization_->zeroTolerance(CoinMin(abcFactorization_->zeroTolerance(), 1.0e-15));
3452 int maxFactor = abcFactorization_->maximumPivots();
3453 if (maxFactor > 10) {
3454 if (forceFactorization_ < 0)
3455 forceFactorization_ = maxFactor;
3456 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
3457 }
3458 // later we may need to unwind more e.g. fake bounds
3459 if (lastGoodIteration_ != numberIterations_) {
3460 clearAll();
3461 pivotRow_ = -1;
3462 sequenceIn_ = -1;
3463 sequenceOut_ = -1;
3464 returnCode = -4;
3465 break;
3466 } else {
3467 // need to reject something
3468 char x = isColumn(sequenceIn_) ? 'C' : 'R';
3469 handler_->message(CLP_SIMPLEX_FLAG, messages_)
3470 << x << sequenceWithin(sequenceIn_)
3471 << CoinMessageEol;
3472 setFlagged(sequenceIn_);
3473 abcProgress_.incrementTimesFlagged();
3474 abcProgress_.clearBadTimes();
3475 lastBadIteration_ = numberIterations_; // say be more cautious
3476 clearAll();
3477 pivotRow_ = -1;
3478 sequenceIn_ = -1;
3479 sequenceOut_ = -1;
3480 returnCode = -5;
3481 break;
3482 }
3483 } else if (updateStatus == 3) {
3484 // out of memory
3485 // increase space if not many iterations
3486 if (abcFactorization_->pivots() < 0.5 * abcFactorization_->maximumPivots() && abcFactorization_->pivots() < 200)
3487 abcFactorization_->areaFactor(
3488 abcFactorization_->areaFactor() * 1.1);
3489 returnCode = -2; // factorize now
3490 } else if (updateStatus == 5) {
3491 problemStatus_ = -2; // factorize now
3492 }
3493 // here do part of steepest - ready for next iteration
3494 if (!ifValuesPass)
3495 abcPrimalColumnPivot_->updateWeights(&usefulArray_[arrayForFtran_]);
3496 } else {
3497 if (pivotRow_ == -1) {
3498 // no outgoing row is valid
3499 if (valueOut_ != COIN_DBL_MAX) {
3500 double objectiveChange = 0.0;
3501 theta_ = valueOut_ - valueIn_;
3502 updatePrimalsInPrimal(&usefulArray_[arrayForFtran_], theta_, objectiveChange, ifValuesPass);
3503 abcSolution_[sequenceIn_] += theta_;
3504 }
3505 #ifdef CLP_USER_DRIVEN1
3506 /* Note if valueOut_ < COIN_DBL_MAX and
3507 theta_ reasonable then this may be a valid sub flip */
3508 if (!userChoiceValid2(this)) {
3509 if (abcFactorization_->pivots() < 5) {
3510 // flag variable
3511 char x = isColumn(sequenceIn_) ? 'C' : 'R';
3512 handler_->message(CLP_SIMPLEX_FLAG, messages_)
3513 << x << sequenceWithin(sequenceIn_)
3514 << CoinMessageEol;
3515 setFlagged(sequenceIn_);
3516 abcProgress_.incrementTimesFlagged();
3517 abcProgress_.clearBadTimes();
3518 roundAgain = true;
3519 continue;
3520 } else {
3521 // try refactorizing first
3522 returnCode = 4; //say looks odd but has iterated
3523 break;
3524 }
3525 }
3526 #endif
3527 if (!abcFactorization_->pivots() && acceptablePivot_ <= 1.0e-8) {
3528 returnCode = 2; //say looks unbounded
3529 // do ray
3530 primalRay(&usefulArray_[arrayForFtran_]);
3531 } else {
3532 acceptablePivot_ = 1.0e-8;
3533 returnCode = 4; //say looks unbounded but has iterated
3534 }
3535 break;
3536 } else {
3537 // flipping from bound to bound
3538 }
3539 }
3540
3541 double oldCost = 0.0;
3542 if (sequenceOut_ >= 0)
3543 oldCost = abcCost_[sequenceOut_];
3544 // update primal solution
3545
3546 double objectiveChange = 0.0;
3547 // after this usefulArray_[arrayForFtran_] is not empty - used to update djs
3548 #ifdef CLP_USER_DRIVEN
3549 if (theta_ < 0.0) {
3550 if (theta_ >= -1.0e-12)
3551 theta_ = 0.0;
3552 //else
3553 //printf("negative theta %g\n",theta_);
3554 }
3555 #endif
3556 updatePrimalsInPrimal(&usefulArray_[arrayForFtran_], theta_, objectiveChange, ifValuesPass);
3557
3558 double oldValue = valueIn_;
3559 if (directionIn_ == -1) {
3560 // as if from upper bound
3561 if (sequenceIn_ != sequenceOut_) {
3562 // variable becoming basic
3563 valueIn_ -= fabs(theta_);
3564 } else {
3565 valueIn_ = lowerIn_;
3566 }
3567 } else {
3568 // as if from lower bound
3569 if (sequenceIn_ != sequenceOut_) {
3570 // variable becoming basic
3571 valueIn_ += fabs(theta_);
3572 } else {
3573 valueIn_ = upperIn_;
3574 }
3575 }
3576 objectiveChange += dualIn_ * (valueIn_ - oldValue);
3577 // outgoing
3578 if (sequenceIn_ != sequenceOut_) {
3579 if (directionOut_ > 0) {
3580 valueOut_ = lowerOut_;
3581 } else {
3582 valueOut_ = upperOut_;
3583 }
3584 if (valueOut_ < abcLower_[sequenceOut_] - primalTolerance_)
3585 valueOut_ = abcLower_[sequenceOut_] - 0.9 * primalTolerance_;
3586 else if (valueOut_ > abcUpper_[sequenceOut_] + primalTolerance_)
3587 valueOut_ = abcUpper_[sequenceOut_] + 0.9 * primalTolerance_;
3588 // may not be exactly at bound and bounds may have changed
3589 // Make sure outgoing looks feasible
3590 directionOut_ = abcNonLinearCost_->setOneOutgoing(pivotRow_, valueOut_);
3591 // May have got inaccurate
3592 //if (oldCost!=abcCost_[sequenceOut_])
3593 //printf("costchange on %d from %g to %g\n",sequenceOut_,
3594 // oldCost,abcCost_[sequenceOut_]);
3595 abcDj_[sequenceOut_] = abcCost_[sequenceOut_] - oldCost; // normally updated next iteration
3596 abcSolution_[sequenceOut_] = valueOut_;
3597 }
3598 // change cost and bounds on incoming if primal
3599 abcNonLinearCost_->setOne(sequenceIn_, valueIn_);
3600 int whatNext = housekeeping(/*objectiveChange*/);
3601 swapPrimalStuff();
3602 if (whatNext == 1) {
3603 returnCode = -2; // refactorize
3604 } else if (whatNext == 2) {
3605 // maximum iterations or equivalent
3606 returnCode = 3;
3607 } else if (numberIterations_ == lastGoodIteration_ + 2 * abcFactorization_->maximumPivots()) {
3608 // done a lot of flips - be safe
3609 returnCode = -2; // refactorize
3610 }
3611 // Check event
3612 {
3613 int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3614 if (status >= 0) {
3615 problemStatus_ = 5;
3616 secondaryStatus_ = ClpEventHandler::endOfIteration;
3617 returnCode = 3;
3618 }
3619 }
3620 }
3621 return returnCode;
3622 }
3623 /*
3624 Reasons to come out (normal mode/user mode):
3625 -1 normal
3626 -2 factorize now - good iteration/ NA
3627 -3 slight inaccuracy - refactorize - iteration done/ same but factor done
3628 -4 inaccuracy - refactorize - no iteration/ NA
3629 -5 something flagged - go round again/ pivot not possible
3630 +2 looks unbounded
3631 +3 max iterations (iteration done)
3632 */
pivotResult4(int ifValuesPass)3633 int AbcSimplexPrimal::pivotResult4(int ifValuesPass)
3634 {
3635
3636 int returnCode = -1;
3637 int numberMinor = 0;
3638 int numberDone = 0;
3639 CoinIndexedVector *vector[16];
3640 pivotStruct stuff[4];
3641 vector[0] = &usefulArray_[arrayForFtran_];
3642 vector[1] = &usefulArray_[arrayForFlipBounds_];
3643 vector[2] = &usefulArray_[arrayForTableauRow_];
3644 vector[3] = &usefulArray_[arrayForBtran_];
3645 /*
3646 For pricing we need to get a vector with difference in costs
3647 later we could modify code so only costBasic changed
3648 and store in first [4*x+1] array to go out the indices of changed
3649 plus at most four for pivots
3650 */
3651 //double * saveCosts=CoinCopyOfArray(costBasic_,numberRows_);
3652 //double saveCostsIn[4];
3653 for (int iMinor = 0; iMinor < 4; iMinor++) {
3654 int sequenceIn = multipleSequenceIn_[iMinor];
3655 if (sequenceIn < 0)
3656 break;
3657 stuff[iMinor].valuesPass_ = ifValuesPass;
3658 stuff[iMinor].lowerIn_ = abcLower_[sequenceIn];
3659 stuff[iMinor].upperIn_ = abcUpper_[sequenceIn];
3660 stuff[iMinor].valueIn_ = abcSolution_[sequenceIn];
3661 stuff[iMinor].sequenceIn_ = sequenceIn;
3662 //saveCostsIn[iMinor]=abcCost_[sequenceIn];
3663 numberMinor++;
3664 if (iMinor) {
3665 vector[4 * iMinor] = rowArray_[2 * iMinor - 2];
3666 vector[4 * iMinor + 1] = rowArray_[2 * iMinor - 1];
3667 vector[4 * iMinor + 2] = columnArray_[2 * iMinor - 2];
3668 vector[4 * iMinor + 3] = columnArray_[2 * iMinor - 1];
3669 }
3670 for (int i = 0; i < 4; i++)
3671 vector[4 * iMinor + i]->checkClear();
3672 unpack(*vector[4 * iMinor], sequenceIn);
3673 }
3674 int numberLeft = numberMinor;
3675 // parallel (with cpu)
3676 for (int iMinor = 1; iMinor < numberMinor; iMinor++) {
3677 // update the incoming columns
3678 cilk_spawn abcFactorization_->updateColumnFT(*vector[4 * iMinor], *vector[4 * iMinor + 3], iMinor);
3679 }
3680 abcFactorization_->updateColumnFT(*vector[0], *vector[+3], 0);
3681 cilk_sync;
3682 for (int iMinor = 0; iMinor < numberMinor; iMinor++) {
3683 // find best (or first if values pass)
3684 int numberDo = 1;
3685 int jMinor = 0;
3686 while (jMinor < numberDo) {
3687 int sequenceIn = stuff[jMinor].sequenceIn_;
3688 double dj = abcDj_[sequenceIn];
3689 bool bad = false;
3690 if (!bad) {
3691 stuff[jMinor].dualIn_ = dj;
3692 stuff[jMinor].saveDualIn_ = dj;
3693 if (dj < 0.0)
3694 stuff[jMinor].directionIn_ = 1;
3695 else
3696 stuff[jMinor].directionIn_ = -1;
3697 jMinor++;
3698 } else {
3699 numberDo--;
3700 numberLeft--;
3701 // throw away
3702 for (int i = 0; i < 4; i++) {
3703 vector[4 * jMinor + i]->clear();
3704 vector[4 * jMinor + i] = vector[4 * numberDo + i];
3705 vector[4 * numberDo + i] = NULL;
3706 }
3707 stuff[jMinor] = stuff[numberDo];
3708 }
3709 }
3710 for (int jMinor = 0; jMinor < numberDo; jMinor++) {
3711 // do ratio test and re-compute dj
3712 primalRow(vector[4 * jMinor], vector[4 * jMinor + 1], vector[4 * jMinor + 2], stuff[jMinor]);
3713 }
3714 // choose best
3715 int iBest = -1;
3716 double bestMovement = -COIN_DBL_MAX;
3717 for (int jMinor = 0; jMinor < numberDo; jMinor++) {
3718 double movement = stuff[jMinor].theta_ * fabs(stuff[jMinor].dualIn_);
3719 if (movement > bestMovement) {
3720 bestMovement = movement;
3721 iBest = jMinor;
3722 }
3723 }
3724 #if 0 //ndef MULTIPLE_PRICE
3725 if (maximumIterations()!=100000)
3726 iBest=0;
3727 #endif
3728 if (iBest >= 0) {
3729 dualIn_ = stuff[iBest].dualIn_;
3730 dualOut_ = stuff[iBest].dualOut_;
3731 lowerOut_ = stuff[iBest].lowerOut_;
3732 upperOut_ = stuff[iBest].upperOut_;
3733 valueOut_ = stuff[iBest].valueOut_;
3734 sequenceOut_ = stuff[iBest].sequenceOut_;
3735 sequenceIn_ = stuff[iBest].sequenceIn_;
3736 lowerIn_ = stuff[iBest].lowerIn_;
3737 upperIn_ = stuff[iBest].upperIn_;
3738 valueIn_ = stuff[iBest].valueIn_;
3739 directionIn_ = stuff[iBest].directionIn_;
3740 directionOut_ = stuff[iBest].directionOut_;
3741 pivotRow_ = stuff[iBest].pivotRow_;
3742 theta_ = stuff[iBest].theta_;
3743 alpha_ = stuff[iBest].alpha_;
3744 #ifdef MULTIPLE_PRICE
3745 for (int i = 0; i < 4 * numberLeft; i++)
3746 vector[i]->clear();
3747 return 0;
3748 #endif
3749 // maybe do this on more?
3750 double theta1 = CoinMax(theta_, 1.0e-12);
3751 double theta2 = numberIterations_ * abcNonLinearCost_->averageTheta();
3752 // Set average theta
3753 abcNonLinearCost_->setAverageTheta((theta1 + theta2) / (static_cast< double >(numberIterations_ + 1)));
3754 if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3755 if (fabs(dualIn_) < 1.0e2 * dualTolerance_) {
3756 // try other way
3757 stuff[iBest].directionIn_ = -directionIn_;
3758 stuff[iBest].valuesPass_ = 0;
3759 primalRow(vector[4 * iBest], vector[4 * iBest + 1], vector[4 * iBest + 2], stuff[iBest]);
3760 dualIn_ = stuff[iBest].dualIn_;
3761 dualOut_ = stuff[iBest].dualOut_;
3762 lowerOut_ = stuff[iBest].lowerOut_;
3763 upperOut_ = stuff[iBest].upperOut_;
3764 valueOut_ = stuff[iBest].valueOut_;
3765 sequenceOut_ = stuff[iBest].sequenceOut_;
3766 sequenceIn_ = stuff[iBest].sequenceIn_;
3767 lowerIn_ = stuff[iBest].lowerIn_;
3768 upperIn_ = stuff[iBest].upperIn_;
3769 valueIn_ = stuff[iBest].valueIn_;
3770 directionIn_ = stuff[iBest].directionIn_;
3771 directionOut_ = stuff[iBest].directionOut_;
3772 pivotRow_ = stuff[iBest].pivotRow_;
3773 theta_ = stuff[iBest].theta_;
3774 alpha_ = stuff[iBest].alpha_;
3775 }
3776 if (pivotRow_ == -1 || (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
3777 // reject it
3778 char x = isColumn(sequenceIn_) ? 'C' : 'R';
3779 handler_->message(CLP_SIMPLEX_FLAG, messages_)
3780 << x << sequenceWithin(sequenceIn_)
3781 << CoinMessageEol;
3782 setFlagged(sequenceIn_);
3783 abcProgress_.incrementTimesFlagged();
3784 abcProgress_.clearBadTimes();
3785 lastBadIteration_ = numberIterations_; // say be more cautious
3786 clearAll();
3787 pivotRow_ = -1;
3788 returnCode = -5;
3789 break;
3790 }
3791 }
3792 CoinIndexedVector *bestUpdate = NULL;
3793 if (pivotRow_ >= 0) {
3794 // Normal pivot
3795 // if stable replace in basis
3796 // check update
3797 CoinIndexedVector *tempVector[4];
3798 for (int i = 0; i < 4; i++)
3799 tempVector[i] = vector[4 * iBest + i];
3800 returnCode = cilk_spawn doFTUpdate(tempVector);
3801 bestUpdate = tempVector[0];
3802 // after this bestUpdate is not empty - used to update djs ??
3803 cilk_spawn updatePrimalsInPrimal(*bestUpdate, theta_, ifValuesPass);
3804 numberDone++;
3805 numberLeft--;
3806 // throw away
3807 for (int i = 0; i < 4; i++) {
3808 vector[4 * iBest + i] = vector[4 * numberLeft + i];
3809 vector[4 * numberLeft + i] = NULL;
3810 }
3811 stuff[iBest] = stuff[numberLeft];
3812 // update pi and other vectors and FT stuff
3813 // parallel (can go 8 way?)
3814 for (int jMinor = 0; jMinor < numberLeft; jMinor++) {
3815 cilk_spawn updatePartialUpdate(*vector[4 * jMinor + 3]);
3816 }
3817 // parallel
3818 for (int jMinor = 0; jMinor < numberLeft; jMinor++) {
3819 stuff[jMinor].dualIn_ = cilk_spawn updateMinorCandidate(*bestUpdate, *vector[4 * jMinor], stuff[jMinor].sequenceIn_);
3820 }
3821 cilk_sync;
3822 // throw away
3823 for (int i = 1; i < 4; i++)
3824 tempVector[i]->clear();
3825 if (returnCode < -3) {
3826 clearAll();
3827 break;
3828 }
3829 // end Normal pivot
3830 } else {
3831 // start Flip
3832 vector[4 * iBest + 3]->clear();
3833 if (pivotRow_ == -1) {
3834 // no outgoing row is valid
3835 if (valueOut_ != COIN_DBL_MAX) {
3836 theta_ = valueOut_ - valueIn_;
3837 updatePrimalsInPrimal(*vector[4 * iBest], theta_, ifValuesPass);
3838 abcSolution_[sequenceIn_] += theta_;
3839 }
3840 if (!abcFactorization_->pivots() && acceptablePivot_ <= 1.0e-8) {
3841 returnCode = 2; //say looks unbounded
3842 // do ray
3843 primalRay(vector[4 * iBest]);
3844 } else {
3845 acceptablePivot_ = 1.0e-8;
3846 returnCode = 4; //say looks unbounded but has iterated
3847 }
3848 break;
3849 } else {
3850 // flipping from bound to bound
3851 bestUpdate = vector[4 * iBest];
3852 // after this bestUpdate is not empty - used to update djs ??
3853 updatePrimalsInPrimal(*bestUpdate, theta_, ifValuesPass);
3854 // throw away best BUT remember we are updating costs somehow
3855 numberDone++;
3856 numberLeft--;
3857 // throw away
3858 for (int i = 0; i < 4; i++) {
3859 if (i)
3860 vector[4 * iBest + i]->clear();
3861 vector[4 * iBest + i] = vector[4 * numberLeft + i];
3862 vector[4 * numberLeft + i] = NULL;
3863 }
3864 stuff[iBest] = stuff[numberLeft];
3865 // parallel
3866 for (int jMinor = 0; jMinor < numberLeft; jMinor++) {
3867 stuff[jMinor].dualIn_ = cilk_spawn updateMinorCandidate(*bestUpdate, *vector[4 * jMinor], stuff[jMinor].sequenceIn_);
3868 }
3869 cilk_sync;
3870 }
3871 // end Flip
3872 }
3873 bestUpdate->clear(); // as only works in values pass
3874
3875 if (directionIn_ == -1) {
3876 // as if from upper bound
3877 if (sequenceIn_ != sequenceOut_) {
3878 // variable becoming basic
3879 valueIn_ -= fabs(theta_);
3880 } else {
3881 valueIn_ = lowerIn_;
3882 }
3883 } else {
3884 // as if from lower bound
3885 if (sequenceIn_ != sequenceOut_) {
3886 // variable becoming basic
3887 valueIn_ += fabs(theta_);
3888 } else {
3889 valueIn_ = upperIn_;
3890 }
3891 }
3892 // outgoing
3893 if (sequenceIn_ != sequenceOut_) {
3894 if (directionOut_ > 0) {
3895 valueOut_ = lowerOut_;
3896 } else {
3897 valueOut_ = upperOut_;
3898 }
3899 if (valueOut_ < abcLower_[sequenceOut_] - primalTolerance_)
3900 valueOut_ = abcLower_[sequenceOut_] - 0.9 * primalTolerance_;
3901 else if (valueOut_ > abcUpper_[sequenceOut_] + primalTolerance_)
3902 valueOut_ = abcUpper_[sequenceOut_] + 0.9 * primalTolerance_;
3903 // may not be exactly at bound and bounds may have changed
3904 // Make sure outgoing looks feasible
3905 directionOut_ = abcNonLinearCost_->setOneOutgoing(pivotRow_, valueOut_);
3906 abcSolution_[sequenceOut_] = valueOut_;
3907 }
3908 // change cost and bounds on incoming if primal
3909 abcNonLinearCost_->setOne(sequenceIn_, valueIn_);
3910 int whatNext = housekeeping(/*objectiveChange*/);
3911 swapPrimalStuff();
3912 if (whatNext == 1) {
3913 returnCode = -2; // refactorize
3914 } else if (whatNext == 2) {
3915 // maximum iterations or equivalent
3916 returnCode = 3;
3917 } else if (numberIterations_ == lastGoodIteration_ + 2 * abcFactorization_->maximumPivots()) {
3918 // done a lot of flips - be safe
3919 returnCode = -2; // refactorize
3920 }
3921 // Check event
3922 {
3923 int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3924 if (status >= 0) {
3925 problemStatus_ = 5;
3926 secondaryStatus_ = ClpEventHandler::endOfIteration;
3927 returnCode = 3;
3928 }
3929 }
3930 }
3931 if (returnCode != -1)
3932 break;
3933 }
3934 for (int i = 0; i < 16; i++) {
3935 if (vector[i])
3936 vector[i]->clear();
3937 if (vector[i])
3938 vector[i]->checkClear();
3939 }
3940 //delete [] saveCosts;
3941 return returnCode;
3942 }
3943 // Do FT update as separate function for minor iterations (nonzero return code on problems)
doFTUpdate(CoinIndexedVector * vector[4])3944 int AbcSimplexPrimal::doFTUpdate(CoinIndexedVector *vector[4])
3945 {
3946 ftAlpha_ = abcFactorization_->checkReplacePart1(&usefulArray_[arrayForReplaceColumn_],
3947 vector[3], pivotRow_);
3948 int updateStatus = abcFactorization_->checkReplacePart2(pivotRow_, btranAlpha_, alpha_, ftAlpha_);
3949 if (!updateStatus)
3950 abcFactorization_->replaceColumnPart3(this,
3951 &usefulArray_[arrayForReplaceColumn_],
3952 vector[0],
3953 vector[3],
3954 pivotRow_,
3955 ftAlpha_ ? ftAlpha_ : alpha_);
3956 else
3957 vector[3]->clear();
3958 int returnCode = 0;
3959 // if no pivots, bad update but reasonable alpha - take and invert
3960 if (updateStatus == 2 && lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
3961 updateStatus = 4;
3962 if (updateStatus == 1 || updateStatus == 4) {
3963 // slight error
3964 if (abcFactorization_->pivots() > 5 || updateStatus == 4) {
3965 returnCode = -3;
3966 }
3967 } else if (updateStatus == 2) {
3968 // major error
3969 // better to have small tolerance even if slower
3970 abcFactorization_->zeroTolerance(CoinMin(abcFactorization_->zeroTolerance(), 1.0e-15));
3971 int maxFactor = abcFactorization_->maximumPivots();
3972 if (maxFactor > 10) {
3973 if (forceFactorization_ < 0)
3974 forceFactorization_ = maxFactor;
3975 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
3976 }
3977 // later we may need to unwind more e.g. fake bounds
3978 if (lastGoodIteration_ != numberIterations_) {
3979 pivotRow_ = -1;
3980 sequenceIn_ = -1;
3981 sequenceOut_ = -1;
3982 returnCode = -4;
3983 } else {
3984 // need to reject something
3985 char x = isColumn(sequenceIn_) ? 'C' : 'R';
3986 handler_->message(CLP_SIMPLEX_FLAG, messages_)
3987 << x << sequenceWithin(sequenceIn_)
3988 << CoinMessageEol;
3989 setFlagged(sequenceIn_);
3990 abcProgress_.incrementTimesFlagged();
3991 abcProgress_.clearBadTimes();
3992 lastBadIteration_ = numberIterations_; // say be more cautious
3993 pivotRow_ = -1;
3994 sequenceIn_ = -1;
3995 sequenceOut_ = -1;
3996 returnCode = -5;
3997 }
3998 } else if (updateStatus == 3) {
3999 // out of memory
4000 // increase space if not many iterations
4001 if (abcFactorization_->pivots() < 0.5 * abcFactorization_->maximumPivots() && abcFactorization_->pivots() < 200)
4002 abcFactorization_->areaFactor(
4003 abcFactorization_->areaFactor() * 1.1);
4004 returnCode = -2; // factorize now
4005 } else if (updateStatus == 5) {
4006 problemStatus_ = -2; // factorize now
4007 returnCode = -2;
4008 }
4009 return returnCode;
4010 }
4011 /* The primals are updated by the given array.
4012 costs are changed
4013 */
updatePrimalsInPrimal(CoinIndexedVector & rowArray,double theta,bool valuesPass)4014 void AbcSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector &rowArray,
4015 double theta, bool valuesPass)
4016 {
4017 // Cost on pivot row may change - may need to change dualIn
4018 //double oldCost = 0.0;
4019 //if (pivotRow_ >= 0)
4020 //oldCost = abcCost_[sequenceOut_];
4021 double *work = rowArray.denseVector();
4022 int number = rowArray.getNumElements();
4023 int *which = rowArray.getIndices();
4024 // allow for case where bound+tolerance == bound
4025 //double tolerance = 0.999999*primalTolerance_;
4026 double relaxedTolerance = 1.001 * primalTolerance_;
4027 int iIndex;
4028 if (!valuesPass) {
4029 for (iIndex = 0; iIndex < number; iIndex++) {
4030
4031 int iRow = which[iIndex];
4032 double alpha = work[iRow];
4033 //work[iRow] = 0.0;
4034 int iPivot = abcPivotVariable_[iRow];
4035 double change = theta * alpha;
4036 double value = abcSolution_[iPivot] - change;
4037 abcSolution_[iPivot] = value;
4038 value = solutionBasic_[iRow] - change;
4039 solutionBasic_[iRow] = value;
4040 if (active(iRow) || theta_ < 0.0) {
4041 clearActive(iRow);
4042 // But make sure one going out is feasible
4043 if (change > 0.0) {
4044 // going down
4045 if (value <= abcLower_[iPivot] + primalTolerance_) {
4046 if (iPivot == sequenceOut_ && value >= abcLower_[iPivot] - relaxedTolerance)
4047 value = abcLower_[iPivot];
4048 abcNonLinearCost_->setOneBasic(iRow, value);
4049 }
4050 } else {
4051 // going up
4052 if (value >= abcUpper_[iPivot] - primalTolerance_) {
4053 if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance)
4054 value = abcUpper_[iPivot];
4055 abcNonLinearCost_->setOneBasic(iRow, value);
4056 }
4057 }
4058 }
4059 }
4060 } else {
4061 // values pass so look at all
4062 for (iIndex = 0; iIndex < number; iIndex++) {
4063
4064 int iRow = which[iIndex];
4065 double alpha = work[iRow];
4066 //work[iRow] = 0.0;
4067 int iPivot = abcPivotVariable_[iRow];
4068 double change = theta * alpha;
4069 double value = abcSolution_[iPivot] - change;
4070 abcSolution_[iPivot] = value;
4071 value = solutionBasic_[iRow] - change;
4072 solutionBasic_[iRow] = value;
4073 clearActive(iRow);
4074 // But make sure one going out is feasible
4075 if (change > 0.0) {
4076 // going down
4077 if (value <= abcLower_[iPivot] + primalTolerance_) {
4078 if (iPivot == sequenceOut_ && value > abcLower_[iPivot] - relaxedTolerance)
4079 value = abcLower_[iPivot];
4080 abcNonLinearCost_->setOneBasic(iRow, value);
4081 }
4082 } else {
4083 // going up
4084 if (value >= abcUpper_[iPivot] - primalTolerance_) {
4085 if (iPivot == sequenceOut_ && value < abcUpper_[iPivot] + relaxedTolerance)
4086 value = abcUpper_[iPivot];
4087 abcNonLinearCost_->setOneBasic(iRow, value);
4088 }
4089 }
4090 }
4091 }
4092 //rowArray.setNumElements(0);
4093 }
4094 /* After rowArray will have cost changes for use next major iteration
4095 */
createUpdateDuals(CoinIndexedVector & rowArray,const double * originalCost,const double extraCost[4],double &,int)4096 void AbcSimplexPrimal::createUpdateDuals(CoinIndexedVector &rowArray,
4097 const double *originalCost,
4098 const double extraCost[4],
4099 double & /*objectiveChange*/,
4100 int /*valuesPass*/)
4101 {
4102 int number = 0;
4103 double *work = rowArray.denseVector();
4104 int *which = rowArray.getIndices();
4105 for (int iRow = 0; iRow < numberRows_; iRow++) {
4106 if (originalCost[iRow] != costBasic_[iRow]) {
4107 work[iRow] = costBasic_[iRow] - originalCost[iRow];
4108 which[number++] = iRow;
4109 }
4110 }
4111 rowArray.setNumElements(number);
4112 }
4113 // Update minor candidate vector
4114 double
updateMinorCandidate(const CoinIndexedVector & updateBy,CoinIndexedVector & candidate,int sequenceIn)4115 AbcSimplexPrimal::updateMinorCandidate(const CoinIndexedVector &updateBy,
4116 CoinIndexedVector &candidate,
4117 int sequenceIn)
4118 {
4119 const double *COIN_RESTRICT regionUpdate = updateBy.denseVector();
4120 const int *COIN_RESTRICT regionIndexUpdate = updateBy.getIndices();
4121 int numberNonZeroUpdate = updateBy.getNumElements();
4122 double *COIN_RESTRICT regionCandidate = candidate.denseVector();
4123 int *COIN_RESTRICT regionIndexCandidate = candidate.getIndices();
4124 int numberNonZeroCandidate = candidate.getNumElements();
4125 assert(fabs(alpha_ - regionUpdate[pivotRow_]) < 1.0e-5);
4126 double value = regionCandidate[pivotRow_];
4127 if (value) {
4128 value /= alpha_;
4129 for (int i = 0; i < numberNonZeroUpdate; i++) {
4130 int iRow = regionIndexUpdate[i];
4131 double oldValue = regionCandidate[iRow];
4132 double newValue = oldValue - value * regionUpdate[iRow];
4133 if (oldValue) {
4134 if (!newValue)
4135 newValue = COIN_INDEXED_REALLY_TINY_ELEMENT;
4136 regionCandidate[iRow] = newValue;
4137 } else {
4138 assert(newValue);
4139 regionCandidate[iRow] = newValue;
4140 regionIndexCandidate[numberNonZeroCandidate++] = iRow;
4141 }
4142 }
4143 double oldValue = regionCandidate[pivotRow_];
4144 double newValue = value;
4145 if (oldValue) {
4146 if (!newValue)
4147 newValue = COIN_INDEXED_REALLY_TINY_ELEMENT;
4148 regionCandidate[pivotRow_] = newValue;
4149 } else {
4150 assert(newValue);
4151 regionCandidate[pivotRow_] = newValue;
4152 regionIndexCandidate[numberNonZeroCandidate++] = pivotRow_;
4153 }
4154 candidate.setNumElements(numberNonZeroCandidate);
4155 }
4156 double newDj = abcCost_[sequenceIn];
4157 for (int i = 0; i < numberNonZeroCandidate; i++) {
4158 int iRow = regionIndexCandidate[i];
4159 double alpha = regionCandidate[iRow];
4160 newDj -= alpha * costBasic_[iRow];
4161 }
4162 return newDj;
4163 }
4164 // Update partial Ftran by R update
updatePartialUpdate(CoinIndexedVector & partialUpdate)4165 void AbcSimplexPrimal::updatePartialUpdate(CoinIndexedVector &partialUpdate)
4166 {
4167 #ifndef ABC_USE_COIN_FACTORIZATION
4168 CoinAbcFactorization *factorization = dynamic_cast< CoinAbcFactorization * >(abcFactorization_->factorization());
4169 if (factorization)
4170 factorization->updatePartialUpdate(partialUpdate);
4171 #else
4172 abcFactorization_->factorization()->updatePartialUpdate(partialUpdate);
4173 #endif
4174 }
4175 // Create primal ray
primalRay(CoinIndexedVector * rowArray)4176 void AbcSimplexPrimal::primalRay(CoinIndexedVector *rowArray)
4177 {
4178 delete[] ray_;
4179 ray_ = new double[numberColumns_];
4180 CoinZeroN(ray_, numberColumns_);
4181 int number = rowArray->getNumElements();
4182 int *index = rowArray->getIndices();
4183 double *array = rowArray->denseVector();
4184 double way = -directionIn_;
4185 int i;
4186 double zeroTolerance = 1.0e-12;
4187 if (sequenceIn_ < numberColumns_)
4188 ray_[sequenceIn_] = directionIn_;
4189 for (i = 0; i < number; i++) {
4190 int iRow = index[i];
4191 int iPivot = abcPivotVariable_[iRow];
4192 double arrayValue = array[iRow];
4193 if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
4194 ray_[iPivot] = way * arrayValue;
4195 }
4196 }
4197 /* Get next superbasic -1 if none,
4198 Normal type is 1
4199 If type is 3 then initializes sorted list
4200 if 2 uses list.
4201 */
nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray)4202 int AbcSimplexPrimal::nextSuperBasic(int superBasicType,
4203 CoinIndexedVector *columnArray)
4204 {
4205 int returnValue = -1;
4206 bool finished = false;
4207 while (!finished) {
4208 returnValue = firstFree_;
4209 int iColumn = firstFree_ + 1;
4210 if (superBasicType > 1) {
4211 if (superBasicType > 2) {
4212 // Initialize list
4213 // Wild guess that lower bound more natural than upper
4214 int number = 0;
4215 double *work = columnArray->denseVector();
4216 int *which = columnArray->getIndices();
4217 for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
4218 if (!flagged(iColumn)) {
4219 if (getInternalStatus(iColumn) == superBasic) {
4220 if (fabs(abcSolution_[iColumn] - abcLower_[iColumn]) <= primalTolerance_) {
4221 abcSolution_[iColumn] = abcLower_[iColumn];
4222 setInternalStatus(iColumn, atLowerBound);
4223 } else if (fabs(abcSolution_[iColumn] - abcUpper_[iColumn])
4224 <= primalTolerance_) {
4225 abcSolution_[iColumn] = abcUpper_[iColumn];
4226 setInternalStatus(iColumn, atUpperBound);
4227 } else if (abcLower_[iColumn] < -1.0e20 && abcUpper_[iColumn] > 1.0e20) {
4228 setInternalStatus(iColumn, isFree);
4229 break;
4230 } else if (!flagged(iColumn)) {
4231 // put ones near bounds at end after sorting
4232 work[number] = -CoinMin(0.1 * (abcSolution_[iColumn] - abcLower_[iColumn]),
4233 abcUpper_[iColumn] - abcSolution_[iColumn]);
4234 which[number++] = iColumn;
4235 }
4236 }
4237 }
4238 }
4239 CoinSort_2(work, work + number, which);
4240 columnArray->setNumElements(number);
4241 CoinZeroN(work, number);
4242 }
4243 int *which = columnArray->getIndices();
4244 int number = columnArray->getNumElements();
4245 if (!number) {
4246 // finished
4247 iColumn = numberRows_ + numberColumns_;
4248 returnValue = -1;
4249 } else {
4250 number--;
4251 returnValue = which[number];
4252 iColumn = returnValue;
4253 columnArray->setNumElements(number);
4254 }
4255 } else {
4256 for (; iColumn < numberRows_ + numberColumns_; iColumn++) {
4257 if (!flagged(iColumn)) {
4258 if (getInternalStatus(iColumn) == superBasic || getInternalStatus(iColumn) == isFree) {
4259 if (fabs(abcSolution_[iColumn] - abcLower_[iColumn]) <= primalTolerance_) {
4260 abcSolution_[iColumn] = abcLower_[iColumn];
4261 setInternalStatus(iColumn, atLowerBound);
4262 } else if (fabs(abcSolution_[iColumn] - abcUpper_[iColumn])
4263 <= primalTolerance_) {
4264 abcSolution_[iColumn] = abcUpper_[iColumn];
4265 setInternalStatus(iColumn, atUpperBound);
4266 } else if (abcLower_[iColumn] < -1.0e20 && abcUpper_[iColumn] > 1.0e20) {
4267 setInternalStatus(iColumn, isFree);
4268 if (fabs(abcDj_[iColumn]) > 10.0 * dualTolerance_)
4269 break;
4270 } else {
4271 break;
4272 }
4273 }
4274 }
4275 }
4276 }
4277 firstFree_ = iColumn;
4278 finished = true;
4279 if (firstFree_ == numberRows_ + numberColumns_)
4280 firstFree_ = -1;
4281 if (returnValue >= 0 && getInternalStatus(returnValue) != superBasic && getInternalStatus(returnValue) != isFree)
4282 finished = false; // somehow picked up odd one
4283 }
4284 return returnValue;
4285 }
clearAll()4286 void AbcSimplexPrimal::clearAll()
4287 {
4288 int number = usefulArray_[arrayForFtran_].getNumElements();
4289 int *which = usefulArray_[arrayForFtran_].getIndices();
4290
4291 int iIndex;
4292 for (iIndex = 0; iIndex < number; iIndex++) {
4293
4294 int iRow = which[iIndex];
4295 clearActive(iRow);
4296 }
4297 usefulArray_[arrayForFtran_].clear();
4298 for (int i = 0; i < 6; i++) {
4299 if (rowArray_[i])
4300 rowArray_[i]->clear();
4301 if (columnArray_[i])
4302 columnArray_[i]->clear();
4303 }
4304 }
4305
4306 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
4307 */
4308