source: trunk/Clp/src/ClpSimplexPrimal.cpp @ 1717

Last change on this file since 1717 was 1717, checked in by forrest, 9 years ago

changes to try and make more stable

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