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

Last change on this file since 1376 was 1376, checked in by forrest, 11 years ago

chnages for speed and alternate factorization

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