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

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

allow infeasibilityray more and prepare

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