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

Last change on this file since 1868 was 1868, checked in by forrest, 7 years ago

allow fix to force to stay in dual

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