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