# source:trunk/Clp/src/ClpSimplexPrimal.cpp@1551

Last change on this file since 1551 was 1551, checked in by mjs, 11 years ago

Rerun formatting after JJF patches (1533 and 1535).

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