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

Last change on this file since 2235 was 2235, checked in by forrest, 2 years ago

stuff for vector matrix

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