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

Last change on this file since 1095 was 1095, checked in by forrest, 14 years ago

event handling and nonlinear

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