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

Last change on this file since 1972 was 1972, checked in by forrest, 6 years ago

changes to allow more options and stop on feasible (and a few other things)

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