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

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

change some ints to CoinBigIndex?

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