1 /* $Id: ClpSimplexDual.cpp 1870 2012-07-22 16:13:48Z stefan $ */
2 // Copyright (C) 2002, International Business Machines
3 // Corporation and others. All Rights Reserved.
4 // This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7 /* Notes on implementation of dual simplex algorithm.
8
9 When dual feasible:
10
11 If primal feasible, we are optimal. Otherwise choose an infeasible
12 basic variable to leave basis (normally going to nearest bound) (B). We
13 now need to find an incoming variable which will leave problem
14 dual feasible so we get the row of the tableau corresponding to
15 the basic variable (with the correct sign depending if basic variable
16 above or below feasibility region - as that affects whether reduced
17 cost on outgoing variable has to be positive or negative).
18
19 We now perform a ratio test to determine which incoming variable will
20 preserve dual feasibility (C). If no variable found then problem
21 is infeasible (in primal sense). If there is a variable, we then
22 perform pivot and repeat. Trivial?
23
24 -------------------------------------------
25
26 A) How do we get dual feasible? If all variables have bounds then
27 it is trivial to get feasible by putting non-basic variables to
28 correct bounds. OSL did not have a phase 1/phase 2 approach but
29 instead effectively put fake bounds on variables and this is the
30 approach here, although I had hoped to make it cleaner.
31
32 If there is a weight of X on getting dual feasible:
33 Non-basic variables with negative reduced costs are put to
34 lesser of their upper bound and their lower bound + X.
35 Similarly, mutatis mutandis, for positive reduced costs.
36
37 Free variables should normally be in basis, otherwise I have
38 coding which may be able to come out (and may not be correct).
39
40 In OSL, this weight was changed heuristically, here at present
41 it is only increased if problem looks finished. If problem is
42 feasible I check for unboundedness. If not unbounded we
43 could play with going into primal. As long as weights increase
44 any algorithm would be finite.
45
46 B) Which outgoing variable to choose is a virtual base class.
47 For difficult problems steepest edge is preferred while for
48 very easy (large) problems we will need partial scan.
49
50 C) Sounds easy, but this is hardest part of algorithm.
51 1) Instead of stopping at first choice, we may be able
52 to flip that variable to other bound and if objective
53 still improving choose again. These mini iterations can
54 increase speed by orders of magnitude but we may need to
55 go to more of a bucket choice of variable rather than looking
56 at them one by one (for speed).
57 2) Accuracy. Reduced costs may be of wrong sign but less than
58 tolerance. Pivoting on these makes objective go backwards.
59 OSL modified cost so a zero move was made, Gill et al
60 (in primal analogue) modified so a strictly positive move was
61 made. It is not quite as neat in dual but that is what we
62 try and do. The two problems are that re-factorizations can
63 change reduced costs above and below tolerances and that when
64 finished we need to reset costs and try again.
65 3) Degeneracy. Gill et al helps but may not be enough. We
66 may need more. Also it can improve speed a lot if we perturb
67 the costs significantly.
68
69 References:
70 Forrest and Goldfarb, Steepest-edge simplex algorithms for
71 linear programming - Mathematical Programming 1992
72 Forrest and Tomlin, Implementing the simplex method for
73 the Optimization Subroutine Library - IBM Systems Journal 1992
74 Gill, Murray, Saunders, Wright A Practical Anti-Cycling
75 Procedure for Linear and Nonlinear Programming SOL report 1988
76
77
78 TODO:
79
80 a) Better recovery procedures. At present I never check on forward
81 progress. There is checkpoint/restart with reducing
82 re-factorization frequency, but this is only on singular
83 factorizations.
84 b) Fast methods for large easy problems (and also the option for
85 the code to automatically choose which method).
86 c) We need to be able to stop in various ways for OSI - this
87 is fairly easy.
88
89 */
90 #ifdef COIN_DEVELOP
91 #undef COIN_DEVELOP
92 #define COIN_DEVELOP 2
93 #endif
94
95 #include "CoinPragma.hpp"
96
97 #include <math.h>
98
99 #include "CoinHelperFunctions.hpp"
100 #include "ClpHelperFunctions.hpp"
101 #include "ClpSimplexDual.hpp"
102 #include "ClpEventHandler.hpp"
103 #include "ClpFactorization.hpp"
104 #include "CoinPackedMatrix.hpp"
105 #include "CoinIndexedVector.hpp"
106 #include "CoinFloatEqual.hpp"
107 #include "ClpDualRowDantzig.hpp"
108 #include "ClpMessage.hpp"
109 #include "ClpLinearObjective.hpp"
110 #include <cfloat>
111 #include <cassert>
112 #include <string>
113 #include <stdio.h>
114 #include <iostream>
115 //#define CLP_DEBUG 1
116 // To force to follow another run put logfile name here and define
117 //#define FORCE_FOLLOW
118 #ifdef FORCE_FOLLOW
119 static FILE * fpFollow = NULL;
120 static char * forceFile = "old.log";
121 static int force_in = -1;
122 static int force_out = -1;
123 static int force_iteration = 0;
124 #endif
125 //#define VUB
126 #ifdef VUB
127 extern int * vub;
128 extern int * toVub;
129 extern int * nextDescendent;
130 #endif
131 #ifdef NDEBUG
132 #define NDEBUG_CLP
133 #endif
134 #ifndef CLP_INVESTIGATE
135 #define NDEBUG_CLP
136 #endif
137 // dual
138
139 /* *** Method
140 This is a vanilla version of dual simplex.
141
142 It tries to be a single phase approach with a weight of 1.0 being
143 given to getting optimal and a weight of dualBound_ being
144 given to getting dual feasible. In this version I have used the
145 idea that this weight can be thought of as a fake bound. If the
146 distance between the lower and upper bounds on a variable is less
147 than the feasibility weight then we are always better off flipping
148 to other bound to make dual feasible. If the distance is greater
149 then we make up a fake bound dualBound_ away from one bound.
150 If we end up optimal or primal infeasible, we check to see if
151 bounds okay. If so we have finished, if not we increase dualBound_
152 and continue (after checking if unbounded). I am undecided about
153 free variables - there is coding but I am not sure about it. At
154 present I put them in basis anyway.
155
156 The code is designed to take advantage of sparsity so arrays are
157 seldom zeroed out from scratch or gone over in their entirety.
158 The only exception is a full scan to find outgoing variable. This
159 will be changed to keep an updated list of infeasibilities (or squares
160 if steepest edge). Also on easy problems we don't need full scan - just
161 pick first reasonable.
162
163 One problem is how to tackle degeneracy and accuracy. At present
164 I am using the modification of costs which I put in OSL and which was
165 extended by Gill et al. I am still not sure of the exact details.
166
167 The flow of dual is three while loops as follows:
168
169 while (not finished) {
170
171 while (not clean solution) {
172
173 Factorize and/or clean up solution by flipping variables so
174 dual feasible. If looks finished check fake dual bounds.
175 Repeat until status is iterating (-1) or finished (0,1,2)
176
177 }
178
179 while (status==-1) {
180
181 Iterate until no pivot in or out or time to re-factorize.
182
183 Flow is:
184
185 choose pivot row (outgoing variable). if none then
186 we are primal feasible so looks as if done but we need to
187 break and check bounds etc.
188
189 Get pivot row in tableau
190
191 Choose incoming column. If we don't find one then we look
192 primal infeasible so break and check bounds etc. (Also the
193 pivot tolerance is larger after any iterations so that may be
194 reason)
195
196 If we do find incoming column, we may have to adjust costs to
197 keep going forwards (anti-degeneracy). Check pivot will be stable
198 and if unstable throw away iteration (we will need to implement
199 flagging of basic variables sometime) and break to re-factorize.
200 If minor error re-factorize after iteration.
201
202 Update everything (this may involve flipping variables to stay
203 dual feasible.
204
205 }
206
207 }
208
209 At present we never check we are going forwards. I overdid that in
210 OSL so will try and make a last resort.
211
212 Needs partial scan pivot out option.
213 Needs dantzig, uninitialized and full steepest edge options (can still
214 use partial scan)
215
216 May need other anti-degeneracy measures, especially if we try and use
217 loose tolerances as a way to solve in fewer iterations.
218
219 I like idea of dynamic scaling. This gives opportunity to decouple
220 different implications of scaling for accuracy, iteration count and
221 feasibility tolerance.
222
223 */
224 #define CLEAN_FIXED 0
225 // Startup part of dual (may be extended to other algorithms)
226 int
startupSolve(int ifValuesPass,double * saveDuals,int startFinishOptions)227 ClpSimplexDual::startupSolve(int ifValuesPass, double * saveDuals, int startFinishOptions)
228 {
229 // If values pass then save given duals round check solution
230 // sanity check
231 // initialize - no values pass and algorithm_ is -1
232 // put in standard form (and make row copy)
233 // create modifiable copies of model rim and do optional scaling
234 // If problem looks okay
235 // Do initial factorization
236 // If user asked for perturbation - do it
237 numberFake_ = 0; // Number of variables at fake bounds
238 numberChanged_ = 0; // Number of variables with changed costs
239 if (!startup(0, startFinishOptions)) {
240 int usePrimal = 0;
241 // looks okay
242 // Superbasic variables not allowed
243 // If values pass then scale pi
244 if (ifValuesPass) {
245 if (problemStatus_ && perturbation_ < 100)
246 usePrimal = perturb();
247 int i;
248 if (scalingFlag_ > 0) {
249 for (i = 0; i < numberRows_; i++) {
250 dual_[i] = saveDuals[i] * inverseRowScale_[i];
251 }
252 } else {
253 CoinMemcpyN(saveDuals, numberRows_, dual_);
254 }
255 // now create my duals
256 for (i = 0; i < numberRows_; i++) {
257 // slack
258 double value = dual_[i];
259 value += rowObjectiveWork_[i];
260 saveDuals[i+numberColumns_] = value;
261 }
262 CoinMemcpyN(objectiveWork_, numberColumns_, saveDuals);
263 transposeTimes(-1.0, dual_, saveDuals);
264 // make reduced costs okay
265 for (i = 0; i < numberColumns_; i++) {
266 if (getStatus(i) == atLowerBound) {
267 if (saveDuals[i] < 0.0) {
268 //if (saveDuals[i]<-1.0e-3)
269 //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
270 saveDuals[i] = 0.0;
271 }
272 } else if (getStatus(i) == atUpperBound) {
273 if (saveDuals[i] > 0.0) {
274 //if (saveDuals[i]>1.0e-3)
275 //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
276 saveDuals[i] = 0.0;
277 }
278 }
279 }
280 CoinMemcpyN(saveDuals, (numberColumns_ + numberRows_), dj_);
281 // set up possible ones
282 for (i = 0; i < numberRows_ + numberColumns_; i++)
283 clearPivoted(i);
284 int iRow;
285 for (iRow = 0; iRow < numberRows_; iRow++) {
286 int iPivot = pivotVariable_[iRow];
287 if (fabs(saveDuals[iPivot]) > dualTolerance_) {
288 if (getStatus(iPivot) != isFree)
289 setPivoted(iPivot);
290 }
291 }
292 } else if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
293 // set up possible ones
294 for (int i = 0; i < numberRows_ + numberColumns_; i++)
295 clearPivoted(i);
296 int iRow;
297 for (iRow = 0; iRow < numberRows_; iRow++) {
298 int iPivot = pivotVariable_[iRow];
299 if (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]) {
300 setPivoted(iPivot);
301 }
302 }
303 }
304
305 double objectiveChange;
306 assert (!numberFake_);
307 assert (numberChanged_ == 0);
308 if (!numberFake_) // if nonzero then adjust
309 changeBounds(1, NULL, objectiveChange);
310
311 if (!ifValuesPass) {
312 // Check optimal
313 if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_)
314 problemStatus_ = 0;
315 }
316 if (problemStatus_ < 0 && perturbation_ < 100) {
317 bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
318 if (!inCbcOrOther)
319 usePrimal = perturb();
320 // Can't get here if values pass
321 gutsOfSolution(NULL, NULL);
322 #ifdef CLP_INVESTIGATE
323 if (numberDualInfeasibilities_)
324 printf("ZZZ %d primal %d dual - sumdinf %g\n",
325 numberPrimalInfeasibilities_,
326 numberDualInfeasibilities_, sumDualInfeasibilities_);
327 #endif
328 if (handler_->logLevel() > 2) {
329 handler_->message(CLP_SIMPLEX_STATUS, messages_)
330 << numberIterations_ << objectiveValue();
331 handler_->printing(sumPrimalInfeasibilities_ > 0.0)
332 << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
333 handler_->printing(sumDualInfeasibilities_ > 0.0)
334 << sumDualInfeasibilities_ << numberDualInfeasibilities_;
335 handler_->printing(numberDualInfeasibilitiesWithoutFree_
336 < numberDualInfeasibilities_)
337 << numberDualInfeasibilitiesWithoutFree_;
338 handler_->message() << CoinMessageEol;
339 }
340 if (inCbcOrOther) {
341 if (numberPrimalInfeasibilities_) {
342 usePrimal = perturb();
343 if (perturbation_ >= 101) {
344 computeDuals(NULL);
345 //gutsOfSolution(NULL,NULL);
346 checkDualSolution(); // recompute objective
347 }
348 } else if (numberDualInfeasibilities_) {
349 problemStatus_ = 10;
350 if ((moreSpecialOptions_ & 32) != 0 && false)
351 problemStatus_ = 0; // say optimal!!
352 #if COIN_DEVELOP>2
353
354 printf("returning at %d\n", __LINE__);
355 #endif
356 return 1; // to primal
357 }
358 }
359 } else if (!ifValuesPass) {
360 gutsOfSolution(NULL, NULL);
361 // double check
362 if (numberDualInfeasibilities_ || numberPrimalInfeasibilities_)
363 problemStatus_ = -1;
364 }
365 if (usePrimal) {
366 problemStatus_ = 10;
367 #if COIN_DEVELOP>2
368 printf("returning to use primal (no obj) at %d\n", __LINE__);
369 #endif
370 }
371 return usePrimal;
372 } else {
373 return 1;
374 }
375 }
376 void
finishSolve(int startFinishOptions)377 ClpSimplexDual::finishSolve(int startFinishOptions)
378 {
379 assert (problemStatus_ || !sumPrimalInfeasibilities_);
380
381 // clean up
382 finish(startFinishOptions);
383 }
384 //#define CLP_REPORT_PROGRESS
385 #ifdef CLP_REPORT_PROGRESS
386 static int ixxxxxx = 0;
387 static int ixxyyyy = 90;
388 #endif
389 #ifdef CLP_INVESTIGATE_SERIAL
390 static int z_reason[7] = {0, 0, 0, 0, 0, 0, 0};
391 static int z_thinks = -1;
392 #endif
393 void
gutsOfDual(int ifValuesPass,double * & saveDuals,int initialStatus,ClpDataSave & data)394 ClpSimplexDual::gutsOfDual(int ifValuesPass, double * & saveDuals, int initialStatus,
395 ClpDataSave & data)
396 {
397 #ifdef CLP_INVESTIGATE_SERIAL
398 z_reason[0]++;
399 z_thinks = -1;
400 int nPivots = 9999;
401 #endif
402 double largestPrimalError = 0.0;
403 double largestDualError = 0.0;
404 // Start can skip some things in transposeTimes
405 specialOptions_ |= 131072;
406 int lastCleaned = 0; // last time objective or bounds cleaned up
407
408 // This says whether to restore things etc
409 // startup will have factorized so can skip
410 int factorType = 0;
411 // Start check for cycles
412 progress_.startCheck();
413 // Say change made on first iteration
414 changeMade_ = 1;
415 // Say last objective infinite
416 //lastObjectiveValue_=-COIN_DBL_MAX;
417 progressFlag_ = 0;
418 /*
419 Status of problem:
420 0 - optimal
421 1 - infeasible
422 2 - unbounded
423 -1 - iterating
424 -2 - factorization wanted
425 -3 - redo checking without factorization
426 -4 - looks infeasible
427 */
428 while (problemStatus_ < 0) {
429 int iRow, iColumn;
430 // clear
431 for (iRow = 0; iRow < 4; iRow++) {
432 rowArray_[iRow]->clear();
433 }
434
435 for (iColumn = 0; iColumn < 2; iColumn++) {
436 columnArray_[iColumn]->clear();
437 }
438
439 // give matrix (and model costs and bounds a chance to be
440 // refreshed (normally null)
441 matrix_->refresh(this);
442 // If getting nowhere - why not give it a kick
443 // does not seem to work too well - do some more work
444 if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_)
445 && initialStatus != 10) {
446 perturb();
447 // Can't get here if values pass
448 gutsOfSolution(NULL, NULL);
449 if (handler_->logLevel() > 2) {
450 handler_->message(CLP_SIMPLEX_STATUS, messages_)
451 << numberIterations_ << objectiveValue();
452 handler_->printing(sumPrimalInfeasibilities_ > 0.0)
453 << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
454 handler_->printing(sumDualInfeasibilities_ > 0.0)
455 << sumDualInfeasibilities_ << numberDualInfeasibilities_;
456 handler_->printing(numberDualInfeasibilitiesWithoutFree_
457 < numberDualInfeasibilities_)
458 << numberDualInfeasibilitiesWithoutFree_;
459 handler_->message() << CoinMessageEol;
460 }
461 }
462 // see if in Cbc etc
463 bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
464 #if 0
465 bool gotoPrimal = false;
466 if (inCbcOrOther && numberIterations_ > disasterArea_ + numberRows_ &&
467 numberDualInfeasibilitiesWithoutFree_ && largestDualError_ > 1.0e-1) {
468 if (!disasterArea_) {
469 printf("trying all slack\n");
470 // try all slack basis
471 allSlackBasis(true);
472 disasterArea_ = 2 * numberRows_;
473 } else {
474 printf("going to primal\n");
475 // go to primal
476 gotoPrimal = true;
477 allSlackBasis(true);
478 }
479 }
480 #endif
481 bool disaster = false;
482 if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
483 disasterArea_->saveInfo();
484 disaster = true;
485 }
486 // may factorize, checks if problem finished
487 statusOfProblemInDual(lastCleaned, factorType, saveDuals, data,
488 ifValuesPass);
489 largestPrimalError = CoinMax(largestPrimalError, largestPrimalError_);
490 largestDualError = CoinMax(largestDualError, largestDualError_);
491 if (disaster)
492 problemStatus_ = 3;
493 // If values pass then do easy ones on first time
494 if (ifValuesPass &&
495 progress_.lastIterationNumber(0) < 0 && saveDuals) {
496 doEasyOnesInValuesPass(saveDuals);
497 }
498
499 // Say good factorization
500 factorType = 1;
501 if (data.sparseThreshold_) {
502 // use default at present
503 factorization_->sparseThreshold(0);
504 factorization_->goSparse();
505 }
506
507 // exit if victory declared
508 if (problemStatus_ >= 0)
509 break;
510
511 // test for maximum iterations
512 if (hitMaximumIterations() || (ifValuesPass == 2 && !saveDuals)) {
513 problemStatus_ = 3;
514 break;
515 }
516 if (ifValuesPass && !saveDuals) {
517 // end of values pass
518 ifValuesPass = 0;
519 int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
520 if (status >= 0) {
521 problemStatus_ = 5;
522 secondaryStatus_ = ClpEventHandler::endOfValuesPass;
523 break;
524 }
525 }
526 // Check event
527 {
528 int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
529 if (status >= 0) {
530 problemStatus_ = 5;
531 secondaryStatus_ = ClpEventHandler::endOfFactorization;
532 break;
533 }
534 }
535 // Do iterations
536 int returnCode = whileIterating(saveDuals, ifValuesPass);
537 if (problemStatus_ == 1 && (progressFlag_&8) != 0 &&
538 fabs(objectiveValue_) > 1.0e10 )
539 problemStatus_ = 10; // infeasible - but has looked feasible
540 #ifdef CLP_INVESTIGATE_SERIAL
541 nPivots = factorization_->pivots();
542 #endif
543 if (!problemStatus_ && factorization_->pivots())
544 computeDuals(NULL); // need to compute duals
545 if (returnCode == -2)
546 factorType = 3;
547 }
548 #ifdef CLP_INVESTIGATE_SERIAL
549 // NOTE - can fail if parallel
550 if (z_thinks != -1) {
551 assert (z_thinks < 4);
552 if ((!factorization_->pivots() && nPivots < 20) && z_thinks >= 0 && z_thinks < 2)
553 z_thinks += 4;
554 z_reason[1+z_thinks]++;
555 }
556 if ((z_reason[0] % 1000) == 0) {
557 printf("Reason");
558 for (int i = 0; i < 7; i++)
559 printf(" %d", z_reason[i]);
560 printf("\n");
561 }
562 #endif
563 // Stop can skip some things in transposeTimes
564 specialOptions_ &= ~131072;
565 largestPrimalError_ = largestPrimalError;
566 largestDualError_ = largestDualError;
567 }
568 int
dual(int ifValuesPass,int startFinishOptions)569 ClpSimplexDual::dual(int ifValuesPass, int startFinishOptions)
570 {
571 bestObjectiveValue_ = -COIN_DBL_MAX;
572 algorithm_ = -1;
573 moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
574 // save data
575 ClpDataSave data = saveData();
576 double * saveDuals = NULL;
577 int saveDont = dontFactorizePivots_;
578 if ((specialOptions_ & 2048) == 0)
579 dontFactorizePivots_ = 0;
580 else if(!dontFactorizePivots_)
581 dontFactorizePivots_ = 20;
582 if (ifValuesPass) {
583 saveDuals = new double [numberRows_+numberColumns_];
584 CoinMemcpyN(dual_, numberRows_, saveDuals);
585 }
586 if (alphaAccuracy_ != -1.0)
587 alphaAccuracy_ = 1.0;
588 int returnCode = startupSolve(ifValuesPass, saveDuals, startFinishOptions);
589 // Save so can see if doing after primal
590 int initialStatus = problemStatus_;
591 if (!returnCode && !numberDualInfeasibilities_ &&
592 !numberPrimalInfeasibilities_ && perturbation_ < 101) {
593 returnCode = 1; // to skip gutsOfDual
594 problemStatus_ = 0;
595 }
596
597 if (!returnCode)
598 gutsOfDual(ifValuesPass, saveDuals, initialStatus, data);
599 if (!problemStatus_) {
600 // see if cutoff reached
601 double limit = 0.0;
602 getDblParam(ClpDualObjectiveLimit, limit);
603 if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
604 limit + 1.0e-7 + 1.0e-8 * fabs(limit)) {
605 // actually infeasible on objective
606 problemStatus_ = 1;
607 secondaryStatus_ = 1;
608 }
609 }
610 // If infeasible but primal errors - try dual
611 if (problemStatus_==1 && numberPrimalInfeasibilities_) {
612 bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
613 double factor = (!inCbcOrOther) ? 1.0 : 0.3;
614 double averageInfeasibility = sumPrimalInfeasibilities_/
615 static_cast<double>(numberPrimalInfeasibilities_);
616 if (averageInfeasibility<factor*largestPrimalError_)
617 problemStatus_= 10;
618 }
619
620 if (problemStatus_ == 10)
621 startFinishOptions |= 1;
622 finishSolve(startFinishOptions);
623 delete [] saveDuals;
624
625 // Restore any saved stuff
626 restoreData(data);
627 dontFactorizePivots_ = saveDont;
628 if (problemStatus_ == 3)
629 objectiveValue_ = CoinMax(bestObjectiveValue_, objectiveValue_ - bestPossibleImprovement_);
630 return problemStatus_;
631 }
632 // old way
633 #if 0
634 int ClpSimplexDual::dual (int ifValuesPass , int startFinishOptions)
635 {
636 algorithm_ = -1;
637
638 // save data
639 ClpDataSave data = saveData();
640 // Save so can see if doing after primal
641 int initialStatus = problemStatus_;
642
643 // If values pass then save given duals round check solution
644 double * saveDuals = NULL;
645 if (ifValuesPass) {
646 saveDuals = new double [numberRows_+numberColumns_];
647 CoinMemcpyN(dual_, numberRows_, saveDuals);
648 }
649 // sanity check
650 // initialize - no values pass and algorithm_ is -1
651 // put in standard form (and make row copy)
652 // create modifiable copies of model rim and do optional scaling
653 // If problem looks okay
654 // Do initial factorization
655 // If user asked for perturbation - do it
656 if (!startup(0, startFinishOptions)) {
657 // looks okay
658 // Superbasic variables not allowed
659 // If values pass then scale pi
660 if (ifValuesPass) {
661 if (problemStatus_ && perturbation_ < 100)
662 perturb();
663 int i;
664 if (scalingFlag_ > 0) {
665 for (i = 0; i < numberRows_; i++) {
666 dual_[i] = saveDuals[i] * inverseRowScale_[i];
667 }
668 } else {
669 CoinMemcpyN(saveDuals, numberRows_, dual_);
670 }
671 // now create my duals
672 for (i = 0; i < numberRows_; i++) {
673 // slack
674 double value = dual_[i];
675 value += rowObjectiveWork_[i];
676 saveDuals[i+numberColumns_] = value;
677 }
678 CoinMemcpyN(objectiveWork_, numberColumns_, saveDuals);
679 transposeTimes(-1.0, dual_, saveDuals);
680 // make reduced costs okay
681 for (i = 0; i < numberColumns_; i++) {
682 if (getStatus(i) == atLowerBound) {
683 if (saveDuals[i] < 0.0) {
684 //if (saveDuals[i]<-1.0e-3)
685 //printf("bad dj at lb %d %g\n",i,saveDuals[i]);
686 saveDuals[i] = 0.0;
687 }
688 } else if (getStatus(i) == atUpperBound) {
689 if (saveDuals[i] > 0.0) {
690 //if (saveDuals[i]>1.0e-3)
691 //printf("bad dj at ub %d %g\n",i,saveDuals[i]);
692 saveDuals[i] = 0.0;
693 }
694 }
695 }
696 CoinMemcpyN(saveDuals, numberColumns_ + numberRows_, dj_);
697 // set up possible ones
698 for (i = 0; i < numberRows_ + numberColumns_; i++)
699 clearPivoted(i);
700 int iRow;
701 for (iRow = 0; iRow < numberRows_; iRow++) {
702 int iPivot = pivotVariable_[iRow];
703 if (fabs(saveDuals[iPivot]) > dualTolerance_) {
704 if (getStatus(iPivot) != isFree)
705 setPivoted(iPivot);
706 }
707 }
708 } else if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
709 // set up possible ones
710 for (int i = 0; i < numberRows_ + numberColumns_; i++)
711 clearPivoted(i);
712 int iRow;
713 for (iRow = 0; iRow < numberRows_; iRow++) {
714 int iPivot = pivotVariable_[iRow];
715 if (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]) {
716 setPivoted(iPivot);
717 }
718 }
719 }
720
721 double objectiveChange;
722 numberFake_ = 0; // Number of variables at fake bounds
723 numberChanged_ = 0; // Number of variables with changed costs
724 changeBounds(1, NULL, objectiveChange);
725
726 int lastCleaned = 0; // last time objective or bounds cleaned up
727
728 if (!ifValuesPass) {
729 // Check optimal
730 if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_)
731 problemStatus_ = 0;
732 }
733 if (problemStatus_ < 0 && perturbation_ < 100) {
734 perturb();
735 // Can't get here if values pass
736 gutsOfSolution(NULL, NULL);
737 if (handler_->logLevel() > 2) {
738 handler_->message(CLP_SIMPLEX_STATUS, messages_)
739 << numberIterations_ << objectiveValue();
740 handler_->printing(sumPrimalInfeasibilities_ > 0.0)
741 << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
742 handler_->printing(sumDualInfeasibilities_ > 0.0)
743 << sumDualInfeasibilities_ << numberDualInfeasibilities_;
744 handler_->printing(numberDualInfeasibilitiesWithoutFree_
745 < numberDualInfeasibilities_)
746 << numberDualInfeasibilitiesWithoutFree_;
747 handler_->message() << CoinMessageEol;
748 }
749 }
750
751 // This says whether to restore things etc
752 // startup will have factorized so can skip
753 int factorType = 0;
754 // Start check for cycles
755 progress_.startCheck();
756 // Say change made on first iteration
757 changeMade_ = 1;
758 /*
759 Status of problem:
760 0 - optimal
761 1 - infeasible
762 2 - unbounded
763 -1 - iterating
764 -2 - factorization wanted
765 -3 - redo checking without factorization
766 -4 - looks infeasible
767 */
768 while (problemStatus_ < 0) {
769 int iRow, iColumn;
770 // clear
771 for (iRow = 0; iRow < 4; iRow++) {
772 rowArray_[iRow]->clear();
773 }
774
775 for (iColumn = 0; iColumn < 2; iColumn++) {
776 columnArray_[iColumn]->clear();
777 }
778
779 // give matrix (and model costs and bounds a chance to be
780 // refreshed (normally null)
781 matrix_->refresh(this);
782 // If getting nowhere - why not give it a kick
783 // does not seem to work too well - do some more work
784 if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_)
785 && initialStatus != 10) {
786 perturb();
787 // Can't get here if values pass
788 gutsOfSolution(NULL, NULL);
789 if (handler_->logLevel() > 2) {
790 handler_->message(CLP_SIMPLEX_STATUS, messages_)
791 << numberIterations_ << objectiveValue();
792 handler_->printing(sumPrimalInfeasibilities_ > 0.0)
793 << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
794 handler_->printing(sumDualInfeasibilities_ > 0.0)
795 << sumDualInfeasibilities_ << numberDualInfeasibilities_;
796 handler_->printing(numberDualInfeasibilitiesWithoutFree_
797 < numberDualInfeasibilities_)
798 << numberDualInfeasibilitiesWithoutFree_;
799 handler_->message() << CoinMessageEol;
800 }
801 }
802 // may factorize, checks if problem finished
803 statusOfProblemInDual(lastCleaned, factorType, saveDuals, data,
804 ifValuesPass);
805 // If values pass then do easy ones on first time
806 if (ifValuesPass &&
807 progress_.lastIterationNumber(0) < 0 && saveDuals) {
808 doEasyOnesInValuesPass(saveDuals);
809 }
810
811 // Say good factorization
812 factorType = 1;
813 if (data.sparseThreshold_) {
814 // use default at present
815 factorization_->sparseThreshold(0);
816 factorization_->goSparse();
817 }
818
819 // exit if victory declared
820 if (problemStatus_ >= 0)
821 break;
822
823 // test for maximum iterations
824 if (hitMaximumIterations() || (ifValuesPass == 2 && !saveDuals)) {
825 problemStatus_ = 3;
826 break;
827 }
828 if (ifValuesPass && !saveDuals) {
829 // end of values pass
830 ifValuesPass = 0;
831 int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
832 if (status >= 0) {
833 problemStatus_ = 5;
834 secondaryStatus_ = ClpEventHandler::endOfValuesPass;
835 break;
836 }
837 }
838 // Check event
839 {
840 int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
841 if (status >= 0) {
842 problemStatus_ = 5;
843 secondaryStatus_ = ClpEventHandler::endOfFactorization;
844 break;
845 }
846 }
847 // Do iterations
848 whileIterating(saveDuals, ifValuesPass);
849 }
850 }
851
852 assert (problemStatus_ || !sumPrimalInfeasibilities_);
853
854 // clean up
855 finish(startFinishOptions);
856 delete [] saveDuals;
857
858 // Restore any saved stuff
859 restoreData(data);
860 return problemStatus_;
861 }
862 #endif
863 //#define CHECK_ACCURACY
864 #ifdef CHECK_ACCURACY
865 static double zzzzzz[100000];
866 #endif
867 /* Reasons to come out:
868 -1 iterations etc
869 -2 inaccuracy
870 -3 slight inaccuracy (and done iterations)
871 +0 looks optimal (might be unbounded - but we will investigate)
872 +1 looks infeasible
873 +3 max iterations
874 */
875 int
whileIterating(double * & givenDuals,int ifValuesPass)876 ClpSimplexDual::whileIterating(double * & givenDuals, int ifValuesPass)
877 {
878 #ifdef CLP_INVESTIGATE_SERIAL
879 z_thinks = -1;
880 #endif
881 #ifdef CLP_DEBUG
882 int debugIteration = -1;
883 #endif
884 {
885 int i;
886 for (i = 0; i < 4; i++) {
887 rowArray_[i]->clear();
888 }
889 for (i = 0; i < 2; i++) {
890 columnArray_[i]->clear();
891 }
892 }
893 #ifdef CLP_REPORT_PROGRESS
894 double * savePSol = new double [numberRows_+numberColumns_];
895 double * saveDj = new double [numberRows_+numberColumns_];
896 double * saveCost = new double [numberRows_+numberColumns_];
897 unsigned char * saveStat = new unsigned char [numberRows_+numberColumns_];
898 #endif
899 // if can't trust much and long way from optimal then relax
900 if (largestPrimalError_ > 10.0)
901 factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
902 else
903 factorization_->relaxAccuracyCheck(1.0);
904 // status stays at -1 while iterating, >=0 finished, -2 to invert
905 // status -3 to go to top without an invert
906 int returnCode = -1;
907 double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
908
909 #if 0
910 // compute average infeasibility for backward test
911 double averagePrimalInfeasibility = sumPrimalInfeasibilities_ /
912 ((double ) (numberPrimalInfeasibilities_ + 1));
913 #endif
914
915 // Get dubious weights
916 CoinBigIndex * dubiousWeights = NULL;
917 #ifdef DUBIOUS_WEIGHTS
918 factorization_->getWeights(rowArray_[0]->getIndices());
919 dubiousWeights = matrix_->dubiousWeights(this, rowArray_[0]->getIndices());
920 #endif
921 // If values pass then get list of candidates
922 int * candidateList = NULL;
923 int numberCandidates = 0;
924 #ifdef CLP_DEBUG
925 bool wasInValuesPass = (givenDuals != NULL);
926 #endif
927 int candidate = -1;
928 if (givenDuals) {
929 assert (ifValuesPass);
930 ifValuesPass = 1;
931 candidateList = new int[numberRows_];
932 // move reduced costs across
933 CoinMemcpyN(givenDuals, numberRows_ + numberColumns_, dj_);
934 int iRow;
935 for (iRow = 0; iRow < numberRows_; iRow++) {
936 int iPivot = pivotVariable_[iRow];
937 if (flagged(iPivot))
938 continue;
939 if (fabs(dj_[iPivot]) > dualTolerance_) {
940 // for now safer to ignore free ones
941 if (lower_[iPivot] > -1.0e50 || upper_[iPivot] < 1.0e50)
942 if (pivoted(iPivot))
943 candidateList[numberCandidates++] = iRow;
944 } else {
945 clearPivoted(iPivot);
946 }
947 }
948 // and set first candidate
949 if (!numberCandidates) {
950 delete [] candidateList;
951 delete [] givenDuals;
952 givenDuals = NULL;
953 candidateList = NULL;
954 int iRow;
955 for (iRow = 0; iRow < numberRows_; iRow++) {
956 int iPivot = pivotVariable_[iRow];
957 clearPivoted(iPivot);
958 }
959 }
960 } else {
961 assert (!ifValuesPass);
962 }
963 #ifdef CHECK_ACCURACY
964 {
965 if (numberIterations_) {
966 int il = -1;
967 double largest = 1.0e-1;
968 int ilnb = -1;
969 double largestnb = 1.0e-8;
970 for (int i = 0; i < numberRows_ + numberColumns_; i++) {
971 double diff = fabs(solution_[i] - zzzzzz[i]);
972 if (diff > largest) {
973 largest = diff;
974 il = i;
975 }
976 if (getColumnStatus(i) != basic) {
977 if (diff > largestnb) {
978 largestnb = diff;
979 ilnb = i;
980 }
981 }
982 }
983 if (il >= 0 && ilnb < 0)
984 printf("largest diff of %g at %d, nonbasic %g at %d\n",
985 largest, il, largestnb, ilnb);
986 }
987 }
988 #endif
989 while (problemStatus_ == -1) {
990 //if (numberIterations_>=101624)
991 //resetFakeBounds(-1);
992 #ifdef CLP_DEBUG
993 if (givenDuals) {
994 double value5 = 0.0;
995 int i;
996 for (i = 0; i < numberRows_ + numberColumns_; i++) {
997 if (dj_[i] < -1.0e-6)
998 if (upper_[i] < 1.0e20)
999 value5 += dj_[i] * upper_[i];
1000 else
1001 printf("bad dj %g on %d with large upper status %d\n",
1002 dj_[i], i, status_[i] & 7);
1003 else if (dj_[i] > 1.0e-6)
1004 if (lower_[i] > -1.0e20)
1005 value5 += dj_[i] * lower_[i];
1006 else
1007 printf("bad dj %g on %d with large lower status %d\n",
1008 dj_[i], i, status_[i] & 7);
1009 }
1010 printf("Values objective Value %g\n", value5);
1011 }
1012 if ((handler_->logLevel() & 32) && wasInValuesPass) {
1013 double value5 = 0.0;
1014 int i;
1015 for (i = 0; i < numberRows_ + numberColumns_; i++) {
1016 if (dj_[i] < -1.0e-6)
1017 if (upper_[i] < 1.0e20)
1018 value5 += dj_[i] * upper_[i];
1019 else if (dj_[i] > 1.0e-6)
1020 if (lower_[i] > -1.0e20)
1021 value5 += dj_[i] * lower_[i];
1022 }
1023 printf("Values objective Value %g\n", value5);
1024 {
1025 int i;
1026 for (i = 0; i < numberRows_ + numberColumns_; i++) {
1027 int iSequence = i;
1028 double oldValue;
1029
1030 switch(getStatus(iSequence)) {
1031
1032 case basic:
1033 case ClpSimplex::isFixed:
1034 break;
1035 case isFree:
1036 case superBasic:
1037 abort();
1038 break;
1039 case atUpperBound:
1040 oldValue = dj_[iSequence];
1041 //assert (oldValue<=tolerance);
1042 assert (fabs(solution_[iSequence] - upper_[iSequence]) < 1.0e-7);
1043 break;
1044 case atLowerBound:
1045 oldValue = dj_[iSequence];
1046 //assert (oldValue>=-tolerance);
1047 assert (fabs(solution_[iSequence] - lower_[iSequence]) < 1.0e-7);
1048 break;
1049 }
1050 }
1051 }
1052 }
1053 #endif
1054 #ifdef CLP_DEBUG
1055 {
1056 int i;
1057 for (i = 0; i < 4; i++) {
1058 rowArray_[i]->checkClear();
1059 }
1060 for (i = 0; i < 2; i++) {
1061 columnArray_[i]->checkClear();
1062 }
1063 }
1064 #endif
1065 #if CLP_DEBUG>2
1066 // very expensive
1067 if (numberIterations_ > 3063 && numberIterations_ < 30700) {
1068 //handler_->setLogLevel(63);
1069 double saveValue = objectiveValue_;
1070 double * saveRow1 = new double[numberRows_];
1071 double * saveRow2 = new double[numberRows_];
1072 CoinMemcpyN(rowReducedCost_, numberRows_, saveRow1);
1073 CoinMemcpyN(rowActivityWork_, numberRows_, saveRow2);
1074 double * saveColumn1 = new double[numberColumns_];
1075 double * saveColumn2 = new double[numberColumns_];
1076 CoinMemcpyN(reducedCostWork_, numberColumns_, saveColumn1);
1077 CoinMemcpyN(columnActivityWork_, numberColumns_, saveColumn2);
1078 gutsOfSolution(NULL, NULL);
1079 printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
1080 numberIterations_,
1081 saveValue, objectiveValue_, sumDualInfeasibilities_);
1082 if (saveValue > objectiveValue_ + 1.0e-2)
1083 printf("**bad**\n");
1084 CoinMemcpyN(saveRow1, numberRows_, rowReducedCost_);
1085 CoinMemcpyN(saveRow2, numberRows_, rowActivityWork_);
1086 CoinMemcpyN(saveColumn1, numberColumns_, reducedCostWork_);
1087 CoinMemcpyN(saveColumn2, numberColumns_, columnActivityWork_);
1088 delete [] saveRow1;
1089 delete [] saveRow2;
1090 delete [] saveColumn1;
1091 delete [] saveColumn2;
1092 objectiveValue_ = saveValue;
1093 }
1094 #endif
1095 #if 0
1096 // if (factorization_->pivots()){
1097 {
1098 int iPivot;
1099 double * array = rowArray_[3]->denseVector();
1100 int i;
1101 for (iPivot = 0; iPivot < numberRows_; iPivot++) {
1102 int iSequence = pivotVariable_[iPivot];
1103 unpack(rowArray_[3], iSequence);
1104 factorization_->updateColumn(rowArray_[2], rowArray_[3]);
1105 assert (fabs(array[iPivot] - 1.0) < 1.0e-4);
1106 array[iPivot] = 0.0;
1107 for (i = 0; i < numberRows_; i++)
1108 assert (fabs(array[i]) < 1.0e-4);
1109 rowArray_[3]->clear();
1110 }
1111 }
1112 #endif
1113 #ifdef CLP_DEBUG
1114 {
1115 int iSequence, number = numberRows_ + numberColumns_;
1116 for (iSequence = 0; iSequence < number; iSequence++) {
1117 double lowerValue = lower_[iSequence];
1118 double upperValue = upper_[iSequence];
1119 double value = solution_[iSequence];
1120 if(getStatus(iSequence) != basic && getStatus(iSequence) != isFree) {
1121 assert(lowerValue > -1.0e20);
1122 assert(upperValue < 1.0e20);
1123 }
1124 switch(getStatus(iSequence)) {
1125
1126 case basic:
1127 break;
1128 case isFree:
1129 case superBasic:
1130 break;
1131 case atUpperBound:
1132 assert (fabs(value - upperValue) <= primalTolerance_) ;
1133 break;
1134 case atLowerBound:
1135 case ClpSimplex::isFixed:
1136 assert (fabs(value - lowerValue) <= primalTolerance_) ;
1137 break;
1138 }
1139 }
1140 }
1141 if(numberIterations_ == debugIteration) {
1142 printf("dodgy iteration coming up\n");
1143 }
1144 #endif
1145 #if 0
1146 printf("checking nz\n");
1147 for (int i = 0; i < 3; i++) {
1148 if (!rowArray_[i]->getNumElements())
1149 rowArray_[i]->checkClear();
1150 }
1151 #endif
1152 // choose row to go out
1153 // dualRow will go to virtual row pivot choice algorithm
1154 // make sure values pass off if it should be
1155 if (numberCandidates)
1156 candidate = candidateList[--numberCandidates];
1157 else
1158 candidate = -1;
1159 dualRow(candidate);
1160 if (pivotRow_ >= 0) {
1161 // we found a pivot row
1162 if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
1163 handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
1164 << pivotRow_
1165 << CoinMessageEol;
1166 }
1167 // check accuracy of weights
1168 dualRowPivot_->checkAccuracy();
1169 // Get good size for pivot
1170 // Allow first few iterations to take tiny
1171 double acceptablePivot = 1.0e-1 * acceptablePivot_;
1172 if (numberIterations_ > 100)
1173 acceptablePivot = acceptablePivot_;
1174 if (factorization_->pivots() > 10 ||
1175 (factorization_->pivots() && saveSumDual))
1176 acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
1177 else if (factorization_->pivots() > 5)
1178 acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
1179 else if (factorization_->pivots())
1180 acceptablePivot = acceptablePivot_; // relax
1181 // But factorizations complain if <1.0e-8
1182 //acceptablePivot=CoinMax(acceptablePivot,1.0e-8);
1183 double bestPossiblePivot = 1.0;
1184 // get sign for finding row of tableau
1185 if (candidate < 0) {
1186 // normal iteration
1187 // create as packed
1188 double direction = directionOut_;
1189 rowArray_[0]->createPacked(1, &pivotRow_, &direction);
1190 factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
1191 // Allow to do dualColumn0
1192 if (numberThreads_ < -1)
1193 spareIntArray_[0] = 1;
1194 spareDoubleArray_[0] = acceptablePivot;
1195 rowArray_[3]->clear();
1196 sequenceIn_ = -1;
1197 // put row of tableau in rowArray[0] and columnArray[0]
1198 assert (!rowArray_[1]->getNumElements());
1199 if (!scaledMatrix_) {
1200 if ((moreSpecialOptions_ & 8) != 0 && !rowScale_)
1201 spareIntArray_[0] = 1;
1202 matrix_->transposeTimes(this, -1.0,
1203 rowArray_[0], rowArray_[1], columnArray_[0]);
1204 } else {
1205 double * saveR = rowScale_;
1206 double * saveC = columnScale_;
1207 rowScale_ = NULL;
1208 columnScale_ = NULL;
1209 if ((moreSpecialOptions_ & 8) != 0)
1210 spareIntArray_[0] = 1;
1211 scaledMatrix_->transposeTimes(this, -1.0,
1212 rowArray_[0], rowArray_[1], columnArray_[0]);
1213 rowScale_ = saveR;
1214 columnScale_ = saveC;
1215 }
1216 #ifdef CLP_REPORT_PROGRESS
1217 memcpy(savePSol, solution_, (numberColumns_ + numberRows_)*sizeof(double));
1218 memcpy(saveDj, dj_, (numberColumns_ + numberRows_)*sizeof(double));
1219 memcpy(saveCost, cost_, (numberColumns_ + numberRows_)*sizeof(double));
1220 memcpy(saveStat, status_, (numberColumns_ + numberRows_)*sizeof(char));
1221 #endif
1222 // do ratio test for normal iteration
1223 bestPossiblePivot = dualColumn(rowArray_[0], columnArray_[0], rowArray_[3],
1224 columnArray_[1], acceptablePivot, dubiousWeights);
1225 } else {
1226 // Make sure direction plausible
1227 CoinAssert (upperOut_ < 1.0e50 || lowerOut_ > -1.0e50);
1228 // If in integer cleanup do direction using duals
1229 // may be wrong way round
1230 if(ifValuesPass == 2) {
1231 if (dual_[pivotRow_] > 0.0) {
1232 // this will give a -1 in pivot row (as slacks are -1.0)
1233 directionOut_ = 1;
1234 } else {
1235 directionOut_ = -1;
1236 }
1237 }
1238 if (directionOut_ < 0 && fabs(valueOut_ - upperOut_) > dualBound_ + primalTolerance_) {
1239 if (fabs(valueOut_ - upperOut_) > fabs(valueOut_ - lowerOut_))
1240 directionOut_ = 1;
1241 } else if (directionOut_ > 0 && fabs(valueOut_ - lowerOut_) > dualBound_ + primalTolerance_) {
1242 if (fabs(valueOut_ - upperOut_) < fabs(valueOut_ - lowerOut_))
1243 directionOut_ = -1;
1244 }
1245 double direction = directionOut_;
1246 rowArray_[0]->createPacked(1, &pivotRow_, &direction);
1247 factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
1248 // put row of tableau in rowArray[0] and columnArray[0]
1249 if (!scaledMatrix_) {
1250 matrix_->transposeTimes(this, -1.0,
1251 rowArray_[0], rowArray_[3], columnArray_[0]);
1252 } else {
1253 double * saveR = rowScale_;
1254 double * saveC = columnScale_;
1255 rowScale_ = NULL;
1256 columnScale_ = NULL;
1257 scaledMatrix_->transposeTimes(this, -1.0,
1258 rowArray_[0], rowArray_[3], columnArray_[0]);
1259 rowScale_ = saveR;
1260 columnScale_ = saveC;
1261 }
1262 acceptablePivot *= 10.0;
1263 // do ratio test
1264 if (ifValuesPass == 1) {
1265 checkPossibleValuesMove(rowArray_[0], columnArray_[0],
1266 acceptablePivot);
1267 } else {
1268 checkPossibleCleanup(rowArray_[0], columnArray_[0],
1269 acceptablePivot);
1270 if (sequenceIn_ < 0) {
1271 rowArray_[0]->clear();
1272 columnArray_[0]->clear();
1273 continue; // can't do anything
1274 }
1275 }
1276
1277 // recompute true dualOut_
1278 if (directionOut_ < 0) {
1279 dualOut_ = valueOut_ - upperOut_;
1280 } else {
1281 dualOut_ = lowerOut_ - valueOut_;
1282 }
1283 // check what happened if was values pass
1284 // may want to move part way i.e. movement
1285 bool normalIteration = (sequenceIn_ != sequenceOut_);
1286
1287 clearPivoted(sequenceOut_); // make sure won't be done again
1288 // see if end of values pass
1289 if (!numberCandidates) {
1290 int iRow;
1291 delete [] candidateList;
1292 delete [] givenDuals;
1293 candidate = -2; // -2 signals end
1294 givenDuals = NULL;
1295 candidateList = NULL;
1296 ifValuesPass = 1;
1297 for (iRow = 0; iRow < numberRows_; iRow++) {
1298 int iPivot = pivotVariable_[iRow];
1299 //assert (fabs(dj_[iPivot]),1.0e-5);
1300 clearPivoted(iPivot);
1301 }
1302 }
1303 if (!normalIteration) {
1304 //rowArray_[0]->cleanAndPackSafe(1.0e-60);
1305 //columnArray_[0]->cleanAndPackSafe(1.0e-60);
1306 updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
1307 if (candidate == -2)
1308 problemStatus_ = -2;
1309 continue; // skip rest of iteration
1310 } else {
1311 // recompute dualOut_
1312 if (directionOut_ < 0) {
1313 dualOut_ = valueOut_ - upperOut_;
1314 } else {
1315 dualOut_ = lowerOut_ - valueOut_;
1316 }
1317 }
1318 }
1319 if (sequenceIn_ >= 0) {
1320 // normal iteration
1321 // update the incoming column
1322 double btranAlpha = -alpha_ * directionOut_; // for check
1323 unpackPacked(rowArray_[1]);
1324 // moved into updateWeights - factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
1325 // and update dual weights (can do in parallel - with extra array)
1326 alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
1327 rowArray_[2],
1328 rowArray_[3],
1329 rowArray_[1]);
1330 // see if update stable
1331 #ifdef CLP_DEBUG
1332 if ((handler_->logLevel() & 32))
1333 printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
1334 #endif
1335 double checkValue = 1.0e-7;
1336 // if can't trust much and long way from optimal then relax
1337 if (largestPrimalError_ > 10.0)
1338 checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
1339 if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
1340 fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
1341 handler_->message(CLP_DUAL_CHECK, messages_)
1342 << btranAlpha
1343 << alpha_
1344 << CoinMessageEol;
1345 if (factorization_->pivots()) {
1346 dualRowPivot_->unrollWeights();
1347 problemStatus_ = -2; // factorize now
1348 rowArray_[0]->clear();
1349 rowArray_[1]->clear();
1350 columnArray_[0]->clear();
1351 returnCode = -2;
1352 break;
1353 } else {
1354 // take on more relaxed criterion
1355 double test;
1356 if (fabs(btranAlpha) < 1.0e-8 || fabs(alpha_) < 1.0e-8)
1357 test = 1.0e-1 * fabs(alpha_);
1358 else
1359 test = 1.0e-4 * (1.0 + fabs(alpha_));
1360 if (fabs(btranAlpha) < 1.0e-12 || fabs(alpha_) < 1.0e-12 ||
1361 fabs(btranAlpha - alpha_) > test) {
1362 dualRowPivot_->unrollWeights();
1363 // need to reject something
1364 char x = isColumn(sequenceOut_) ? 'C' : 'R';
1365 handler_->message(CLP_SIMPLEX_FLAG, messages_)
1366 << x << sequenceWithin(sequenceOut_)
1367 << CoinMessageEol;
1368 #ifdef COIN_DEVELOP
1369 printf("flag a %g %g\n", btranAlpha, alpha_);
1370 #endif
1371 //#define FEB_TRY
1372 #if 1 //def FEB_TRY
1373 // Make safer?
1374 factorization_->saferTolerances (-0.99, -1.03);
1375 #endif
1376 setFlagged(sequenceOut_);
1377 progress_.clearBadTimes();
1378 lastBadIteration_ = numberIterations_; // say be more cautious
1379 rowArray_[0]->clear();
1380 rowArray_[1]->clear();
1381 columnArray_[0]->clear();
1382 if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
1383 //printf("I think should declare infeasible\n");
1384 problemStatus_ = 1;
1385 returnCode = 1;
1386 break;
1387 }
1388 continue;
1389 }
1390 }
1391 }
1392 // update duals BEFORE replaceColumn so can do updateColumn
1393 double objectiveChange = 0.0;
1394 // do duals first as variables may flip bounds
1395 // rowArray_[0] and columnArray_[0] may have flips
1396 // so use rowArray_[3] for work array from here on
1397 int nswapped = 0;
1398 //rowArray_[0]->cleanAndPackSafe(1.0e-60);
1399 //columnArray_[0]->cleanAndPackSafe(1.0e-60);
1400 if (candidate == -1) {
1401 // make sure incoming doesn't count
1402 Status saveStatus = getStatus(sequenceIn_);
1403 setStatus(sequenceIn_, basic);
1404 nswapped = updateDualsInDual(rowArray_[0], columnArray_[0],
1405 rowArray_[2], theta_,
1406 objectiveChange, false);
1407 setStatus(sequenceIn_, saveStatus);
1408 } else {
1409 updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
1410 }
1411 double oldDualOut = dualOut_;
1412 // which will change basic solution
1413 if (nswapped) {
1414 if (rowArray_[2]->getNumElements()) {
1415 factorization_->updateColumn(rowArray_[3], rowArray_[2]);
1416 dualRowPivot_->updatePrimalSolution(rowArray_[2],
1417 1.0, objectiveChange);
1418 }
1419 // recompute dualOut_
1420 valueOut_ = solution_[sequenceOut_];
1421 if (directionOut_ < 0) {
1422 dualOut_ = valueOut_ - upperOut_;
1423 } else {
1424 dualOut_ = lowerOut_ - valueOut_;
1425 }
1426 #if 0
1427 if (dualOut_ < 0.0) {
1428 #ifdef CLP_DEBUG
1429 if (handler_->logLevel() & 32) {
1430 printf(" dualOut_ %g %g save %g\n", dualOut_, averagePrimalInfeasibility, saveDualOut);
1431 printf("values %g %g %g %g %g %g %g\n", lowerOut_, valueOut_, upperOut_,
1432 objectiveChange,);
1433 }
1434 #endif
1435 if (upperOut_ == lowerOut_)
1436 dualOut_ = 0.0;
1437 }
1438 if(dualOut_ < -CoinMax(1.0e-12 * averagePrimalInfeasibility, 1.0e-8)
1439 && factorization_->pivots() > 100 &&
1440 getStatus(sequenceIn_) != isFree) {
1441 // going backwards - factorize
1442 dualRowPivot_->unrollWeights();
1443 problemStatus_ = -2; // factorize now
1444 returnCode = -2;
1445 break;
1446 }
1447 #endif
1448 }
1449 // amount primal will move
1450 double movement = -dualOut_ * directionOut_ / alpha_;
1451 double movementOld = oldDualOut * directionOut_ / alpha_;
1452 // so objective should increase by fabs(dj)*movement
1453 // but we already have objective change - so check will be good
1454 if (objectiveChange + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(objectiveValue_))) {
1455 #ifdef CLP_DEBUG
1456 if (handler_->logLevel() & 32)
1457 printf("movement %g, swap change %g, rest %g * %g\n",
1458 objectiveChange + fabs(movement * dualIn_),
1459 objectiveChange, movement, dualIn_);
1460 #endif
1461 if(factorization_->pivots()) {
1462 // going backwards - factorize
1463 dualRowPivot_->unrollWeights();
1464 problemStatus_ = -2; // factorize now
1465 returnCode = -2;
1466 break;
1467 }
1468 }
1469 // if stable replace in basis
1470 int updateStatus = factorization_->replaceColumn(this,
1471 rowArray_[2],
1472 rowArray_[1],
1473 pivotRow_,
1474 alpha_,
1475 (moreSpecialOptions_ & 16) != 0,
1476 acceptablePivot);
1477 // If looks like bad pivot - refactorize
1478 if (fabs(dualOut_) > 1.0e50)
1479 updateStatus = 2;
1480 // if no pivots, bad update but reasonable alpha - take and invert
1481 if (updateStatus == 2 &&
1482 !factorization_->pivots() && fabs(alpha_) > 1.0e-5)
1483 updateStatus = 4;
1484 if (updateStatus == 1 || updateStatus == 4) {
1485 // slight error
1486 if (factorization_->pivots() > 5 || updateStatus == 4) {
1487 problemStatus_ = -2; // factorize now
1488 returnCode = -3;
1489 }
1490 } else if (updateStatus == 2) {
1491 // major error
1492 dualRowPivot_->unrollWeights();
1493 // later we may need to unwind more e.g. fake bounds
1494 if (factorization_->pivots() &&
1495 ((moreSpecialOptions_ & 16) == 0 || factorization_->pivots() > 4)) {
1496 problemStatus_ = -2; // factorize now
1497 returnCode = -2;
1498 moreSpecialOptions_ |= 16;
1499 break;
1500 } else {
1501 // need to reject something
1502 char x = isColumn(sequenceOut_) ? 'C' : 'R';
1503 handler_->message(CLP_SIMPLEX_FLAG, messages_)
1504 << x << sequenceWithin(sequenceOut_)
1505 << CoinMessageEol;
1506 #ifdef COIN_DEVELOP
1507 printf("flag b %g\n", alpha_);
1508 #endif
1509 setFlagged(sequenceOut_);
1510 progress_.clearBadTimes();
1511 lastBadIteration_ = numberIterations_; // say be more cautious
1512 rowArray_[0]->clear();
1513 rowArray_[1]->clear();
1514 columnArray_[0]->clear();
1515 // make sure dual feasible
1516 // look at all rows and columns
1517 double objectiveChange = 0.0;
1518 updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
1519 0.0, objectiveChange, true);
1520 rowArray_[1]->clear();
1521 columnArray_[0]->clear();
1522 continue;
1523 }
1524 } else if (updateStatus == 3) {
1525 // out of memory
1526 // increase space if not many iterations
1527 if (factorization_->pivots() <
1528 0.5 * factorization_->maximumPivots() &&
1529 factorization_->pivots() < 200)
1530 factorization_->areaFactor(
1531 factorization_->areaFactor() * 1.1);
1532 problemStatus_ = -2; // factorize now
1533 } else if (updateStatus == 5) {
1534 problemStatus_ = -2; // factorize now
1535 }
1536 // update primal solution
1537 if (theta_ < 0.0 && candidate == -1) {
1538 #ifdef CLP_DEBUG
1539 if (handler_->logLevel() & 32)
1540 printf("negative theta %g\n", theta_);
1541 #endif
1542 theta_ = 0.0;
1543 }
1544 // do actual flips
1545 flipBounds(rowArray_[0], columnArray_[0]);
1546 //rowArray_[1]->expand();
1547 dualRowPivot_->updatePrimalSolution(rowArray_[1],
1548 movement,
1549 objectiveChange);
1550 #ifdef CLP_DEBUG
1551 double oldobj = objectiveValue_;
1552 #endif
1553 // modify dualout
1554 dualOut_ /= alpha_;
1555 dualOut_ *= -directionOut_;
1556 //setStatus(sequenceIn_,basic);
1557 dj_[sequenceIn_] = 0.0;
1558 double oldValue = valueIn_;
1559 if (directionIn_ == -1) {
1560 // as if from upper bound
1561 valueIn_ = upperIn_ + dualOut_;
1562 } else {
1563 // as if from lower bound
1564 valueIn_ = lowerIn_ + dualOut_;
1565 }
1566 objectiveChange += cost_[sequenceIn_] * (valueIn_ - oldValue);
1567 // outgoing
1568 // set dj to zero unless values pass
1569 if (directionOut_ > 0) {
1570 valueOut_ = lowerOut_;
1571 if (candidate == -1)
1572 dj_[sequenceOut_] = theta_;
1573 } else {
1574 valueOut_ = upperOut_;
1575 if (candidate == -1)
1576 dj_[sequenceOut_] = -theta_;
1577 }
1578 solution_[sequenceOut_] = valueOut_;
1579 int whatNext = housekeeping(objectiveChange);
1580 #ifdef CLP_REPORT_PROGRESS
1581 if (ixxxxxx > ixxyyyy - 5) {
1582 handler_->setLogLevel(63);
1583 int nTotal = numberColumns_ + numberRows_;
1584 double oldObj = 0.0;
1585 double newObj = 0.0;
1586 for (int i = 0; i < nTotal; i++) {
1587 if (savePSol[i])
1588 oldObj += savePSol[i] * saveCost[i];
1589 if (solution_[i])
1590 newObj += solution_[i] * cost_[i];
1591 bool printIt = false;
1592 if (cost_[i] != saveCost[i])
1593 printIt = true;
1594 if (status_[i] != saveStat[i])
1595 printIt = true;
1596 if (printIt)
1597 printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
1598 i, saveStat[i], saveCost[i], savePSol[i],
1599 status_[i], cost_[i], solution_[i]);
1600 // difference
1601 savePSol[i] = solution_[i] - savePSol[i];
1602 }
1603 printf("pivots %d, old obj %g new %g\n",
1604 factorization_->pivots(),
1605 oldObj, newObj);
1606 memset(saveDj, 0, numberRows_ * sizeof(double));
1607 times(1.0, savePSol, saveDj);
1608 double largest = 1.0e-6;
1609 int k = -1;
1610 for (int i = 0; i < numberRows_; i++) {
1611 saveDj[i] -= savePSol[i+numberColumns_];
1612 if (fabs(saveDj[i]) > largest) {
1613 largest = fabs(saveDj[i]);
1614 k = i;
1615 }
1616 }
1617 if (k >= 0)
1618 printf("Not null %d %g\n", k, largest);
1619 }
1620 #endif
1621 #ifdef VUB
1622 {
1623 if ((sequenceIn_ < numberColumns_ && vub[sequenceIn_] >= 0) || toVub[sequenceIn_] >= 0 ||
1624 (sequenceOut_ < numberColumns_ && vub[sequenceOut_] >= 0) || toVub[sequenceOut_] >= 0) {
1625 int inSequence = sequenceIn_;
1626 int inVub = -1;
1627 if (sequenceIn_ < numberColumns_)
1628 inVub = vub[sequenceIn_];
1629 int inBack = toVub[inSequence];
1630 int inSlack = -1;
1631 if (inSequence >= numberColumns_ && inBack >= 0) {
1632 inSlack = inSequence - numberColumns_;
1633 inSequence = inBack;
1634 inBack = toVub[inSequence];
1635 }
1636 if (inVub >= 0)
1637 printf("Vub %d in ", inSequence);
1638 if (inBack >= 0 && inSlack < 0)
1639 printf("%d (descendent of %d) in ", inSequence, inBack);
1640 if (inSlack >= 0)
1641 printf("slack for row %d -> %d (descendent of %d) in ", inSlack, inSequence, inBack);
1642 int outSequence = sequenceOut_;
1643 int outVub = -1;
1644 if (sequenceOut_ < numberColumns_)
1645 outVub = vub[sequenceOut_];
1646 int outBack = toVub[outSequence];
1647 int outSlack = -1;
1648 if (outSequence >= numberColumns_ && outBack >= 0) {
1649 outSlack = outSequence - numberColumns_;
1650 outSequence = outBack;
1651 outBack = toVub[outSequence];
1652 }
1653 if (outVub >= 0)
1654 printf("Vub %d out ", outSequence);
1655 if (outBack >= 0 && outSlack < 0)
1656 printf("%d (descendent of %d) out ", outSequence, outBack);
1657 if (outSlack >= 0)
1658 printf("slack for row %d -> %d (descendent of %d) out ", outSlack, outSequence, outBack);
1659 printf("\n");
1660 }
1661 }
1662 #endif
1663 #if 0
1664 if (numberIterations_ > 206033)
1665 handler_->setLogLevel(63);
1666 if (numberIterations_ > 210567)
1667 exit(77);
1668 #endif
1669 if (!givenDuals && ifValuesPass && ifValuesPass != 2) {
1670 handler_->message(CLP_END_VALUES_PASS, messages_)
1671 << numberIterations_;
1672 whatNext = 1;
1673 }
1674 #ifdef CHECK_ACCURACY
1675 if (whatNext) {
1676 CoinMemcpyN(solution_, (numberRows_ + numberColumns_), zzzzzz);
1677 }
1678 #endif
1679 //if (numberIterations_==1890)
1680 //whatNext=1;
1681 //if (numberIterations_>2000)
1682 //exit(77);
1683 // and set bounds correctly
1684 originalBound(sequenceIn_);
1685 changeBound(sequenceOut_);
1686 #ifdef CLP_DEBUG
1687 if (objectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16))
1688 printf("obj backwards %g %g\n", objectiveValue_, oldobj);
1689 #endif
1690 #if 0
1691 {
1692 for (int i = 0; i < numberRows_ + numberColumns_; i++) {
1693 FakeBound bound = getFakeBound(i);
1694 if (bound == ClpSimplexDual::upperFake) {
1695 assert (upper_[i] < 1.0e20);
1696 } else if (bound == ClpSimplexDual::lowerFake) {
1697 assert (lower_[i] > -1.0e20);
1698 } else if (bound == ClpSimplexDual::bothFake) {
1699 assert (upper_[i] < 1.0e20);
1700 assert (lower_[i] > -1.0e20);
1701 }
1702 }
1703 }
1704 #endif
1705 if (whatNext == 1 || candidate == -2) {
1706 problemStatus_ = -2; // refactorize
1707 } else if (whatNext == 2) {
1708 // maximum iterations or equivalent
1709 problemStatus_ = 3;
1710 returnCode = 3;
1711 break;
1712 }
1713 // Check event
1714 {
1715 int status = eventHandler_->event(ClpEventHandler::endOfIteration);
1716 if (status >= 0) {
1717 problemStatus_ = 5;
1718 secondaryStatus_ = ClpEventHandler::endOfIteration;
1719 returnCode = 4;
1720 break;
1721 }
1722 }
1723 } else {
1724 #ifdef CLP_INVESTIGATE_SERIAL
1725 z_thinks = 1;
1726 #endif
1727 // no incoming column is valid
1728 pivotRow_ = -1;
1729 #ifdef CLP_DEBUG
1730 if (handler_->logLevel() & 32)
1731 printf("** no column pivot\n");
1732 #endif
1733 if (factorization_->pivots() < 2 && acceptablePivot_ <= 1.0e-8) {
1734 //&&goodAccuracy()) {
1735 // If not in branch and bound etc save ray
1736 delete [] ray_;
1737 if ((specialOptions_&(1024 | 4096)) == 0 || (specialOptions_ & 32) != 0) {
1738 // create ray anyway
1739 ray_ = new double [ numberRows_];
1740 rowArray_[0]->expand(); // in case packed
1741 CoinMemcpyN(rowArray_[0]->denseVector(), numberRows_, ray_);
1742 } else {
1743 ray_ = NULL;
1744 }
1745 // If we have just factorized and infeasibility reasonable say infeas
1746 double dualTest = ((specialOptions_ & 4096) != 0) ? 1.0e8 : 1.0e13;
1747 if (((specialOptions_ & 4096) != 0 || bestPossiblePivot < 1.0e-11) && dualBound_ > dualTest) {
1748 double testValue = 1.0e-4;
1749 if (!factorization_->pivots() && numberPrimalInfeasibilities_ == 1)
1750 testValue = 1.0e-6;
1751 if (valueOut_ > upperOut_ + testValue || valueOut_ < lowerOut_ - testValue
1752 || (specialOptions_ & 64) == 0) {
1753 // say infeasible
1754 problemStatus_ = 1;
1755 // unless primal feasible!!!!
1756 //printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
1757 // numberDualInfeasibilities_,sumDualInfeasibilities_);
1758 //#define TEST_CLP_NODE
1759 #ifndef TEST_CLP_NODE
1760 // Should be correct - but ...
1761 int numberFake = numberAtFakeBound();
1762 double sumPrimal = (!numberFake) ? 2.0e5 : sumPrimalInfeasibilities_;
1763 if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-5 ||
1764 (sumPrimal < 1.0e5 && (specialOptions_ & 1024) != 0 && factorization_->pivots())) {
1765 if (sumPrimal > 50.0 && factorization_->pivots() > 2) {
1766 problemStatus_ = -4;
1767 #ifdef COIN_DEVELOP
1768 printf("status to -4 at %d - primalinf %g pivots %d\n",
1769 __LINE__, sumPrimalInfeasibilities_,
1770 factorization_->pivots());
1771 #endif
1772 } else {
1773 problemStatus_ = 10;
1774 #if COIN_DEVELOP>1
1775 printf("returning at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d - options (1024-16384) %d %d %d %d %d\n",
1776 __LINE__, numberPrimalInfeasibilities_,
1777 sumPrimalInfeasibilities_,
1778 numberDualInfeasibilities_, sumDualInfeasibilities_,
1779 numberFake_, dualBound_, factorization_->pivots(),
1780 (specialOptions_ & 1024) != 0 ? 1 : 0,
1781 (specialOptions_ & 2048) != 0 ? 1 : 0,
1782 (specialOptions_ & 4096) != 0 ? 1 : 0,
1783 (specialOptions_ & 8192) != 0 ? 1 : 0,
1784 (specialOptions_ & 16384) != 0 ? 1 : 0
1785 );
1786 #endif
1787 // Get rid of objective
1788 if ((specialOptions_ & 16384) == 0)
1789 objective_ = new ClpLinearObjective(NULL, numberColumns_);
1790 }
1791 }
1792 #else
1793 if (sumPrimalInfeasibilities_ < 1.0e-3 || sumDualInfeasibilities_ > 1.0e-6) {
1794 #ifdef COIN_DEVELOP
1795 printf("at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d\n",
1796 __LINE__, numberPrimalInfeasibilities_,
1797 sumPrimalInfeasibilities_,
1798 numberDualInfeasibilities_, sumDualInfeasibilities_,
1799 numberFake_, dualBound_, factorization_->pivots());
1800 #endif
1801 if ((specialOptions_ & 1024) != 0 && factorization_->pivots()) {
1802 problemStatus_ = 10;
1803 #if COIN_DEVELOP>1
1804 printf("returning at %d\n", __LINE__);
1805 #endif
1806 // Get rid of objective
1807 if ((specialOptions_ & 16384) == 0)
1808 objective_ = new ClpLinearObjective(NULL, numberColumns_);
1809 }
1810 }
1811 #endif
1812 rowArray_[0]->clear();
1813 columnArray_[0]->clear();
1814 returnCode = 1;
1815 break;
1816 }
1817 }
1818 // If special option set - put off as long as possible
1819 if ((specialOptions_ & 64) == 0 || (moreSpecialOptions_ & 64) != 0) {
1820 if (factorization_->pivots() == 0)
1821 problemStatus_ = -4; //say looks infeasible
1822 } else {
1823 // flag
1824 char x = isColumn(sequenceOut_) ? 'C' : 'R';
1825 handler_->message(CLP_SIMPLEX_FLAG, messages_)
1826 << x << sequenceWithin(sequenceOut_)
1827 << CoinMessageEol;
1828 #ifdef COIN_DEVELOP
1829 printf("flag c\n");
1830 #endif
1831 setFlagged(sequenceOut_);
1832 if (!factorization_->pivots()) {
1833 rowArray_[0]->clear();
1834 columnArray_[0]->clear();
1835 continue;
1836 }
1837 }
1838 }
1839 if (factorization_->pivots() < 5 && acceptablePivot_ > 1.0e-8)
1840 acceptablePivot_ = 1.0e-8;
1841 rowArray_[0]->clear();
1842 columnArray_[0]->clear();
1843 returnCode = 1;
1844 break;
1845 }
1846 } else {
1847 #ifdef CLP_INVESTIGATE_SERIAL
1848 z_thinks = 0;
1849 #endif
1850 // no pivot row
1851 #ifdef CLP_DEBUG
1852 if (handler_->logLevel() & 32)
1853 printf("** no row pivot\n");
1854 #endif
1855 // If in branch and bound try and get rid of fixed variables
1856 if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
1857 assert (!candidateList);
1858 candidateList = new int[numberRows_];
1859 int iRow;
1860 for (iRow = 0; iRow < numberRows_; iRow++) {
1861 int iPivot = pivotVariable_[iRow];
1862 if (flagged(iPivot) || !pivoted(iPivot))
1863 continue;
1864 assert (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]);
1865 candidateList[numberCandidates++] = iRow;
1866 }
1867 // and set first candidate
1868 if (!numberCandidates) {
1869 delete [] candidateList;
1870 candidateList = NULL;
1871 int iRow;
1872 for (iRow = 0; iRow < numberRows_; iRow++) {
1873 int iPivot = pivotVariable_[iRow];
1874 clearPivoted(iPivot);
1875 }
1876 } else {
1877 ifValuesPass = 2;
1878 continue;
1879 }
1880 }
1881 int numberPivots = factorization_->pivots();
1882 bool specialCase;
1883 int useNumberFake;
1884 returnCode = 0;
1885 if (numberPivots <= CoinMax(dontFactorizePivots_, 20) &&
1886 (specialOptions_ & 2048) != 0 && (true || !numberChanged_ || perturbation_ == 101)
1887 && dualBound_ >= 1.0e8) {
1888 specialCase = true;
1889 // as dual bound high - should be okay
1890 useNumberFake = 0;
1891 } else {
1892 specialCase = false;
1893 useNumberFake = numberFake_;
1894 }
1895 if (!numberPivots || specialCase) {
1896 // may have crept through - so may be optimal
1897 // check any flagged variables
1898 int iRow;
1899 for (iRow = 0; iRow < numberRows_; iRow++) {
1900 int iPivot = pivotVariable_[iRow];
1901 if (flagged(iPivot))
1902 break;
1903 }
1904 if (iRow < numberRows_ && numberPivots) {
1905 // try factorization
1906 returnCode = -2;
1907 }
1908
1909 if (useNumberFake || numberDualInfeasibilities_) {
1910 // may be dual infeasible
1911 if ((specialOptions_ & 1024) == 0)
1912 problemStatus_ = -5;
1913 else if (!useNumberFake && numberPrimalInfeasibilities_
1914 && !numberPivots)
1915 problemStatus_ = 1;
1916 } else {
1917 if (iRow < numberRows_) {
1918 #ifdef COIN_DEVELOP
1919 std::cout << "Flagged variables at end - infeasible?" << std::endl;
1920 printf("Probably infeasible - pivot was %g\n", alpha_);
1921 #endif
1922 //if (fabs(alpha_)<1.0e-4) {
1923 //problemStatus_=1;
1924 //} else {
1925 #ifdef CLP_DEBUG
1926 abort();
1927 #endif
1928 //}
1929 problemStatus_ = -5;
1930 } else {
1931 problemStatus_ = 0;
1932 #ifndef CLP_CHECK_NUMBER_PIVOTS
1933 #define CLP_CHECK_NUMBER_PIVOTS 10
1934 #endif
1935 #if CLP_CHECK_NUMBER_PIVOTS < 20
1936 if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
1937 #ifndef NDEBUG_CLP
1938 int nTotal = numberRows_ + numberColumns_;
1939 double * comp = CoinCopyOfArray(solution_, nTotal);
1940 #endif
1941 computePrimals(rowActivityWork_, columnActivityWork_);
1942 #ifndef NDEBUG_CLP
1943 double largest = 1.0e-5;
1944 int bad = -1;
1945 for (int i = 0; i < nTotal; i++) {
1946 double value = solution_[i];
1947 double larger = CoinMax(fabs(value), fabs(comp[i]));
1948 double tol = 1.0e-5 + 1.0e-5 * larger;
1949 double diff = fabs(value - comp[i]);
1950 if (diff - tol > largest) {
1951 bad = i;
1952 largest = diff - tol;
1953 }
1954 }
1955 if (bad >= 0)
1956 COIN_DETAIL_PRINT(printf("bad %d old %g new %g\n", bad, comp[bad], solution_[bad]));
1957 #endif
1958 checkPrimalSolution(rowActivityWork_, columnActivityWork_);
1959 if (numberPrimalInfeasibilities_) {
1960 #ifdef CLP_INVESTIGATE
1961 printf("XXX Infeas ? %d inf summing to %g\n", numberPrimalInfeasibilities_,
1962 sumPrimalInfeasibilities_);
1963 #endif
1964 problemStatus_ = -1;
1965 returnCode = -2;
1966 }
1967 #ifndef NDEBUG_CLP
1968 memcpy(solution_, comp, nTotal * sizeof(double));
1969 delete [] comp;
1970 #endif
1971 }
1972 #endif
1973 if (!problemStatus_) {
1974 // make it look OK
1975 numberPrimalInfeasibilities_ = 0;
1976 sumPrimalInfeasibilities_ = 0.0;
1977 numberDualInfeasibilities_ = 0;
1978 sumDualInfeasibilities_ = 0.0;
1979 // May be perturbed
1980 if (perturbation_ == 101 || numberChanged_) {
1981 numberChanged_ = 0; // Number of variables with changed costs
1982 perturbation_ = 102; // stop any perturbations
1983 //double changeCost;
1984 //changeBounds(1,NULL,changeCost);
1985 createRim4(false);
1986 // make sure duals are current
1987 computeDuals(givenDuals);
1988 checkDualSolution();
1989 progress_.modifyObjective(-COIN_DBL_MAX);
1990 if (numberDualInfeasibilities_) {
1991 problemStatus_ = 10; // was -3;
1992 } else {
1993 computeObjectiveValue(true);
1994 }
1995 } else if (numberPivots) {
1996 computeObjectiveValue(true);
1997 }
1998 if (numberPivots < -1000) {
1999 // objective may be wrong
2000 objectiveValue_ = innerProduct(cost_, numberColumns_ + numberRows_, solution_);
2001 objectiveValue_ += objective_->nonlinearOffset();
2002 objectiveValue_ /= (objectiveScale_ * rhsScale_);
2003 if ((specialOptions_ & 16384) == 0) {
2004 // and dual_ may be wrong (i.e. for fixed or basic)
2005 CoinIndexedVector * arrayVector = rowArray_[1];
2006 arrayVector->clear();
2007 int iRow;
2008 double * array = arrayVector->denseVector();
2009 /* Use dual_ instead of array
2010 Even though dual_ is only numberRows_ long this is
2011 okay as gets permuted to longer rowArray_[2]
2012 */
2013 arrayVector->setDenseVector(dual_);
2014 int * index = arrayVector->getIndices();
2015 int number = 0;
2016 for (iRow = 0; iRow < numberRows_; iRow++) {
2017 int iPivot = pivotVariable_[iRow];
2018 double value = cost_[iPivot];
2019 dual_[iRow] = value;
2020 if (value) {
2021 index[number++] = iRow;
2022 }
2023 }
2024 arrayVector->setNumElements(number);
2025 // Extended duals before "updateTranspose"
2026 matrix_->dualExpanded(this, arrayVector, NULL, 0);
2027 // Btran basic costs
2028 rowArray_[2]->clear();
2029 factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
2030 // and return vector
2031 arrayVector->setDenseVector(array);
2032 }
2033 }
2034 sumPrimalInfeasibilities_ = 0.0;
2035 }
2036 if ((specialOptions_&(1024 + 16384)) != 0 && !problemStatus_) {
2037 CoinIndexedVector * arrayVector = rowArray_[1];
2038 arrayVector->clear();
2039 double * rhs = arrayVector->denseVector();
2040 times(1.0, solution_, rhs);
2041 #ifdef CHECK_ACCURACY
2042 bool bad = false;
2043 #endif
2044 bool bad2 = false;
2045 int i;
2046 for ( i = 0; i < numberRows_; i++) {
2047 if (rhs[i] < rowLowerWork_[i] - primalTolerance_ ||
2048 rhs[i] > rowUpperWork_[i] + primalTolerance_) {
2049 bad2 = true;
2050 #ifdef CHECK_ACCURACY
2051 printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
2052 rowLowerWork_[i], rowUpperWork_[i],
2053 rhs[i], rowActivityWork_[i]);
2054 #endif
2055 } else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
2056 #ifdef CHECK_ACCURACY
2057 bad = true;
2058 printf("row %d correct %g bad %g\n", i, rhs[i], rowActivityWork_[i]);
2059 #endif
2060 }
2061 rhs[i] = 0.0;
2062 }
2063 for ( i = 0; i < numberColumns_; i++) {
2064 if (solution_[i] < columnLowerWork_[i] - primalTolerance_ ||
2065 solution_[i] > columnUpperWork_[i] + primalTolerance_) {
2066 bad2 = true;
2067 #ifdef CHECK_ACCURACY
2068 printf("column %d out of bounds %g, %g correct %g bad %g\n", i,
2069 columnLowerWork_[i], columnUpperWork_[i],
2070 solution_[i], columnActivityWork_[i]);
2071 #endif
2072 }
2073 }
2074 if (bad2) {
2075 problemStatus_ = -3;
2076 returnCode = -2;
2077 // Force to re-factorize early next time
2078 int numberPivots = factorization_->pivots();
2079 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2080 }
2081 }
2082 }
2083 }
2084 } else {
2085 problemStatus_ = -3;
2086 returnCode = -2;
2087 // Force to re-factorize early next time
2088 int numberPivots = factorization_->pivots();
2089 forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2090 }
2091 break;
2092 }
2093 }
2094 if (givenDuals) {
2095 CoinMemcpyN(dj_, numberRows_ + numberColumns_, givenDuals);
2096 // get rid of any values pass array
2097 delete [] candidateList;
2098 }
2099 delete [] dubiousWeights;
2100 #ifdef CLP_REPORT_PROGRESS
2101 if (ixxxxxx > ixxyyyy - 5) {
2102 int nTotal = numberColumns_ + numberRows_;
2103 double oldObj = 0.0;
2104 double newObj = 0.0;
2105 for (int i = 0; i < nTotal; i++) {
2106 if (savePSol[i])
2107 oldObj += savePSol[i] * saveCost[i];
2108 if (solution_[i])
2109 newObj += solution_[i] * cost_[i];
2110 bool printIt = false;
2111 if (cost_[i] != saveCost[i])
2112 printIt = true;
2113 if (status_[i] != saveStat[i])
2114 printIt = true;
2115 if (printIt)
2116 printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
2117 i, saveStat[i], saveCost[i], savePSol[i],
2118 status_[i], cost_[i], solution_[i]);
2119 // difference
2120 savePSol[i] = solution_[i] - savePSol[i];
2121 }
2122 printf("exit pivots %d, old obj %g new %g\n",
2123 factorization_->pivots(),
2124 oldObj, newObj);
2125 memset(saveDj, 0, numberRows_ * sizeof(double));
2126 times(1.0, savePSol, saveDj);
2127 double largest = 1.0e-6;
2128 int k = -1;
2129 for (int i = 0; i < numberRows_; i++) {
2130 saveDj[i] -= savePSol[i+numberColumns_];
2131 if (fabs(saveDj[i]) > largest) {
2132 largest = fabs(saveDj[i]);
2133 k = i;
2134 }
2135 }
2136 if (k >= 0)
2137 printf("Not null %d %g\n", k, largest);
2138 }
2139 delete [] savePSol ;
2140 delete [] saveDj ;
2141 delete [] saveCost ;
2142 delete [] saveStat ;
2143 #endif
2144 return returnCode;
2145 }
2146 /* The duals are updated by the given arrays.
2147 Returns number of infeasibilities.
2148 rowArray and columnarray will have flipped
2149 The output vector has movement (row length array) */
2150 int
updateDualsInDual(CoinIndexedVector * rowArray,CoinIndexedVector * columnArray,CoinIndexedVector * outputArray,double theta,double & objectiveChange,bool fullRecompute)2151 ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
2152 CoinIndexedVector * columnArray,
2153 CoinIndexedVector * outputArray,
2154 double theta,
2155 double & objectiveChange,
2156 bool fullRecompute)
2157 {
2158
2159 outputArray->clear();
2160
2161
2162 int numberInfeasibilities = 0;
2163 int numberRowInfeasibilities = 0;
2164
2165 // get a tolerance
2166 double tolerance = dualTolerance_;
2167 // we can't really trust infeasibilities if there is dual error
2168 double error = CoinMin(1.0e-2, largestDualError_);
2169 // allow tolerance at least slightly bigger than standard
2170 tolerance = tolerance + error;
2171
2172 double changeObj = 0.0;
2173
2174 // Coding is very similar but we can save a bit by splitting
2175 // Do rows
2176 if (!fullRecompute) {
2177 int i;
2178 double * COIN_RESTRICT reducedCost = djRegion(0);
2179 const double * COIN_RESTRICT lower = lowerRegion(0);
2180 const double * COIN_RESTRICT upper = upperRegion(0);
2181 const double * COIN_RESTRICT cost = costRegion(0);
2182 double * COIN_RESTRICT work;
2183 int number;
2184 int * COIN_RESTRICT which;
2185 const unsigned char * COIN_RESTRICT statusArray = status_ + numberColumns_;
2186 assert(rowArray->packedMode());
2187 work = rowArray->denseVector();
2188 number = rowArray->getNumElements();
2189 which = rowArray->getIndices();
2190 double multiplier[] = { -1.0, 1.0};
2191 for (i = 0; i < number; i++) {
2192 int iSequence = which[i];
2193 double alphaI = work[i];
2194 work[i] = 0.0;
2195 int iStatus = (statusArray[iSequence] & 3) - 1;
2196 if (iStatus) {
2197 double value = reducedCost[iSequence] - theta * alphaI;
2198 reducedCost[iSequence] = value;
2199 double mult = multiplier[iStatus-1];
2200 value *= mult;
2201 if (value < -tolerance) {
2202 // flipping bounds
2203 double movement = mult * (lower[iSequence] - upper[iSequence]);
2204 which[numberInfeasibilities++] = iSequence;
2205 #ifndef NDEBUG
2206 if (fabs(movement) >= 1.0e30)
2207 resetFakeBounds(-1000 - iSequence);
2208 #endif
2209 #ifdef CLP_DEBUG
2210 if ((handler_->logLevel() & 32))
2211 printf("%d %d, new dj %g, alpha %g, movement %g\n",
2212 0, iSequence, value, alphaI, movement);
2213 #endif
2214 changeObj -= movement * cost[iSequence];
2215 outputArray->quickAdd(iSequence, movement);
2216 }
2217 }
2218 }
2219 // Do columns
2220 reducedCost = djRegion(1);
2221 lower = lowerRegion(1);
2222 upper = upperRegion(1);
2223 cost = costRegion(1);
2224 // set number of infeasibilities in row array
2225 numberRowInfeasibilities = numberInfeasibilities;
2226 rowArray->setNumElements(numberInfeasibilities);
2227 numberInfeasibilities = 0;
2228 work = columnArray->denseVector();
2229 number = columnArray->getNumElements();
2230 which = columnArray->getIndices();
2231 if ((moreSpecialOptions_ & 8) != 0) {
2232 const unsigned char * COIN_RESTRICT statusArray = status_;
2233 for (i = 0; i < number; i++) {
2234 int iSequence = which[i];
2235 double alphaI = work[i];
2236 work[i] = 0.0;
2237
2238 int iStatus = (statusArray[iSequence] & 3) - 1;
2239 if (iStatus) {
2240 double value = reducedCost[iSequence] - theta * alphaI;
2241 reducedCost[iSequence] = value;
2242 double mult = multiplier[iStatus-1];
2243 value *= mult;
2244 if (value < -tolerance) {
2245 // flipping bounds
2246 double movement = mult * (upper[iSequence] - lower[iSequence]);
2247 which[numberInfeasibilities++] = iSequence;
2248 #ifndef NDEBUG
2249 if (fabs(movement) >= 1.0e30)
2250 resetFakeBounds(-1000 - iSequence);
2251 #endif
2252 #ifdef CLP_DEBUG
2253 if ((handler_->logLevel() & 32))
2254 printf("%d %d, new dj %g, alpha %g, movement %g\n",
2255 1, iSequence, value, alphaI, movement);
2256 #endif
2257 changeObj += movement * cost[iSequence];
2258 matrix_->add(this, outputArray, iSequence, movement);
2259 }
2260 }
2261 }
2262 } else {
2263 for (i = 0; i < number; i++) {
2264 int iSequence = which[i];
2265 double alphaI = work[i];
2266 work[i] = 0.0;
2267
2268 Status status = getStatus(iSequence);
2269 if (status == atLowerBound) {
2270 double value = reducedCost[iSequence] - theta * alphaI;
2271 reducedCost[iSequence] = value;
2272 double movement = 0.0;
2273
2274 if (value < -tolerance) {
2275 // to upper bound
2276 which[numberInfeasibilities++] = iSequence;
2277 movement = upper[iSequence] - lower[iSequence];
2278 #ifndef NDEBUG
2279 if (fabs(movement) >= 1.0e30)
2280 resetFakeBounds(-1000 - iSequence);
2281 #endif
2282 #ifdef CLP_DEBUG
2283 if ((handler_->logLevel() & 32))
2284 printf("%d %d, new dj %g, alpha %g, movement %g\n",
2285 1, iSequence, value, alphaI, movement);
2286 #endif
2287 changeObj += movement * cost[iSequence];
2288 matrix_->add(this, outputArray, iSequence, movement);
2289 }
2290 } else if (status == atUpperBound) {
2291 double value = reducedCost[iSequence] - theta * alphaI;
2292 reducedCost[iSequence] = value;
2293 double movement = 0.0;
2294
2295 if (value > tolerance) {
2296 // to lower bound (if swap)
2297 which[numberInfeasibilities++] = iSequence;
2298 movement = lower[iSequence] - upper[iSequence];
2299 #ifndef NDEBUG
2300 if (fabs(movement) >= 1.0e30)
2301 resetFakeBounds(-1000 - iSequence);
2302 #endif
2303 #ifdef CLP_DEBUG
2304 if ((handler_->logLevel() & 32))
2305 printf("%d %d, new dj %g, alpha %g, movement %g\n",
2306 1, iSequence, value, alphaI, movement);
2307 #endif
2308 changeObj += movement * cost[iSequence];
2309 matrix_->add(this, outputArray, iSequence, movement);
2310 }
2311 } else if (status == isFree) {
2312 double value = reducedCost[iSequence] - theta * alphaI;
2313 reducedCost[iSequence] = value;
2314 }
2315 }
2316 }
2317 } else {
2318 double * COIN_RESTRICT solution = solutionRegion(0);
2319 double * COIN_RESTRICT reducedCost = djRegion(0);
2320 const double * COIN_RESTRICT lower = lowerRegion(0);
2321 const double * COIN_RESTRICT upper = upperRegion(0);
2322 const double * COIN_RESTRICT cost = costRegion(0);
2323 int * COIN_RESTRICT which;
2324 which = rowArray->getIndices();
2325 int iSequence;
2326 for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2327 double value = reducedCost[iSequence];
2328
2329 Status status = getStatus(iSequence + numberColumns_);
2330 // more likely to be at upper bound ?
2331 if (status == atUpperBound) {
2332 double movement = 0.0;
2333 //#define NO_SWAP7
2334 if (value > tolerance) {
2335 // to lower bound (if swap)
2336 // put back alpha
2337 which[numberInfeasibilities++] = iSequence;
2338 movement = lower[iSequence] - upper[iSequence];
2339 changeObj += movement * cost[iSequence];
2340 outputArray->quickAdd(iSequence, -movement);
2341 #ifndef NDEBUG
2342 if (fabs(movement) >= 1.0e30)
2343 resetFakeBounds(-1000 - iSequence);
2344 #endif
2345 #ifndef NO_SWAP7
2346 } else if (value > -tolerance) {
2347 // at correct bound but may swap
2348 FakeBound bound = getFakeBound(iSequence + numberColumns_);
2349 if (bound == ClpSimplexDual::upperFake) {
2350 movement = lower[iSequence] - upper[iSequence];
2351 #ifndef NDEBUG
2352 if (fabs(movement) >= 1.0e30)
2353 resetFakeBounds(-1000 - iSequence);
2354 #endif
2355 setStatus(iSequence + numberColumns_, atLowerBound);
2356 solution[iSequence] = lower[iSequence];
2357 changeObj += movement * cost[iSequence];
2358 //numberFake_--;
2359 //setFakeBound(iSequence+numberColumns_,noFake);
2360 }
2361 #endif
2362 }
2363 } else if (status == atLowerBound) {
2364 double movement = 0.0;
2365
2366 if (value < -tolerance) {
2367 // to upper bound
2368 // put back alpha
2369 which[numberInfeasibilities++] = iSequence;
2370 movement = upper[iSequence] - lower[iSequence];
2371 #ifndef NDEBUG
2372 if (fabs(movement) >= 1.0e30)
2373 resetFakeBounds(-1000 - iSequence);
2374 #endif
2375 changeObj += movement * cost[iSequence];
2376 outputArray->quickAdd(iSequence, -movement);
2377 #ifndef NO_SWAP7
2378 } else if (value < tolerance) {
2379 // at correct bound but may swap
2380 FakeBound bound = getFakeBound(iSequence + numberColumns_);
2381 if (bound == ClpSimplexDual::lowerFake) {
2382 movement = upper[iSequence] - lower[iSequence];
2383 #ifndef NDEBUG
2384 if (fabs(movement) >= 1.0e30)
2385 resetFakeBounds(-1000 - iSequence);
2386 #endif
2387 setStatus(iSequence + numberColumns_, atUpperBound);
2388 solution[iSequence] = upper[iSequence];
2389 changeObj += movement * cost[iSequence];
2390 //numberFake_--;
2391 //setFakeBound(iSequence+numberColumns_,noFake);
2392 }
2393 #endif
2394 }
2395 }
2396 }
2397 // Do columns
2398 solution = solutionRegion(1);
2399 reducedCost = djRegion(1);
2400 lower = lowerRegion(1);
2401 upper = upperRegion(1);
2402 cost = costRegion(1);
2403 // set number of infeasibilities in row array
2404 numberRowInfeasibilities = numberInfeasibilities;
2405 rowArray->setNumElements(numberInfeasibilities);
2406 numberInfeasibilities = 0;
2407 which = columnArray->getIndices();
2408 for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2409 double value = reducedCost[iSequence];
2410
2411 Status status = getStatus(iSequence);
2412 if (status == atLowerBound) {
2413 double movement = 0.0;
2414
2415 if (value < -tolerance) {
2416 // to upper bound
2417 // put back alpha
2418 which[numberInfeasibilities++] = iSequence;
2419 movement = upper[iSequence] - lower[iSequence];
2420 #ifndef NDEBUG
2421 if (fabs(movement) >= 1.0e30)
2422 resetFakeBounds(-1000 - iSequence);
2423 #endif
2424 changeObj += movement * cost[iSequence];
2425 matrix_->add(this, outputArray, iSequence, movement);
2426 #ifndef NO_SWAP7
2427 } else if (value < tolerance) {
2428 // at correct bound but may swap
2429 FakeBound bound = getFakeBound(iSequence);
2430 if (bound == ClpSimplexDual::lowerFake) {
2431 movement = upper[iSequence] - lower[iSequence];
2432 #ifndef NDEBUG
2433 if (fabs(movement) >= 1.0e30)
2434 resetFakeBounds(-1000 - iSequence);
2435 #endif
2436 setStatus(iSequence, atUpperBound);
2437 solution[iSequence] = upper[iSequence];
2438 changeObj += movement * cost[iSequence];
2439 //numberFake_--;
2440 //setFakeBound(iSequence,noFake);
2441 }
2442 #endif
2443 }
2444 } else if (status == atUpperBound) {
2445 double movement = 0.0;
2446
2447 if (value > tolerance) {
2448 // to lower bound (if swap)
2449 // put back alpha
2450 which[numberInfeasibilities++] = iSequence;
2451 movement = lower[iSequence] - upper[iSequence];
2452 #ifndef NDEBUG
2453 if (fabs(movement) >= 1.0e30)
2454 resetFakeBounds(-1000 - iSequence);
2455 #endif
2456 changeObj += movement * cost[iSequence];
2457 matrix_->add(this, outputArray, iSequence, movement);
2458 #ifndef NO_SWAP7
2459 } else if (value > -tolerance) {
2460 // at correct bound but may swap
2461 FakeBound bound = getFakeBound(iSequence);
2462 if (bound == ClpSimplexDual::upperFake) {
2463 movement = lower[iSequence] - upper[iSequence];
2464 #ifndef NDEBUG
2465 if (fabs(movement) >= 1.0e30)
2466 resetFakeBounds(-1000 - iSequence);
2467 #endif
2468 setStatus(iSequence, atLowerBound);
2469 solution[iSequence] = lower[iSequence];
2470 changeObj += movement * cost[iSequence];
2471 //numberFake_--;
2472 //setFakeBound(iSequence,noFake);
2473 }
2474 #endif
2475 }
2476 }
2477 }
2478 }
2479
2480 #ifdef CLP_DEBUG
2481 if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0)
2482 printf("%d fake after full update\n", numberFake_);
2483 #endif
2484 // set number of infeasibilities
2485 columnArray->setNumElements(numberInfeasibilities);
2486 numberInfeasibilities += numberRowInfeasibilities;
2487 if (fullRecompute) {
2488 // do actual flips
2489 flipBounds(rowArray, columnArray);
2490 }
2491 objectiveChange += changeObj;
2492 return numberInfeasibilities;
2493 }
2494 void
updateDualsInValuesPass(CoinIndexedVector * rowArray,CoinIndexedVector * columnArray,double theta)2495 ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
2496 CoinIndexedVector * columnArray,
2497 double theta)
2498 {
2499
2500 // use a tighter tolerance except for all being okay
2501 double tolerance = dualTolerance_;
2502
2503 // Coding is very similar but we can save a bit by splitting
2504 // Do rows
2505 {
2506 int i;
2507 double * reducedCost = djRegion(0);
2508 double * work;
2509 int number;
2510 int * which;
2511 work = rowArray->denseVector();
2512 number = rowArray->getNumElements();
2513 which = rowArray->getIndices();
2514 for (i = 0; i < number; i++) {
2515 int iSequence = which[i];
2516 double alphaI = work[i];
2517 double value = reducedCost[iSequence] - theta * alphaI;
2518 work[i] = 0.0;
2519 reducedCost[iSequence] = value;
2520
2521 Status status = getStatus(iSequence + numberColumns_);
2522 // more likely to be at upper bound ?
2523 if (status == atUpperBound) {
2524
2525 if (value > tolerance)
2526 reducedCost[iSequence] = 0.0;
2527 } else if (status == atLowerBound) {
2528
2529 if (value < -tolerance) {
2530 reducedCost[iSequence] = 0.0;
2531 }
2532 }
2533 }
2534 }
2535 rowArray->setNumElements(0);
2536
2537 // Do columns
2538 {
2539 int i;
2540 double * reducedCost = djRegion(1);
2541 double * work;
2542 int number;
2543 int * which;
2544 work = columnArray->denseVector();
2545 number = columnArray->getNumElements();
2546 which = columnArray->getIndices();
2547
2548 for (i = 0; i < number; i++) {
2549 int iSequence = which[i];
2550 double alphaI = work[i];
2551 double value = reducedCost[iSequence] - theta * alphaI;
2552 work[i] = 0.0;
2553 reducedCost[iSequence] = value;
2554
2555 Status status = getStatus(iSequence);
2556 if (status == atLowerBound) {
2557 if (value < -tolerance)
2558 reducedCost[iSequence] = 0.0;
2559 } else if (status == atUpperBound) {
2560 if (value > tolerance)
2561 reducedCost[iSequence] = 0.0;
2562 }
2563 }
2564 }
2565 columnArray->setNumElements(0);
2566 }
2567 /*
2568 Chooses dual pivot row
2569 Would be faster with separate region to scan
2570 and will have this (with square of infeasibility) when steepest
2571 For easy problems we can just choose one of the first rows we look at
2572 */
2573 void
dualRow(int alreadyChosen)2574 ClpSimplexDual::dualRow(int alreadyChosen)
2575 {
2576 // get pivot row using whichever method it is
2577 int chosenRow = -1;
2578 #ifdef FORCE_FOLLOW
2579 bool forceThis = false;
2580 if (!fpFollow && strlen(forceFile)) {
2581 fpFollow = fopen(forceFile, "r");
2582 assert (fpFollow);
2583 }
2584 if (fpFollow) {
2585 if (numberIterations_ <= force_iteration) {
2586 // read to next Clp0102
2587 char temp[300];
2588 while (fgets(temp, 250, fpFollow)) {
2589 if (strncmp(temp, "Clp0102", 7))
2590 continue;
2591 char cin, cout;
2592 sscanf(temp + 9, "%d%*f%*s%*c%c%d%*s%*c%c%d",
2593 &force_iteration, &cin, &force_in, &cout, &force_out);
2594 if (cin == 'R')
2595 force_in += numberColumns_;
2596 if (cout == 'R')
2597 force_out += numberColumns_;
2598 forceThis = true;
2599 assert (numberIterations_ == force_iteration - 1);
2600 printf("Iteration %d will force %d out and %d in\n",
2601 force_iteration, force_out, force_in);
2602 alreadyChosen = force_out;
2603 break;
2604 }
2605 } else {
2606 // use old
2607 forceThis = true;
2608 }
2609 if (!forceThis) {
2610 fclose(fpFollow);
2611 fpFollow = NULL;
2612 forceFile = "";
2613 }
2614 }
2615 #endif
2616 double freeAlpha = 0.0;
2617 if (alreadyChosen < 0) {
2618 // first see if any free variables and put them in basis
2619 int nextFree = nextSuperBasic();
2620 //nextFree=-1; //off
2621 if (nextFree >= 0) {
2622 // unpack vector and find a good pivot
2623 unpack(rowArray_[1], nextFree);
2624 factorization_->updateColumn(rowArray_[2], rowArray_[1]);
2625
2626 double * work = rowArray_[1]->denseVector();
2627 int number = rowArray_[1]->getNumElements();
2628 int * which = rowArray_[1]->getIndices();
2629 double bestFeasibleAlpha = 0.0;
2630 int bestFeasibleRow = -1;
2631 double bestInfeasibleAlpha = 0.0;
2632 int bestInfeasibleRow = -1;
2633 int i;
2634
2635 for (i = 0; i < number; i++) {
2636 int iRow = which[i];
2637 double alpha = fabs(work[iRow]);
2638 if (alpha > 1.0e-3) {
2639 int iSequence = pivotVariable_[iRow];
2640 double value = solution_[iSequence];
2641 double lower = lower_[iSequence];
2642 double upper = upper_[iSequence];
2643 double infeasibility = 0.0;
2644 if (value > upper)
2645 infeasibility = value - upper;
2646 else if (value < lower)
2647 infeasibility = lower - value;
2648 if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
2649 if (!flagged(iSequence)) {
2650 bestInfeasibleAlpha = infeasibility * alpha;
2651 bestInfeasibleRow = iRow;
2652 }
2653 }
2654 if (alpha > bestFeasibleAlpha && (lower > -1.0e20 || upper < 1.0e20)) {
2655 bestFeasibleAlpha = alpha;
2656 bestFeasibleRow = iRow;
2657 }
2658 }
2659 }
2660 if (bestInfeasibleRow >= 0)
2661 chosenRow = bestInfeasibleRow;
2662 else if (bestFeasibleAlpha > 1.0e-2)
2663 chosenRow = bestFeasibleRow;
2664 if (chosenRow >= 0) {
2665 pivotRow_ = chosenRow;
2666 freeAlpha = work[chosenRow];
2667 }
2668 rowArray_[1]->clear();
2669 }
2670 } else {
2671 // in values pass
2672 chosenRow = alreadyChosen;
2673 #ifdef FORCE_FOLLOW
2674 if(forceThis) {
2675 alreadyChosen = -1;
2676 chosenRow = -1;
2677 for (int i = 0; i < numberRows_; i++) {
2678 if (pivotVariable_[i] == force_out) {
2679 chosenRow = i;
2680 break;
2681 }
2682 }
2683 assert (chosenRow >= 0);
2684 }
2685 #endif
2686 pivotRow_ = chosenRow;
2687 }
2688 if (chosenRow < 0)
2689 pivotRow_ = dualRowPivot_->pivotRow();
2690
2691 if (pivotRow_ >= 0) {
2692 sequenceOut_ = pivotVariable_[pivotRow_];
2693 valueOut_ = solution_[sequenceOut_];
2694 lowerOut_ = lower_[sequenceOut_];
2695 upperOut_ = upper_[sequenceOut_];
2696 if (alreadyChosen < 0) {
2697 // if we have problems we could try other way and hope we get a
2698 // zero pivot?
2699 if (valueOut_ > upperOut_) {
2700 directionOut_ = -1;
2701 dualOut_ = valueOut_ - upperOut_;
2702 } else if (valueOut_ < lowerOut_) {
2703 directionOut_ = 1;
2704 dualOut_ = lowerOut_ - valueOut_;
2705 } else {
2706 #if 1
2707 // odd (could be free) - it's feasible - go to nearest
2708 if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
2709 directionOut_ = 1;
2710 dualOut_ = lowerOut_ - valueOut_;
2711 } else {
2712 directionOut_ = -1;
2713 dualOut_ = valueOut_ - upperOut_;
2714 }
2715 #else
2716 // odd (could be free) - it's feasible - improve obj
2717 printf("direction from alpha of %g is %d\n",
2718 freeAlpha, freeAlpha > 0.0 ? 1 : -1);
2719 if (valueOut_ - lowerOut_ > 1.0e20)
2720 freeAlpha = 1.0;
2721 else if(upperOut_ - valueOut_ > 1.0e20)
2722 freeAlpha = -1.0;
2723 //if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2724 if (freeAlpha < 0.0) {
2725 directionOut_ = 1;
2726 dualOut_ = lowerOut_ - valueOut_;
2727 } else {
2728 directionOut_ = -1;
2729 dualOut_ = valueOut_ - upperOut_;
2730 }
2731 printf("direction taken %d - bounds %g %g %g\n",
2732 directionOut_, lowerOut_, valueOut_, upperOut_);
2733 #endif
2734 }
2735 #ifdef CLP_DEBUG
2736 assert(dualOut_ >= 0.0);
2737 #endif
2738 } else {
2739 // in values pass so just use sign of dj
2740 // We don't want to go through any barriers so set dualOut low
2741 // free variables will never be here
2742 dualOut_ = 1.0e-6;
2743 if (dj_[sequenceOut_] > 0.0) {
2744 // this will give a -1 in pivot row (as slacks are -1.0)
2745 directionOut_ = 1;
2746 } else {
2747 directionOut_ = -1;
2748 }
2749 }
2750 }
2751 return ;
2752 }
2753 // Checks if any fake bounds active - if so returns number and modifies
2754 // dualBound_ and everything.
2755 // Free variables will be left as free
2756 // Returns number of bounds changed if >=0
2757 // Returns -1 if not initialize and no effect
2758 // Fills in changeVector which can be used to see if unbounded
2759 // and cost of change vector
2760 int
changeBounds(int initialize,CoinIndexedVector * outputArray,double & changeCost)2761 ClpSimplexDual::changeBounds(int initialize,
2762 CoinIndexedVector * outputArray,
2763 double & changeCost)
2764 {
2765 numberFake_ = 0;
2766 if (!initialize) {
2767 int numberInfeasibilities;
2768 double newBound;
2769 newBound = 5.0 * dualBound_;
2770 numberInfeasibilities = 0;
2771 changeCost = 0.0;
2772 // put back original bounds and then check
2773 createRim1(false);
2774 int iSequence;
2775 // bounds will get bigger - just look at ones at bounds
2776 for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2777 double lowerValue = lower_[iSequence];
2778 double upperValue = upper_[iSequence];
2779 double value = solution_[iSequence];
2780 setFakeBound(iSequence, ClpSimplexDual::noFake);
2781 switch(getStatus(iSequence)) {
2782
2783 case basic:
2784 case ClpSimplex::isFixed:
2785 break;
2786 case isFree:
2787 case superBasic:
2788 break;
2789 case atUpperBound:
2790 if (fabs(value - upperValue) > primalTolerance_)
2791 numberInfeasibilities++;
2792 break;
2793 case atLowerBound:
2794 if (fabs(value - lowerValue) > primalTolerance_)
2795 numberInfeasibilities++;
2796 break;
2797 }
2798 }
2799 // If dual infeasible then carry on
2800 if (numberInfeasibilities) {
2801 handler_->message(CLP_DUAL_CHECKB, messages_)
2802 << newBound
2803 << CoinMessageEol;
2804 int iSequence;
2805 for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2806 double lowerValue = lower_[iSequence];
2807 double upperValue = upper_[iSequence];
2808 double newLowerValue;
2809 double newUpperValue;
2810 Status status = getStatus(iSequence);
2811 if (status == atUpperBound ||
2812 status == atLowerBound) {
2813 double value = solution_[iSequence];
2814 if (value - lowerValue <= upperValue - value) {
2815 newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound);
2816 newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
2817 } else {
2818 newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound);
2819 newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
2820 }
2821 lower_[iSequence] = newLowerValue;
2822 upper_[iSequence] = newUpperValue;
2823 if (newLowerValue > lowerValue) {
2824 if (newUpperValue < upperValue) {
2825 setFakeBound(iSequence, ClpSimplexDual::bothFake);
2826 #ifdef CLP_INVESTIGATE
2827 abort(); // No idea what should happen here - I have never got here
2828 #endif
2829 numberFake_++;
2830 } else {
2831 setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2832 numberFake_++;
2833 }
2834 } else {
2835 if (newUpperValue < upperValue) {
2836 setFakeBound(iSequence, ClpSimplexDual::upperFake);
2837 numberFake_++;
2838 }
2839 }
2840 if (status == atUpperBound)
2841 solution_[iSequence] = newUpperValue;
2842 else
2843 solution_[iSequence] = newLowerValue;
2844 double movement = solution_[iSequence] - value;
2845 if (movement && outputArray) {
2846 if (iSequence >= numberColumns_) {
2847 outputArray->quickAdd(iSequence, -movement);
2848 changeCost += movement * cost_[iSequence];
2849 } else {
2850 matrix_->add(this, outputArray, iSequence, movement);
2851 changeCost += movement * cost_[iSequence];
2852 }
2853 }
2854 }
2855 }
2856 dualBound_ = newBound;
2857 } else {
2858 numberInfeasibilities = -1;
2859 }
2860 return numberInfeasibilities;
2861 } else if (initialize == 1 || initialize == 3) {
2862 int iSequence;
2863 if (initialize == 3) {
2864 for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2865 setFakeBound(iSequence, ClpSimplexDual::noFake);
2866 }
2867 }
2868 double testBound = 0.999999 * dualBound_;
2869 for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2870 Status status = getStatus(iSequence);
2871 if (status == atUpperBound ||
2872 status == atLowerBound) {
2873 double lowerValue = lower_[iSequence];
2874 double upperValue = upper_[iSequence];
2875 double value = solution_[iSequence];
2876 if (lowerValue > -largeValue_ || upperValue < largeValue_) {
2877 if (true || lowerValue - value > -0.5 * dualBound_ ||
2878 upperValue - value < 0.5 * dualBound_) {
2879 if (fabs(lowerValue - value) <= fabs(upperValue - value)) {
2880 if (upperValue > lowerValue + testBound) {
2881 if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2882 numberFake_++;
2883 upper_[iSequence] = lowerValue + dualBound_;
2884 setFakeBound(iSequence, ClpSimplexDual::upperFake);
2885 }
2886 } else {
2887 if (lowerValue < upperValue - testBound) {
2888 if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2889 numberFake_++;
2890 lower_[iSequence] = upperValue - dualBound_;
2891 setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2892 }
2893 }
2894 } else {
2895 if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2896 numberFake_++;
2897 lower_[iSequence] = -0.5 * dualBound_;
2898 upper_[iSequence] = 0.5 * dualBound_;
2899 setFakeBound(iSequence, ClpSimplexDual::bothFake);
2900 abort();
2901 }
2902 if (status == atUpperBound)
2903 solution_[iSequence] = upper_[iSequence];
2904 else
2905 solution_[iSequence] = lower_[iSequence];
2906 } else {
2907 // set non basic free variables to fake bounds
2908 // I don't think we should ever get here
2909 CoinAssert(!("should not be here"));
2910 lower_[iSequence] = -0.5 * dualBound_;
2911 upper_[iSequence] = 0.5 * dualBound_;
2912 setFakeBound(iSequence, ClpSimplexDual::bothFake);
2913 numberFake_++;
2914 setStatus(iSequence, atUpperBound);
2915 solution_[iSequence] = 0.5 * dualBound_;
2916 }
2917 } else if (status == basic) {
2918 // make sure not at fake bound and bounds correct
2919 setFakeBound(iSequence, ClpSimplexDual::noFake);
2920 double gap = upper_[iSequence] - lower_[iSequence];
2921 if (gap > 0.5 * dualBound_ && gap < 2.0 * dualBound_) {
2922 if (iSequence < numberColumns_) {
2923 if (columnScale_) {
2924 double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2925 // lower
2926 double value = columnLower_[iSequence];
2927 if (value > -1.0e30) {
2928 value *= multiplier;
2929 }
2930 lower_[iSequence] = value;
2931 // upper
2932 value = columnUpper_[iSequence];
2933 if (value < 1.0e30) {
2934 value *= multiplier;
2935 }
2936 upper_[iSequence] = value;
2937 } else {
2938 lower_[iSequence] = columnLower_[iSequence];
2939 upper_[iSequence] = columnUpper_[iSequence];
2940 }
2941 } else {
2942 int iRow = iSequence - numberColumns_;
2943 if (rowScale_) {
2944 // lower
2945 double multiplier = rhsScale_ * rowScale_[iRow];
2946 double value = rowLower_[iRow];
2947 if (value > -1.0e30) {
2948 value *= multiplier;
2949 }
2950 lower_[iSequence] = value;
2951 // upper
2952 value = rowUpper_[iRow];
2953 if (value < 1.0e30) {
2954 value *= multiplier;
2955 }
2956 upper_[iSequence] = value;
2957 } else {
2958 lower_[iSequence] = rowLower_[iRow];
2959 upper_[iSequence] = rowUpper_[iRow];
2960 }
2961 }
2962 }
2963 }
2964 }
2965
2966 return 1;
2967 } else {
2968 // just reset changed ones
2969 if (columnScale_) {
2970 int iSequence;
2971 for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2972 FakeBound fakeStatus = getFakeBound(iSequence);
2973 if (fakeStatus != noFake) {
2974 if ((static_cast<int> (fakeStatus) & 1) != 0) {
2975 // lower
2976 double value = columnLower_[iSequence];
2977 if (value > -1.0e30) {
2978 double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2979 value *= multiplier;
2980 }
2981 columnLowerWork_[iSequence] = value;
2982 }
2983 if ((static_cast<int> (fakeStatus) & 2) != 0) {
2984 // upper
2985 double value = columnUpper_[iSequence];
2986 if (value < 1.0e30) {
2987 double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2988 value *= multiplier;
2989 }
2990 columnUpperWork_[iSequence] = value;
2991 }
2992 }
2993 }
2994 for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2995 FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
2996 if (fakeStatus != noFake) {
2997 if ((static_cast<int> (fakeStatus) & 1) != 0) {
2998 // lower
2999 double value = rowLower_[iSequence];
3000 if (value > -1.0e30) {
3001 double multiplier = rhsScale_ * rowScale_[iSequence];
3002 value *= multiplier;
3003 }
3004 rowLowerWork_[iSequence] = value;
3005 }
3006 if ((static_cast<int> (fakeStatus) & 2) != 0) {
3007 // upper
3008 double value = rowUpper_[iSequence];
3009 if (value < 1.0e30) {
3010 double multiplier = rhsScale_ * rowScale_[iSequence];
3011 value *= multiplier;
3012 }
3013 rowUpperWork_[iSequence] = value;
3014 }
3015 }
3016 }
3017 } else {
3018 int iSequence;
3019 for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3020 FakeBound fakeStatus = getFakeBound(iSequence);
3021 if ((static_cast<int> (fakeStatus) & 1) != 0) {
3022 // lower
3023 columnLowerWork_[iSequence] = columnLower_[iSequence];
3024 }
3025 if ((static_cast<int> (fakeStatus) & 2) != 0) {
3026 // upper
3027 columnUpperWork_[iSequence] = columnUpper_[iSequence];
3028 }
3029 }
3030 for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3031 FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3032 if ((static_cast<int> (fakeStatus) & 1) != 0) {
3033 // lower
3034 rowLowerWork_[iSequence] = rowLower_[iSequence];
3035 }
3036 if ((static_cast<int> (fakeStatus) & 2) != 0) {
3037 // upper
3038 rowUpperWork_[iSequence] = rowUpper_[iSequence];
3039 }
3040 }
3041 }
3042 return 0;
3043 }
3044 }
3045 int
dualColumn0(const CoinIndexedVector * rowArray,const CoinIndexedVector * columnArray,CoinIndexedVector * spareArray,double acceptablePivot,double & upperReturn,double & bestReturn,double & badFree)3046 ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
3047 const CoinIndexedVector * columnArray,
3048 CoinIndexedVector * spareArray,
3049 double acceptablePivot,
3050 double & upperReturn, double &bestReturn, double & badFree)
3051 {
3052 // do first pass to get possibles
3053 double * spare = spareArray->denseVector();
3054 int * index = spareArray->getIndices();
3055 const double * work;
3056 int number;
3057 const int * which;
3058 const double * reducedCost;
3059 // We can also see if infeasible or pivoting on free
3060 double tentativeTheta = 1.0e15;
3061 double upperTheta = 1.0e31;
3062 double freePivot = acceptablePivot;
3063 double bestPossible = 0.0;
3064 int numberRemaining = 0;
3065 int i;
3066 badFree = 0.0;
3067 if ((moreSpecialOptions_ & 8) != 0) {
3068 // No free or super basic
3069 double multiplier[] = { -1.0, 1.0};
3070 double dualT = - dualTolerance_;
3071 for (int iSection = 0; iSection < 2; iSection++) {
3072
3073 int addSequence;
3074 unsigned char * statusArray;
3075 if (!iSection) {
3076 work = rowArray->denseVector();
3077 number = rowArray->getNumElements();
3078 which = rowArray->getIndices();
3079 reducedCost = rowReducedCost_;
3080 addSequence = numberColumns_;
3081 statusArray = status_ + numberColumns_;
3082 } else {
3083 work = columnArray->denseVector();
3084 number = columnArray->getNumElements();
3085 which = columnArray->getIndices();
3086 reducedCost = reducedCostWork_;
3087 addSequence = 0;
3088 statusArray = status_;
3089 }
3090
3091 for (i = 0; i < number; i++) {
3092 int iSequence = which[i];
3093 double alpha;
3094 double oldValue;
3095 double value;
3096
3097 assert (getStatus(iSequence + addSequence) != isFree
3098 && getStatus(iSequence + addSequence) != superBasic);
3099 int iStatus = (statusArray[iSequence] & 3) - 1;
3100 if (iStatus) {
3101 double mult = multiplier[iStatus-1];
3102 alpha = work[i] * mult;
3103 if (alpha > 0.0) {
3104 oldValue = reducedCost[iSequence] * mult;
3105 value = oldValue - tentativeTheta * alpha;
3106 if (value < dualT) {
3107 bestPossible = CoinMax(bestPossible, alpha);
3108 value = oldValue - upperTheta * alpha;
3109 if (value < dualT && alpha >= acceptablePivot) {
3110 upperTheta = (oldValue - dualT) / alpha;
3111 //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3112 }
3113 // add to list
3114 spare[numberRemaining] = alpha * mult;
3115 index[numberRemaining++] = iSequence + addSequence;
3116 }
3117 }
3118 }
3119 }
3120 }
3121 } else {
3122 // some free or super basic
3123 for (int iSection = 0; iSection < 2; iSection++) {
3124
3125 int addSequence;
3126
3127 if (!iSection) {
3128 work = rowArray->denseVector();
3129 number = rowArray->getNumElements();
3130 which = rowArray->getIndices();
3131 reducedCost = rowReducedCost_;
3132 addSequence = numberColumns_;
3133 } else {
3134 work = columnArray->denseVector();
3135 number = columnArray->getNumElements();
3136 which = columnArray->getIndices();
3137 reducedCost = reducedCostWork_;
3138 addSequence = 0;
3139 }
3140
3141 for (i = 0; i < number; i++) {
3142 int iSequence = which[i];
3143 double alpha;
3144 double oldValue;
3145 double value;
3146 bool keep;
3147
3148 switch(getStatus(iSequence + addSequence)) {
3149
3150 case basic:
3151 case ClpSimplex::isFixed:
3152 break;
3153 case isFree:
3154 case superBasic:
3155 alpha = work[i];
3156 bestPossible = CoinMax(bestPossible, fabs(alpha));
3157 oldValue = reducedCost[iSequence];
3158 // If free has to be very large - should come in via dualRow
3159 //if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
3160 //break;
3161 if (oldValue > dualTolerance_) {
3162 keep = true;
3163 } else if (oldValue < -dualTolerance_) {
3164 keep = true;
3165 } else {
3166 if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
3167 keep = true;
3168 } else {
3169 keep = false;
3170 badFree = CoinMax(badFree, fabs(alpha));
3171 }
3172 }
3173 if (keep) {
3174 // free - choose largest
3175 if (fabs(alpha) > freePivot) {
3176 freePivot = fabs(alpha);
3177 sequenceIn_ = iSequence + addSequence;
3178 theta_ = oldValue / alpha;
3179 alpha_ = alpha;
3180 }
3181 }
3182 break;
3183 case atUpperBound:
3184 alpha = work[i];
3185 oldValue = reducedCost[iSequence];
3186 value = oldValue - tentativeTheta * alpha;
3187 //assert (oldValue<=dualTolerance_*1.0001);
3188 if (value > dualTolerance_) {
3189 bestPossible = CoinMax(bestPossible, -alpha);
3190 value = oldValue - upperTheta * alpha;
3191 if (value > dualTolerance_ && -alpha >= acceptablePivot) {
3192 upperTheta = (oldValue - dualTolerance_) / alpha;
3193 //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3194 }
3195 // add to list
3196 spare[numberRemaining] = alpha;
3197 index[numberRemaining++] = iSequence + addSequence;
3198 }
3199 break;
3200 case atLowerBound:
3201 alpha = work[i];
3202 oldValue = reducedCost[iSequence];
3203 value = oldValue - tentativeTheta * alpha;
3204 //assert (oldValue>=-dualTolerance_*1.0001);
3205 if (value < -dualTolerance_) {
3206 bestPossible = CoinMax(bestPossible, alpha);
3207 value = oldValue - upperTheta * alpha;
3208 if (value < -dualTolerance_ && alpha >= acceptablePivot) {
3209 upperTheta = (oldValue + dualTolerance_) / alpha;
3210 //tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3211 }
3212 // add to list
3213 spare[numberRemaining] = alpha;
3214 index[numberRemaining++] = iSequence + addSequence;
3215 }
3216 break;
3217 }
3218 }
3219 }
3220 }
3221 upperReturn = upperTheta;
3222 bestReturn = bestPossible;
3223 return numberRemaining;
3224 }
3225 /*
3226 Row array has row part of pivot row (as duals so sign may be switched)
3227 Column array has column part.
3228 This chooses pivot column.
3229 Spare array will be needed when we start getting clever.
3230 We will check for basic so spare array will never overflow.
3231 If necessary will modify costs
3232 */
3233 double
dualColumn(CoinIndexedVector * rowArray,CoinIndexedVector * columnArray,CoinIndexedVector * spareArray,CoinIndexedVector * spareArray2,double acceptablePivot,CoinBigIndex *)3234 ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
3235 CoinIndexedVector * columnArray,
3236 CoinIndexedVector * spareArray,
3237 CoinIndexedVector * spareArray2,
3238 double acceptablePivot,
3239 CoinBigIndex * /*dubiousWeights*/)
3240 {
3241 int numberPossiblySwapped = 0;
3242 int numberRemaining = 0;
3243
3244 double totalThru = 0.0; // for when variables flip
3245 //double saveAcceptable=acceptablePivot;
3246 //acceptablePivot=1.0e-9;
3247
3248 double bestEverPivot = acceptablePivot;
3249 int lastSequence = -1;
3250 double lastPivot = 0.0;
3251 double upperTheta;
3252 double newTolerance = dualTolerance_;
3253 //newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
3254 // will we need to increase tolerance
3255 bool thisIncrease = false;
3256 // If we think we need to modify costs (not if something from broad sweep)
3257 bool modifyCosts = false;
3258 // Increase in objective due to swapping bounds (may be negative)
3259 double increaseInObjective = 0.0;
3260
3261 // use spareArrays to put ones looked at in
3262 // we are going to flip flop between
3263 int iFlip = 0;
3264 // Possible list of pivots
3265 int interesting[2];
3266 // where possible swapped ones are
3267 int swapped[2];
3268 // for zeroing out arrays after
3269 int marker[2][2];
3270 // pivot elements
3271 double * array[2], * spare, * spare2;
3272 // indices
3273 int * indices[2], * index, * index2;
3274 spareArray2->clear();
3275 array[0] = spareArray->denseVector();
3276 indices[0] = spareArray->getIndices();
3277 spare = array[0];
3278 index = indices[0];
3279 array[1] = spareArray2->denseVector();
3280 indices[1] = spareArray2->getIndices();
3281 int i;
3282
3283 // initialize lists
3284 for (i = 0; i < 2; i++) {
3285 interesting[i] = 0;
3286 swapped[i] = numberColumns_;
3287 marker[i][0] = 0;
3288 marker[i][1] = numberColumns_;
3289 }
3290 /*
3291 First we get a list of possible pivots. We can also see if the
3292 problem looks infeasible or whether we want to pivot in free variable.
3293 This may make objective go backwards but can only happen a finite
3294 number of times and I do want free variables basic.
3295
3296 Then we flip back and forth. At the start of each iteration
3297 interesting[iFlip] should have possible candidates and swapped[iFlip]
3298 will have pivots if we decide to take a previous pivot.
3299 At end of each iteration interesting[1-iFlip] should have
3300 candidates if we go through this theta and swapped[1-iFlip]
3301 pivots if we don't go through.
3302
3303 At first we increase theta and see what happens. We start
3304 theta at a reasonable guess. If in right area then we do bit by bit.
3305
3306 */
3307
3308 // do first pass to get possibles
3309 upperTheta = 1.0e31;
3310 double bestPossible = 0.0;
3311 double badFree = 0.0;
3312 alpha_ = 0.0;
3313 if (spareIntArray_[0] >= 0) {
3314 numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3315 acceptablePivot, upperTheta, bestPossible, badFree);
3316 } else {
3317 // already done
3318 numberRemaining = spareArray->getNumElements();
3319 spareArray->setNumElements(0);
3320 upperTheta = spareDoubleArray_[0];
3321 bestPossible = spareDoubleArray_[1];
3322 if (spareIntArray_[0] == -1) {
3323 theta_ = spareDoubleArray_[2];
3324 alpha_ = spareDoubleArray_[3];
3325 sequenceIn_ = spareIntArray_[1];
3326 } else {
3327 #if 0
3328 int n = numberRemaining;
3329 double u = upperTheta;
3330 double b = bestPossible;
3331 upperTheta = 1.0e31;
3332 bestPossible = 0.0;
3333 numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3334 acceptablePivot, upperTheta, bestPossible, badFree);
3335 assert (n == numberRemaining);
3336 assert (fabs(b - bestPossible) < 1.0e-7);
3337 assert (fabs(u - upperTheta) < 1.0e-7);
3338 #endif
3339 }
3340 }
3341 // switch off
3342 spareIntArray_[0] = 0;
3343 // We can also see if infeasible or pivoting on free
3344 double tentativeTheta = 1.0e25;
3345 interesting[0] = numberRemaining;
3346 marker[0][0] = numberRemaining;
3347
3348 if (!numberRemaining && sequenceIn_ < 0)
3349 return 0.0; // Looks infeasible
3350
3351 // If sum of bad small pivots too much
3352 #define MORE_CAREFUL
3353 #ifdef MORE_CAREFUL
3354 bool badSumPivots = false;
3355 #endif
3356 if (sequenceIn_ >= 0) {
3357 // free variable - always choose
3358 } else {
3359
3360 theta_ = 1.0e50;
3361 // now flip flop between spare arrays until reasonable theta
3362 tentativeTheta = CoinMax(10.0 * upperTheta, 1.0e-7);
3363
3364 // loops increasing tentative theta until can't go through
3365
3366 while (tentativeTheta < 1.0e22) {
3367 double thruThis = 0.0;
3368
3369 double bestPivot = acceptablePivot;
3370 int bestSequence = -1;
3371
3372 numberPossiblySwapped = numberColumns_;
3373 numberRemaining = 0;
3374
3375 upperTheta = 1.0e50;
3376
3377 spare = array[iFlip];
3378 index = indices[iFlip];
3379 spare2 = array[1-iFlip];
3380 index2 = indices[1-iFlip];
3381
3382 // try 3 different ways
3383 // 1 bias increase by ones with slightly wrong djs
3384 // 2 bias by all
3385 // 3 bias by all - tolerance
3386 #define TRYBIAS 3
3387
3388
3389 double increaseInThis = 0.0; //objective increase in this loop
3390
3391 for (i = 0; i < interesting[iFlip]; i++) {
3392 int iSequence = index[i];
3393 double alpha = spare[i];
3394 double oldValue = dj_[iSequence];
3395 double value = oldValue - tentativeTheta * alpha;
3396
3397 if (alpha < 0.0) {
3398 //at upper bound
3399 if (value > newTolerance) {
3400 double range = upper_[iSequence] - lower_[iSequence];
3401 thruThis -= range * alpha;
3402 #if TRYBIAS==1
3403 if (oldValue > 0.0)
3404 increaseInThis -= oldValue * range;
3405 #elif TRYBIAS==2
3406 increaseInThis -= oldValue * range;
3407 #else
3408 increaseInThis -= (oldValue + dualTolerance_) * range;
3409 #endif
3410 // goes on swapped list (also means candidates if too many)
3411 spare2[--numberPossiblySwapped] = alpha;
3412 index2[numberPossiblySwapped] = iSequence;
3413 if (fabs(alpha) > bestPivot) {
3414 bestPivot = fabs(alpha);
3415 bestSequence = numberPossiblySwapped;
3416 }
3417 } else {
3418 value = oldValue - upperTheta * alpha;
3419 if (value > newTolerance && -alpha >= acceptablePivot)
3420 upperTheta = (oldValue - newTolerance) / alpha;
3421 spare2[numberRemaining] = alpha;
3422 index2[numberRemaining++] = iSequence;
3423 }
3424 } else {
3425 // at lower bound
3426 if (value < -newTolerance) {
3427 double range = upper_[iSequence] - lower_[iSequence];
3428 thruThis += range * alpha;
3429 //?? is this correct - and should we look at good ones
3430 #if TRYBIAS==1
3431 if (oldValue < 0.0)
3432 increaseInThis += oldValue * range;
3433 #elif TRYBIAS==2
3434 increaseInThis += oldValue * range;
3435 #else
3436 increaseInThis += (oldValue - dualTolerance_) * range;
3437 #endif
3438 // goes on swapped list (also means candidates if too many)
3439 spare2[--numberPossiblySwapped] = alpha;
3440 index2[numberPossiblySwapped] = iSequence;
3441 if (fabs(alpha) > bestPivot) {
3442 bestPivot = fabs(alpha);
3443 bestSequence = numberPossiblySwapped;
3444 }
3445 } else {
3446 value = oldValue - upperTheta * alpha;
3447 if (value < -newTolerance && alpha >= acceptablePivot)
3448 upperTheta = (oldValue + newTolerance) / alpha;
3449 spare2[numberRemaining] = alpha;
3450 index2[numberRemaining++] = iSequence;
3451 }
3452 }
3453 }
3454 swapped[1-iFlip] = numberPossiblySwapped;
3455 interesting[1-iFlip] = numberRemaining;
3456 marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3457 marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3458
3459 if (totalThru + thruThis >= fabs(dualOut_) ||
3460 increaseInObjective + increaseInThis < 0.0) {
3461 // We should be pivoting in this batch
3462 // so compress down to this lot
3463 numberRemaining = 0;
3464 for (i = numberColumns_ - 1; i >= swapped[1-iFlip]; i--) {
3465 spare[numberRemaining] = spare2[i];
3466 index[numberRemaining++] = index2[i];
3467 }
3468 interesting[iFlip] = numberRemaining;
3469 int iTry;
3470 #define MAXTRY 100
3471 // first get ratio with tolerance
3472 for (iTry = 0; iTry < MAXTRY; iTry++) {
3473
3474 upperTheta = 1.0e50;
3475 numberPossiblySwapped = numberColumns_;
3476 numberRemaining = 0;
3477
3478 increaseInThis = 0.0; //objective increase in this loop
3479
3480 thruThis = 0.0;
3481
3482 spare = array[iFlip];
3483 index = indices[iFlip];
3484 spare2 = array[1-iFlip];
3485 index2 = indices[1-iFlip];
3486 for (i = 0; i < interesting[iFlip]; i++) {
3487 int iSequence = index[i];
3488 double alpha = spare[i];
3489 double oldValue = dj_[iSequence];
3490 double value = oldValue - upperTheta * alpha;
3491
3492 if (alpha < 0.0) {
3493 //at upper bound
3494 if (value > newTolerance) {
3495 if (-alpha >= acceptablePivot) {
3496 upperTheta = (oldValue - newTolerance) / alpha;
3497 }
3498 }
3499 } else {
3500 // at lower bound
3501 if (value < -newTolerance) {
3502 if (alpha >= acceptablePivot) {
3503 upperTheta = (oldValue + newTolerance) / alpha;
3504 }
3505 }
3506 }
3507 }
3508 bestPivot = acceptablePivot;
3509 sequenceIn_ = -1;
3510 #ifdef DUBIOUS_WEIGHTS
3511 double bestWeight = COIN_DBL_MAX;
3512 #endif
3513 double largestPivot = acceptablePivot;
3514 // now choose largest and sum all ones which will go through
3515 //printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
3516 // Sum of bad small pivots
3517 #ifdef MORE_CAREFUL
3518 double sumBadPivots = 0.0;
3519 badSumPivots = false;
3520 #endif
3521 // Make sure upperTheta will work (-O2 and above gives problems)
3522 upperTheta *= 1.0000000001;
3523 for (i = 0; i < interesting[iFlip]; i++) {
3524 int iSequence = index[i];
3525 double alpha = spare[i];
3526 double value = dj_[iSequence] - upperTheta * alpha;
3527 double badDj = 0.0;
3528
3529 bool addToSwapped = false;
3530
3531 if (alpha < 0.0) {
3532 //at upper bound
3533 if (value >= 0.0) {
3534 addToSwapped = true;
3535 #if TRYBIAS==1
3536 badDj = -CoinMax(dj_[iSequence], 0.0);
3537 #elif TRYBIAS==2
3538 badDj = -dj_[iSequence];
3539 #else
3540 badDj = -dj_[iSequence] - dualTolerance_;
3541 #endif
3542 }
3543 } else {
3544 // at lower bound
3545 if (value <= 0.0) {
3546 addToSwapped = true;
3547 #if TRYBIAS==1
3548 badDj = CoinMin(dj_[iSequence], 0.0);
3549 #elif TRYBIAS==2
3550 badDj = dj_[iSequence];
3551 #else
3552 badDj = dj_[iSequence] - dualTolerance_;
3553 #endif
3554 }
3555 }
3556 if (!addToSwapped) {
3557 // add to list of remaining
3558 spare2[numberRemaining] = alpha;
3559 index2[numberRemaining++] = iSequence;
3560 } else {
3561 // add to list of swapped
3562 spare2[--numberPossiblySwapped] = alpha;
3563 index2[numberPossiblySwapped] = iSequence;
3564 // select if largest pivot
3565 bool take = false;
3566 double absAlpha = fabs(alpha);
3567 #ifdef DUBIOUS_WEIGHTS
3568 // User could do anything to break ties here
3569 double weight;
3570 if (dubiousWeights)
3571 weight = dubiousWeights[iSequence];
3572 else
3573 weight = 1.0;
3574 weight += randomNumberGenerator_.randomDouble() * 1.0e-2;
3575 if (absAlpha > 2.0 * bestPivot) {
3576 take = true;
3577 } else if (absAlpha > largestPivot) {
3578 // could multiply absAlpha and weight
3579 if (weight * bestPivot < bestWeight * absAlpha)
3580 take = true;
3581 }
3582 #else
3583 if (absAlpha > bestPivot)
3584 take = true;
3585 #endif
3586 #ifdef MORE_CAREFUL
3587 if (absAlpha < acceptablePivot && upperTheta < 1.0e20) {
3588 if (alpha < 0.0) {
3589 //at upper bound
3590 if (value > dualTolerance_) {
3591 double gap = upper_[iSequence] - lower_[iSequence];
3592 if (gap < 1.0e20)
3593 sumBadPivots += value * gap;
3594 else
3595 sumBadPivots += 1.0e20;
3596 //printf("bad %d alpha %g dj at upper %g\n",
3597 // iSequence,alpha,value);
3598 }
3599 } else {
3600 //at lower bound
3601 if (value < -dualTolerance_) {
3602 double gap = upper_[iSequence] - lower_[iSequence];
3603 if (gap < 1.0e20)
3604 sumBadPivots -= value * gap;
3605 else
3606 sumBadPivots += 1.0e20;
3607 //printf("bad %d alpha %g dj at lower %g\n",
3608 // iSequence,alpha,value);
3609 }
3610 }
3611 }
3612 #endif
3613 #ifdef FORCE_FOLLOW
3614 if (iSequence == force_in) {
3615 printf("taking %d - alpha %g best %g\n", force_in, absAlpha, largestPivot);
3616 take = true;
3617 }
3618 #endif
3619 if (take) {
3620 sequenceIn_ = numberPossiblySwapped;
3621 bestPivot = absAlpha;
3622 theta_ = dj_[iSequence] / alpha;
3623 largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
3624 #ifdef DUBIOUS_WEIGHTS
3625 bestWeight = weight;
3626 #endif
3627 //printf(" taken seq %d alpha %g weight %d\n",
3628 // iSequence,absAlpha,dubiousWeights[iSequence]);
3629 } else {
3630 //printf(" not taken seq %d alpha %g weight %d\n",
3631 // iSequence,absAlpha,dubiousWeights[iSequence]);
3632 }
3633 double range = upper_[iSequence] - lower_[iSequence];
3634 thruThis += range * fabs(alpha);
3635 increaseInThis += badDj * range;
3636 }
3637 }
3638 marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3639 marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3640 #ifdef MORE_CAREFUL
3641 // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3642 if (sumBadPivots > 1.0e4) {
3643 if (handler_->logLevel() > 1)
3644 printf("maybe forcing re-factorization - sum %g %d pivots\n", sumBadPivots,
3645 factorization_->pivots());
3646 if(factorization_->pivots() > 3) {
3647 badSumPivots = true;
3648 break;
3649 }
3650 }
3651 #endif
3652 swapped[1-iFlip] = numberPossiblySwapped;
3653 interesting[1-iFlip] = numberRemaining;
3654 // If we stop now this will be increase in objective (I think)
3655 double increase = (fabs(dualOut_) - totalThru) * theta_;
3656 increase += increaseInObjective;
3657 if (theta_ < 0.0)
3658 thruThis += fabs(dualOut_); // force using this one
3659 if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
3660 // back
3661 // We may need to be more careful - we could do by
3662 // switch so we always do fine grained?
3663 bestPivot = 0.0;
3664 } else {
3665 // add in
3666 totalThru += thruThis;
3667 increaseInObjective += increaseInThis;
3668 }
3669 if (bestPivot < 0.1 * bestEverPivot &&
3670 bestEverPivot > 1.0e-6 &&
3671 (bestPivot < 1.0e-3 || totalThru * 2.0 > fabs(dualOut_))) {
3672 // back to previous one
3673 sequenceIn_ = lastSequence;
3674 // swap regions
3675 iFlip = 1 - iFlip;
3676 break;
3677 } else if (sequenceIn_ == -1 && upperTheta > largeValue_) {
3678 if (lastPivot > acceptablePivot) {
3679 // back to previous one
3680 sequenceIn_ = lastSequence;
3681 // swap regions
3682 iFlip = 1 - iFlip;
3683 } else {
3684 // can only get here if all pivots too small
3685 }
3686 break;
3687 } else if (totalThru >= fabs(dualOut_)) {
3688 modifyCosts = true; // fine grain - we can modify costs
3689 break; // no point trying another loop
3690 } else {
3691 lastSequence = sequenceIn_;
3692 if (bestPivot > bestEverPivot)
3693 bestEverPivot = bestPivot;
3694 iFlip = 1 - iFlip;
3695 modifyCosts = true; // fine grain - we can modify costs
3696 }
3697 }
3698 if (iTry == MAXTRY)
3699 iFlip = 1 - iFlip; // flip back
3700 break;
3701 } else {
3702 // skip this lot
3703 if (bestPivot > 1.0e-3 || bestPivot > bestEverPivot) {
3704 bestEverPivot = bestPivot;
3705 lastSequence = bestSequence;
3706 } else {
3707 // keep old swapped
3708 CoinMemcpyN(array[iFlip] + swapped[iFlip],
3709 numberColumns_ - swapped[iFlip], array[1-iFlip] + swapped[iFlip]);
3710 CoinMemcpyN(indices[iFlip] + swapped[iFlip],
3711 numberColumns_ - swapped[iFlip], indices[1-iFlip] + swapped[iFlip]);
3712 marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], swapped[iFlip]);
3713 swapped[1-iFlip] = swapped[iFlip];
3714 }
3715 increaseInObjective += increaseInThis;
3716 iFlip = 1 - iFlip; // swap regions
3717 tentativeTheta = 2.0 * upperTheta;
3718 totalThru += thruThis;
3719 }
3720 }
3721
3722 // can get here without sequenceIn_ set but with lastSequence
3723 if (sequenceIn_ < 0 && lastSequence >= 0) {
3724 // back to previous one
3725 sequenceIn_ = lastSequence;
3726 // swap regions
3727 iFlip = 1 - iFlip;
3728 }
3729
3730 #define MINIMUMTHETA 1.0e-18
3731 // Movement should be minimum for anti-degeneracy - unless
3732 // fixed variable out
3733 double minimumTheta;
3734 if (upperOut_ > lowerOut_)
3735 minimumTheta = MINIMUMTHETA;
3736 else
3737 minimumTheta = 0.0;
3738 if (sequenceIn_ >= 0) {
3739 // at this stage sequenceIn_ is just pointer into index array
3740 // flip just so we can use iFlip
3741 iFlip = 1 - iFlip;
3742 spare = array[iFlip];
3743 index = indices[iFlip];
3744 double oldValue;
3745 alpha_ = spare[sequenceIn_];
3746 sequenceIn_ = indices[iFlip][sequenceIn_];
3747 oldValue = dj_[sequenceIn_];
3748 theta_ = CoinMax(oldValue / alpha_, 0.0);
3749 if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
3750 // can't pivot to zero
3751 #if 0
3752 if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) {
3753 theta_ = minimumTheta;
3754 } else if (oldValue - minimumTheta*alpha_ >= -newTolerance) {
3755 theta_ = minimumTheta;
3756 thisIncrease = true;
3757 } else {
3758 theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0);
3759 thisIncrease = true;
3760 }
3761 #else
3762 theta_ = minimumTheta;
3763 #endif
3764 }
3765 // may need to adjust costs so all dual feasible AND pivoted is exactly 0
3766 //int costOffset = numberRows_+numberColumns_;
3767 if (modifyCosts && !badSumPivots) {
3768 int i;
3769 for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) {
3770 int iSequence = index[i];
3771 double alpha = spare[i];
3772 double value = dj_[iSequence] - theta_ * alpha;
3773
3774 // can't be free here
3775
3776 if (alpha < 0.0) {
3777 //at upper bound
3778 if (value > dualTolerance_) {
3779 thisIncrease = true;
3780 #define MODIFYCOST 2
3781 #if MODIFYCOST
3782 // modify cost to hit new tolerance
3783 double modification = alpha * theta_ - dj_[iSequence]
3784 + newTolerance;
3785 if ((specialOptions_&(2048 + 4096 + 16384)) != 0) {
3786 if ((specialOptions_ & 16384) != 0) {
3787 if (fabs(modification) < 1.0e-8)
3788 modification = 0.0;
3789 } else if ((specialOptions_ & 2048) != 0) {
3790 if (fabs(modification) < 1.0e-10)
3791 modification = 0.0;
3792 } else {
3793 if (fabs(modification) < 1.0e-12)
3794 modification = 0.0;
3795 }
3796 }
3797 dj_[iSequence] += modification;
3798 cost_[iSequence] += modification;
3799 if (modification)
3800 numberChanged_ ++; // Say changed costs
3801 //cost_[iSequence+costOffset] += modification; // save change
3802 #endif
3803 }
3804 } else {
3805 // at lower bound
3806 if (-value > dualTolerance_) {
3807 thisIncrease = true;
3808 #if MODIFYCOST
3809 // modify cost to hit new tolerance
3810 double modification = alpha * theta_ - dj_[iSequence]
3811 - newTolerance;
3812 //modification = CoinMax(modification,-dualTolerance_);
3813 //assert (fabs(modification)<1.0e-7);
3814 if ((specialOptions_&(2048 + 4096)) != 0) {
3815 if ((specialOptions_ & 2048) != 0) {
3816 if (fabs(modification) < 1.0e-10)
3817 modification = 0.0;
3818 } else {
3819 if (fabs(modification) < 1.0e-12)
3820 modification = 0.0;
3821 }
3822 }
3823 dj_[iSequence] += modification;
3824 cost_[iSequence] += modification;
3825 if (modification)
3826 numberChanged_ ++; // Say changed costs
3827 //cost_[iSequence+costOffset] += modification; // save change
3828 #endif
3829 }
3830 }
3831 }
3832 }
3833 }
3834 }
3835
3836 #ifdef MORE_CAREFUL
3837 // If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3838 if ((badSumPivots ||
3839 fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) {
3840 if (handler_->logLevel() > 1)
3841 printf("forcing re-factorization\n");
3842 sequenceIn_ = -1;
3843 }
3844 #endif
3845 if (sequenceIn_ >= 0) {
3846 lowerIn_ = lower_[sequenceIn_];
3847 upperIn_ = upper_[sequenceIn_];
3848 valueIn_ = solution_[sequenceIn_];
3849 dualIn_ = dj_[sequenceIn_];
3850
3851 if (numberTimesOptimal_) {
3852 // can we adjust cost back closer to original
3853 //*** add coding
3854 }
3855 #if MODIFYCOST>1
3856 // modify cost to hit zero exactly
3857 // so (dualIn_+modification)==theta_*alpha_
3858 double modification = theta_ * alpha_ - dualIn_;
3859 // But should not move objective too much ??
3860 #define DONT_MOVE_OBJECTIVE
3861 #ifdef DONT_MOVE_OBJECTIVE
3862 double moveObjective = fabs(modification * solution_[sequenceIn_]);
3863 double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3);
3864 if (moveObjective > smallMove) {
3865 if (handler_->logLevel() > 1)
3866 printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
3867 modification, solution_[sequenceIn_]);
3868 modification *= smallMove / moveObjective;
3869 }
3870 #endif
3871 if (badSumPivots)
3872 modification = 0.0;
3873 if ((specialOptions_&(2048 + 4096)) != 0) {
3874 if ((specialOptions_ & 16384) != 0) {
3875 // in fast dual
3876 if (fabs(modification) < 1.0e-7)
3877 modification = 0.0;
3878 } else if ((specialOptions_ & 2048) != 0) {
3879 if (fabs(modification) < 1.0e-10)
3880 modification = 0.0;
3881 } else {
3882 if (fabs(modification) < 1.0e-12)
3883 modification = 0.0;
3884 }
3885 }
3886 dualIn_ += modification;
3887 dj_[sequenceIn_] = dualIn_;
3888 cost_[sequenceIn_] += modification;
3889 if (modification)
3890 numberChanged_ ++; // Say changed costs
3891 //int costOffset = numberRows_+numberColumns_;
3892 //cost_[sequenceIn_+costOffset] += modification; // save change
3893 //assert (fabs(modification)<1.0e-6);
3894 #ifdef CLP_DEBUG
3895 if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
3896 printf("exact %d new cost %g, change %g\n", sequenceIn_,
3897 cost_[sequenceIn_], modification);
3898 #endif
3899 #endif
3900
3901 if (alpha_ < 0.0) {
3902 // as if from upper bound
3903 directionIn_ = -1;
3904 upperIn_ = valueIn_;
3905 } else {
3906 // as if from lower bound
3907 directionIn_ = 1;
3908 lowerIn_ = valueIn_;
3909 }
3910 } else {
3911 // no pivot
3912 bestPossible = 0.0;
3913 alpha_ = 0.0;
3914 }
3915 //if (thisIncrease)
3916 //dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3917
3918 // clear arrays
3919
3920 for (i = 0; i < 2; i++) {
3921 CoinZeroN(array[i], marker[i][0]);
3922 CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]);
3923 }
3924 return bestPossible;
3925 }
3926 #ifdef CLP_ALL_ONE_FILE
3927 #undef MAXTRY
3928 #endif
3929 /* Checks if tentative optimal actually means unbounded
3930 Returns -3 if not, 2 if is unbounded */
3931 int
checkUnbounded(CoinIndexedVector * ray,CoinIndexedVector * spare,double changeCost)3932 ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3933 CoinIndexedVector * spare,
3934 double changeCost)
3935 {
3936 int status = 2; // say unbounded
3937 factorization_->updateColumn(spare, ray);
3938 // get reduced cost
3939 int i;
3940 int number = ray->getNumElements();
3941 int * index = ray->getIndices();
3942 double * array = ray->denseVector();
3943 for (i = 0; i < number; i++) {
3944 int iRow = index[i];
3945 int iPivot = pivotVariable_[iRow];
3946 changeCost -= cost(iPivot) * array[iRow];
3947 }
3948 double way;
3949 if (changeCost > 0.0) {
3950 //try going down
3951 way = 1.0;
3952 } else if (changeCost < 0.0) {
3953 //try going up
3954 way = -1.0;
3955 } else {
3956 #ifdef CLP_DEBUG
3957 printf("can't decide on up or down\n");
3958 #endif
3959 way = 0.0;
3960 status = -3;
3961 }
3962 double movement = 1.0e10 * way; // some largish number
3963 double zeroTolerance = 1.0e-14 * dualBound_;
3964 for (i = 0; i < number; i++) {
3965 int iRow = index[i];
3966 int iPivot = pivotVariable_[iRow];
3967 double arrayValue = array[iRow];
3968 if (fabs(arrayValue) < zeroTolerance)
3969 arrayValue = 0.0;
3970 double newValue = solution(iPivot) + movement * arrayValue;
3971 if (newValue > upper(iPivot) + primalTolerance_ ||
3972 newValue < lower(iPivot) - primalTolerance_)
3973 status = -3; // not unbounded
3974 }
3975 if (status == 2) {
3976 // create ray
3977 delete [] ray_;
3978 ray_ = new double [numberColumns_];
3979 CoinZeroN(ray_, numberColumns_);
3980 for (i = 0; i < number; i++) {
3981 int iRow = index[i];
3982 int iPivot = pivotVariable_[iRow];
3983 double arrayValue = array[iRow];
3984 if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
3985 ray_[iPivot] = way * array[iRow];
3986 }
3987 }
3988 ray->clear();
3989 return status;
3990 }
3991 //static int count_alpha=0;
3992 /* Checks if finished. Updates status */
3993 void
statusOfProblemInDual(int & lastCleaned,int type,double * givenDuals,ClpDataSave & saveData,int ifValuesPass)3994 ClpSimplexDual::statusOfProblemInDual(int & lastCleaned, int type,
3995 double * givenDuals, ClpDataSave & saveData,
3996 int ifValuesPass)
3997 {
3998 #ifdef CLP_INVESTIGATE_SERIAL
3999 if (z_thinks > 0 && z_thinks < 2)
4000 z_thinks += 2;
4001 #endif
4002 bool arraysNotCreated = (type==0);
4003 // If lots of iterations then adjust costs if large ones
4004 if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) {
4005 double largest = 0.0;
4006 for (int i = 0; i < numberRows_; i++) {
4007 int iColumn = pivotVariable_[i];
4008 largest = CoinMax(largest, fabs(cost_[iColumn]));
4009 }
4010 if (largest > 1.0e6) {
4011 objectiveScale_ = 1.0e6 / largest;
4012 for (int i = 0; i < numberRows_ + numberColumns_; i++)
4013 cost_[i] *= objectiveScale_;
4014 }
4015 }
4016 int numberPivots = factorization_->pivots();
4017 double realDualInfeasibilities = 0.0;
4018 if (type == 2) {
4019 if (alphaAccuracy_ != -1.0)
4020 alphaAccuracy_ = -2.0;
4021 // trouble - restore solution
4022 CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4023 CoinMemcpyN(savedSolution_ + numberColumns_ ,
4024 numberRows_, rowActivityWork_);
4025 CoinMemcpyN(savedSolution_ ,
4026 numberColumns_, columnActivityWork_);
4027 // restore extra stuff
4028 int dummy;
4029 matrix_->generalExpanded(this, 6, dummy);
4030 forceFactorization_ = 1; // a bit drastic but ..
4031 changeMade_++; // say something changed
4032 // get correct bounds on all variables
4033 resetFakeBounds(0);
4034 }
4035 int tentativeStatus = problemStatus_;
4036 double changeCost;
4037 bool unflagVariables = true;
4038 bool weightsSaved = false;
4039 bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_;
4040 int dontFactorizePivots = dontFactorizePivots_;
4041 if (type == 3) {
4042 type = 1;
4043 dontFactorizePivots = 1;
4044 }
4045 if (alphaAccuracy_ < 0.0 || !numberPivots || alphaAccuracy_ > 1.0e4 || numberPivots > 20) {
4046 if (problemStatus_ > -3 || numberPivots > dontFactorizePivots) {
4047 // factorize
4048 // later on we will need to recover from singularities
4049 // also we could skip if first time
4050 // save dual weights
4051 dualRowPivot_->saveWeights(this, 1);
4052 weightsSaved = true;
4053 if (type) {
4054 // is factorization okay?
4055 if (internalFactorize(1)) {
4056 // no - restore previous basis
4057 unflagVariables = false;
4058 assert (type == 1);
4059 changeMade_++; // say something changed
4060 // Keep any flagged variables
4061 int i;
4062 for (i = 0; i < numberRows_ + numberColumns_; i++) {
4063 if (flagged(i))
4064 saveStatus_[i] |= 64; //say flagged
4065 }
4066 CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4067 CoinMemcpyN(savedSolution_ + numberColumns_ ,
4068 numberRows_, rowActivityWork_);
4069 CoinMemcpyN(savedSolution_ ,
4070 numberColumns_, columnActivityWork_);
4071 // restore extra stuff
4072 int dummy;
4073 matrix_->generalExpanded(this, 6, dummy);
4074 // get correct bounds on all variables
4075 resetFakeBounds(1);
4076 // need to reject something
4077 char x = isColumn(sequenceOut_) ? 'C' : 'R';
4078 handler_->message(CLP_SIMPLEX_FLAG, messages_)
4079 << x << sequenceWithin(sequenceOut_)
4080 << CoinMessageEol;
4081 #ifdef COIN_DEVELOP
4082 printf("flag d\n");
4083 #endif
4084 setFlagged(sequenceOut_);
4085 progress_.clearBadTimes();
4086
4087 // Go to safe
4088 factorization_->pivotTolerance(0.99);
4089 forceFactorization_ = 1; // a bit drastic but ..
4090 type = 2;
4091 //assert (internalFactorize(1)==0);
4092 if (internalFactorize(1)) {
4093 CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4094 CoinMemcpyN(savedSolution_ + numberColumns_ ,
4095 numberRows_, rowActivityWork_);
4096 CoinMemcpyN(savedSolution_ ,
4097 numberColumns_, columnActivityWork_);
4098 // restore extra stuff
4099 int dummy;
4100 matrix_->generalExpanded(this, 6, dummy);
4101 // debug
4102 int returnCode = internalFactorize(1);
4103 while (returnCode) {
4104 // ouch
4105 // switch off dense
4106 int saveDense = factorization_->denseThreshold();
4107 factorization_->setDenseThreshold(0);
4108 // Go to safe
4109 factorization_->pivotTolerance(0.99);
4110 // make sure will do safe factorization
4111 pivotVariable_[0] = -1;
4112 returnCode = internalFactorize(2);
4113 factorization_->setDenseThreshold(saveDense);
4114 }
4115 // get correct bounds on all variables
4116 resetFakeBounds(1);
4117 }
4118 }
4119 }
4120 if (problemStatus_ != -4 || numberPivots > 10)
4121 problemStatus_ = -3;
4122 }
4123 } else {
4124 //printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
4125 //count_alpha++;
4126 //if ((count_alpha%5000)==0)
4127 //printf("count alpha %d\n",count_alpha);
4128 }
4129 // at this stage status is -3 or -4 if looks infeasible
4130 // get primal and dual solutions
4131 #if 0
4132 {
4133 int numberTotal = numberRows_ + numberColumns_;
4134 double * saveSol = CoinCopyOfArray(solution_, numberTotal);
4135 double * saveDj = CoinCopyOfArray(dj_, numberTotal);
4136 double tolerance = type ? 1.0e-4 : 1.0e-8;
4137 // always if values pass
4138 double saveObj = objectiveValue_;
4139 double sumPrimal = sumPrimalInfeasibilities_;
4140 int numberPrimal = numberPrimalInfeasibilities_;
4141 double sumDual = sumDualInfeasibilities_;
4142 int numberDual = numberDualInfeasibilities_;
4143 gutsOfSolution(givenDuals, NULL);
4144 int j;
4145 double largestPrimal = tolerance;
4146 int iPrimal = -1;
4147 for (j = 0; j < numberTotal; j++) {
4148 double difference = solution_[j] - saveSol[j];
4149 if (fabs(difference) > largestPrimal) {
4150 iPrimal = j;
4151 largestPrimal = fabs(difference);
4152 }
4153 }
4154 double largestDual = tolerance;
4155 int iDual = -1;
4156 for (j = 0; j < numberTotal; j++) {
4157 double difference = dj_[j] - saveDj[j];
4158 if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
4159 iDual = j;
4160 largestDual = fabs(difference);
4161 }
4162 }
4163 if (!type) {
4164 if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4165 numberPrimal != numberPrimalInfeasibilities_ || numberPrimal != 1 ||
4166 fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 || iPrimal >= 0 ||
4167 numberDual != numberDualInfeasibilities_ || numberDual != 0 ||
4168 fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 || iDual >= 0)
4169 printf("type %d its %d pivots %d primal n(%d,%d) s(%g,%g) diff(%g,%d) dual n(%d,%d) s(%g,%g) diff(%g,%d) obj(%g,%g)\n",
4170 type, numberIterations_, numberPivots,
4171 numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4172 largestPrimal, iPrimal,
4173 numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4174 largestDual, iDual,
4175 saveObj, objectiveValue_);
4176 } else {
4177 if (fabs(saveObj - objectiveValue_) > 1.0e-5 ||
4178 numberPrimalInfeasibilities_ || iPrimal >= 0 ||
4179 numberDualInfeasibilities_ || iDual >= 0)
4180 printf("type %d its %d pivots %d primal n(%d,%d) s(%g,%g) diff(%g,%d) dual n(%d,%d) s(%g,%g) diff(%g,%d) obj(%g,%g)\n",
4181 type, numberIterations_, numberPivots,
4182 numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4183 largestPrimal, iPrimal,
4184 numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4185 largestDual, iDual,
4186 saveObj, objectiveValue_);
4187 }
4188 delete [] saveSol;
4189 delete [] saveDj;
4190 }
4191 #else
4192 if (type || ifValuesPass)
4193 gutsOfSolution(givenDuals, NULL);
4194 #endif
4195 // If bad accuracy treat as singular
4196 if ((largestPrimalError_ > 1.0e15 || largestDualError_ > 1.0e15) && numberIterations_) {
4197 // restore previous basis
4198 unflagVariables = false;
4199 changeMade_++; // say something changed
4200 // Keep any flagged variables
4201 int i;
4202 for (i = 0; i < numberRows_ + numberColumns_; i++) {
4203 if (flagged(i))
4204 saveStatus_[i] |= 64; //say flagged
4205 }
4206 CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4207 CoinMemcpyN(savedSolution_ + numberColumns_ ,
4208 numberRows_, rowActivityWork_);
4209 CoinMemcpyN(savedSolution_ ,
4210 numberColumns_, columnActivityWork_);
4211 // restore extra stuff
4212 int dummy;
4213 matrix_->generalExpanded(this, 6, dummy);
4214 // get correct bounds on all variables
4215 resetFakeBounds(1);
4216 // need to reject something
4217 char x = isColumn(sequenceOut_) ? 'C' : 'R';
4218 handler_->message(CLP_SIMPLEX_FLAG, messages_)
4219 << x << sequenceWithin(sequenceOut_)
4220 << CoinMessageEol;
4221 #ifdef COIN_DEVELOP
4222 printf("flag e\n");
4223 #endif
4224 setFlagged(sequenceOut_);
4225 progress_.clearBadTimes();
4226
4227 // Go to safer
4228 double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99);
4229 factorization_->pivotTolerance(newTolerance);
4230 forceFactorization_ = 1; // a bit drastic but ..
4231 if (alphaAccuracy_ != -1.0)
4232 alphaAccuracy_ = -2.0;
4233 type = 2;
4234 //assert (internalFactorize(1)==0);
4235 if (internalFactorize(1)) {
4236 CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4237 CoinMemcpyN(savedSolution_ + numberColumns_ ,
4238 numberRows_, rowActivityWork_);
4239 CoinMemcpyN(savedSolution_ ,
4240 numberColumns_, columnActivityWork_);
4241 // restore extra stuff
4242 int dummy;
4243 matrix_->generalExpanded(this, 6, dummy);
4244 // debug
4245 int returnCode = internalFactorize(1);
4246 while (returnCode) {
4247 // ouch
4248 // switch off dense
4249 int saveDense = factorization_->denseThreshold();
4250 factorization_->setDenseThreshold(0);
4251 // Go to safe
4252 factorization_->pivotTolerance(0.99);
4253 // make sure will do safe factorization
4254 pivotVariable_[0] = -1;
4255 returnCode = internalFactorize(2);
4256 factorization_->setDenseThreshold(saveDense);
4257 }
4258 // get correct bounds on all variables
4259 resetFakeBounds(1);
4260 }
4261 // get primal and dual solutions
4262 gutsOfSolution(givenDuals, NULL);
4263 } else if (goodAccuracy()) {
4264 // Can reduce tolerance
4265 double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
4266 factorization_->pivotTolerance(newTolerance);
4267 }
4268 bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
4269 objectiveValue_ - bestPossibleImprovement_);
4270 bool reallyBadProblems = false;
4271 // Double check infeasibility if no action
4272 if (progress_.lastIterationNumber(0) == numberIterations_) {
4273 if (dualRowPivot_->looksOptimal()) {
4274 numberPrimalInfeasibilities_ = 0;
4275 sumPrimalInfeasibilities_ = 0.0;
4276 }
4277 #if 1
4278 } else {
4279 double thisObj = objectiveValue_ - bestPossibleImprovement_;
4280 #ifdef CLP_INVESTIGATE
4281 assert (bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100);
4282 if (bestPossibleImprovement_)
4283 printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_,
4284 thisObj);
4285 #endif
4286 double lastObj = progress_.lastObjective(0);
4287 #ifndef NDEBUG
4288 #ifdef COIN_DEVELOP
4289 resetFakeBounds(-1);
4290 #endif
4291 #endif
4292 #ifdef CLP_REPORT_PROGRESS
4293 ixxxxxx++;
4294 if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4295 char temp[20];
4296 sprintf(temp, "sol%d.out", ixxxxxx);
4297 printf("sol%d.out\n", ixxxxxx);
4298 FILE * fp = fopen(temp, "w");
4299 int nTotal = numberRows_ + numberColumns_;
4300 for (int i = 0; i < nTotal; i++)
4301 fprintf(fp, "%d %d %g %g %g %g %g\n",
4302 i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4303 fclose(fp);
4304 }
4305 #endif
4306 if(!ifValuesPass && firstFree_ < 0) {
4307 double testTol = 5.0e-3;
4308 if (progress_.timesFlagged() > 10) {
4309 testTol *= pow(2.0, progress_.timesFlagged() - 8);
4310 } else if (progress_.timesFlagged() > 5) {
4311 testTol *= 5.0;
4312 }
4313 if (lastObj > thisObj +
4314 testTol*(fabs(thisObj) + fabs(lastObj)) + testTol) {
4315 int maxFactor = factorization_->maximumPivots();
4316 if ((specialOptions_ & 1048576) == 0) {
4317 if (progress_.timesFlagged() > 10)
4318 progress_.incrementReallyBadTimes();
4319 if (maxFactor > 10 - 9) {
4320 #ifdef COIN_DEVELOP
4321 printf("lastobj %g thisobj %g\n", lastObj, thisObj);
4322 #endif
4323 //if (forceFactorization_<0)
4324 //forceFactorization_= maxFactor;
4325 //forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4326 if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 &&
4327 largestPrimalError_ < 1.0e2) {
4328 // Just save costs
4329 // save extra copy of cost_
4330 int nTotal = numberRows_ + numberColumns_;
4331 double * temp = new double [2*nTotal];
4332 memcpy(temp, cost_, nTotal * sizeof(double));
4333 memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4334 delete [] cost_;
4335 cost_ = temp;
4336 objectiveWork_ = cost_;
4337 rowObjectiveWork_ = cost_ + numberColumns_;
4338 progressFlag_ |= 4;
4339 } else {
4340 forceFactorization_ = 1;
4341 #ifdef COIN_DEVELOP
4342 printf("Reducing factorization frequency - bad backwards\n");
4343 #endif
4344 #if 1
4345 unflagVariables = false;
4346 changeMade_++; // say something changed
4347 int nTotal = numberRows_ + numberColumns_;
4348 CoinMemcpyN(saveStatus_, nTotal, status_);
4349 CoinMemcpyN(savedSolution_ + numberColumns_ ,
4350 numberRows_, rowActivityWork_);
4351 CoinMemcpyN(savedSolution_ ,
4352 numberColumns_, columnActivityWork_);
4353 if ((progressFlag_ & 4) == 0) {
4354 // save extra copy of cost_
4355 double * temp = new double [2*nTotal];
4356 memcpy(temp, cost_, nTotal * sizeof(double));
4357 memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4358 delete [] cost_;
4359 cost_ = temp;
4360 objectiveWork_ = cost_;
4361 rowObjectiveWork_ = cost_ + numberColumns_;
4362 progressFlag_ |= 4;
4363 } else {
4364 memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double));
4365 }
4366 // restore extra stuff
4367 int dummy;
4368 matrix_->generalExpanded(this, 6, dummy);
4369 double pivotTolerance = factorization_->pivotTolerance();
4370 if(pivotTolerance < 0.2)
4371 factorization_->pivotTolerance(0.2);
4372 else if(progress_.timesFlagged() > 2)
4373 factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99));
4374 if (alphaAccuracy_ != -1.0)
4375 alphaAccuracy_ = -2.0;
4376 if (internalFactorize(1)) {
4377 CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4378 CoinMemcpyN(savedSolution_ + numberColumns_ ,
4379 numberRows_, rowActivityWork_);
4380 CoinMemcpyN(savedSolution_ ,
4381 numberColumns_, columnActivityWork_);
4382 // restore extra stuff
4383 int dummy;
4384 matrix_->generalExpanded(this, 6, dummy);
4385 // debug
4386 int returnCode = internalFactorize(1);
4387 while (returnCode) {
4388 // ouch
4389 // switch off dense
4390 int saveDense = factorization_->denseThreshold();
4391 factorization_->setDenseThreshold(0);
4392 // Go to safe
4393 factorization_->pivotTolerance(0.99);
4394 // make sure will do safe factorization
4395 pivotVariable_[0] = -1;
4396 returnCode = internalFactorize(2);
4397 factorization_->setDenseThreshold(saveDense);
4398 }
4399 }
4400 resetFakeBounds(0);
4401 type = 2; // so will restore weights
4402 // get primal and dual solutions
4403 gutsOfSolution(givenDuals, NULL);
4404 if (numberPivots < 2) {
4405 // need to reject something
4406 char x = isColumn(sequenceOut_) ? 'C' : 'R';
4407 handler_->message(CLP_SIMPLEX_FLAG, messages_)
4408 << x << sequenceWithin(sequenceOut_)
4409 << CoinMessageEol;
4410 #ifdef COIN_DEVELOP
4411 printf("flag d\n");
4412 #endif
4413 setFlagged(sequenceOut_);
4414 progress_.clearBadTimes();
4415 progress_.incrementTimesFlagged();
4416 }
4417 if (numberPivots < 10)
4418 reallyBadProblems = true;
4419 #ifdef COIN_DEVELOP
4420 printf("obj now %g\n", objectiveValue_);
4421 #endif
4422 progress_.modifyObjective(objectiveValue_
4423 - bestPossibleImprovement_);
4424 #endif
4425 }
4426 }
4427 } else {
4428 // in fast dual give up
4429 #ifdef COIN_DEVELOP
4430 printf("In fast dual?\n");
4431 #endif
4432 problemStatus_ = 3;
4433 }
4434 } else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) {
4435 numberTimesOptimal_ = 0;
4436 }
4437 }
4438 #endif
4439 }
4440 // Up tolerance if looks a bit odd
4441 if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
4442 if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
4443 int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1);
4444 if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
4445 if (factorization_->pivotTolerance() < 0.9) {
4446 // up tolerance
4447 factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
4448 //printf("tol now %g\n",factorization_->pivotTolerance());
4449 progress_.clearIterationNumbers();
4450 }
4451 }
4452 }
4453 }
4454 // Check if looping
4455 int loop;
4456 if (!givenDuals && type != 2)
4457 loop = progress_.looping();
4458 else
4459 loop = -1;
4460 if (progress_.reallyBadTimes() > 10) {
4461 problemStatus_ = 10; // instead - try other algorithm
4462 #if COIN_DEVELOP>2
4463 printf("returning at %d\n", __LINE__);
4464 #endif
4465 }
4466 int situationChanged = 0;
4467 if (loop >= 0) {
4468 problemStatus_ = loop; //exit if in loop
4469 if (!problemStatus_) {
4470 // declaring victory
4471 numberPrimalInfeasibilities_ = 0;
4472 sumPrimalInfeasibilities_ = 0.0;
4473 } else {
4474 problemStatus_ = 10; // instead - try other algorithm
4475 #if COIN_DEVELOP>2
4476 printf("returning at %d\n", __LINE__);
4477 #endif
4478 }
4479 return;
4480 } else if (loop < -1) {
4481 // something may have changed
4482 gutsOfSolution(NULL, NULL);
4483 situationChanged = 1;
4484 }
4485 // really for free variables in
4486 if((progressFlag_ & 2) != 0) {
4487 situationChanged = 2;
4488 }
4489 progressFlag_ &= (~3); //reset progress flag
4490 if ((progressFlag_ & 4) != 0) {
4491 // save copy of cost_
4492 int nTotal = numberRows_ + numberColumns_;
4493 memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double));
4494 }
4495 /*if (!numberIterations_&&sumDualInfeasibilities_)
4496 printf("OBJ %g sumPinf %g sumDinf %g\n",
4497 objectiveValue(),sumPrimalInfeasibilities_,
4498 sumDualInfeasibilities_);*/
4499 // mark as having gone optimal if looks like it
4500 if (!numberPrimalInfeasibilities_&&
4501 !numberDualInfeasibilities_)
4502 progressFlag_ |= 8;
4503 if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
4504 handler_->message(CLP_SIMPLEX_STATUS, messages_)
4505 << numberIterations_ << objectiveValue();
4506 handler_->printing(sumPrimalInfeasibilities_ > 0.0)
4507 << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
4508 handler_->printing(sumDualInfeasibilities_ > 0.0)
4509 << sumDualInfeasibilities_ << numberDualInfeasibilities_;
4510 handler_->printing(numberDualInfeasibilitiesWithoutFree_
4511 < numberDualInfeasibilities_)
4512 << numberDualInfeasibilitiesWithoutFree_;
4513 handler_->message() << CoinMessageEol;
4514 }
4515 #if 0
4516 printf("IT %d %g %g(%d) %g(%d)\n",
4517 numberIterations_, objectiveValue(),
4518 sumPrimalInfeasibilities_, numberPrimalInfeasibilities_,
4519 sumDualInfeasibilities_, numberDualInfeasibilities_);
4520 #endif
4521 double approximateObjective = objectiveValue_;
4522 #ifdef CLP_REPORT_PROGRESS
4523 if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4524 char temp[20];
4525 sprintf(temp, "x_sol%d.out", ixxxxxx);
4526 FILE * fp = fopen(temp, "w");
4527 int nTotal = numberRows_ + numberColumns_;
4528 for (int i = 0; i < nTotal; i++)
4529 fprintf(fp, "%d %d %g %g %g %g %g\n",
4530 i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4531 fclose(fp);
4532 if (ixxxxxx == ixxyyyy)
4533 exit(6);
4534 }
4535 #endif
4536 realDualInfeasibilities = sumDualInfeasibilities_;
4537 double saveTolerance = dualTolerance_;
4538 // If we need to carry on cleaning variables
4539 if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
4540 for (int iRow = 0; iRow < numberRows_; iRow++) {
4541 int iPivot = pivotVariable_[iRow];
4542 if (!flagged(iPivot) && pivoted(iPivot)) {
4543 // carry on
4544 numberPrimalInfeasibilities_ = -1;
4545 sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4546 sumPrimalInfeasibilities_ = 1.0;
4547 break;
4548 }
4549 }
4550 }
4551 /* If we are primal feasible and any dual infeasibilities are on
4552 free variables then it is better to go to primal */
4553 if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
4554 numberDualInfeasibilities_)
4555 problemStatus_ = 10;
4556 // dual bound coming in
4557 //double saveDualBound = dualBound_;
4558 bool needCleanFake = false;
4559 while (problemStatus_ <= -3) {
4560 int cleanDuals = 0;
4561 if (situationChanged != 0)
4562 cleanDuals = 1;
4563 int numberChangedBounds = 0;
4564 int doOriginalTolerance = 0;
4565 if ( lastCleaned == numberIterations_)
4566 doOriginalTolerance = 1;
4567 // check optimal
4568 // give code benefit of doubt
4569 if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
4570 sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4571 // say optimal (with these bounds etc)
4572 numberDualInfeasibilities_ = 0;
4573 sumDualInfeasibilities_ = 0.0;
4574 numberPrimalInfeasibilities_ = 0;
4575 sumPrimalInfeasibilities_ = 0.0;
4576 }
4577 //if (dualFeasible()||problemStatus_==-4||(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4578 if (dualFeasible() || problemStatus_ == -4) {
4579 progress_.modifyObjective(objectiveValue_
4580 - bestPossibleImprovement_);
4581 #ifdef COIN_DEVELOP
4582 if (sumDualInfeasibilities_ || bestPossibleImprovement_)
4583 printf("improve %g dualinf %g -> %g\n",
4584 bestPossibleImprovement_, sumDualInfeasibilities_,
4585 sumDualInfeasibilities_ * dualBound_);
4586 #endif
4587 // see if cutoff reached
4588 double limit = 0.0;
4589 getDblParam(ClpDualObjectiveLimit, limit);
4590 #if 0
4591 if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
4592 limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) {
4593 //looks infeasible on objective
4594 if (perturbation_ == 101) {
4595 cleanDuals = 1;
4596 // Save costs
4597 int numberTotal = numberRows_ + numberColumns_;
4598 double * saveCost = CoinCopyOfArray(cost_, numberTotal);
4599 // make sure fake bounds are back
4600 changeBounds(1, NULL, changeCost);
4601 createRim4(false);
4602 // make sure duals are current
4603 computeDuals(givenDuals);
4604 checkDualSolution();
4605 if(objectiveValue()*optimizationDirection_ >
4606 limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) {
4607 perturbation_ = 102; // stop any perturbations
4608 printf("cutoff test succeeded\n");
4609 } else {
4610 printf("cutoff test failed\n");
4611 // put back
4612 memcpy(cost_, saveCost, numberTotal * sizeof(double));
4613 // make sure duals are current
4614 computeDuals(givenDuals);
4615 checkDualSolution();
4616 progress_.modifyObjective(-COIN_DBL_MAX);
4617 problemStatus_ = -1;
4618 }
4619 delete [] saveCost;
4620 }
4621 }
4622 #endif
4623 if (primalFeasible() && !givenDuals) {
4624 // may be optimal - or may be bounds are wrong
4625 handler_->message(CLP_DUAL_BOUNDS, messages_)
4626 << dualBound_
4627 << CoinMessageEol;
4628 // save solution in case unbounded
4629 double * saveColumnSolution = NULL;
4630 double * saveRowSolution = NULL;
4631 bool inCbc = (specialOptions_ & (0x01000000 | 16384)) != 0;
4632 if (!inCbc) {
4633 saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_);
4634 saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_);
4635 }
4636 numberChangedBounds = changeBounds(0, rowArray_[3], changeCost);
4637 if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
4638 //looks optimal - do we need to reset tolerance
4639 if (perturbation_ == 101) {
4640 perturbation_ = 102; // stop any perturbations
4641 cleanDuals = 1;
4642 // make sure fake bounds are back
4643 //computeObjectiveValue();
4644 changeBounds(1, NULL, changeCost);
4645 //computeObjectiveValue();
4646 createRim4(false);
4647 // make sure duals are current
4648 computeDuals(givenDuals);
4649 checkDualSolution();
4650 progress_.modifyObjective(-COIN_DBL_MAX);
4651 #define DUAL_TRY_FASTER
4652 #ifdef DUAL_TRY_FASTER
4653 if (numberDualInfeasibilities_) {
4654 #endif
4655 numberChanged_ = 1; // force something to happen
4656 lastCleaned = numberIterations_ - 1;
4657 #ifdef DUAL_TRY_FASTER
4658 } else {
4659 //double value = objectiveValue_;
4660 computeObjectiveValue(true);
4661 //printf("old %g new %g\n",value,objectiveValue_);
4662 //numberChanged_=1;
4663 }
4664 #endif
4665 }
4666 if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 &&
4667 (numberChanged_ || (specialOptions_ & 4096) == 0)) {
4668 doOriginalTolerance = 2;
4669 numberTimesOptimal_++;
4670 changeMade_++; // say something changed
4671 if (numberTimesOptimal_ == 1) {
4672 dualTolerance_ = dblParam_[ClpDualTolerance];
4673 } else {
4674 if (numberTimesOptimal_ == 2) {
4675 // better to have small tolerance even if slower
4676 factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
4677 }
4678 dualTolerance_ = dblParam_[ClpDualTolerance];
4679 dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1);
4680 }
4681 cleanDuals = 2; // If nothing changed optimal else primal
4682 } else {
4683 problemStatus_ = 0; // optimal
4684 if (lastCleaned < numberIterations_ && numberChanged_) {
4685 handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
4686 << CoinMessageEol;
4687 }
4688 }
4689 } else {
4690 cleanDuals = 1;
4691 if (doOriginalTolerance == 1) {
4692 // check unbounded
4693 // find a variable with bad dj
4694 int iSequence;
4695 int iChosen = -1;
4696 if (!inCbc) {
4697 double largest = 100.0 * primalTolerance_;
4698 for (iSequence = 0; iSequence < numberRows_ + numberColumns_;
4699 iSequence++) {
4700 double djValue = dj_[iSequence];
4701 double originalLo = originalLower(iSequence);
4702 double originalUp = originalUpper(iSequence);
4703 if (fabs(djValue) > fabs(largest)) {
4704 if (getStatus(iSequence) != basic) {
4705 if (djValue > 0 && originalLo < -1.0e20) {
4706 if (djValue > fabs(largest)) {
4707 largest = djValue;
4708 iChosen = iSequence;
4709 }
4710 } else if (djValue < 0 && originalUp > 1.0e20) {
4711 if (-djValue > fabs(largest)) {
4712 largest = djValue;
4713 iChosen = iSequence;
4714 }
4715 }
4716 }
4717 }
4718 }
4719 }
4720 if (iChosen >= 0) {
4721 int iSave = sequenceIn_;
4722 sequenceIn_ = iChosen;
4723 unpack(rowArray_[1]);
4724 sequenceIn_ = iSave;
4725 // if dual infeasibilities then must be free vector so add in dual
4726 if (numberDualInfeasibilities_) {
4727 if (fabs(changeCost) > 1.0e-5)
4728 COIN_DETAIL_PRINT(printf("Odd free/unbounded combo\n"));
4729 changeCost += cost_[iChosen];
4730 }
4731 problemStatus_ = checkUnbounded(rowArray_[1], rowArray_[0],
4732 changeCost);
4733 rowArray_[1]->clear();
4734 } else {
4735 problemStatus_ = -3;
4736 }
4737 if (problemStatus_ == 2 && perturbation_ == 101) {
4738 perturbation_ = 102; // stop any perturbations
4739 cleanDuals = 1;
4740 createRim4(false);
4741 progress_.modifyObjective(-COIN_DBL_MAX);
4742 problemStatus_ = -1;
4743 }
4744 if (problemStatus_ == 2) {
4745 // it is unbounded - restore solution
4746 // but first add in changes to non-basic
4747 int iColumn;
4748 double * original = columnArray_[0]->denseVector();
4749 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4750 if(getColumnStatus(iColumn) != basic)
4751 ray_[iColumn] +=
4752 saveColumnSolution[iColumn] - original[iColumn];
4753 columnActivityWork_[iColumn] = original[iColumn];
4754 }
4755 CoinMemcpyN(saveRowSolution, numberRows_,
4756 rowActivityWork_);
4757 }
4758 } else {
4759 doOriginalTolerance = 2;
4760 rowArray_[0]->clear();
4761 }
4762 }
4763 delete [] saveColumnSolution;
4764 delete [] saveRowSolution;
4765 }
4766 if (problemStatus_ == -4 || problemStatus_ == -5) {
4767 // may be infeasible - or may be bounds are wrong
4768 numberChangedBounds = changeBounds(0, NULL, changeCost);
4769 needCleanFake = true;
4770 /* Should this be here as makes no difference to being feasible.
4771 But seems to make a difference to run times. */
4772 if (perturbation_ == 101 && 0) {
4773 perturbation_ = 102; // stop any perturbations
4774 cleanDuals = 1;
4775 numberChangedBounds = 1;
4776 // make sure fake bounds are back
4777 changeBounds(1, NULL, changeCost);
4778 needCleanFake = true;
4779 createRim4(false);
4780 progress_.modifyObjective(-COIN_DBL_MAX);
4781 }
4782 if ((numberChangedBounds <= 0 || dualBound_ > 1.0e20 ||
4783 (largestPrimalError_ > 1.0 && dualBound_ > 1.0e17)) &&
4784 (numberPivots < 4 || sumPrimalInfeasibilities_ > 1.0e-6)) {
4785 problemStatus_ = 1; // infeasible
4786 if (perturbation_ == 101) {
4787 perturbation_ = 102; // stop any perturbations
4788 //cleanDuals=1;
4789 //numberChangedBounds=1;
4790 //createRim4(false);
4791 }
4792 } else {
4793 problemStatus_ = -1; //iterate
4794 cleanDuals = 1;
4795 if (numberChangedBounds <= 0)
4796 doOriginalTolerance = 2;
4797 // and delete ray which has been created
4798 delete [] ray_;
4799 ray_ = NULL;
4800 }
4801
4802 }
4803 } else {
4804 cleanDuals = 1;
4805 }
4806 if (problemStatus_ < 0) {
4807 if (doOriginalTolerance == 2) {
4808 // put back original tolerance
4809 lastCleaned = numberIterations_;
4810 numberChanged_ = 0; // Number of variables with changed costs
4811 handler_->message(CLP_DUAL_ORIGINAL, messages_)
4812 << CoinMessageEol;
4813 perturbation_ = 102; // stop any perturbations
4814 #if 0
4815 double * xcost = new double[numberRows_+numberColumns_];
4816 double * xlower = new double[numberRows_+numberColumns_];
4817 double * xupper = new double[numberRows_+numberColumns_];
4818 double * xdj = new double[numberRows_+numberColumns_];
4819 double * xsolution = new double[numberRows_+numberColumns_];
4820 CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4821 CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4822 CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4823 CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4824 CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4825 #endif
4826 createRim4(false);
4827 progress_.modifyObjective(-COIN_DBL_MAX);
4828 // make sure duals are current
4829 computeDuals(givenDuals);
4830 checkDualSolution();
4831 #if 0
4832 int i;
4833 for (i = 0; i < numberRows_ + numberColumns_; i++) {
4834 if (cost_[i] != xcost[i])
4835 printf("** %d old cost %g new %g sol %g\n",
4836 i, xcost[i], cost_[i], solution_[i]);
4837 if (lower_[i] != xlower[i])
4838 printf("** %d old lower %g new %g sol %g\n",
4839 i, xlower[i], lower_[i], solution_[i]);
4840 if (upper_[i] != xupper[i])
4841 printf("** %d old upper %g new %g sol %g\n",
4842 i, xupper[i], upper_[i], solution_[i]);
4843 if (dj_[i] != xdj[i])
4844 printf("** %d old dj %g new %g sol %g\n",
4845 i, xdj[i], dj_[i], solution_[i]);
4846 if (solution_[i] != xsolution[i])
4847 printf("** %d old solution %g new %g sol %g\n",
4848 i, xsolution[i], solution_[i], solution_[i]);
4849 }
4850 //delete [] xcost;
4851 //delete [] xupper;
4852 //delete [] xlower;
4853 //delete [] xdj;
4854 //delete [] xsolution;
4855 #endif
4856 // put back bounds as they were if was optimal
4857 if (doOriginalTolerance == 2 && cleanDuals != 2) {
4858 changeMade_++; // say something changed
4859 /* We may have already changed some bounds in this function
4860 so save numberFake_ and add in.
4861
4862 Worst that can happen is that we waste a bit of time - but it must be finite.
4863 */
4864 //int saveNumberFake = numberFake_;
4865 //resetFakeBounds(-1);
4866 changeBounds(3, NULL, changeCost);
4867 needCleanFake = true;
4868 //numberFake_ += saveNumberFake;
4869 //resetFakeBounds(-1);
4870 cleanDuals = 2;
4871 //cleanDuals=1;
4872 }
4873 #if 0
4874 //int i;
4875 for (i = 0; i < numberRows_ + numberColumns_; i++) {
4876 if (cost_[i] != xcost[i])
4877 printf("** %d old cost %g new %g sol %g\n",
4878 i, xcost[i], cost_[i], solution_[i]);
4879 if (lower_[i] != xlower[i])
4880 printf("** %d old lower %g new %g sol %g\n",
4881 i, xlower[i], lower_[i], solution_[i]);
4882 if (upper_[i] != xupper[i])
4883 printf("** %d old upper %g new %g sol %g\n",
4884 i, xupper[i], upper_[i], solution_[i]);
4885 if (dj_[i] != xdj[i])
4886 printf("** %d old dj %g new %g sol %g\n",
4887 i, xdj[i], dj_[i], solution_[i]);
4888 if (solution_[i] != xsolution[i])
4889 printf("** %d old solution %g new %g sol %g\n",
4890 i, xsolution[i], solution_[i], solution_[i]);
4891 }
4892 delete [] xcost;
4893 delete [] xupper;
4894 delete [] xlower;
4895 delete [] xdj;
4896 delete [] xsolution;
4897 #endif
4898 }
4899 if (cleanDuals == 1 || (cleanDuals == 2 && !numberDualInfeasibilities_)) {
4900 // make sure dual feasible
4901 // look at all rows and columns
4902 rowArray_[0]->clear();
4903 columnArray_[0]->clear();
4904 double objectiveChange = 0.0;
4905 double savePrimalInfeasibilities = sumPrimalInfeasibilities_;
4906 if (!numberIterations_) {
4907 int nTotal = numberRows_ + numberColumns_;
4908 if (arraysNotCreated) {
4909 // create save arrays
4910 delete [] saveStatus_;
4911 delete [] savedSolution_;
4912 saveStatus_ = new unsigned char [nTotal];
4913 savedSolution_ = new double [nTotal];
4914 arraysNotCreated = false;
4915 }
4916 // save arrays
4917 CoinMemcpyN(status_, nTotal, saveStatus_);
4918 CoinMemcpyN(rowActivityWork_,
4919 numberRows_, savedSolution_ + numberColumns_);
4920 CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
4921 }
4922 #if 0
4923 double * xcost = new double[numberRows_+numberColumns_];
4924 double * xlower = new double[numberRows_+numberColumns_];
4925 double * xupper = new double[numberRows_+numberColumns_];
4926 double * xdj = new double[numberRows_+numberColumns_];
4927 double * xsolution = new double[numberRows_+numberColumns_];
4928 CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4929 CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4930 CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4931 CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4932 CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4933 #endif
4934 if (givenDuals)
4935 dualTolerance_ = 1.0e50;
4936 updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
4937 0.0, objectiveChange, true);
4938 dualTolerance_ = saveTolerance;
4939 #if 0
4940 int i;
4941 for (i = 0; i < numberRows_ + numberColumns_; i++) {
4942 if (cost_[i] != xcost[i])
4943 printf("** %d old cost %g new %g sol %g\n",
4944 i, xcost[i], cost_[i], solution_[i]);
4945 if (lower_[i] != xlower[i])
4946 printf("** %d old lower %g new %g sol %g\n",
4947 i, xlower[i], lower_[i], solution_[i]);
4948 if (upper_[i] != xupper[i])
4949 printf("** %d old upper %g new %g sol %g\n",
4950 i, xupper[i], upper_[i], solution_[i]);
4951 if (dj_[i] != xdj[i])
4952 printf("** %d old dj %g new %g sol %g\n",
4953 i, xdj[i], dj_[i], solution_[i]);
4954 if (solution_[i] != xsolution[i])
4955 printf("** %d old solution %g new %g sol %g\n",
4956 i, xsolution[i], solution_[i], solution_[i]);
4957 }
4958 delete [] xcost;
4959 delete [] xupper;
4960 delete [] xlower;
4961 delete [] xdj;
4962 delete [] xsolution;
4963 #endif
4964 // for now - recompute all
4965 gutsOfSolution(NULL, NULL);
4966 if (givenDuals)
4967 dualTolerance_ = 1.0e50;
4968 updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
4969 0.0, objectiveChange, true);
4970 dualTolerance_ = saveTolerance;
4971 if (!numberIterations_ && sumPrimalInfeasibilities_ >
4972 1.0e5*(savePrimalInfeasibilities+1.0e3) &&
4973 (moreSpecialOptions_ & (256|8192)) == 0) {
4974 // Use primal
4975 int nTotal = numberRows_ + numberColumns_;
4976 CoinMemcpyN(saveStatus_, nTotal, status_);
4977 CoinMemcpyN(savedSolution_ + numberColumns_ ,
4978 numberRows_, rowActivityWork_);
4979 CoinMemcpyN(savedSolution_ ,
4980 numberColumns_, columnActivityWork_);
4981 problemStatus_ = 10;
4982 situationChanged = 0;
4983 }
4984 //assert(numberDualInfeasibilitiesWithoutFree_==0);
4985 if (numberDualInfeasibilities_) {
4986 if ((numberPrimalInfeasibilities_ || numberPivots)
4987 && problemStatus_!=10) {
4988 problemStatus_ = -1; // carry on as normal
4989 } else {
4990 problemStatus_ = 10; // try primal
4991 #if COIN_DEVELOP>1
4992 printf("returning at %d\n", __LINE__);
4993 #endif
4994 }
4995 } else if (situationChanged == 2) {
4996 problemStatus_ = -1; // carry on as normal
4997 // need to reset bounds
4998 changeBounds(3, NULL, changeCost);
4999 }
5000 situationChanged = 0;
5001 } else {
5002 // iterate
5003 if (cleanDuals != 2) {
5004 problemStatus_ = -1;
5005 } else {
5006 problemStatus_ = 10; // try primal
5007 #if COIN_DEVELOP>2
5008 printf("returning at %d\n", __LINE__);
5009 #endif
5010 }
5011 }
5012 }
5013 }
5014 // unflag all variables (we may want to wait a bit?)
5015 if ((tentativeStatus != -2 && tentativeStatus != -1) && unflagVariables) {
5016 int iRow;
5017 int numberFlagged = 0;
5018 for (iRow = 0; iRow < numberRows_; iRow++) {
5019 int iPivot = pivotVariable_[iRow];
5020 if (flagged(iPivot)) {
5021 numberFlagged++;
5022 clearFlagged(iPivot);
5023 }
5024 }
5025 #ifdef COIN_DEVELOP
5026 if (numberFlagged) {
5027 printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus,
5028 problemStatus_, numberPrimalInfeasibilities_,
5029 numberTimesOptimal_);
5030 }
5031 #endif
5032 unflagVariables = numberFlagged > 0;
5033 if (numberFlagged && !numberPivots) {
5034 /* looks like trouble as we have not done any iterations.
5035 Try changing pivot tolerance then give it a few goes and give up */
5036 if (factorization_->pivotTolerance() < 0.9) {
5037 factorization_->pivotTolerance(0.99);
5038 problemStatus_ = -1;
5039 } else if (numberTimesOptimal_ < 3) {
5040 numberTimesOptimal_++;
5041 problemStatus_ = -1;
5042 } else {
5043 unflagVariables = false;
5044 //secondaryStatus_ = 1; // and say probably infeasible
5045 if ((moreSpecialOptions_ & 256) == 0) {
5046 // try primal
5047 problemStatus_ = 10;
5048 } else {
5049 // almost certainly infeasible
5050 problemStatus_ = 1;
5051 }
5052 #if COIN_DEVELOP>1
5053 printf("returning at %d\n", __LINE__);
5054 #endif
5055 }
5056 }
5057 }
5058 if (problemStatus_ < 0) {
5059 if (needCleanFake) {
5060 double dummyChangeCost = 0.0;
5061 changeBounds(3, NULL, dummyChangeCost);
5062 }
5063 #if 0
5064 if (objectiveValue_ < lastObjectiveValue_ - 1.0e-8 *
5065 CoinMax(fabs(objectivevalue_), fabs(lastObjectiveValue_))) {
5066 } else {
5067 lastObjectiveValue_ = objectiveValue_;
5068 }
5069 #endif
5070 if (type == 0 || type == 1) {
5071 if (!type && arraysNotCreated) {
5072 // create save arrays
5073 delete [] saveStatus_;
5074 delete [] savedSolution_;
5075 saveStatus_ = new unsigned char [numberRows_+numberColumns_];
5076 savedSolution_ = new double [numberRows_+numberColumns_];
5077 }
5078 // save arrays
5079 CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
5080 CoinMemcpyN(rowActivityWork_,
5081 numberRows_, savedSolution_ + numberColumns_);
5082 CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
5083 // save extra stuff
5084 int dummy;
5085 matrix_->generalExpanded(this, 5, dummy);
5086 }
5087 if (weightsSaved) {
5088 // restore weights (if saved) - also recompute infeasibility list
5089 if (!reallyBadProblems && (largestPrimalError_ < 100.0 || numberPivots > 10)) {
5090 if (tentativeStatus > -3)
5091 dualRowPivot_->saveWeights(this, (type < 2) ? 2 : 4);
5092 else
5093 dualRowPivot_->saveWeights(this, 3);
5094 } else {
5095 // reset weights or scale back
5096 dualRowPivot_->saveWeights(this, 6);
5097 }
5098 } else if (weightsSaved2 && numberPrimalInfeasibilities_) {
5099 dualRowPivot_->saveWeights(this, 3);
5100 }
5101 }
5102 // see if cutoff reached
5103 double limit = 0.0;
5104 getDblParam(ClpDualObjectiveLimit, limit);
5105 #if 0
5106 if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5107 limit + 100.0) {
5108 printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
5109 limit, objectiveValue_, objectiveValue(), computeInternalObjectiveValue(), sumDualInfeasibilities_);
5110 }
5111 #endif
5112 if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5113 limit && !numberAtFakeBound()) {
5114 bool looksInfeasible = !numberDualInfeasibilities_;
5115 if (objectiveValue()*optimizationDirection_ > limit + fabs(0.1 * limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 &&
5116 sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000)
5117 looksInfeasible = true;
5118 if (looksInfeasible) {
5119 // Even if not perturbed internal costs may have changed
5120 // be careful
5121 if (true || numberIterations_) {
5122 if(computeInternalObjectiveValue() > limit) {
5123 problemStatus_ = 1;
5124 secondaryStatus_ = 1; // and say was on cutoff
5125 }
5126 } else {
5127 problemStatus_ = 1;
5128 secondaryStatus_ = 1; // and say was on cutoff
5129 }
5130 }
5131 }
5132 // If we are in trouble and in branch and bound give up
5133 if ((specialOptions_ & 1024) != 0) {
5134 int looksBad = 0;
5135 if (largestPrimalError_ * largestDualError_ > 1.0e2) {
5136 looksBad = 1;
5137 } else if (largestPrimalError_ > 1.0e-2
5138 && objectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) {
5139 looksBad = 2;
5140 }
5141 if (looksBad) {
5142 if (factorization_->pivotTolerance() < 0.9) {
5143 // up tolerance
5144 factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
5145 } else if (numberIterations_ > 10000) {
5146 if (handler_->logLevel() > 2)
5147 printf("bad dual - saying infeasible %d\n", looksBad);
5148 problemStatus_ = 1;
5149 secondaryStatus_ = 1; // and say was on cutoff
5150 } else if (largestPrimalError_ > 1.0e5) {
5151 {
5152 int iBigB = -1;
5153 double bigB = 0.0;
5154 int iBigN = -1;
5155 double bigN = 0.0;
5156 for (int i = 0; i < numberRows_ + numberColumns_; i++) {
5157 double value = fabs(solution_[i]);
5158 if (getStatus(i) == basic) {
5159 if (value > bigB) {
5160 bigB = value;
5161 iBigB = i;
5162 }
5163 } else {
5164 if (value > bigN) {
5165 bigN = value;
5166 iBigN = i;
5167 }
5168 }
5169 }
5170 #ifdef CLP_INVESTIGATE
5171 if (bigB > 1.0e8 || bigN > 1.0e8) {
5172 if (handler_->logLevel() > 0)
5173 printf("it %d - basic %d %g, nonbasic %d %g\n",
5174 numberIterations_, iBigB, bigB, iBigN, bigN);
5175 }
5176 #endif
5177 }
5178 #if COIN_DEVELOP!=2
5179 if (handler_->logLevel() > 2)
5180 #endif
5181 printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_);
5182 allSlackBasis(true);
5183 problemStatus_ = 10;
5184 }
5185 }
5186 }
5187 if (problemStatus_ < 0 && !changeMade_) {
5188 problemStatus_ = 4; // unknown
5189 }
5190 lastGoodIteration_ = numberIterations_;
5191 if (numberIterations_ > lastBadIteration_ + 100)
5192 moreSpecialOptions_ &= ~16; // clear check accuracy flag
5193 if (problemStatus_ < 0) {
5194 sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
5195 if (sumDualInfeasibilities_)
5196 numberDualInfeasibilities_ = 1;
5197 }
5198 #ifdef CLP_REPORT_PROGRESS
5199 if (ixxxxxx > ixxyyyy - 3) {
5200 printf("objectiveValue_ %g\n", objectiveValue_);
5201 handler_->setLogLevel(63);
5202 int nTotal = numberColumns_ + numberRows_;
5203 double newObj = 0.0;
5204 for (int i = 0; i < nTotal; i++) {
5205 if (solution_[i])
5206 newObj += solution_[i] * cost_[i];
5207 }
5208 printf("xxx obj %g\n", newObj);
5209 // for now - recompute all
5210 gutsOfSolution(NULL, NULL);
5211 newObj = 0.0;
5212 for (int i = 0; i < nTotal; i++) {
5213 if (solution_[i])
5214 newObj += solution_[i] * cost_[i];
5215 }
5216 printf("yyy obj %g %g\n", newObj, objectiveValue_);
5217 progress_.modifyObjective(objectiveValue_
5218 - bestPossibleImprovement_);
5219 }
5220 #endif
5221 #if 1
5222 double thisObj = progress_.lastObjective(0);
5223 double lastObj = progress_.lastObjective(1);
5224 if (lastObj > thisObj + 1.0e-4 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4
5225 && givenDuals == NULL && firstFree_ < 0) {
5226 int maxFactor = factorization_->maximumPivots();
5227 if (maxFactor > 10) {
5228 if (forceFactorization_ < 0)
5229 forceFactorization_ = maxFactor;
5230 forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
5231 //printf("Reducing factorization frequency\n");
5232 }
5233 }
5234 #endif
5235 // Allow matrices to be sorted etc
5236 int fake = -999; // signal sort
5237 matrix_->correctSequence(this, fake, fake);
5238 if (alphaAccuracy_ > 0.0)
5239 alphaAccuracy_ = 1.0;
5240 // If we are stopping - use plausible objective
5241 // Maybe only in fast dual
5242 if (problemStatus_ > 2)
5243 objectiveValue_ = approximateObjective;
5244 if (problemStatus_ == 1 && (progressFlag_&8) != 0 &&
5245 fabs(objectiveValue_) > 1.0e10 )
5246 problemStatus_ = 10; // infeasible - but has looked feasible
5247 }
5248 /* While updateDualsInDual sees what effect is of flip
5249 this does actual flipping.
5250 If change >0.0 then value in array >0.0 => from lower to upper
5251 */
5252 void
flipBounds(CoinIndexedVector * rowArray,CoinIndexedVector * columnArray)5253 ClpSimplexDual::flipBounds(CoinIndexedVector * rowArray,
5254 CoinIndexedVector * columnArray)
5255 {
5256 int number;
5257 int * which;
5258
5259 int iSection;
5260
5261 for (iSection = 0; iSection < 2; iSection++) {
5262 int i;
5263 double * solution = solutionRegion(iSection);
5264 double * lower = lowerRegion(iSection);
5265 double * upper = upperRegion(iSection);
5266 int addSequence;
5267 if (!iSection) {
5268 number = rowArray->getNumElements();
5269 which = rowArray->getIndices();
5270 addSequence = numberColumns_;
5271 } else {
5272 number = columnArray->getNumElements();
5273 which = columnArray->getIndices();
5274 addSequence = 0;
5275 }
5276
5277 for (i = 0; i < number; i++) {
5278 int iSequence = which[i];
5279 Status status = getStatus(iSequence + addSequence);
5280
5281 switch(status) {
5282
5283 case basic:
5284 case isFree:
5285 case superBasic:
5286 case ClpSimplex::isFixed:
5287 break;
5288 case atUpperBound:
5289 // to lower bound
5290 setStatus(iSequence + addSequence, atLowerBound);
5291 solution[iSequence] = lower[iSequence];
5292 break;
5293 case atLowerBound:
5294 // to upper bound
5295 setStatus(iSequence + addSequence, atUpperBound);
5296 solution[iSequence] = upper[iSequence];
5297 break;
5298 }
5299 }
5300 }
5301 rowArray->setNumElements(0);
5302 columnArray->setNumElements(0);
5303 }
5304 // Restores bound to original bound
5305 void
originalBound(int iSequence)5306 ClpSimplexDual::originalBound( int iSequence)
5307 {
5308 if (getFakeBound(iSequence) != noFake) {
5309 numberFake_--;
5310 setFakeBound(iSequence, noFake);
5311 if (iSequence >= numberColumns_) {
5312 // rows
5313 int iRow = iSequence - numberColumns_;
5314 rowLowerWork_[iRow] = rowLower_[iRow];
5315 rowUpperWork_[iRow] = rowUpper_[iRow];
5316 if (rowScale_) {
5317 if (rowLowerWork_[iRow] > -1.0e50)
5318 rowLowerWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5319 if (rowUpperWork_[iRow] < 1.0e50)
5320 rowUpperWork_[iRow] *= rowScale_[iRow] * rhsScale_;
5321 } else if (rhsScale_ != 1.0) {
5322 if (rowLowerWork_[iRow] > -1.0e50)
5323 rowLowerWork_[iRow] *= rhsScale_;
5324 if (rowUpperWork_[iRow] < 1.0e50)
5325 rowUpperWork_[iRow] *= rhsScale_;
5326 }
5327 } else {
5328 // columns
5329 columnLowerWork_[iSequence] = columnLower_[iSequence];
5330 columnUpperWork_[iSequence] = columnUpper_[iSequence];
5331 if (rowScale_) {
5332 double multiplier = 1.0 * inverseColumnScale_[iSequence];
5333 if (columnLowerWork_[iSequence] > -1.0e50)
5334 columnLowerWork_[iSequence] *= multiplier * rhsScale_;
5335 if (columnUpperWork_[iSequence] < 1.0e50)
5336 columnUpperWork_[iSequence] *= multiplier * rhsScale_;
5337 } else if (rhsScale_ != 1.0) {
5338 if (columnLowerWork_[iSequence] > -1.0e50)
5339 columnLowerWork_[iSequence] *= rhsScale_;
5340 if (columnUpperWork_[iSequence] < 1.0e50)
5341 columnUpperWork_[iSequence] *= rhsScale_;
5342 }
5343 }
5344 }
5345 }
5346 /* As changeBounds but just changes new bounds for a single variable.
5347 Returns true if change */
5348 bool
changeBound(int iSequence)5349 ClpSimplexDual::changeBound( int iSequence)
5350 {
5351 // old values
5352 double oldLower = lower_[iSequence];
5353 double oldUpper = upper_[iSequence];
5354 double value = solution_[iSequence];
5355 bool modified = false;
5356 originalBound(iSequence);
5357 // original values
5358 double lowerValue = lower_[iSequence];
5359 double upperValue = upper_[iSequence];
5360 // back to altered values
5361 lower_[iSequence] = oldLower;
5362 upper_[iSequence] = oldUpper;
5363 assert (getFakeBound(iSequence) == noFake);
5364 //if (getFakeBound(iSequence)!=noFake)
5365 //numberFake_--;
5366 if (value == oldLower) {
5367 if (upperValue > oldLower + dualBound_) {
5368 upper_[iSequence] = oldLower + dualBound_;
5369 setFakeBound(iSequence, upperFake);
5370 modified = true;
5371 numberFake_++;
5372 }
5373 } else if (value == oldUpper) {
5374 if (lowerValue < oldUpper - dualBound_) {
5375 lower_[iSequence] = oldUpper - dualBound_;
5376 setFakeBound(iSequence, lowerFake);
5377 modified = true;
5378 numberFake_++;
5379 }
5380 } else {
5381 assert(value == oldLower || value == oldUpper);
5382 }
5383 return modified;
5384 }
5385 // Perturbs problem
5386 int
perturb()5387 ClpSimplexDual::perturb()
5388 {
5389 if (perturbation_ > 100)
5390 return 0; //perturbed already
5391 if (perturbation_ == 100)
5392 perturbation_ = 50; // treat as normal
5393 int savePerturbation = perturbation_;
5394 bool modifyRowCosts = false;
5395 // dual perturbation
5396 double perturbation = 1.0e-20;
5397 // maximum fraction of cost to perturb
5398 double maximumFraction = 1.0e-5;
5399 double constantPerturbation = 100.0 * dualTolerance_;
5400 int maxLength = 0;
5401 int minLength = numberRows_;
5402 double averageCost = 0.0;
5403 #if 0
5404 // look at element range
5405 double smallestNegative;
5406 double largestNegative;
5407 double smallestPositive;
5408 double largestPositive;
5409 matrix_->rangeOfElements(smallestNegative, largestNegative,
5410 smallestPositive, largestPositive);
5411 smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
5412 largestPositive = CoinMax(fabs(largestNegative), largestPositive);
5413 double elementRatio = largestPositive / smallestPositive;
5414 #endif
5415 int numberNonZero = 0;
5416 if (!numberIterations_ && perturbation_ >= 50) {
5417 // See if we need to perturb
5418 double * sort = new double[numberColumns_];
5419 // Use objective BEFORE scaling
5420 const double * obj = ((moreSpecialOptions_ & 128) == 0) ? objective() : cost_;
5421 int i;
5422 for (i = 0; i < numberColumns_; i++) {
5423 double value = fabs(obj[i]);
5424 sort[i] = value;
5425 averageCost += value;
5426 if (value)
5427 numberNonZero++;
5428 }
5429 if (numberNonZero)
5430 averageCost /= static_cast<double> (numberNonZero);
5431 else
5432 averageCost = 1.0;
5433 std::sort(sort, sort + numberColumns_);
5434 int number = 1;
5435 double last = sort[0];
5436 for (i = 1; i < numberColumns_; i++) {
5437 if (last != sort[i])
5438 number++;
5439 last = sort[i];
5440 }
5441 delete [] sort;
5442 if (!numberNonZero && perturbation_ < 55)
5443 return 1; // safer to use primal
5444 #if 0
5445 printf("nnz %d percent %d", number, (number * 100) / numberColumns_);
5446 if (number * 4 > numberColumns_)
5447 printf(" - Would not perturb\n");
5448 else
5449 printf(" - Would perturb\n");
5450 //exit(0);
5451 #endif
5452 //printf("ratio number diff costs %g, element ratio %g\n",((double)number)/((double) numberColumns_),
5453 // elementRatio);
5454 //number=0;
5455 //if (number*4>numberColumns_||elementRatio>1.0e12) {
5456 if (number * 4 > numberColumns_) {
5457 perturbation_ = 100;
5458 return 0; // good enough
5459 }
5460 }
5461 int iColumn;
5462 const int * columnLength = matrix_->getVectorLengths();
5463 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5464 if (columnLowerWork_[iColumn] < columnUpperWork_[iColumn]) {
5465 int length = columnLength[iColumn];
5466 if (length > 2) {
5467 maxLength = CoinMax(maxLength, length);
5468 minLength = CoinMin(minLength, length);
5469 }
5470 }
5471 }
5472 // If > 70 then do rows
5473 if (perturbation_ >= 70) {
5474 modifyRowCosts = true;
5475 perturbation_ -= 20;
5476 printf("Row costs modified, ");
5477 }
5478 bool uniformChange = false;
5479 bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
5480 if (perturbation_ > 50) {
5481 // Experiment
5482 // maximumFraction could be 1.0e-10 to 1.0
5483 double m[] = {1.0e-10, 1.0e-9, 1.0e-8, 1.0e-7, 1.0e-6, 1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1, 1.0};
5484 int whichOne = perturbation_ - 51;
5485 //if (inCbcOrOther&&whichOne>0)
5486 //whichOne--;
5487 maximumFraction = m[CoinMin(whichOne, 10)];
5488 } else if (inCbcOrOther) {
5489 //maximumFraction = 1.0e-6;
5490 }
5491 int iRow;
5492 double smallestNonZero = 1.0e100;
5493 numberNonZero = 0;
5494 if (perturbation_ >= 50) {
5495 perturbation = 1.0e-8;
5496 if (perturbation_ > 50 && perturbation_ < 60)
5497 perturbation = CoinMax(1.0e-8,maximumFraction);
5498 bool allSame = true;
5499 double lastValue = 0.0;
5500 for (iRow = 0; iRow < numberRows_; iRow++) {
5501 double lo = rowLowerWork_[iRow];
5502 double up = rowUpperWork_[iRow];
5503 if (lo < up) {
5504 double value = fabs(rowObjectiveWork_[iRow]);
5505 perturbation = CoinMax(perturbation, value);
5506 if (value) {
5507 modifyRowCosts = true;
5508 smallestNonZero = CoinMin(smallestNonZero, value);
5509 }
5510 }
5511 if (lo && lo > -1.0e10) {
5512 numberNonZero++;
5513 lo = fabs(lo);
5514 if (!lastValue)
5515 lastValue = lo;
5516 else if (fabs(lo - lastValue) > 1.0e-7)
5517 allSame = false;
5518 }
5519 if (up && up < 1.0e10) {
5520 numberNonZero++;
5521 up = fabs(up);
5522 if (!lastValue)
5523 lastValue = up;
5524 else if (fabs(up - lastValue) > 1.0e-7)
5525 allSame = false;
5526 }
5527 }
5528 double lastValue2 = 0.0;
5529 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5530 double lo = columnLowerWork_[iColumn];
5531 double up = columnUpperWork_[iColumn];
5532 if (lo < up) {
5533 double value =
5534 fabs(objectiveWork_[iColumn]);
5535 perturbation = CoinMax(perturbation, value);
5536 if (value) {
5537 smallestNonZero = CoinMin(smallestNonZero, value);
5538 }
5539 }
5540 if (lo && lo > -1.0e10) {
5541 //numberNonZero++;
5542 lo = fabs(lo);
5543 if (!lastValue2)
5544 lastValue2 = lo;
5545 else if (fabs(lo - lastValue2) > 1.0e-7)
5546 allSame = false;
5547 }
5548 if (up && up < 1.0e10) {
5549 //numberNonZero++;
5550 up = fabs(up);
5551 if (!lastValue2)
5552 lastValue2 = up;
5553 else if (fabs(up - lastValue2) > 1.0e-7)
5554 allSame = false;
5555 }
5556 }
5557 if (allSame) {
5558 // Check elements
5559 double smallestNegative;
5560 double largestNegative;
5561 double smallestPositive;
5562 double largestPositive;
5563 matrix_->rangeOfElements(smallestNegative, largestNegative,
5564 smallestPositive, largestPositive);
5565 if (smallestNegative == largestNegative &&
5566 smallestPositive == largestPositive) {
5567 // Really hit perturbation
5568 double adjust = CoinMin(100.0 * maximumFraction, 1.0e-3 * CoinMax(lastValue, lastValue2));
5569 maximumFraction = CoinMax(adjust, maximumFraction);
5570 }
5571 }
5572 perturbation = CoinMin(perturbation, smallestNonZero / maximumFraction);
5573 } else {
5574 // user is in charge
5575 maximumFraction = 1.0e-1;
5576 // but some experiments
5577 if (perturbation_ <= -900) {
5578 modifyRowCosts = true;
5579 perturbation_ += 1000;
5580 printf("Row costs modified, ");
5581 }
5582 if (perturbation_ <= -10) {
5583 perturbation_ += 10;
5584 maximumFraction = 1.0;
5585 if ((-perturbation_) % 100 >= 10) {
5586 uniformChange = true;
5587 perturbation_ += 20;
5588 }
5589 while (perturbation_ < -10) {
5590 perturbation_ += 100;
5591 maximumFraction *= 1.0e-1;
5592 }
5593 }
5594 perturbation = pow(10.0, perturbation_);
5595 }
5596 double largestZero = 0.0;
5597 double largest = 0.0;
5598 double largestPerCent = 0.0;
5599 // modify costs
5600 bool printOut = (handler_->logLevel() == 63);
5601 printOut = false;
5602 //assert (!modifyRowCosts);
5603 modifyRowCosts = false;
5604 if (modifyRowCosts) {
5605 for (iRow = 0; iRow < numberRows_; iRow++) {
5606 if (rowLowerWork_[iRow] < rowUpperWork_[iRow]) {
5607 double value = perturbation;
5608 double currentValue = rowObjectiveWork_[iRow];
5609 value = CoinMin(value, maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-3));
5610 if (rowLowerWork_[iRow] > -largeValue_) {
5611 if (fabs(rowLowerWork_[iRow]) < fabs(rowUpperWork_[iRow]))
5612 value *= randomNumberGenerator_.randomDouble();
5613 else
5614 value *= -randomNumberGenerator_.randomDouble();
5615 } else if (rowUpperWork_[iRow] < largeValue_) {
5616 value *= -randomNumberGenerator_.randomDouble();
5617 } else {
5618 value = 0.0;
5619 }
5620 if (currentValue) {
5621 largest = CoinMax(largest, fabs(value));
5622 if (fabs(value) > fabs(currentValue)*largestPerCent)
5623 largestPerCent = fabs(value / currentValue);
5624 } else {
5625 largestZero = CoinMax(largestZero, fabs(value));
5626 }
5627 if (printOut)
5628 printf("row %d cost %g change %g\n", iRow, rowObjectiveWork_[iRow], value);
5629 rowObjectiveWork_[iRow] += value;
5630 }
5631 }
5632 }
5633 // more its but faster double weight[]={1.0e-4,1.0e-2,1.0e-1,1.0,2.0,10.0,100.0,200.0,400.0,600.0,1000.0};
5634 // good its double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,2.0,5.0,10.0,20.0,30.0,40.0,100.0};
5635 double weight[] = {1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 5.0, 10.0, 20.0, 30.0, 40.0, 100.0};
5636 //double weight[]={1.0e-4,1.0e-2,5.0e-1,1.0,20.0,50.0,100.0,120.0,130.0,140.0,200.0};
5637 double extraWeight = 10.0;
5638 // Scale back if wanted
5639 double weight2[] = {1.0e-4, 1.0e-2, 5.0e-1, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0};
5640 if (constantPerturbation < 99.0 * dualTolerance_) {
5641 perturbation *= 0.1;
5642 extraWeight = 0.5;
5643 memcpy(weight, weight2, sizeof(weight2));
5644 }
5645 // adjust weights if all columns long
5646 double factor = 1.0;
5647 if (maxLength) {
5648 factor = 3.0 / static_cast<double> (minLength);
5649 }
5650 // Make variables with more elements more expensive
5651 const double m1 = 0.5;
5652 double smallestAllowed = CoinMin(1.0e-2 * dualTolerance_, maximumFraction);
5653 double largestAllowed = CoinMax(1.0e3 * dualTolerance_, maximumFraction * averageCost);
5654 // smaller if in BAB
5655 //if (inCbcOrOther)
5656 //largestAllowed=CoinMin(largestAllowed,1.0e-5);
5657 //smallestAllowed = CoinMin(smallestAllowed,0.1*largestAllowed);
5658 #define SAVE_PERT
5659 #ifdef SAVE_PERT
5660 if (2 * numberColumns_ > maximumPerturbationSize_) {
5661 delete [] perturbationArray_;
5662 maximumPerturbationSize_ = 2 * numberColumns_;
5663 perturbationArray_ = new double [maximumPerturbationSize_];
5664 for (iColumn = 0; iColumn < maximumPerturbationSize_; iColumn++) {
5665 perturbationArray_[iColumn] = randomNumberGenerator_.randomDouble();
5666 }
5667 }
5668 #endif
5669 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5670 if (columnLowerWork_[iColumn] < columnUpperWork_[iColumn] && getStatus(iColumn) != basic) {
5671 double value = perturbation;
5672 double currentValue = objectiveWork_[iColumn];
5673 value = CoinMin(value, constantPerturbation + maximumFraction * (fabs(currentValue) + 1.0e-1 * perturbation + 1.0e-8));
5674 //value = CoinMin(value,constantPerturbation;+maximumFraction*fabs(currentValue));
5675 double value2 = constantPerturbation + 1.0e-1 * smallestNonZero;
5676 if (uniformChange) {
5677 value = maximumFraction;
5678 value2 = maximumFraction;
5679 }
5680 if (columnLowerWork_[iColumn] > -largeValue_) {
5681 if (fabs(columnLowerWork_[iColumn]) <
5682 fabs(columnUpperWork_[iColumn])) {
5683 #ifndef SAVE_PERT
5684 value *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
5685 value2 *= (1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
5686 #else
5687 value *= (1.0 - m1 + m1 * perturbationArray_[2*iColumn]);
5688 value2 *= (1.0 - m1 + m1 * perturbationArray_[2*iColumn+1]);
5689 #endif
5690 } else {
5691 //value *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5692 //value2 *= -(1.0-m1+m1*randomNumberGenerator_.randomDouble());
5693 value = 0.0;
5694 }
5695 } else if (columnUpperWork_[iColumn] < largeValue_) {
5696 #ifndef SAVE_PERT
5697 value *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
5698 value2 *= -(1.0 - m1 + m1 * randomNumberGenerator_.randomDouble());
5699 #else
5700 value *= -(1.0 - m1 + m1 * perturbationArray_[2*iColumn]);
5701 value2 *= -(1.0 - m1 + m1 * perturbationArray_[2*iColumn+1]);
5702 #endif
5703 } else {
5704 value = 0.0;
5705 }
5706 if (value) {
5707 int length = columnLength[iColumn];
5708 if (length > 3) {
5709 length = static_cast<int> (static_cast<double> (length) * factor);
5710 length = CoinMax(3, length);
5711 }
5712 double multiplier;
5713 #if 1
5714 if (length < 10)
5715 multiplier = weight[length];
5716 else
5717 multiplier = weight[10];
5718 #else
5719 if (length < 10)
5720 multiplier = weight[length];
5721 else
5722 multiplier = weight[10] + extraWeight * (length - 10);
5723 multiplier *= 0.5;
5724 #endif
5725 value *= multiplier;
5726 value = CoinMin(value, value2);
5727 if (savePerturbation < 50 || savePerturbation > 60) {
5728 if (fabs(value) <= dualTolerance_)
5729 value = 0.0;
5730 } else if (value) {
5731 // get in range
5732 if (fabs(value) <= smallestAllowed) {
5733 value *= 10.0;
5734 while (fabs(value) <= smallestAllowed)
5735 value *= 10.0;
5736 } else if (fabs(value) > largestAllowed) {
5737 value *= 0.1;
5738 while (fabs(value) > largestAllowed)
5739 value *= 0.1;
5740 }
5741 }
5742 if (currentValue) {
5743 largest = CoinMax(largest, fabs(value));
5744 if (fabs(value) > fabs(currentValue)*largestPerCent)
5745 largestPerCent = fabs(value / currentValue);
5746 } else {
5747 largestZero = CoinMax(largestZero, fabs(value));
5748 }
5749 // but negative if at ub
5750 if (getStatus(iColumn) == atUpperBound)
5751 value = -value;
5752 if (printOut)
5753 printf("col %d cost %g change %g\n", iColumn, objectiveWork_[iColumn], value);
5754 objectiveWork_[iColumn] += value;
5755 }
5756 }
5757 }
5758 handler_->message(CLP_SIMPLEX_PERTURB, messages_)
5759 << 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
5760 << CoinMessageEol;
5761 // and zero changes
5762 //int nTotal = numberRows_+numberColumns_;
5763 //CoinZeroN(cost_+nTotal,nTotal);
5764 // say perturbed
5765 perturbation_ = 101;
5766 return 0;
5767 }
5768 /* For strong branching. On input lower and upper are new bounds
5769 while on output they are change in objective function values
5770 (>1.0e50 infeasible).
5771 Return code is 0 if nothing interesting, -1 if infeasible both
5772 ways and +1 if infeasible one way (check values to see which one(s))
5773 Returns -2 if bad factorization
5774 */
strongBranching(int numberVariables,const int * variables,double * newLower,double * newUpper,double ** outputSolution,int * outputStatus,int * outputIterations,bool stopOnFirstInfeasible,bool alwaysFinish,int startFinishOptions)5775 int ClpSimplexDual::strongBranching(int numberVariables, const int * variables,
5776 double * newLower, double * newUpper,
5777 double ** outputSolution,
5778 int * outputStatus, int * outputIterations,
5779 bool stopOnFirstInfeasible,
5780 bool alwaysFinish,
5781 int startFinishOptions)
5782 {
5783 int i;
5784 int returnCode = 0;
5785 double saveObjectiveValue = objectiveValue_;
5786 algorithm_ = -1;
5787
5788 //scaling(false);
5789
5790 // put in standard form (and make row copy)
5791 // create modifiable copies of model rim and do optional scaling
5792 createRim(7 + 8 + 16 + 32, true, startFinishOptions);
5793
5794 // change newLower and newUpper if scaled
5795
5796 // Do initial factorization
5797 // and set certain stuff
5798 // We can either set increasing rows so ...IsBasic gives pivot row
5799 // or we can just increment iBasic one by one
5800 // for now let ...iBasic give pivot row
5801 int useFactorization = false;
5802 if ((startFinishOptions & 2) != 0 && (whatsChanged_&(2 + 512)) == 2 + 512)
5803 useFactorization = true; // Keep factorization if possible
5804 // switch off factorization if bad
5805 if (pivotVariable_[0] < 0)
5806 useFactorization = false;
5807 if (!useFactorization || factorization_->numberRows() != numberRows_) {
5808 useFactorization = false;
5809 factorization_->setDefaultValues();
5810
5811 int factorizationStatus = internalFactorize(0);
5812 if (factorizationStatus < 0) {
5813 // some error
5814 // we should either debug or ignore
5815 #ifndef NDEBUG
5816 printf("***** ClpDual strong branching factorization error - debug\n");
5817 #endif
5818 return -2;
5819 } else if (factorizationStatus && factorizationStatus <= numberRows_) {
5820 handler_->message(CLP_SINGULARITIES, messages_)
5821 << factorizationStatus
5822 << CoinMessageEol;
5823 }
5824 }
5825 // save stuff
5826 ClpFactorization saveFactorization(*factorization_);
5827 // Get fake bounds correctly
5828 double changeCost;
5829 changeBounds(3, NULL, changeCost);
5830 int saveNumberFake = numberFake_;
5831 // save basis and solution
5832 double * saveSolution = new double[numberRows_+numberColumns_];
5833 CoinMemcpyN(solution_,
5834 numberRows_ + numberColumns_, saveSolution);
5835 unsigned char * saveStatus =
5836 new unsigned char [numberRows_+numberColumns_];
5837 CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus);
5838 // save bounds as createRim makes clean copies
5839 double * saveLower = new double[numberRows_+numberColumns_];
5840 CoinMemcpyN(lower_,
5841 numberRows_ + numberColumns_, saveLower);
5842 double * saveUpper = new double[numberRows_+numberColumns_];
5843 CoinMemcpyN(upper_,
5844 numberRows_ + numberColumns_, saveUpper);
5845 double * saveObjective = new double[numberRows_+numberColumns_];
5846 CoinMemcpyN(cost_,
5847 numberRows_ + numberColumns_, saveObjective);
5848 int * savePivot = new int [numberRows_];
5849 CoinMemcpyN(pivotVariable_, numberRows_, savePivot);
5850 // need to save/restore weights.
5851
5852 int iSolution = 0;
5853 for (i = 0; i < numberVariables; i++) {
5854 int iColumn = variables[i];
5855 double objectiveChange;
5856 double saveBound;
5857
5858 // try down
5859
5860 saveBound = columnUpper_[iColumn];
5861 // external view - in case really getting optimal
5862 columnUpper_[iColumn] = newUpper[i];
5863 assert (inverseColumnScale_ || scalingFlag_ <= 0);
5864 if (scalingFlag_ <= 0)
5865 upper_[iColumn] = newUpper[i] * rhsScale_;
5866 else
5867 upper_[iColumn] = (newUpper[i] * inverseColumnScale_[iColumn]) * rhsScale_; // scale
5868 // Start of fast iterations
5869 int status = fastDual(alwaysFinish);
5870 CoinAssert (problemStatus_ || objectiveValue_ < 1.0e50);
5871 #ifdef CLP_DEBUG
5872 printf("first status %d obj %g\n",problemStatus_,objectiveValue_);
5873 #endif
5874 if(problemStatus_==10)
5875 problemStatus_=3;
5876 // make sure plausible
5877 double obj = CoinMax(objectiveValue_, saveObjectiveValue);
5878 if (status && problemStatus_ != 3) {
5879 // not finished - might be optimal
5880 checkPrimalSolution(rowActivityWork_, columnActivityWork_);
5881 double limit = 0.0;
5882 getDblParam(ClpDualObjectiveLimit, limit);
5883 if (!numberPrimalInfeasibilities_ && obj < limit) {
5884 problemStatus_ = 0;
5885 }
5886 status = problemStatus_;
5887 }
5888 if (problemStatus_ == 3)
5889 status = 2;
5890 if (status || (problemStatus_ == 0 && !isDualObjectiveLimitReached())) {
5891 objectiveChange = obj - saveObjectiveValue;
5892 } else {
5893 objectiveChange = 1.0e100;
5894 status = 1;
5895 }
5896
5897 if (scalingFlag_ <= 0) {
5898 CoinMemcpyN(solution_, numberColumns_, outputSolution[iSolution]);
5899 } else {
5900 int j;
5901 double * sol = outputSolution[iSolution];
5902 for (j = 0; j < numberColumns_; j++)
5903 sol[j] = solution_[j] * columnScale_[j];
5904 }
5905 outputStatus[iSolution] = status;
5906 outputIterations[iSolution] = numberIterations_;
5907 iSolution++;
5908 // restore
5909 numberFake_ = saveNumberFake;
5910 CoinMemcpyN(saveSolution,
5911 numberRows_ + numberColumns_, solution_);
5912 CoinMemcpyN(saveStatus, numberColumns_ + numberRows_, status_);
5913 CoinMemcpyN(saveLower,
5914 numberRows_ + numberColumns_, lower_);
5915 CoinMemcpyN(saveUpper,
5916 numberRows_ + numberColumns_, upper_);
5917 CoinMemcpyN(saveObjective,
5918 numberRows_ + numberColumns_, cost_);
5919 columnUpper_[iColumn] = saveBound;
5920 CoinMemcpyN(savePivot, numberRows_, pivotVariable_);
5921 //delete factorization_;
5922 //factorization_ = new ClpFactorization(saveFactorization,numberRows_);
5923 setFactorization(saveFactorization);
5924 newUpper[i] = objectiveChange;
5925 #ifdef CLP_DEBUG
5926 printf("down on %d costs %g\n", iColumn, objectiveChange);
5927 #endif
5928
5929 // try up
5930
5931 saveBound = columnLower_[iColumn];
5932 // external view - in case really getting optimal
5933 columnLower_[iColumn] = newLower[i];
5934 assert (inverseColumnScale_ || scalingFlag_ <= 0);
5935 if (scalingFlag_ <= 0)
5936 lower_[iColumn] = newLower[i] * rhsScale_;
5937 else
5938 lower_[iColumn] = (newLower[i] * inverseColumnScale_[iColumn]) * rhsScale_; // scale
5939 // Start of fast iterations
5940 status = fastDual(alwaysFinish);
5941 CoinAssert (problemStatus_||objectiveValue_<1.0e50);
5942 #ifdef CLP_DEBUG
5943 printf("second status %d obj %g\n",problemStatus_,objectiveValue_);
5944 #endif
5945 if(problemStatus_==10)
5946 problemStatus_=3;
5947 // make sure plausible
5948 obj = CoinMax(objectiveValue_, saveObjectiveValue);
5949 if (status && problemStatus_ != 3) {
5950 // not finished - might be optimal
5951 checkPrimalSolution(rowActivityWork_, columnActivityWork_);
5952 double limit = 0.0;
5953 getDblParam(ClpDualObjectiveLimit, limit);
5954 if (!numberPrimalInfeasibilities_ && obj < limit) {
5955 problemStatus_ = 0;
5956 }
5957 status = problemStatus_;
5958 }
5959 if (problemStatus_ == 3)
5960 status = 2;
5961 if (status || (problemStatus_ == 0 && !isDualObjectiveLimitReached())) {
5962 objectiveChange = obj - saveObjectiveValue;
5963 } else {
5964 objectiveChange = 1.0e100;
5965 status = 1;
5966 }
5967 if (scalingFlag_ <= 0) {
5968 CoinMemcpyN(solution_, numberColumns_, outputSolution[iSolution]);
5969 } else {
5970 int j;
5971 double * sol = outputSolution[iSolution];
5972 for (j = 0; j < numberColumns_; j++)
5973 sol[j] = solution_[j] * columnScale_[j];
5974 }
5975 outputStatus[iSolution] = status;
5976 outputIterations[iSolution] = numberIterations_;
5977 iSolution++;
5978
5979 // restore
5980 numberFake_ = saveNumberFake;
5981 CoinMemcpyN(saveSolution,
5982 numberRows_ + numberColumns_, solution_);
5983 CoinMemcpyN(saveStatus, numberColumns_ + numberRows_, status_);
5984 CoinMemcpyN(saveLower,
5985 numberRows_ + numberColumns_, lower_);
5986 CoinMemcpyN(saveUpper,
5987 numberRows_ + numberColumns_, upper_);
5988 CoinMemcpyN(saveObjective,
5989 numberRows_ + numberColumns_, cost_);
5990 columnLower_[iColumn] = saveBound;
5991 CoinMemcpyN(savePivot, numberRows_, pivotVariable_);
5992 //delete factorization_;
5993 //factorization_ = new ClpFactorization(saveFactorization,numberRows_);
5994 setFactorization(saveFactorization);
5995
5996 newLower[i] = objectiveChange;
5997 #ifdef CLP_DEBUG
5998 printf("up on %d costs %g\n", iColumn, objectiveChange);
5999 #endif
6000
6001 /* Possibilities are:
6002 Both sides feasible - store
6003 Neither side feasible - set objective high and exit if desired
6004 One side feasible - change bounds and resolve
6005 */
6006 if (newUpper[i] < 1.0e100) {
6007 if(newLower[i] < 1.0e100) {
6008 // feasible - no action
6009 } else {
6010 // up feasible, down infeasible
6011 returnCode = 1;
6012 if (stopOnFirstInfeasible)
6013 break;
6014 }
6015 } else {
6016 if(newLower[i] < 1.0e100) {
6017 // down feasible, up infeasible
6018 returnCode = 1;
6019 if (stopOnFirstInfeasible)
6020 break;
6021 } else {
6022 // neither side feasible
6023 returnCode = -1;
6024 break;
6025 }
6026 }
6027 }
6028 delete [] saveSolution;
6029 delete [] saveLower;
6030 delete [] saveUpper;
6031 delete [] saveObjective;
6032 delete [] saveStatus;
6033 delete [] savePivot;
6034 if ((startFinishOptions & 1) == 0) {
6035 deleteRim(1);
6036 whatsChanged_ &= ~0xffff;
6037 } else {
6038 // Original factorization will have been put back by last loop
6039 //delete factorization_;
6040 //factorization_ = new ClpFactorization(saveFactorization);
6041 deleteRim(0);
6042 // mark all as current
6043 whatsChanged_ = 0x3ffffff;
6044 }
6045 objectiveValue_ = saveObjectiveValue;
6046 return returnCode;
6047 }
6048 // treat no pivot as finished (unless interesting)
fastDual(bool alwaysFinish)6049 int ClpSimplexDual::fastDual(bool alwaysFinish)
6050 {
6051 progressFlag_ = 0;
6052 bestObjectiveValue_ = objectiveValue_;
6053 algorithm_ = -1;
6054 secondaryStatus_ = 0;
6055 // Say in fast dual
6056 if (!alwaysFinish)
6057 specialOptions_ |= 1048576;
6058 specialOptions_ |= 16384;
6059 int saveDont = dontFactorizePivots_;
6060 if ((specialOptions_ & 2048) == 0)
6061 dontFactorizePivots_ = 0;
6062 else if(!dontFactorizePivots_)
6063 dontFactorizePivots_ = 20;
6064 //handler_->setLogLevel(63);
6065 // save data
6066 ClpDataSave data = saveData();
6067 dualTolerance_ = dblParam_[ClpDualTolerance];
6068 primalTolerance_ = dblParam_[ClpPrimalTolerance];
6069
6070 // save dual bound
6071 double saveDualBound = dualBound_;
6072
6073 // Start can skip some things in transposeTimes
6074 specialOptions_ |= 131072;
6075 if (alphaAccuracy_ != -1.0)
6076 alphaAccuracy_ = 1.0;
6077 // for dual we will change bounds using dualBound_
6078 // for this we need clean basis so it is after factorize
6079 #if 0
6080 {
6081 int numberTotal = numberRows_ + numberColumns_;
6082 double * saveSol = CoinCopyOfArray(solution_, numberTotal);
6083 double * saveDj = CoinCopyOfArray(dj_, numberTotal);
6084 double tolerance = 1.0e-8;
6085 gutsOfSolution(NULL, NULL);
6086 int j;
6087 double largestPrimal = tolerance;
6088 int iPrimal = -1;
6089 for (j = 0; j < numberTotal; j++) {
6090 double difference = solution_[j] - saveSol[j];
6091 if (fabs(difference) > largestPrimal) {
6092 iPrimal = j;
6093 largestPrimal = fabs(difference);
6094 }
6095 }
6096 double largestDual = tolerance;
6097 int iDual = -1;
6098 for (j = 0; j < numberTotal; j++) {
6099 double difference = dj_[j] - saveDj[j];
6100 if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
6101 iDual = j;
6102 largestDual = fabs(difference);
6103 }
6104 }
6105 if (iPrimal >= 0 || iDual >= 0)
6106 printf("pivots %d primal diff(%g,%d) dual diff(%g,%d)\n",
6107 factorization_->pivots(),
6108 largestPrimal, iPrimal,
6109 largestDual, iDual);
6110 delete [] saveSol;
6111 delete [] saveDj;
6112 }
6113 #else
6114 if ((specialOptions_ & 524288) == 0)
6115 gutsOfSolution(NULL, NULL);
6116 #endif
6117 #if 0
6118 if (numberPrimalInfeasibilities_ != 1 ||
6119 numberDualInfeasibilities_)
6120 printf("dual %g (%d) primal %g (%d)\n",
6121 sumDualInfeasibilities_, numberDualInfeasibilities_,
6122 sumPrimalInfeasibilities_, numberPrimalInfeasibilities_);
6123 #endif
6124 #ifndef NDEBUG
6125 #ifdef COIN_DEVELOP
6126 resetFakeBounds(-1);
6127 #endif
6128 #endif
6129 //numberFake_ =0; // Number of variables at fake bounds
6130 numberChanged_ = 0; // Number of variables with changed costs
6131 //changeBounds(1,NULL,objectiveChange);
6132
6133 problemStatus_ = -1;
6134 numberIterations_ = 0;
6135 if ((specialOptions_ & 524288) == 0) {
6136 factorization_->sparseThreshold(0);
6137 factorization_->goSparse();
6138 }
6139
6140 int lastCleaned = 0; // last time objective or bounds cleaned up
6141
6142 // number of times we have declared optimality
6143 numberTimesOptimal_ = 0;
6144
6145 // This says whether to restore things etc
6146 int factorType = 0;
6147 /*
6148 Status of problem:
6149 0 - optimal
6150 1 - infeasible
6151 2 - unbounded
6152 -1 - iterating
6153 -2 - factorization wanted
6154 -3 - redo checking without factorization
6155 -4 - looks infeasible
6156
6157 BUT also from whileIterating return code is:
6158
6159 -1 iterations etc
6160 -2 inaccuracy
6161 -3 slight inaccuracy (and done iterations)
6162 +0 looks optimal (might be unbounded - but we will investigate)
6163 +1 looks infeasible
6164 +3 max iterations
6165
6166 */
6167
6168 int returnCode = 0;
6169
6170 int iRow, iColumn;
6171 int maxPass = maximumIterations();
6172 while (problemStatus_ < 0) {
6173 // clear
6174 for (iRow = 0; iRow < 4; iRow++) {
6175 rowArray_[iRow]->clear();
6176 }
6177
6178 for (iColumn = 0; iColumn < 2; iColumn++) {
6179 columnArray_[iColumn]->clear();
6180 }
6181
6182 // give matrix (and model costs and bounds a chance to be
6183 // refreshed (normally null)
6184 matrix_->refresh(this);
6185 // If getting nowhere - why not give it a kick
6186 // does not seem to work too well - do some more work
6187 if ((specialOptions_ & 524288) != 0 && (moreSpecialOptions_&2048) == 0 &&
6188 perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_)) {
6189 perturb();
6190 // Can't get here if values pass
6191 gutsOfSolution(NULL, NULL);
6192 if (handler_->logLevel() > 2) {
6193 handler_->message(CLP_SIMPLEX_STATUS, messages_)
6194 << numberIterations_ << objectiveValue();
6195 handler_->printing(sumPrimalInfeasibilities_ > 0.0)
6196 << sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
6197 handler_->printing(sumDualInfeasibilities_ > 0.0)
6198 << sumDualInfeasibilities_ << numberDualInfeasibilities_;
6199 handler_->printing(numberDualInfeasibilitiesWithoutFree_
6200 < numberDualInfeasibilities_)
6201 << numberDualInfeasibilitiesWithoutFree_;
6202 handler_->message() << CoinMessageEol;
6203 }
6204 }
6205 // may factorize, checks if problem finished
6206 // should be able to speed this up on first time
6207 statusOfProblemInDual(lastCleaned, factorType, NULL, data, 0);
6208
6209 // Say good factorization
6210 factorType = 1;
6211 maxPass--;
6212 if (maxPass < -10) {
6213 // odd
6214 returnCode = 1;
6215 problemStatus_ = 3;
6216 // can't say anything interesting - might as well return
6217 #ifdef CLP_DEBUG
6218 printf("returning from fastDual after %d iterations with code %d because of loop\n",
6219 numberIterations_, returnCode);
6220 #endif
6221 break;
6222 }
6223
6224 // Do iterations
6225 if (problemStatus_ < 0) {
6226 double * givenPi = NULL;
6227 returnCode = whileIterating(givenPi, 0);
6228 if ((!alwaysFinish && returnCode < 0) || returnCode == 3) {
6229 if (returnCode != 3)
6230 assert (problemStatus_ < 0);
6231 returnCode = 1;
6232 problemStatus_ = 3;
6233 // can't say anything interesting - might as well return
6234 #ifdef CLP_DEBUG
6235 printf("returning from fastDual after %d iterations with code %d\n",
6236 numberIterations_, returnCode);
6237 #endif
6238 break;
6239 }
6240 if (returnCode == -2)
6241 factorType = 3;
6242 returnCode = 0;
6243 }
6244 }
6245
6246 // clear
6247 for (iRow = 0; iRow < 4; iRow++) {
6248 rowArray_[iRow]->clear();
6249 }
6250
6251 for (iColumn = 0; iColumn < 2; iColumn++) {
6252 columnArray_[iColumn]->clear();
6253 }
6254 // Say not in fast dual
6255 specialOptions_ &= ~(16384 | 1048576);
6256 assert(!numberFake_ || ((specialOptions_&(2048 | 4096)) != 0 && dualBound_ >= 1.0e8)
6257 || returnCode || problemStatus_); // all bounds should be okay
6258 if (numberFake_ > 0 && false) {
6259 // Set back
6260 double dummy;
6261 changeBounds(2, NULL, dummy);
6262 }
6263 // Restore any saved stuff
6264 restoreData(data);
6265 dontFactorizePivots_ = saveDont;
6266 dualBound_ = saveDualBound;
6267 // Stop can skip some things in transposeTimes
6268 specialOptions_ &= ~131072;
6269 if (!problemStatus_) {
6270 // see if cutoff reached
6271 double limit = 0.0;
6272 getDblParam(ClpDualObjectiveLimit, limit);
6273 if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
6274 limit + 1.0e-7 + 1.0e-8 * fabs(limit)) {
6275 // actually infeasible on objective
6276 problemStatus_ = 1;
6277 secondaryStatus_ = 1;
6278 }
6279 }
6280 if (problemStatus_ == 3)
6281 objectiveValue_ = CoinMax(bestObjectiveValue_, objectiveValue_ - bestPossibleImprovement_);
6282 return returnCode;
6283 }
6284 // This does first part of StrongBranching
6285 ClpFactorization *
setupForStrongBranching(char * arrays,int numberRows,int numberColumns,bool solveLp)6286 ClpSimplexDual::setupForStrongBranching(char * arrays, int numberRows,
6287 int numberColumns, bool solveLp)
6288 {
6289 if (solveLp) {
6290 // make sure won't create fake objective
6291 int saveOptions = specialOptions_;
6292 specialOptions_ |= 16384;
6293 // solve
6294 int saveMaximumIterations = intParam_[ClpMaxNumIteration];
6295 intParam_[ClpMaxNumIteration] = 100+numberRows_+numberColumns_;
6296 dual(0, 7);
6297 if (problemStatus_ == 10) {
6298 ClpSimplex::dual(0, 0);
6299 assert (problemStatus_ != 10);
6300 if (problemStatus_ == 0) {
6301 dual(0, 7);
6302 //assert (problemStatus_!=10);
6303 }
6304 }
6305 intParam_[ClpMaxNumIteration] = saveMaximumIterations;
6306 specialOptions_ = saveOptions;
6307 if (problemStatus_ != 0 )
6308 return NULL; // say infeasible or odd
6309 // May be empty
6310 solveLp = (solution_ != NULL && problemStatus_ == 0);
6311 }
6312 problemStatus_ = 0;
6313 if (!solveLp) {
6314 algorithm_ = -1;
6315 // put in standard form (and make row copy)
6316 // create modifiable copies of model rim and do optional scaling
6317 int startFinishOptions;
6318 if((specialOptions_ & 4096) == 0) {
6319 startFinishOptions = 0;
6320 } else {
6321 startFinishOptions = 1 + 2 + 4;
6322 }
6323 createRim(7 + 8 + 16 + 32, true, startFinishOptions);
6324 // Do initial factorization
6325 // and set certain stuff
6326 // We can either set increasing rows so ...IsBasic gives pivot row
6327 // or we can just increment iBasic one by one
6328 // for now let ...iBasic give pivot row
6329 bool useFactorization = false;
6330 if ((startFinishOptions & 2) != 0 && (whatsChanged_&(2 + 512)) == 2 + 512) {
6331 useFactorization = true; // Keep factorization if possible
6332 // switch off factorization if bad
6333 if (pivotVariable_[0] < 0 || factorization_->numberRows() != numberRows_)
6334 useFactorization = false;
6335 }
6336 if (!useFactorization) {
6337 factorization_->setDefaultValues();
6338
6339 int factorizationStatus = internalFactorize(0);
6340 if (factorizationStatus < 0) {
6341 // some error
6342 // we should either debug or ignore
6343 #ifndef NDEBUG
6344 printf("***** ClpDual strong branching factorization error - debug\n");
6345 #endif
6346 } else if (factorizationStatus && factorizationStatus <= numberRows_) {
6347 handler_->message(CLP_SINGULARITIES, messages_)
6348 << factorizationStatus
6349 << CoinMessageEol;
6350 }
6351 }
6352 }
6353 // Get fake bounds correctly
6354 double dummyChangeCost;
6355 changeBounds(3, NULL, dummyChangeCost);
6356 double * arrayD = reinterpret_cast<double *> (arrays);
6357 arrayD[0] = objectiveValue() * optimizationDirection_;
6358 double * saveSolution = arrayD + 1;
6359 double * saveLower = saveSolution + (numberRows + numberColumns);
6360 double * saveUpper = saveLower + (numberRows + numberColumns);
6361 double * saveObjective = saveUpper + (numberRows + numberColumns);
6362 double * saveLowerOriginal = saveObjective + (numberRows + numberColumns);
6363 double * saveUpperOriginal = saveLowerOriginal + numberColumns;
6364 arrayD = saveUpperOriginal + numberColumns;
6365 int * savePivot = reinterpret_cast<int *> (arrayD);
6366 int * whichRow = savePivot + numberRows;
6367 int * whichColumn = whichRow + 3 * numberRows;
6368 int * arrayI = whichColumn + 2 * numberColumns;
6369 unsigned char * saveStatus = reinterpret_cast<unsigned char *> (arrayI + 1);
6370 // save stuff
6371 // save basis and solution
6372 CoinMemcpyN(solution_,
6373 numberRows_ + numberColumns_, saveSolution);
6374 CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus);
6375 CoinMemcpyN(lower_,
6376 numberRows_ + numberColumns_, saveLower);
6377 CoinMemcpyN(upper_,
6378 numberRows_ + numberColumns_, saveUpper);
6379 CoinMemcpyN(cost_,
6380 numberRows_ + numberColumns_, saveObjective);
6381 CoinMemcpyN(pivotVariable_, numberRows_, savePivot);
6382 ClpFactorization * factorization = factorization_;
6383 factorization_ = NULL;
6384 return factorization;
6385 }
6386 // This cleans up after strong branching
6387 void
cleanupAfterStrongBranching(ClpFactorization * factorization)6388 ClpSimplexDual::cleanupAfterStrongBranching(ClpFactorization * factorization)
6389 {
6390 int startFinishOptions;
6391 /* COIN_CLP_VETTED
6392 Looks safe for Cbc
6393 */
6394 if((specialOptions_ & 4096) == 0) {
6395 startFinishOptions = 0;
6396 } else {
6397 startFinishOptions = 1 + 2 + 4;
6398 }
6399 if ((startFinishOptions & 1) == 0 && cost_) {
6400 deleteRim(1);
6401 } else {
6402 // Original factorization will have been put back by last loop
6403 delete factorization_;
6404 factorization_ = factorization;
6405 //deleteRim(0);
6406 // mark all as current
6407 }
6408 whatsChanged_ &= ~0xffff;
6409 }
6410 /* Checks number of variables at fake bounds. This is used by fastDual
6411 so can exit gracefully before end */
6412 int
numberAtFakeBound()6413 ClpSimplexDual::numberAtFakeBound()
6414 {
6415 int iSequence;
6416 int numberFake = 0;
6417
6418 for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
6419 FakeBound bound = getFakeBound(iSequence);
6420 switch(getStatus(iSequence)) {
6421
6422 case basic:
6423 break;
6424 case isFree:
6425 case superBasic:
6426 case ClpSimplex::isFixed:
6427 //setFakeBound (iSequence, noFake);
6428 break;
6429 case atUpperBound:
6430 if (bound == upperFake || bound == bothFake)
6431 numberFake++;
6432 break;
6433 case atLowerBound:
6434 if (bound == lowerFake || bound == bothFake)
6435 numberFake++;
6436 break;
6437 }
6438 }
6439 //numberFake_ = numberFake;
6440 return numberFake;
6441 }
6442 /* Pivot out a variable and choose an incoing one. Assumes dual
6443 feasible - will not go through a reduced cost.
6444 Returns step length in theta
6445 Returns ray in ray_ (or NULL if no pivot)
6446 Return codes as before but -1 means no acceptable pivot
6447 */
6448 int
pivotResult()6449 ClpSimplexDual::pivotResult()
6450 {
6451 abort();
6452 return 0;
6453 }
6454 /*
6455 Row array has row part of pivot row
6456 Column array has column part.
6457 This is used in dual values pass
6458 */
6459 void
checkPossibleValuesMove(CoinIndexedVector * rowArray,CoinIndexedVector * columnArray,double acceptablePivot)6460 ClpSimplexDual::checkPossibleValuesMove(CoinIndexedVector * rowArray,
6461 CoinIndexedVector * columnArray,
6462 double acceptablePivot)
6463 {
6464 double * work;
6465 int number;
6466 int * which;
6467 int iSection;
6468
6469 double tolerance = dualTolerance_ * 1.001;
6470
6471 double thetaDown = 1.0e31;
6472 double changeDown ;
6473 double thetaUp = 1.0e31;
6474 double bestAlphaDown = acceptablePivot * 0.99999;
6475 double bestAlphaUp = acceptablePivot * 0.99999;
6476 int sequenceDown = -1;
6477 int sequenceUp = sequenceOut_;
6478
6479 double djBasic = dj_[sequenceOut_];
6480 if (djBasic > 0.0) {
6481 // basic at lower bound so directionOut_ 1 and -1 in pivot row
6482 // dj will go to zero on other way
6483 thetaUp = djBasic;
6484 changeDown = -lower_[sequenceOut_];
6485 } else {
6486 // basic at upper bound so directionOut_ -1 and 1 in pivot row
6487 // dj will go to zero on other way
6488 thetaUp = -djBasic;
6489 changeDown = upper_[sequenceOut_];
6490 }
6491 bestAlphaUp = 1.0;
6492 int addSequence;
6493
6494 double alphaUp = 0.0;
6495 double alphaDown = 0.0;
6496
6497 for (iSection = 0; iSection < 2; iSection++) {
6498
6499 int i;
6500 if (!iSection) {
6501 work = rowArray->denseVector();
6502 number = rowArray->getNumElements();
6503 which = rowArray->getIndices();
6504 addSequence = numberColumns_;
6505 } else {
6506 work = columnArray->denseVector();
6507 number = columnArray->getNumElements();
6508 which = columnArray->getIndices();
6509 addSequence = 0;
6510 }
6511
6512 for (i = 0; i < number; i++) {
6513 int iSequence = which[i];
6514 int iSequence2 = iSequence + addSequence;
6515 double alpha;
6516 double oldValue;
6517 double value;
6518
6519 switch(getStatus(iSequence2)) {
6520
6521 case basic:
6522 break;
6523 case ClpSimplex::isFixed:
6524 alpha = work[i];
6525 changeDown += alpha * upper_[iSequence2];
6526 break;
6527 case isFree:
6528 case superBasic:
6529 alpha = work[i];
6530 // dj must be effectively zero as dual feasible
6531 if (fabs(alpha) > bestAlphaUp) {
6532 thetaDown = 0.0;
6533 thetaUp = 0.0;
6534 bestAlphaDown = fabs(alpha);
6535 bestAlphaUp = bestAlphaDown;
6536 sequenceDown = iSequence2;
6537 sequenceUp = sequenceDown;
6538 alphaUp = alpha;
6539 alphaDown = alpha;
6540 }
6541 break;
6542 case atUpperBound:
6543 alpha = work[i];
6544 oldValue = dj_[iSequence2];
6545 changeDown += alpha * upper_[iSequence2];
6546 if (alpha >= acceptablePivot) {
6547 // might do other way
6548 value = oldValue + thetaUp * alpha;
6549 if (value > -tolerance) {
6550 if (value > tolerance || fabs(alpha) > bestAlphaUp) {
6551 thetaUp = -oldValue / alpha;
6552 bestAlphaUp = fabs(alpha);
6553 sequenceUp = iSequence2;
6554 alphaUp = alpha;
6555 }
6556 }
6557 } else if (alpha <= -acceptablePivot) {
6558 // might do this way
6559 value = oldValue - thetaDown * alpha;
6560 if (value > -tolerance) {
6561 if (value > tolerance || fabs(alpha) > bestAlphaDown) {
6562 thetaDown = oldValue / alpha;
6563 bestAlphaDown = fabs(alpha);
6564 sequenceDown = iSequence2;
6565 alphaDown = alpha;
6566 }
6567 }
6568 }
6569 break;
6570 case atLowerBound:
6571 alpha = work[i];
6572 oldValue = dj_[iSequence2];
6573 changeDown += alpha * lower_[iSequence2];
6574 if (alpha <= -acceptablePivot) {
6575 // might do other way
6576 value = oldValue + thetaUp * alpha;
6577 if (value < tolerance) {
6578 if (value < -tolerance || fabs(alpha) > bestAlphaUp) {
6579 thetaUp = -oldValue / alpha;
6580 bestAlphaUp = fabs(alpha);
6581 sequenceUp = iSequence2;
6582 alphaUp = alpha;
6583 }
6584 }
6585 } else if (alpha >= acceptablePivot) {
6586 // might do this way
6587 value = oldValue - thetaDown * alpha;
6588 if (value < tolerance) {
6589 if (value < -tolerance || fabs(alpha) > bestAlphaDown) {
6590 thetaDown = oldValue / alpha;
6591 bestAlphaDown = fabs(alpha);
6592 sequenceDown = iSequence2;
6593 alphaDown = alpha;
6594 }
6595 }
6596 }
6597 break;
6598 }
6599 }
6600 }
6601 thetaUp *= -1.0;
6602 double changeUp = -thetaUp * changeDown;
6603 changeDown = -thetaDown * changeDown;
6604 if (CoinMax(fabs(thetaDown), fabs(thetaUp)) < 1.0e-8) {
6605 // largest
6606 if (fabs(alphaDown) < fabs(alphaUp)) {
6607 sequenceDown = -1;
6608 }
6609 }
6610 // choose
6611 sequenceIn_ = -1;
6612 if (changeDown > changeUp && sequenceDown >= 0) {
6613 theta_ = thetaDown;
6614 if (fabs(changeDown) < 1.0e30)
6615 sequenceIn_ = sequenceDown;
6616 alpha_ = alphaDown;
6617 #ifdef CLP_DEBUG
6618 if ((handler_->logLevel() & 32))
6619 printf("predicted way - dirout %d, change %g,%g theta %g\n",
6620 directionOut_, changeDown, changeUp, theta_);
6621 #endif
6622 } else {
6623 theta_ = thetaUp;
6624 if (fabs(changeUp) < 1.0e30)
6625 sequenceIn_ = sequenceUp;
6626 alpha_ = alphaUp;
6627 if (sequenceIn_ != sequenceOut_) {
6628 #ifdef CLP_DEBUG
6629 if ((handler_->logLevel() & 32))
6630 printf("opposite way - dirout %d, change %g,%g theta %g\n",
6631 directionOut_, changeDown, changeUp, theta_);
6632 #endif
6633 } else {
6634 #ifdef CLP_DEBUG
6635 if ((handler_->logLevel() & 32))
6636 printf("opposite way to zero dj - dirout %d, change %g,%g theta %g\n",
6637 directionOut_, changeDown, changeUp, theta_);
6638 #endif
6639 }
6640 }
6641 if (sequenceIn_ >= 0) {
6642 lowerIn_ = lower_[sequenceIn_];
6643 upperIn_ = upper_[sequenceIn_];
6644 valueIn_ = solution_[sequenceIn_];
6645 dualIn_ = dj_[sequenceIn_];
6646
6647 if (alpha_ < 0.0) {
6648 // as if from upper bound
6649 directionIn_ = -1;
6650 upperIn_ = valueIn_;
6651 } else {
6652 // as if from lower bound
6653 directionIn_ = 1;
6654 lowerIn_ = valueIn_;
6655 }
6656 }
6657 }
6658 /*
6659 Row array has row part of pivot row
6660 Column array has column part.
6661 This is used in cleanup
6662 */
6663 void
checkPossibleCleanup(CoinIndexedVector * rowArray,CoinIndexedVector * columnArray,double acceptablePivot)6664 ClpSimplexDual::checkPossibleCleanup(CoinIndexedVector * rowArray,
6665 CoinIndexedVector * columnArray,
6666 double acceptablePivot)
6667 {
6668 double * work;
6669 int number;
6670 int * which;
6671 int iSection;
6672
6673 double tolerance = dualTolerance_ * 1.001;
6674
6675 double thetaDown = 1.0e31;
6676 double thetaUp = 1.0e31;
6677 double bestAlphaDown = acceptablePivot * 10.0;
6678 double bestAlphaUp = acceptablePivot * 10.0;
6679 int sequenceDown = -1;
6680 int sequenceUp = -1;
6681
6682 double djSlack = dj_[pivotRow_];
6683 if (getRowStatus(pivotRow_) == basic)
6684 djSlack = COIN_DBL_MAX;
6685 if (fabs(djSlack) < tolerance)
6686 djSlack = 0.0;
6687 int addSequence;
6688
6689 double alphaUp = 0.0;
6690 double alphaDown = 0.0;
6691 for (iSection = 0; iSection < 2; iSection++) {
6692
6693 int i;
6694 if (!iSection) {
6695 work = rowArray->denseVector();
6696 number = rowArray->getNumElements();
6697 which = rowArray->getIndices();
6698 addSequence = numberColumns_;
6699 } else {
6700 work = columnArray->denseVector();
6701 number = columnArray->getNumElements();
6702 which = columnArray->getIndices();
6703 addSequence = 0;
6704 }
6705
6706 for (i = 0; i < number; i++) {
6707 int iSequence = which[i];
6708 int iSequence2 = iSequence + addSequence;
6709 double alpha;
6710 double oldValue;
6711 double value;
6712
6713 switch(getStatus(iSequence2)) {
6714
6715 case basic:
6716 break;
6717 case ClpSimplex::isFixed:
6718 alpha = work[i];
6719 if (addSequence) {
6720 COIN_DETAIL_PRINT(printf("possible - pivot row %d this %d\n", pivotRow_, iSequence));
6721 oldValue = dj_[iSequence2];
6722 if (alpha <= -acceptablePivot) {
6723 // might do other way
6724 value = oldValue + thetaUp * alpha;
6725 if (value < tolerance) {
6726 if (value < -tolerance || fabs(alpha) > bestAlphaUp) {
6727 thetaUp = -oldValue / alpha;
6728 bestAlphaUp = fabs(alpha);
6729 sequenceUp = iSequence2;
6730 alphaUp = alpha;
6731 }
6732 }
6733 } else if (alpha >= acceptablePivot) {
6734 // might do this way
6735 value = oldValue - thetaDown * alpha;
6736 if (value < tolerance) {
6737 if (value < -tolerance || fabs(alpha) > bestAlphaDown) {
6738 thetaDown = oldValue / alpha;
6739 bestAlphaDown = fabs(alpha);
6740 sequenceDown = iSequence2;
6741 alphaDown = alpha;
6742 }
6743 }
6744 }
6745 }
6746 break;
6747 case isFree:
6748 case superBasic:
6749 alpha = work[i];
6750 // dj must be effectively zero as dual feasible
6751 if (fabs(alpha) > bestAlphaUp) {
6752 thetaDown = 0.0;
6753 thetaUp = 0.0;
6754 bestAlphaDown = fabs(alpha);
6755 bestAlphaUp = bestAlphaDown;
6756 sequenceDown = iSequence2;
6757 sequenceUp = sequenceDown;
6758 alphaUp = alpha;
6759 alphaDown = alpha;
6760 }
6761 break;
6762 case atUpperBound:
6763 alpha = work[i];
6764 oldValue = dj_[iSequence2];
6765 if (alpha >= acceptablePivot) {
6766 // might do other way
6767 value = oldValue + thetaUp * alpha;
6768 if (value > -tolerance) {
6769 if (value > tolerance || fabs(alpha) > bestAlphaUp) {
6770 thetaUp = -oldValue / alpha;
6771 bestAlphaUp = fabs(alpha);
6772 sequenceUp = iSequence2;
6773 alphaUp = alpha;
6774 }
6775 }
6776 } else if (alpha <= -acceptablePivot) {
6777 // might do this way
6778 value = oldValue - thetaDown * alpha;
6779 if (value > -tolerance) {
6780 if (value > tolerance || fabs(alpha) > bestAlphaDown) {
6781 thetaDown = oldValue / alpha;
6782 bestAlphaDown = fabs(alpha);
6783 sequenceDown = iSequence2;
6784 alphaDown = alpha;
6785 }
6786 }
6787 }
6788 break;
6789 case atLowerBound:
6790 alpha = work[i];
6791 oldValue = dj_[iSequence2];
6792 if (alpha <= -acceptablePivot) {
6793 // might do other way
6794 value = oldValue + thetaUp * alpha;
6795 if (value < tolerance) {
6796 if (value < -tolerance || fabs(alpha) > bestAlphaUp) {
6797 thetaUp = -oldValue / alpha;
6798 bestAlphaUp = fabs(alpha);
6799 sequenceUp = iSequence2;
6800 alphaUp = alpha;
6801 }
6802 }
6803 } else if (alpha >= acceptablePivot) {
6804 // might do this way
6805 value = oldValue - thetaDown * alpha;
6806 if (value < tolerance) {
6807 if (value < -tolerance || fabs(alpha) > bestAlphaDown) {
6808 thetaDown = oldValue / alpha;
6809 bestAlphaDown = fabs(alpha);
6810 sequenceDown = iSequence2;
6811 alphaDown = alpha;
6812 }
6813 }
6814 }
6815 break;
6816 }
6817 }
6818 }
6819 thetaUp *= -1.0;
6820 // largest
6821 if (bestAlphaDown < bestAlphaUp)
6822 sequenceDown = -1;
6823 else
6824 sequenceUp = -1;
6825
6826 sequenceIn_ = -1;
6827
6828 if (sequenceDown >= 0) {
6829 theta_ = thetaDown;
6830 sequenceIn_ = sequenceDown;
6831 alpha_ = alphaDown;
6832 #ifdef CLP_DEBUG
6833 if ((handler_->logLevel() & 32))
6834 printf("predicted way - dirout %d, theta %g\n",
6835 directionOut_, theta_);
6836 #endif
6837 } else if (sequenceUp >= 0) {
6838 theta_ = thetaUp;
6839 sequenceIn_ = sequenceUp;
6840 alpha_ = alphaUp;
6841 #ifdef CLP_DEBUG
6842 if ((handler_->logLevel() & 32))
6843 printf("opposite way - dirout %d,theta %g\n",
6844 directionOut_, theta_);
6845 #endif
6846 }
6847 if (sequenceIn_ >= 0) {
6848 lowerIn_ = lower_[sequenceIn_];
6849 upperIn_ = upper_[sequenceIn_];
6850 valueIn_ = solution_[sequenceIn_];
6851 dualIn_ = dj_[sequenceIn_];
6852
6853 if (alpha_ < 0.0) {
6854 // as if from upper bound
6855 directionIn_ = -1;
6856 upperIn_ = valueIn_;
6857 } else {
6858 // as if from lower bound
6859 directionIn_ = 1;
6860 lowerIn_ = valueIn_;
6861 }
6862 }
6863 }
6864 /*
6865 This sees if we can move duals in dual values pass.
6866 This is done before any pivoting
6867 */
doEasyOnesInValuesPass(double * dj)6868 void ClpSimplexDual::doEasyOnesInValuesPass(double * dj)
6869 {
6870 // Get column copy
6871 CoinPackedMatrix * columnCopy = matrix();
6872 // Get a row copy in standard format
6873 CoinPackedMatrix copy;
6874 copy.setExtraGap(0.0);
6875 copy.setExtraMajor(0.0);
6876 copy.reverseOrderedCopyOf(*columnCopy);
6877 // get matrix data pointers
6878 const int * column = copy.getIndices();
6879 const CoinBigIndex * rowStart = copy.getVectorStarts();
6880 const int * rowLength = copy.getVectorLengths();
6881 const double * elementByRow = copy.getElements();
6882 double tolerance = dualTolerance_ * 1.001;
6883
6884 int iRow;
6885 #ifdef CLP_DEBUG
6886 {
6887 double value5 = 0.0;
6888 int i;
6889 for (i = 0; i < numberRows_ + numberColumns_; i++) {
6890 if (dj[i] < -1.0e-6)
6891 value5 += dj[i] * upper_[i];
6892 else if (dj[i] > 1.0e-6)
6893 value5 += dj[i] * lower_[i];
6894 }
6895 printf("Values objective Value before %g\n", value5);
6896 }
6897 #endif
6898 // for scaled row
6899 double * scaled = NULL;
6900 if (rowScale_)
6901 scaled = new double[numberColumns_];
6902 for (iRow = 0; iRow < numberRows_; iRow++) {
6903
6904 int iSequence = iRow + numberColumns_;
6905 double djBasic = dj[iSequence];
6906 if (getRowStatus(iRow) == basic && fabs(djBasic) > tolerance) {
6907
6908 double changeUp ;
6909 // always -1 in pivot row
6910 if (djBasic > 0.0) {
6911 // basic at lower bound
6912 changeUp = -lower_[iSequence];
6913 } else {
6914 // basic at upper bound
6915 changeUp = upper_[iSequence];
6916 }
6917 bool canMove = true;
6918 int i;
6919 const double * thisElements = elementByRow + rowStart[iRow];
6920 const int * thisIndices = column + rowStart[iRow];
6921 if (rowScale_) {
6922 // scale row
6923 double scale = rowScale_[iRow];
6924 for (i = 0; i < rowLength[iRow]; i++) {
6925 int iColumn = thisIndices[i];
6926 double alpha = thisElements[i];
6927 scaled[i] = scale * alpha * columnScale_[iColumn];
6928 }
6929 thisElements = scaled;
6930 }
6931 for (i = 0; i < rowLength[iRow]; i++) {
6932 int iColumn = thisIndices[i];
6933 double alpha = thisElements[i];
6934 double oldValue = dj[iColumn];
6935 double value;
6936
6937 switch(getStatus(iColumn)) {
6938
6939 case basic:
6940 if (dj[iColumn] < -tolerance &&
6941 fabs(solution_[iColumn] - upper_[iColumn]) < 1.0e-8) {
6942 // at ub
6943 changeUp += alpha * upper_[iColumn];
6944 // might do other way
6945 value = oldValue + djBasic * alpha;
6946 if (value > tolerance)
6947 canMove = false;
6948 } else if (dj[iColumn] > tolerance &&
6949 fabs(solution_[iColumn] - lower_[iColumn]) < 1.0e-8) {
6950 changeUp += alpha * lower_[iColumn];
6951 // might do other way
6952 value = oldValue + djBasic * alpha;
6953 if (value < -tolerance)
6954 canMove = false;
6955 } else {
6956 canMove = false;
6957 }
6958 break;
6959 case ClpSimplex::isFixed:
6960 changeUp += alpha * upper_[iColumn];
6961 break;
6962 case isFree:
6963 case superBasic:
6964 canMove = false;
6965 break;
6966 case atUpperBound:
6967 changeUp += alpha * upper_[iColumn];
6968 // might do other way
6969 value = oldValue + djBasic * alpha;
6970 if (value > tolerance)
6971 canMove = false;
6972 break;
6973 case atLowerBound:
6974 changeUp += alpha * lower_[iColumn];
6975 // might do other way
6976 value = oldValue + djBasic * alpha;
6977 if (value < -tolerance)
6978 canMove = false;
6979 break;
6980 }
6981 }
6982 if (canMove) {
6983 if (changeUp * djBasic > 1.0e-12 || fabs(changeUp) < 1.0e-8) {
6984 // move
6985 for (i = 0; i < rowLength[iRow]; i++) {
6986 int iColumn = thisIndices[i];
6987 double alpha = thisElements[i];
6988 dj[iColumn] += djBasic * alpha;
6989 }
6990 dj[iSequence] = 0.0;
6991 #ifdef CLP_DEBUG
6992 {
6993 double value5 = 0.0;
6994 int i;
6995 for (i = 0; i < numberRows_ + numberColumns_; i++) {
6996 if (dj[i] < -1.0e-6)
6997 value5 += dj[i] * upper_[i];
6998 else if (dj[i] > 1.0e-6)
6999 value5 += dj[i] * lower_[i];
7000 }
7001 printf("Values objective Value after row %d old dj %g %g\n",
7002 iRow, djBasic, value5);
7003 }
7004 #endif
7005 }
7006 }
7007 }
7008 }
7009 delete [] scaled;
7010 }
7011 int
nextSuperBasic()7012 ClpSimplexDual::nextSuperBasic()
7013 {
7014 if (firstFree_ >= 0) {
7015 int returnValue = firstFree_;
7016 int iColumn = firstFree_ + 1;
7017 for (; iColumn < numberRows_ + numberColumns_; iColumn++) {
7018 if (getStatus(iColumn) == isFree)
7019 if (fabs(dj_[iColumn]) > 1.0e2 * dualTolerance_)
7020 break;
7021 }
7022 firstFree_ = iColumn;
7023 if (firstFree_ == numberRows_ + numberColumns_)
7024 firstFree_ = -1;
7025 return returnValue;
7026 } else {
7027 return -1;
7028 }
7029 }
7030 void
resetFakeBounds(int type)7031 ClpSimplexDual::resetFakeBounds(int type)
7032 {
7033 if (type == 0) {
7034 // put back original bounds and then check
7035 createRim1(false);
7036 double dummyChangeCost = 0.0;
7037 changeBounds(3, NULL, dummyChangeCost);
7038 } else if (type < 0) {
7039 #ifndef NDEBUG
7040 // just check
7041 int nTotal = numberRows_ + numberColumns_;
7042 double * tempLower = CoinCopyOfArray(lower_, nTotal);
7043 double * tempUpper = CoinCopyOfArray(upper_, nTotal);
7044 int iSequence;
7045 // Get scaled true bounds
7046 if (columnScale_) {
7047 for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
7048 // lower
7049 double value = columnLower_[iSequence];
7050 if (value > -1.0e30) {
7051 double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
7052 value *= multiplier;
7053 }
7054 tempLower[iSequence] = value;
7055 // upper
7056 value = columnUpper_[iSequence];
7057 if (value < 1.0e30) {
7058 double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
7059 value *= multiplier;
7060 }
7061 tempUpper[iSequence] = value;
7062 }
7063 for (iSequence = 0; iSequence < numberRows_; iSequence++) {
7064 // lower
7065 double value = rowLower_[iSequence];
7066 if (value > -1.0e30) {
7067 double multiplier = rhsScale_ * rowScale_[iSequence];
7068 value *= multiplier;
7069 }
7070 tempLower[iSequence+numberColumns_] = value;
7071 // upper
7072 value = rowUpper_[iSequence];
7073 if (value < 1.0e30) {
7074 double multiplier = rhsScale_ * rowScale_[iSequence];
7075 value *= multiplier;
7076 }
7077 tempUpper[iSequence+numberColumns_] = value;
7078 }
7079 } else {
7080 for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
7081 // lower
7082 tempLower[iSequence] = columnLower_[iSequence];
7083 // upper
7084 tempUpper[iSequence] = columnUpper_[iSequence];
7085 }
7086 for (iSequence = 0; iSequence < numberRows_; iSequence++) {
7087 // lower
7088 tempLower[iSequence+numberColumns_] = rowLower_[iSequence];
7089 // upper
7090 tempUpper[iSequence+numberColumns_] = rowUpper_[iSequence];
7091 }
7092 }
7093 int nFake = 0;
7094 int nErrors = 0;
7095 int nSuperBasic = 0;
7096 int nWarnings = 0;
7097 for (iSequence = 0; iSequence < nTotal; iSequence++) {
7098 FakeBound fakeStatus = getFakeBound(iSequence);
7099 Status status = getStatus(iSequence);
7100 bool isFake = false;
7101 char RC = 'C';
7102 int jSequence = iSequence;
7103 if (jSequence >= numberColumns_) {
7104 RC = 'R';
7105 jSequence -= numberColumns_;
7106 }
7107 double lowerValue = tempLower[iSequence];
7108 double upperValue = tempUpper[iSequence];
7109 double value = solution_[iSequence];
7110 CoinRelFltEq equal;
7111 if (status == atUpperBound ||
7112 status == atLowerBound) {
7113 if (fakeStatus == ClpSimplexDual::upperFake) {
7114 if(!equal(upper_[iSequence], (lowerValue + dualBound_)) ||
7115 !(equal(upper_[iSequence], value) ||
7116 equal(lower_[iSequence], value))) {
7117 nErrors++;
7118 #ifdef CLP_INVESTIGATE
7119 printf("** upperFake %c%d %g <= %g <= %g true %g, %g\n",
7120 RC, jSequence, lower_[iSequence], solution_[iSequence],
7121 upper_[iSequence], lowerValue, upperValue);
7122 #endif
7123 }
7124 isFake = true;
7125 } else if (fakeStatus == ClpSimplexDual::lowerFake) {
7126 if(!equal(lower_[iSequence], (upperValue - dualBound_)) ||
7127 !(equal(upper_[iSequence], value) ||
7128 equal(lower_[iSequence], value))) {
7129 nErrors++;
7130 #ifdef CLP_INVESTIGATE
7131 printf("** lowerFake %c%d %g <= %g <= %g true %g, %g\n",
7132 RC, jSequence, lower_[iSequence], solution_[iSequence],
7133 upper_[iSequence], lowerValue, upperValue);
7134 #endif
7135 }
7136 isFake = true;
7137 } else if (fakeStatus == ClpSimplexDual::bothFake) {
7138 nWarnings++;
7139 #ifdef CLP_INVESTIGATE
7140 printf("** %d at bothFake?\n", iSequence);
7141 #endif
7142 } else if (upper_[iSequence] - lower_[iSequence] > 2.0 * dualBound_) {
7143 nErrors++;
7144 #ifdef CLP_INVESTIGATE
7145 printf("** noFake! %c%d %g <= %g <= %g true %g, %g\n",
7146 RC, jSequence, lower_[iSequence], solution_[iSequence],
7147 upper_[iSequence], lowerValue, upperValue);
7148 #endif
7149 }
7150 } else if (status == superBasic || status == isFree) {
7151 nSuperBasic++;
7152 //printf("** free or superbasic %c%d %g <= %g <= %g true %g, %g - status %d\n",
7153 // RC,jSequence,lower_[iSequence],solution_[iSequence],
7154 // upper_[iSequence],lowerValue,upperValue,status);
7155 } else if (status == basic) {
7156 bool odd = false;
7157 if (!equal(lower_[iSequence], lowerValue))
7158 odd = true;
7159 if (!equal(upper_[iSequence], upperValue))
7160 odd = true;
7161 if (odd) {
7162 #ifdef CLP_INVESTIGATE
7163 printf("** basic %c%d %g <= %g <= %g true %g, %g\n",
7164 RC, jSequence, lower_[iSequence], solution_[iSequence],
7165 upper_[iSequence], lowerValue, upperValue);
7166 #endif
7167 nWarnings++;
7168 }
7169 } else if (status == isFixed) {
7170 if (!equal(upper_[iSequence], lower_[iSequence])) {
7171 nErrors++;
7172 #ifdef CLP_INVESTIGATE
7173 printf("** fixed! %c%d %g <= %g <= %g true %g, %g\n",
7174 RC, jSequence, lower_[iSequence], solution_[iSequence],
7175 upper_[iSequence], lowerValue, upperValue);
7176 #endif
7177 }
7178 }
7179 if (isFake) {
7180 nFake++;
7181 } else {
7182 if (fakeStatus != ClpSimplexDual::noFake) {
7183 nErrors++;
7184 #ifdef CLP_INVESTIGATE
7185 printf("** bad fake status %c%d %d\n",
7186 RC, jSequence, fakeStatus);
7187 #endif
7188 }
7189 }
7190 }
7191 if (nFake != numberFake_) {
7192 #ifdef CLP_INVESTIGATE
7193 printf("nfake %d numberFake %d\n", nFake, numberFake_);
7194 #endif
7195 nErrors++;
7196 }
7197 if (nErrors || type <= -1000) {
7198 #ifdef CLP_INVESTIGATE
7199 printf("%d errors, %d warnings, %d free/superbasic, %d fake\n",
7200 nErrors, nWarnings, nSuperBasic, numberFake_);
7201 printf("dualBound %g\n",
7202 dualBound_);
7203 #endif
7204 if (type <= -1000) {
7205 iSequence = -type;
7206 iSequence -= 1000;
7207 char RC = 'C';
7208 int jSequence = iSequence;
7209 if (jSequence >= numberColumns_) {
7210 RC = 'R';
7211 jSequence -= numberColumns_;
7212 }
7213 #ifdef CLP_INVESTIGATE
7214 double lowerValue = tempLower[iSequence];
7215 double upperValue = tempUpper[iSequence];
7216 printf("*** movement>1.0e30 for %c%d %g <= %g <= %g true %g, %g - status %d\n",
7217 RC, jSequence, lower_[iSequence], solution_[iSequence],
7218 upper_[iSequence], lowerValue, upperValue, status_[iSequence]);
7219 #endif
7220 assert (nErrors); // should have been picked up
7221 }
7222 assert (!nErrors);
7223 }
7224 delete [] tempLower;
7225 delete [] tempUpper;
7226 #endif
7227 } else if (lower_) {
7228 // reset using status
7229 int nTotal = numberRows_ + numberColumns_;
7230 int iSequence;
7231 if (columnScale_) {
7232 for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
7233 double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
7234 // lower
7235 double value = columnLower_[iSequence];
7236 if (value > -1.0e30) {
7237 value *= multiplier;
7238 }
7239 lower_[iSequence] = value;
7240 // upper
7241 value = columnUpper_[iSequence];
7242 if (value < 1.0e30) {
7243 value *= multiplier;
7244 }
7245 upper_[iSequence] = value;
7246 }
7247 for (iSequence = 0; iSequence < numberRows_; iSequence++) {
7248 // lower
7249 double multiplier = rhsScale_ * rowScale_[iSequence];
7250 double value = rowLower_[iSequence];
7251 if (value > -1.0e30) {
7252 value *= multiplier;
7253 }
7254 lower_[iSequence+numberColumns_] = value;
7255 // upper
7256 value = rowUpper_[iSequence];
7257 if (value < 1.0e30) {
7258 value *= multiplier;
7259 }
7260 upper_[iSequence+numberColumns_] = value;
7261 }
7262 } else {
7263 memcpy(lower_, columnLower_, numberColumns_ * sizeof(double));
7264 memcpy(upper_, columnUpper_, numberColumns_ * sizeof(double));
7265 memcpy(lower_ + numberColumns_, rowLower_, numberRows_ * sizeof(double));
7266 memcpy(upper_ + numberColumns_, rowUpper_, numberRows_ * sizeof(double));
7267 }
7268 numberFake_ = 0;
7269 for (iSequence = 0; iSequence < nTotal; iSequence++) {
7270 FakeBound fakeStatus = getFakeBound(iSequence);
7271 if (fakeStatus != ClpSimplexDual::noFake) {
7272 Status status = getStatus(iSequence);
7273 if (status == basic) {
7274 setFakeBound(iSequence, ClpSimplexDual::noFake);
7275 continue;
7276 }
7277 double lowerValue = lower_[iSequence];
7278 double upperValue = upper_[iSequence];
7279 double value = solution_[iSequence];
7280 numberFake_++;
7281 if (fakeStatus == ClpSimplexDual::upperFake) {
7282 upper_[iSequence] = lowerValue + dualBound_;
7283 if (status == ClpSimplex::atLowerBound) {
7284 solution_[iSequence] = lowerValue;
7285 } else if (status == ClpSimplex::atUpperBound) {
7286 solution_[iSequence] = upper_[iSequence];
7287 } else {
7288 abort();
7289 }
7290 } else if (fakeStatus == ClpSimplexDual::lowerFake) {
7291 lower_[iSequence] = upperValue - dualBound_;
7292 if (status == ClpSimplex::atLowerBound) {
7293 solution_[iSequence] = lower_[iSequence];
7294 } else if (status == ClpSimplex::atUpperBound) {
7295 solution_[iSequence] = upperValue;
7296 } else {
7297 abort();
7298 }
7299 } else {
7300 assert (fakeStatus == ClpSimplexDual::bothFake);
7301 if (status == ClpSimplex::atLowerBound) {
7302 lower_[iSequence] = value;
7303 upper_[iSequence] = value + dualBound_;
7304 } else if (status == ClpSimplex::atUpperBound) {
7305 upper_[iSequence] = value;
7306 lower_[iSequence] = value - dualBound_;
7307 } else if (status == ClpSimplex::isFree ||
7308 status == ClpSimplex::superBasic) {
7309 lower_[iSequence] = value - 0.5 * dualBound_;
7310 upper_[iSequence] = value + 0.5 * dualBound_;
7311 } else {
7312 abort();
7313 }
7314 }
7315 }
7316 }
7317 #ifndef NDEBUG
7318 } else {
7319 COIN_DETAIL_PRINT(printf("NULL lower\n"));
7320 #endif
7321 }
7322 }
7323