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

Last change on this file since 1088 was 1088, checked in by andreasw, 14 years ago

merging changes from Bug Squashing Party Aug 2007 to regular trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 96.8 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        problemStatus_=3;
665      }
666    } else {
667      // no pivot column
668#ifdef CLP_DEBUG
669      if (handler_->logLevel()&32)
670        printf("** no column pivot\n");
671#endif
672      if (nonLinearCost_->numberInfeasibilities())
673        problemStatus_=-4; // might be infeasible
674      // Force to re-factorize early next time
675      int numberPivots = factorization_->pivots();
676      forceFactorization_=CoinMin(forceFactorization_,(numberPivots+1)>>1);
677      returnCode=0;
678      break;
679    }
680  }
681  if (valuesOption>1) 
682    columnArray_[0]->setNumElements(0);
683  return returnCode;
684}
685/* Checks if finished.  Updates status */
686void 
687ClpSimplexPrimal::statusOfProblemInPrimal(int & lastCleaned,int type,
688                                          ClpSimplexProgress * progress,
689                                          bool doFactorization,
690                                          int ifValuesPass,
691                                          ClpSimplex * originalModel)
692{
693  int dummy; // for use in generalExpanded
694  int saveFirstFree=firstFree_;
695  // number of pivots done
696  int numberPivots = factorization_->pivots();
697  if (type==2) {
698    // trouble - restore solution
699    CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
700    CoinMemcpyN(savedSolution_+numberColumns_ ,
701           numberRows_,rowActivityWork_);
702    CoinMemcpyN(savedSolution_ ,
703           numberColumns_,columnActivityWork_);
704    // restore extra stuff
705    matrix_->generalExpanded(this,6,dummy);
706    forceFactorization_=1; // a bit drastic but ..
707    pivotRow_=-1; // say no weights update
708    changeMade_++; // say change made
709  }
710  int saveThreshold = factorization_->sparseThreshold();
711  int tentativeStatus = problemStatus_;
712  int numberThrownOut=1; // to loop round on bad factorization in values pass
713  double lastSumInfeasibility=COIN_DBL_MAX;
714  if (numberIterations_)
715    lastSumInfeasibility=nonLinearCost_->sumInfeasibilities();
716  int nPass=0;
717  while (numberThrownOut) {
718    int nSlackBasic=0;
719    if (nPass) {
720      for (int i=0;i<numberRows_;i++) {
721        if (getRowStatus(i)==basic)
722          nSlackBasic++;
723      }
724    }
725    nPass++;
726    if (problemStatus_>-3||problemStatus_==-4) {
727      // factorize
728      // later on we will need to recover from singularities
729      // also we could skip if first time
730      // do weights
731      // This may save pivotRow_ for use
732      if (doFactorization)
733        primalColumnPivot_->saveWeights(this,1);
734     
735      if ((type&&doFactorization)||nSlackBasic==numberRows_) {
736        // is factorization okay?
737        int factorStatus = internalFactorize(1);
738        if (factorStatus) {
739          if (solveType_==2+8) {
740            // say odd
741            problemStatus_=5;
742            return;
743          }
744          if (type!=1||largestPrimalError_>1.0e3
745              ||largestDualError_>1.0e3) {
746            // switch off dense
747            int saveDense = factorization_->denseThreshold();
748            factorization_->setDenseThreshold(0);
749            // Go to safe
750            factorization_->pivotTolerance(0.99);
751            // make sure will do safe factorization
752            pivotVariable_[0]=-1;
753            internalFactorize(2);
754            factorization_->setDenseThreshold(saveDense);
755            // restore extra stuff
756            matrix_->generalExpanded(this,6,dummy);
757          } else {
758            // no - restore previous basis
759            // Keep any flagged variables
760            int i;
761            for (i=0;i<numberRows_+numberColumns_;i++) {
762              if (flagged(i))
763                saveStatus_[i] |= 64; //say flagged
764            }
765            CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
766            if (numberPivots<=1) {
767              // throw out something
768              if (sequenceIn_>=0&&getStatus(sequenceIn_)!=basic) {
769                setFlagged(sequenceIn_);
770              } else if (sequenceOut_>=0&&getStatus(sequenceOut_)!=basic) {
771                setFlagged(sequenceOut_);
772              }
773              double newTolerance = CoinMax(0.5 + 0.499*CoinDrand48(),factorization_->pivotTolerance());
774              factorization_->pivotTolerance(newTolerance);
775            } else {
776              // Go to safe
777              factorization_->pivotTolerance(0.99);
778            }
779            CoinMemcpyN(savedSolution_+numberColumns_ ,
780                   numberRows_,rowActivityWork_);
781            CoinMemcpyN(savedSolution_ ,
782                   numberColumns_,columnActivityWork_);
783            // restore extra stuff
784            matrix_->generalExpanded(this,6,dummy);
785            matrix_->generalExpanded(this,5,dummy);
786            forceFactorization_=1; // a bit drastic but ..
787            type = 2;
788            if (internalFactorize(1)!=0)
789               largestPrimalError_=1.0e4; // force other type
790          }
791          changeMade_++; // say change made
792        }
793      }
794      if (problemStatus_!=-4)
795        problemStatus_=-3;
796    }
797    // at this stage status is -3 or -5 if looks unbounded
798    // get primal and dual solutions
799    // put back original costs and then check
800    // createRim(4); // costs do not change
801    // May need to do more if column generation
802    dummy=4;
803    matrix_->generalExpanded(this,9,dummy);
804    numberThrownOut=gutsOfSolution(NULL,NULL,(firstFree_>=0));
805    double sumInfeasibility =  nonLinearCost_->sumInfeasibilities();
806    if (numberThrownOut||
807        (sumInfeasibility>1.0e7&&sumInfeasibility>100.0*lastSumInfeasibility
808         &&factorization_->pivotTolerance()<0.11)||(largestPrimalError_>1.0e10&&largestDualError_>1.0e10)) {
809      problemStatus_=tentativeStatus;
810      doFactorization=true;
811      if (numberPivots) {
812        // go back
813        numberThrownOut=-1;
814        // trouble - restore solution
815        CoinMemcpyN(saveStatus_,numberColumns_+numberRows_,status_);
816        CoinMemcpyN(savedSolution_+numberColumns_ ,
817                    numberRows_,rowActivityWork_);
818        CoinMemcpyN(savedSolution_ ,
819           numberColumns_,columnActivityWork_);
820        // restore extra stuff
821        matrix_->generalExpanded(this,6,dummy);
822        forceFactorization_=1; // a bit drastic but ..
823        // Go to safe
824        factorization_->pivotTolerance(0.99);
825        pivotRow_=-1; // say no weights update
826        changeMade_++; // say change made
827        if (numberPivots==1) {
828          // throw out something
829          if (sequenceIn_>=0&&getStatus(sequenceIn_)!=basic) {
830            setFlagged(sequenceIn_);
831          } else if (sequenceOut_>=0&&getStatus(sequenceOut_)!=basic) {
832            setFlagged(sequenceOut_);
833          }
834        }
835        numberPivots=0;
836      }
837    }
838  }
839  // Double check reduced costs if no action
840  if (progress->lastIterationNumber(0)==numberIterations_) {
841    if (primalColumnPivot_->looksOptimal()) {
842      numberDualInfeasibilities_ = 0;
843      sumDualInfeasibilities_ = 0.0;
844    }
845  }
846  // If in primal and small dj give up
847  if ((specialOptions_&1024)!=0&&!numberPrimalInfeasibilities_&&numberDualInfeasibilities_) {
848    double average = sumDualInfeasibilities_/((double) numberDualInfeasibilities_);
849    if (numberIterations_>300&&average<1.0e-4) {
850      numberDualInfeasibilities_ = 0;
851      sumDualInfeasibilities_ = 0.0;
852    }
853  }
854  // Check if looping
855  int loop;
856  if (type!=2&&!ifValuesPass) 
857    loop = progress->looping();
858  else
859    loop=-1;
860  if (loop>=0) {
861    if (!problemStatus_) {
862      // declaring victory
863      numberPrimalInfeasibilities_ = 0;
864      sumPrimalInfeasibilities_ = 0.0;
865    } else {
866      problemStatus_ = loop; //exit if in loop
867      problemStatus_ = 10; // instead - try other algorithm
868      numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
869    }
870    problemStatus_ = 10; // instead - try other algorithm
871    return ;
872  } else if (loop<-1) {
873    // Is it time for drastic measures
874    if (nonLinearCost_->numberInfeasibilities()&&progress->badTimes()>5&&
875        progress->oddState()<10&&progress->oddState()>=0) {
876      progress->newOddState();
877      nonLinearCost_->zapCosts();
878    }
879    // something may have changed
880    gutsOfSolution(NULL,NULL,ifValuesPass!=0);
881  }
882  // If progress then reset costs
883  if (loop==-1&&!nonLinearCost_->numberInfeasibilities()&&progress->oddState()<0) {
884    createRim(4,false); // costs back
885    delete nonLinearCost_;
886    nonLinearCost_ = new ClpNonLinearCost(this);
887    progress->endOddState();
888    gutsOfSolution(NULL,NULL,ifValuesPass!=0);
889  }
890  // Flag to say whether to go to dual to clean up
891  bool goToDual=false;
892  // really for free variables in
893  //if((progressFlag_&2)!=0)
894  //problemStatus_=-1;;
895  progressFlag_ = 0; //reset progress flag
896
897  handler_->message(CLP_SIMPLEX_STATUS,messages_)
898    <<numberIterations_<<nonLinearCost_->feasibleReportCost();
899  handler_->printing(nonLinearCost_->numberInfeasibilities()>0)
900    <<nonLinearCost_->sumInfeasibilities()<<nonLinearCost_->numberInfeasibilities();
901  handler_->printing(sumDualInfeasibilities_>0.0)
902    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
903  handler_->printing(numberDualInfeasibilitiesWithoutFree_
904                     <numberDualInfeasibilities_)
905                       <<numberDualInfeasibilitiesWithoutFree_;
906  handler_->message()<<CoinMessageEol;
907  if (!primalFeasible()) {
908    nonLinearCost_->checkInfeasibilities(primalTolerance_);
909    gutsOfSolution(NULL,NULL,ifValuesPass!=0);
910    nonLinearCost_->checkInfeasibilities(primalTolerance_);
911  }
912  if (nonLinearCost_->numberInfeasibilities()>0&&!progress->initialWeight_&&!ifValuesPass&&infeasibilityCost_==1.0e10) {
913    // first time infeasible - start up weight computation
914    double * oldDj = dj_;
915    double * oldCost = cost_;
916    int numberRows2 = numberRows_+numberExtraRows_;
917    int numberTotal = numberRows2+numberColumns_;
918    dj_ = new double[numberTotal];
919    cost_ = new double[numberTotal];
920    reducedCostWork_ = dj_;
921    rowReducedCost_ = dj_+numberColumns_;
922    objectiveWork_ = cost_;
923    rowObjectiveWork_ = cost_+numberColumns_;
924    double direction = optimizationDirection_*objectiveScale_;
925    const double * obj = objective();
926    memset(rowObjectiveWork_,0,numberRows_*sizeof(double));
927    int iSequence;
928    if (columnScale_)
929      for (iSequence=0;iSequence<numberColumns_;iSequence++) 
930        cost_[iSequence] = obj[iSequence]*direction*columnScale_[iSequence];
931    else
932      for (iSequence=0;iSequence<numberColumns_;iSequence++) 
933        cost_[iSequence] = obj[iSequence]*direction;
934    computeDuals(NULL);
935    int numberSame=0;
936    int numberDifferent=0;
937    int numberZero=0;
938    int numberFreeSame=0;
939    int numberFreeDifferent=0;
940    int numberFreeZero=0;
941    int n=0;
942    for (iSequence=0;iSequence<numberTotal;iSequence++) {
943      if (getStatus(iSequence) != basic&&!flagged(iSequence)) {
944        // not basic
945        double distanceUp = upper_[iSequence]-solution_[iSequence];
946        double distanceDown = solution_[iSequence]-lower_[iSequence];
947        double feasibleDj = dj_[iSequence];
948        double infeasibleDj = oldDj[iSequence]-feasibleDj;
949        double value = feasibleDj*infeasibleDj;
950        if (distanceUp>primalTolerance_) {
951          // Check if "free"
952          if (distanceDown>primalTolerance_) {
953            // free
954            if (value>dualTolerance_) {
955              numberFreeSame++;
956            } else if(value<-dualTolerance_) {
957              numberFreeDifferent++;
958              dj_[n++] = feasibleDj/infeasibleDj;
959            } else {
960              numberFreeZero++;
961            }
962          } else {
963            // should not be negative
964            if (value>dualTolerance_) {
965              numberSame++;
966            } else if(value<-dualTolerance_) {
967              numberDifferent++;
968              dj_[n++] = feasibleDj/infeasibleDj;
969            } else {
970              numberZero++;
971            }
972          }
973        } else if (distanceDown>primalTolerance_) {
974          // should not be positive
975          if (value>dualTolerance_) {
976              numberSame++;
977            } else if(value<-dualTolerance_) {
978              numberDifferent++;
979              dj_[n++] = feasibleDj/infeasibleDj;
980            } else {
981              numberZero++;
982            }
983        }
984      }
985      progress->initialWeight_=-1.0;
986    }
987    //printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
988    //   numberSame,numberDifferent,numberZero,
989    //   numberFreeSame,numberFreeDifferent,numberFreeZero);
990    // we want most to be same
991    if (n) {
992      double most = 0.95;
993      std::sort(dj_,dj_+n);
994      int which = (int) ((1.0-most)*((double) n));
995      double take = -dj_[which]*infeasibilityCost_;
996      //printf("XXXXZ inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take,-dj_[0]*infeasibilityCost_,-dj_[n-1]*infeasibilityCost_);
997      take = -dj_[0]*infeasibilityCost_;
998      infeasibilityCost_ = CoinMin(CoinMax(1000.0*take,1.0e8),1.0000001e10);;
999      //printf("XXXX increasing weight to %g\n",infeasibilityCost_);
1000    }
1001    delete [] dj_;
1002    delete [] cost_;
1003    dj_= oldDj;
1004    cost_ = oldCost;
1005    reducedCostWork_ = dj_;
1006    rowReducedCost_ = dj_+numberColumns_;
1007    objectiveWork_ = cost_;
1008    rowObjectiveWork_ = cost_+numberColumns_;
1009    if (n)
1010      gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1011  }
1012  double trueInfeasibility =nonLinearCost_->sumInfeasibilities();
1013  if (!nonLinearCost_->numberInfeasibilities()&&infeasibilityCost_==1.0e10&&!ifValuesPass&&true) {
1014    // relax if default
1015    infeasibilityCost_ = CoinMin(CoinMax(100.0*sumDualInfeasibilities_,1.0e8),1.00000001e10);
1016    // reset looping criterion
1017    progress->reset();
1018    trueInfeasibility = 1.123456e10;
1019  }
1020  if (trueInfeasibility>1.0) {
1021    // If infeasibility going up may change weights
1022    double testValue = trueInfeasibility-1.0e-4*(10.0+trueInfeasibility);
1023    double lastInf = progress->lastInfeasibility(1);
1024    double lastInf3 = progress->lastInfeasibility(3);
1025    double thisObj = progress->lastObjective(0);
1026    double thisInf = progress->lastInfeasibility(0);
1027    thisObj += infeasibilityCost_*thisInf;
1028    double lastObj = progress->lastObjective(1);
1029    lastObj += infeasibilityCost_*lastInf;
1030    double lastObj3 = progress->lastObjective(3);
1031    lastObj3 += infeasibilityCost_*lastInf3;
1032    if (lastObj<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj))-1.0e-7
1033        &&firstFree_<0) {
1034      if (handler_->logLevel()==63)
1035        printf("lastobj %g this %g force %d ",lastObj,thisObj,forceFactorization_);
1036      int maxFactor = factorization_->maximumPivots();
1037      if (maxFactor>10) {
1038        if (forceFactorization_<0)
1039          forceFactorization_= maxFactor;
1040        forceFactorization_ = CoinMax(1,(forceFactorization_>>2));
1041        if (handler_->logLevel()==63)
1042          printf("Reducing factorization frequency to %d\n",forceFactorization_);
1043      }
1044    } else if (lastObj3<thisObj-1.0e-5*CoinMax(fabs(thisObj),fabs(lastObj3))-1.0e-7
1045        &&firstFree_<0) {
1046      if (handler_->logLevel()==63)
1047        printf("lastobj3 %g this3 %g `force %d ",lastObj3,thisObj,forceFactorization_);
1048      int maxFactor = factorization_->maximumPivots();
1049      if (maxFactor>10) {
1050        if (forceFactorization_<0)
1051          forceFactorization_= maxFactor;
1052        forceFactorization_ = CoinMax(1,(forceFactorization_*2)/3);
1053        if (handler_->logLevel()==63)
1054        printf("Reducing factorization frequency to %d\n",forceFactorization_);
1055      }
1056    } else if(lastInf<testValue||trueInfeasibility==1.123456e10) {
1057      if (infeasibilityCost_<1.0e14) {
1058        infeasibilityCost_ *= 1.5;
1059        // reset looping criterion
1060        progress->reset();
1061        if (handler_->logLevel()==63)
1062          printf("increasing weight to %g\n",infeasibilityCost_);
1063        gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1064      }
1065    }
1066  }
1067  // we may wish to say it is optimal even if infeasible
1068  bool alwaysOptimal = (specialOptions_&1)!=0;
1069  // give code benefit of doubt
1070  if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
1071      sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1072    // say optimal (with these bounds etc)
1073    numberDualInfeasibilities_ = 0;
1074    sumDualInfeasibilities_ = 0.0;
1075    numberPrimalInfeasibilities_ = 0;
1076    sumPrimalInfeasibilities_ = 0.0;
1077    // But check if in sprint
1078    if (originalModel) {
1079      // Carry on and re-do
1080      numberDualInfeasibilities_ = -776;
1081    }
1082    // But if real primal infeasibilities nonzero carry on
1083    if (nonLinearCost_->numberInfeasibilities()) {
1084      // most likely to happen if infeasible
1085      double relaxedToleranceP=primalTolerance_;
1086      // we can't really trust infeasibilities if there is primal error
1087      double error = CoinMin(1.0e-2,largestPrimalError_);
1088      // allow tolerance at least slightly bigger than standard
1089      relaxedToleranceP = relaxedToleranceP +  error;
1090      int ninfeas = nonLinearCost_->numberInfeasibilities();
1091      double sum = nonLinearCost_->sumInfeasibilities();
1092      double average = sum/ ((double) ninfeas);
1093#ifdef COIN_DEVELOP
1094      if (handler_->logLevel()>0)
1095        printf("nonLinearCost says infeasible %d summing to %g\n",
1096               ninfeas,sum);
1097#endif
1098      if (average>relaxedToleranceP) {
1099        sumOfRelaxedPrimalInfeasibilities_ = sum;
1100        numberPrimalInfeasibilities_ = ninfeas;
1101        sumPrimalInfeasibilities_ = sum;
1102#ifdef COIN_DEVELOP
1103        bool unflagged = 
1104#endif
1105        unflag();
1106#ifdef COIN_DEVELOP
1107        if (unflagged&&handler_->logLevel()>0)
1108          printf(" - but flagged variables\n");
1109#endif
1110      }
1111    }
1112  } 
1113  // had ||(type==3&&problemStatus_!=-5) -- ??? why ????
1114  if ((dualFeasible()||problemStatus_==-4)&&!ifValuesPass) {
1115    // see if extra helps
1116    if (nonLinearCost_->numberInfeasibilities()&&
1117         (nonLinearCost_->sumInfeasibilities()>1.0e-3||sumOfRelaxedPrimalInfeasibilities_)
1118        &&!alwaysOptimal) {
1119      //may need infeasiblity cost changed
1120      // we can see if we can construct a ray
1121      // make up a new objective
1122      double saveWeight = infeasibilityCost_;
1123      // save nonlinear cost as we are going to switch off costs
1124      ClpNonLinearCost * nonLinear = nonLinearCost_;
1125      // do twice to make sure Primal solution has settled
1126      // put non-basics to bounds in case tolerance moved
1127      // put back original costs
1128      createRim(4);
1129      nonLinearCost_->checkInfeasibilities(0.0);
1130      gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1131
1132      infeasibilityCost_=1.0e100;
1133      // put back original costs
1134      createRim(4);
1135      nonLinearCost_->checkInfeasibilities(primalTolerance_);
1136      // may have fixed infeasibilities - double check
1137      if (nonLinearCost_->numberInfeasibilities()==0) {
1138        // carry on
1139        problemStatus_ = -1;
1140        infeasibilityCost_=saveWeight;
1141        nonLinearCost_->checkInfeasibilities(primalTolerance_);
1142      } else {
1143        nonLinearCost_=NULL;
1144        // scale
1145        int i;
1146        for (i=0;i<numberRows_+numberColumns_;i++) 
1147          cost_[i] *= 1.0e-95;
1148        gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1149        nonLinearCost_=nonLinear;
1150        infeasibilityCost_=saveWeight;
1151        if ((infeasibilityCost_>=1.0e18||
1152             numberDualInfeasibilities_==0)&&perturbation_==101) {
1153          goToDual=unPerturb(); // stop any further perturbation
1154          if (nonLinearCost_->sumInfeasibilities()>1.0e-1)
1155            goToDual=false;
1156          nonLinearCost_->checkInfeasibilities(primalTolerance_);
1157          numberDualInfeasibilities_=1; // carry on
1158          problemStatus_=-1;
1159        }
1160        if (infeasibilityCost_>=1.0e20||
1161            numberDualInfeasibilities_==0) {
1162          // we are infeasible - use as ray
1163          delete [] ray_;
1164          ray_ = new double [numberRows_];
1165          CoinMemcpyN(dual_,numberRows_,ray_);
1166          // and get feasible duals
1167          infeasibilityCost_=0.0;
1168          createRim(4);
1169          nonLinearCost_->checkInfeasibilities(primalTolerance_);
1170          gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1171          // so will exit
1172          infeasibilityCost_=1.0e30;
1173          // reset infeasibilities
1174          sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();;
1175          numberPrimalInfeasibilities_=
1176            nonLinearCost_->numberInfeasibilities();
1177        }
1178        if (infeasibilityCost_<1.0e20) {
1179          infeasibilityCost_ *= 5.0;
1180          // reset looping criterion
1181          progress->reset();
1182          changeMade_++; // say change made
1183          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
1184            <<infeasibilityCost_
1185            <<CoinMessageEol;
1186          // put back original costs and then check
1187          createRim(4);
1188          nonLinearCost_->checkInfeasibilities(0.0);
1189          gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1190          problemStatus_=-1; //continue
1191          goToDual=false;
1192        } else {
1193          // say infeasible
1194          problemStatus_ = 1;
1195        }
1196      }
1197    } else {
1198      // may be optimal
1199      if (perturbation_==101) {
1200        goToDual=unPerturb(); // stop any further perturbation
1201        if (numberRows_>20000&&!numberTimesOptimal_)
1202          goToDual=0; // Better to carry on a bit longer
1203        lastCleaned=-1; // carry on
1204      }
1205      bool unflagged = (unflag()!=0);
1206      if ( lastCleaned!=numberIterations_||unflagged) {
1207        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
1208          <<primalTolerance_
1209          <<CoinMessageEol;
1210        if (numberTimesOptimal_<4) {
1211          numberTimesOptimal_++;
1212          changeMade_++; // say change made
1213          if (numberTimesOptimal_==1) {
1214            // better to have small tolerance even if slower
1215            factorization_->zeroTolerance(1.0e-15);
1216          }
1217          lastCleaned=numberIterations_;
1218          if (primalTolerance_!=dblParam_[ClpPrimalTolerance])
1219            handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
1220              <<CoinMessageEol;
1221          double oldTolerance = primalTolerance_;
1222          primalTolerance_=dblParam_[ClpPrimalTolerance];
1223#if 0
1224          double * xcost = new double[numberRows_+numberColumns_];
1225          double * xlower = new double[numberRows_+numberColumns_];
1226          double * xupper = new double[numberRows_+numberColumns_];
1227          double * xdj = new double[numberRows_+numberColumns_];
1228          double * xsolution = new double[numberRows_+numberColumns_];
1229          memcpy(xcost,cost_,(numberRows_+numberColumns_)*sizeof(double));
1230          memcpy(xlower,lower_,(numberRows_+numberColumns_)*sizeof(double));
1231          memcpy(xupper,upper_,(numberRows_+numberColumns_)*sizeof(double));
1232          memcpy(xdj,dj_,(numberRows_+numberColumns_)*sizeof(double));
1233          memcpy(xsolution,solution_,(numberRows_+numberColumns_)*sizeof(double));
1234#endif
1235          // put back original costs and then check
1236          createRim(4);
1237          nonLinearCost_->checkInfeasibilities(oldTolerance);
1238#if 0
1239          int i;
1240          for (i=0;i<numberRows_+numberColumns_;i++) {
1241            if (cost_[i]!=xcost[i])
1242              printf("** %d old cost %g new %g sol %g\n",
1243                     i,xcost[i],cost_[i],solution_[i]);
1244            if (lower_[i]!=xlower[i])
1245              printf("** %d old lower %g new %g sol %g\n",
1246                     i,xlower[i],lower_[i],solution_[i]);
1247            if (upper_[i]!=xupper[i])
1248              printf("** %d old upper %g new %g sol %g\n",
1249                     i,xupper[i],upper_[i],solution_[i]);
1250            if (dj_[i]!=xdj[i])
1251              printf("** %d old dj %g new %g sol %g\n",
1252                     i,xdj[i],dj_[i],solution_[i]);
1253            if (solution_[i]!=xsolution[i])
1254              printf("** %d old solution %g new %g sol %g\n",
1255                     i,xsolution[i],solution_[i],solution_[i]);
1256          }
1257          delete [] xcost;
1258          delete [] xupper;
1259          delete [] xlower;
1260          delete [] xdj;
1261          delete [] xsolution;
1262#endif
1263          gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1264          if (sumOfRelaxedDualInfeasibilities_ == 0.0&&
1265              sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1266            // say optimal (with these bounds etc)
1267            numberDualInfeasibilities_ = 0;
1268            sumDualInfeasibilities_ = 0.0;
1269            numberPrimalInfeasibilities_ = 0;
1270            sumPrimalInfeasibilities_ = 0.0;
1271          }
1272          if (dualFeasible()&&!nonLinearCost_->numberInfeasibilities()&&lastCleaned>=0)
1273            problemStatus_=0;
1274          else
1275            problemStatus_ = -1;
1276        } else {
1277          problemStatus_=0; // optimal
1278          if (lastCleaned<numberIterations_) {
1279            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
1280              <<CoinMessageEol;
1281          }
1282        }
1283      } else {
1284        problemStatus_=0; // optimal
1285      }
1286    }
1287  } else {
1288    // see if looks unbounded
1289    if (problemStatus_==-5) {
1290      if (nonLinearCost_->numberInfeasibilities()) {
1291        if (infeasibilityCost_>1.0e18&&perturbation_==101) {
1292          // back off weight
1293          infeasibilityCost_ = 1.0e13;
1294          // reset looping criterion
1295          progress->reset();
1296          unPerturb(); // stop any further perturbation
1297        }
1298        //we need infeasiblity cost changed
1299        if (infeasibilityCost_<1.0e20) {
1300          infeasibilityCost_ *= 5.0;
1301          // reset looping criterion
1302          progress->reset();
1303          changeMade_++; // say change made
1304          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
1305            <<infeasibilityCost_
1306            <<CoinMessageEol;
1307          // put back original costs and then check
1308          createRim(4);
1309          gutsOfSolution(NULL,NULL,ifValuesPass!=0);
1310          problemStatus_=-1; //continue
1311        } else {
1312          // say infeasible
1313          problemStatus_ = 1;
1314          // we are infeasible - use as ray
1315          delete [] ray_;
1316          ray_ = new double [numberRows_];
1317          CoinMemcpyN(dual_,numberRows_,ray_);
1318        }
1319      } else {
1320        // say unbounded
1321        problemStatus_ = 2;
1322      } 
1323    } else {
1324      // carry on
1325      problemStatus_ = -1;
1326      if(type==3&&problemStatus_!=-5) {
1327        //bool unflagged =
1328        unflag();
1329        if (sumDualInfeasibilities_<1.0e-3||
1330            (sumDualInfeasibilities_/(double) numberDualInfeasibilities_)<1.0e-5||
1331            progress->lastIterationNumber(0)==numberIterations_) {
1332          if (!numberPrimalInfeasibilities_) {
1333            if (numberTimesOptimal_<4) {
1334              numberTimesOptimal_++;
1335              changeMade_++; // say change made
1336            } else {
1337              problemStatus_=0;
1338              secondaryStatus_=5;
1339            }
1340          }
1341        }
1342      }
1343    }
1344  }
1345  // save extra stuff
1346  matrix_->generalExpanded(this,5,dummy);
1347  if (type==0||type==1) {
1348    if (type!=1||!saveStatus_) {
1349      // create save arrays
1350      delete [] saveStatus_;
1351      delete [] savedSolution_;
1352      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
1353      savedSolution_ = new double [numberRows_+numberColumns_];
1354    }
1355    // save arrays
1356    CoinMemcpyN(status_,numberColumns_+numberRows_,saveStatus_);
1357    CoinMemcpyN(rowActivityWork_,
1358           numberRows_,savedSolution_+numberColumns_);
1359    CoinMemcpyN(columnActivityWork_,numberColumns_,savedSolution_);
1360  }
1361  // see if in Cbc etc
1362  bool inCbcOrOther = (specialOptions_&0x03000000)!=0;
1363  bool disaster=false;
1364  if (disasterArea_&&inCbcOrOther&&disasterArea_->check()) {
1365    disasterArea_->saveInfo();
1366    disaster=true;
1367  }
1368  if (disaster)
1369    problemStatus_=3;
1370  if (doFactorization) {
1371    // restore weights (if saved) - also recompute infeasibility list
1372    if (tentativeStatus>-3) 
1373      primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
1374    else
1375      primalColumnPivot_->saveWeights(this,3);
1376    if (saveThreshold) {
1377      // use default at present
1378      factorization_->sparseThreshold(0);
1379      factorization_->goSparse();
1380    }
1381  }
1382  if (problemStatus_<0&&!changeMade_) {
1383    problemStatus_=4; // unknown
1384  }
1385  lastGoodIteration_ = numberIterations_;
1386  if (goToDual) 
1387    problemStatus_=10; // try dual
1388  // make sure first free monotonic
1389  if (firstFree_>=0&&saveFirstFree>=0) {
1390    firstFree_=saveFirstFree;
1391    nextSuperBasic(1,NULL);
1392  }
1393  // Allow matrices to be sorted etc
1394  int fake=-999; // signal sort
1395  matrix_->correctSequence(this,fake,fake);
1396}
1397/*
1398   Row array has pivot column
1399   This chooses pivot row.
1400   For speed, we may need to go to a bucket approach when many
1401   variables go through bounds
1402   On exit rhsArray will have changes in costs of basic variables
1403*/
1404void 
1405ClpSimplexPrimal::primalRow(CoinIndexedVector * rowArray,
1406                            CoinIndexedVector * rhsArray,
1407                            CoinIndexedVector * spareArray,
1408                            CoinIndexedVector * spareArray2,
1409                            int valuesPass)
1410{
1411  double saveDj = dualIn_;
1412  if (valuesPass&&objective_->type()<2) {
1413    dualIn_ = cost_[sequenceIn_];
1414
1415    double * work=rowArray->denseVector();
1416    int number=rowArray->getNumElements();
1417    int * which=rowArray->getIndices();
1418
1419    int iIndex;
1420    for (iIndex=0;iIndex<number;iIndex++) {
1421     
1422      int iRow = which[iIndex];
1423      double alpha = work[iIndex];
1424      int iPivot=pivotVariable_[iRow];
1425      dualIn_ -= alpha*cost_[iPivot];
1426    }
1427    // determine direction here
1428    if (dualIn_<-dualTolerance_) {
1429      directionIn_=1;
1430    } else if (dualIn_>dualTolerance_) {
1431      directionIn_=-1;
1432    } else {
1433      // towards nearest bound
1434      if (valueIn_-lowerIn_<upperIn_-valueIn_) {
1435        directionIn_=-1;
1436        dualIn_=dualTolerance_;
1437      } else {
1438        directionIn_=1;
1439        dualIn_=-dualTolerance_;
1440      }
1441    }
1442  }
1443
1444  // sequence stays as row number until end
1445  pivotRow_=-1;
1446  int numberRemaining=0;
1447
1448  double totalThru=0.0; // for when variables flip
1449  // Allow first few iterations to take tiny
1450  double acceptablePivot=1.0e-1*acceptablePivot_;
1451  if (numberIterations_>100)
1452    acceptablePivot=acceptablePivot_;
1453  if (factorization_->pivots()>10)
1454    acceptablePivot=1.0e+3*acceptablePivot_; // if we have iterated be more strict
1455  else if (factorization_->pivots()>5)
1456    acceptablePivot=1.0e+2*acceptablePivot_; // if we have iterated be slightly more strict
1457  else if (factorization_->pivots())
1458    acceptablePivot=acceptablePivot_; // relax
1459  double bestEverPivot=acceptablePivot;
1460  int lastPivotRow = -1;
1461  double lastPivot=0.0;
1462  double lastTheta=1.0e50;
1463
1464  // use spareArrays to put ones looked at in
1465  // First one is list of candidates
1466  // We could compress if we really know we won't need any more
1467  // Second array has current set of pivot candidates
1468  // with a backup list saved in double * part of indexed vector
1469
1470  // pivot elements
1471  double * spare;
1472  // indices
1473  int * index;
1474  spareArray->clear();
1475  spare = spareArray->denseVector();
1476  index = spareArray->getIndices();
1477
1478  // we also need somewhere for effective rhs
1479  double * rhs=rhsArray->denseVector();
1480  // and we can use indices to point to alpha
1481  // that way we can store fabs(alpha)
1482  int * indexPoint = rhsArray->getIndices();
1483  //int numberFlip=0; // Those which may change if flips
1484
1485  /*
1486    First we get a list of possible pivots.  We can also see if the
1487    problem looks unbounded.
1488
1489    At first we increase theta and see what happens.  We start
1490    theta at a reasonable guess.  If in right area then we do bit by bit.
1491    We save possible pivot candidates
1492
1493   */
1494
1495  // do first pass to get possibles
1496  // We can also see if unbounded
1497
1498  double * work=rowArray->denseVector();
1499  int number=rowArray->getNumElements();
1500  int * which=rowArray->getIndices();
1501
1502  // we need to swap sign if coming in from ub
1503  double way = directionIn_;
1504  double maximumMovement;
1505  if (way>0.0) 
1506    maximumMovement = CoinMin(1.0e30,upperIn_-valueIn_);
1507  else
1508    maximumMovement = CoinMin(1.0e30,valueIn_-lowerIn_);
1509
1510  double averageTheta = nonLinearCost_->averageTheta();
1511  double tentativeTheta = CoinMin(10.0*averageTheta,maximumMovement);
1512  double upperTheta = maximumMovement;
1513  if (tentativeTheta>0.5*maximumMovement)
1514    tentativeTheta=maximumMovement;
1515  bool thetaAtMaximum=tentativeTheta==maximumMovement;
1516  // In case tiny bounds increase
1517  if (maximumMovement<1.0)
1518    tentativeTheta *= 1.1;
1519  double dualCheck = fabs(dualIn_);
1520  // but make a bit more pessimistic
1521  dualCheck=CoinMax(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
1522
1523  int iIndex;
1524  int pivotOne=-1;
1525  //#define CLP_DEBUG
1526#ifdef CLP_DEBUG
1527  if (numberIterations_==-3839||numberIterations_==-3840) {
1528    double dj=cost_[sequenceIn_];
1529    printf("cost in on %d is %g, dual in %g\n",sequenceIn_,dj,dualIn_);
1530    for (iIndex=0;iIndex<number;iIndex++) {
1531
1532      int iRow = which[iIndex];
1533      double alpha = work[iIndex];
1534      int iPivot=pivotVariable_[iRow];
1535      dj -= alpha*cost_[iPivot];
1536      printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n",
1537             iRow,iPivot,lower_[iPivot],solution_[iPivot],upper_[iPivot],
1538             alpha, solution_[iPivot]-1.0e9*alpha,cost_[iPivot],dj);
1539    }
1540  }
1541#endif
1542  while (true) {
1543    pivotOne=-1;
1544    totalThru=0.0;
1545    // We also re-compute reduced cost
1546    numberRemaining=0;
1547    dualIn_ = cost_[sequenceIn_];
1548#ifndef NDEBUG
1549    double tolerance = primalTolerance_*1.002;
1550#endif
1551    for (iIndex=0;iIndex<number;iIndex++) {
1552
1553      int iRow = which[iIndex];
1554      double alpha = work[iIndex];
1555      int iPivot=pivotVariable_[iRow];
1556      if (cost_[iPivot])
1557        dualIn_ -= alpha*cost_[iPivot];
1558      alpha *= way;
1559      double oldValue = solution_[iPivot];
1560      // get where in bound sequence
1561      // note that after this alpha is actually fabs(alpha)
1562      bool possible;
1563      // do computation same way as later on in primal
1564      if (alpha>0.0) {
1565        // basic variable going towards lower bound
1566        double bound = lower_[iPivot];
1567        possible = (oldValue-tentativeTheta*alpha)<=bound+primalTolerance_;
1568        oldValue -= bound;
1569      } else {
1570        // basic variable going towards upper bound
1571        double bound = upper_[iPivot];
1572        possible = (oldValue-tentativeTheta*alpha)>=bound-primalTolerance_;;
1573        oldValue = bound-oldValue;
1574        alpha = - alpha;
1575      }
1576      double value;
1577      assert (oldValue>=-tolerance);
1578      if (possible) {
1579        value=oldValue-upperTheta*alpha;
1580        if (value<-primalTolerance_&&alpha>=acceptablePivot) {
1581          upperTheta = (oldValue+primalTolerance_)/alpha;
1582          pivotOne=numberRemaining;
1583        }
1584        // add to list
1585        spare[numberRemaining]=alpha;
1586        rhs[numberRemaining]=oldValue;
1587        indexPoint[numberRemaining]=iIndex;
1588        index[numberRemaining++]=iRow;
1589        totalThru += alpha;
1590        setActive(iRow);
1591        //} else if (value<primalTolerance_*1.002) {
1592        // May change if is a flip
1593        //indexRhs[numberFlip++]=iRow;
1594      }
1595    }
1596    if (upperTheta<maximumMovement&&totalThru*infeasibilityCost_>=1.0001*dualCheck) {
1597      // Can pivot here
1598      break;
1599    } else if (!thetaAtMaximum) {
1600      //printf("Going round with average theta of %g\n",averageTheta);
1601      tentativeTheta=maximumMovement;
1602      thetaAtMaximum=true; // seems to be odd compiler error
1603    } else {
1604      break;
1605    }
1606  }
1607  totalThru=0.0;
1608
1609  theta_=maximumMovement;
1610
1611  bool goBackOne = false;
1612  if (objective_->type()>1) 
1613    dualIn_=saveDj;
1614
1615  //printf("%d remain out of %d\n",numberRemaining,number);
1616  int iTry=0;
1617#define MAXTRY 1000
1618  if (numberRemaining&&upperTheta<maximumMovement) {
1619    // First check if previously chosen one will work
1620    if (pivotOne>=0&&0) {
1621      double thruCost = infeasibilityCost_*spare[pivotOne];
1622      if (thruCost>=0.99*fabs(dualIn_))
1623        printf("Could pivot on %d as change %g dj %g\n",
1624               index[pivotOne],thruCost,dualIn_);
1625      double alpha = spare[pivotOne];
1626      double oldValue = rhs[pivotOne];
1627      theta_ = oldValue/alpha;
1628      pivotRow_=pivotOne;
1629      // Stop loop
1630      iTry=MAXTRY;
1631    }
1632
1633    // first get ratio with tolerance
1634    for ( ;iTry<MAXTRY;iTry++) {
1635     
1636      upperTheta=maximumMovement;
1637      int iBest=-1;
1638      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1639       
1640        double alpha = spare[iIndex];
1641        double oldValue = rhs[iIndex];
1642        double value = oldValue-upperTheta*alpha;
1643       
1644        if (value<-primalTolerance_ && alpha>=acceptablePivot) {
1645          upperTheta = (oldValue+primalTolerance_)/alpha;
1646          iBest=iIndex; // just in case weird numbers
1647        }
1648      }
1649     
1650      // now look at best in this lot
1651      // But also see how infeasible small pivots will make
1652      double sumInfeasibilities=0.0;
1653      double bestPivot=acceptablePivot;
1654      pivotRow_=-1;
1655      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1656       
1657        int iRow = index[iIndex];
1658        double alpha = spare[iIndex];
1659        double oldValue = rhs[iIndex];
1660        double value = oldValue-upperTheta*alpha;
1661       
1662        if (value<=0||iBest==iIndex) {
1663          // how much would it cost to go thru and modify bound
1664          double trueAlpha=way*work[indexPoint[iIndex]];
1665          totalThru += nonLinearCost_->changeInCost(pivotVariable_[iRow],trueAlpha,rhs[iIndex]);
1666          setActive(iRow);
1667          if (alpha>bestPivot) {
1668            bestPivot=alpha;
1669            theta_ = oldValue/bestPivot;
1670            pivotRow_=iIndex;
1671          } else if (alpha<acceptablePivot) {
1672            if (value<-primalTolerance_)
1673              sumInfeasibilities += -value-primalTolerance_;
1674          }
1675        }
1676      }
1677      if (bestPivot<0.1*bestEverPivot&&
1678          bestEverPivot>1.0e-6&& bestPivot<1.0e-3) {
1679        // back to previous one
1680        goBackOne = true;
1681        break;
1682      } else if (pivotRow_==-1&&upperTheta>largeValue_) {
1683        if (lastPivot>acceptablePivot) {
1684          // back to previous one
1685          goBackOne = true;
1686        } else {
1687          // can only get here if all pivots so far too small
1688        }
1689        break;
1690      } else if (totalThru>=dualCheck) {
1691        if (sumInfeasibilities>primalTolerance_&&!nonLinearCost_->numberInfeasibilities()) {
1692          // Looks a bad choice
1693          if (lastPivot>acceptablePivot) {
1694            goBackOne=true;
1695          } else {
1696            // say no good
1697            dualIn_=0.0;
1698          }
1699        } 
1700        break; // no point trying another loop
1701      } else {
1702        lastPivotRow=pivotRow_;
1703        lastTheta = theta_;
1704        if (bestPivot>bestEverPivot)
1705          bestEverPivot=bestPivot;
1706      }   
1707    }
1708    // can get here without pivotRow_ set but with lastPivotRow
1709    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
1710      // back to previous one
1711      pivotRow_=lastPivotRow;
1712      theta_ = lastTheta;
1713    }
1714  } else if (pivotRow_<0&&maximumMovement>1.0e20) {
1715    // looks unbounded
1716    valueOut_=COIN_DBL_MAX; // say odd
1717    if (nonLinearCost_->numberInfeasibilities()) {
1718      // but infeasible??
1719      // move variable but don't pivot
1720      tentativeTheta=1.0e50;
1721      for (iIndex=0;iIndex<number;iIndex++) {
1722        int iRow = which[iIndex];
1723        double alpha = work[iIndex];
1724        int iPivot=pivotVariable_[iRow];
1725        alpha *= way;
1726        double oldValue = solution_[iPivot];
1727        // get where in bound sequence
1728        // note that after this alpha is actually fabs(alpha)
1729        if (alpha>0.0) {
1730          // basic variable going towards lower bound
1731          double bound = lower_[iPivot];
1732          oldValue -= bound;
1733        } else {
1734          // basic variable going towards upper bound
1735          double bound = upper_[iPivot];
1736          oldValue = bound-oldValue;
1737          alpha = - alpha;
1738        }
1739        if (oldValue-tentativeTheta*alpha<0.0) {
1740          tentativeTheta = oldValue/alpha;
1741        }
1742      }
1743      // If free in then see if we can get to 0.0
1744      if (lowerIn_<-1.0e20&&upperIn_>1.0e20) {
1745        if (dualIn_*valueIn_>0.0) {
1746          if (fabs(valueIn_)<1.0e-2&&(tentativeTheta<fabs(valueIn_)||tentativeTheta>1.0e20)) {
1747            tentativeTheta = fabs(valueIn_);
1748          }
1749        }
1750      }
1751      if (tentativeTheta<1.0e10)
1752        valueOut_=valueIn_+way*tentativeTheta;
1753    }
1754  }
1755  //if (iTry>50)
1756  //printf("** %d tries\n",iTry);
1757  if (pivotRow_>=0) {
1758    int position=pivotRow_; // position in list
1759    pivotRow_=index[position];
1760    alpha_=work[indexPoint[position]];
1761    // translate to sequence
1762    sequenceOut_ = pivotVariable_[pivotRow_];
1763    valueOut_ = solution(sequenceOut_);
1764    lowerOut_=lower_[sequenceOut_];
1765    upperOut_=upper_[sequenceOut_];
1766#define MINIMUMTHETA 1.0e-12
1767    // Movement should be minimum for anti-degeneracy - unless
1768    // fixed variable out
1769    double minimumTheta;
1770    if (upperOut_>lowerOut_)
1771      minimumTheta=MINIMUMTHETA;
1772    else
1773      minimumTheta=0.0;
1774    // But can't go infeasible
1775    double distance;
1776    if (alpha_*way>0.0) 
1777      distance=valueOut_-lowerOut_;
1778    else
1779      distance=upperOut_-valueOut_;
1780    if (distance-minimumTheta*fabs(alpha_)<-primalTolerance_)
1781      minimumTheta = CoinMax(0.0,(distance+0.5*primalTolerance_)/fabs(alpha_));
1782    // will we need to increase tolerance
1783    //#define CLP_DEBUG
1784    double largestInfeasibility = primalTolerance_;
1785    if (theta_<minimumTheta&&(specialOptions_&4)==0&&!valuesPass) {
1786      theta_=minimumTheta;
1787      for (iIndex=0;iIndex<numberRemaining-numberRemaining;iIndex++) {
1788        largestInfeasibility = CoinMax(largestInfeasibility,
1789                                    -(rhs[iIndex]-spare[iIndex]*theta_));
1790      }
1791//#define CLP_DEBUG
1792#ifdef CLP_DEBUG
1793      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32)>-1)
1794        printf("Primal tolerance increased from %g to %g\n",
1795               primalTolerance_,largestInfeasibility);
1796#endif
1797//#undef CLP_DEBUG
1798      primalTolerance_ = CoinMax(primalTolerance_,largestInfeasibility);
1799    }
1800    // Need to look at all in some cases
1801    if (theta_>tentativeTheta) {
1802      for (iIndex=0;iIndex<number;iIndex++) 
1803        setActive(which[iIndex]);
1804    }
1805    if (way<0.0) 
1806      theta_ = - theta_;
1807    double newValue = valueOut_ - theta_*alpha_;
1808    // If  4 bit set - Force outgoing variables to exact bound (primal)
1809    if (alpha_*way<0.0) {
1810      directionOut_=-1;      // to upper bound
1811      if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
1812        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1813      } else {
1814        upperOut_ = newValue;
1815      }
1816    } else {
1817      directionOut_=1;      // to lower bound
1818      if (fabs(theta_)>1.0e-6||(specialOptions_&4)!=0) {
1819        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1820      } else {
1821        lowerOut_ = newValue;
1822      }
1823    }
1824    dualOut_ = reducedCost(sequenceOut_);
1825  } else if (maximumMovement<1.0e20) {
1826    // flip
1827    pivotRow_ = -2; // so we can tell its a flip
1828    sequenceOut_ = sequenceIn_;
1829    valueOut_ = valueIn_;
1830    dualOut_ = dualIn_;
1831    lowerOut_ = lowerIn_;
1832    upperOut_ = upperIn_;
1833    alpha_ = 0.0;
1834    if (way<0.0) {
1835      directionOut_=1;      // to lower bound
1836      theta_ = lowerOut_ - valueOut_;
1837    } else {
1838      directionOut_=-1;      // to upper bound
1839      theta_ = upperOut_ - valueOut_;
1840    }
1841  }
1842
1843  double theta1 = CoinMax(theta_,1.0e-12);
1844  double theta2 = numberIterations_*nonLinearCost_->averageTheta();
1845  // Set average theta
1846  nonLinearCost_->setAverageTheta((theta1+theta2)/((double) (numberIterations_+1)));
1847  //if (numberIterations_%1000==0)
1848  //printf("average theta is %g\n",nonLinearCost_->averageTheta());
1849
1850  // clear arrays
1851
1852  CoinZeroN(spare,numberRemaining);
1853
1854  // put back original bounds etc
1855  CoinMemcpyN(index,numberRemaining,rhsArray->getIndices());
1856  rhsArray->setNumElements(numberRemaining);
1857  rhsArray->setPacked();
1858  nonLinearCost_->goBackAll(rhsArray);
1859  rhsArray->clear();
1860
1861}
1862/*
1863   Chooses primal pivot column
1864   updateArray has cost updates (also use pivotRow_ from last iteration)
1865   Would be faster with separate region to scan
1866   and will have this (with square of infeasibility) when steepest
1867   For easy problems we can just choose one of the first columns we look at
1868*/
1869void 
1870ClpSimplexPrimal::primalColumn(CoinIndexedVector * updates,
1871                               CoinIndexedVector * spareRow1,
1872                               CoinIndexedVector * spareRow2,
1873                               CoinIndexedVector * spareColumn1,
1874                               CoinIndexedVector * spareColumn2)
1875{
1876  sequenceIn_ = primalColumnPivot_->pivotColumn(updates,spareRow1,
1877                                               spareRow2,spareColumn1,
1878                                               spareColumn2);
1879  if (sequenceIn_>=0) {
1880    valueIn_=solution_[sequenceIn_];
1881    dualIn_=dj_[sequenceIn_];
1882    if (nonLinearCost_->lookBothWays()) {
1883      // double check
1884      ClpSimplex::Status status = getStatus(sequenceIn_);
1885     
1886      switch(status) {
1887      case ClpSimplex::atUpperBound:
1888        if (dualIn_<0.0) {
1889          // move to other side
1890          printf("For %d U (%g, %g, %g) dj changed from %g",
1891                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
1892                 upper_[sequenceIn_],dualIn_);
1893          dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
1894          printf(" to %g\n",dualIn_);
1895          nonLinearCost_->setOne(sequenceIn_,upper_[sequenceIn_]+2.0*currentPrimalTolerance());
1896          setStatus(sequenceIn_,ClpSimplex::atLowerBound);
1897        }
1898        break;
1899      case ClpSimplex::atLowerBound:
1900        if (dualIn_>0.0) {
1901          // move to other side
1902          printf("For %d L (%g, %g, %g) dj changed from %g",
1903                 sequenceIn_,lower_[sequenceIn_],solution_[sequenceIn_],
1904                 upper_[sequenceIn_],dualIn_);
1905          dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
1906          printf(" to %g\n",dualIn_);
1907          nonLinearCost_->setOne(sequenceIn_,lower_[sequenceIn_]-2.0*currentPrimalTolerance());
1908          setStatus(sequenceIn_,ClpSimplex::atUpperBound);
1909        }
1910        break;
1911      default:
1912        break;
1913      }
1914    }
1915    lowerIn_=lower_[sequenceIn_];
1916    upperIn_=upper_[sequenceIn_];
1917    if (dualIn_>0.0)
1918      directionIn_ = -1;
1919    else 
1920      directionIn_ = 1;
1921  } else {
1922    sequenceIn_ = -1;
1923  }
1924}
1925/* The primals are updated by the given array.
1926   Returns number of infeasibilities.
1927   After rowArray will have list of cost changes
1928*/
1929int 
1930ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
1931                                        double theta,
1932                                        double & objectiveChange,
1933                                        int valuesPass)
1934{
1935  // Cost on pivot row may change - may need to change dualIn
1936  double oldCost=0.0;
1937  if (pivotRow_>=0)
1938    oldCost = cost_[sequenceOut_];
1939  //rowArray->scanAndPack();
1940  double * work=rowArray->denseVector();
1941  int number=rowArray->getNumElements();
1942  int * which=rowArray->getIndices();
1943
1944  int newNumber = 0;
1945  int pivotPosition = -1;
1946  nonLinearCost_->setChangeInCost(0.0);
1947  //printf("XX 4138 sol %g lower %g upper %g cost %g status %d\n",
1948  //   solution_[4138],lower_[4138],upper_[4138],cost_[4138],status_[4138]);
1949  // allow for case where bound+tolerance == bound
1950  //double tolerance = 0.999999*primalTolerance_;
1951  double relaxedTolerance = 1.001*primalTolerance_;
1952  int iIndex;
1953  if (!valuesPass) {
1954    for (iIndex=0;iIndex<number;iIndex++) {
1955     
1956      int iRow = which[iIndex];
1957      double alpha = work[iIndex];
1958      work[iIndex]=0.0;
1959      int iPivot=pivotVariable_[iRow];
1960      double change = theta*alpha;
1961      double value = solution_[iPivot] - change;
1962      solution_[iPivot]=value;
1963#ifndef NDEBUG
1964      // check if not active then okay
1965      if (!active(iRow)&&(specialOptions_&4)==0&&pivotRow_!=-1) {
1966        // But make sure one going out is feasible
1967        if (change>0.0) {
1968          // going down
1969          if (value<=lower_[iPivot]+primalTolerance_) {
1970            if (iPivot==sequenceOut_&&value>lower_[iPivot]-relaxedTolerance)
1971              value=lower_[iPivot];
1972            double difference = nonLinearCost_->setOne(iPivot,value);
1973            assert (!difference);
1974          }
1975        } else {
1976          // going up
1977          if (value>=upper_[iPivot]-primalTolerance_) {
1978            if (iPivot==sequenceOut_&&value<upper_[iPivot]+relaxedTolerance)
1979              value=upper_[iPivot];
1980            double difference = nonLinearCost_->setOne(iPivot,value);
1981            assert (!difference);
1982          }
1983        }
1984      }
1985#endif   
1986      if (active(iRow)||theta_<0.0) {
1987        clearActive(iRow);
1988        // But make sure one going out is feasible
1989        if (change>0.0) {
1990          // going down
1991          if (value<=lower_[iPivot]+primalTolerance_) {
1992            if (iPivot==sequenceOut_&&value>=lower_[iPivot]-relaxedTolerance)
1993              value=lower_[iPivot];
1994            double difference = nonLinearCost_->setOne(iPivot,value);
1995            if (difference) {
1996              if (iRow==pivotRow_)
1997                pivotPosition=newNumber;
1998              work[newNumber] = difference;
1999              //change reduced cost on this
2000              dj_[iPivot] = -difference;
2001              which[newNumber++]=iRow;
2002            }
2003          }
2004        } else {
2005          // going up
2006          if (value>=upper_[iPivot]-primalTolerance_) {
2007            if (iPivot==sequenceOut_&&value<upper_[iPivot]+relaxedTolerance)
2008              value=upper_[iPivot];
2009            double difference = nonLinearCost_->setOne(iPivot,value);
2010            if (difference) {
2011              if (iRow==pivotRow_)
2012                pivotPosition=newNumber;
2013              work[newNumber] = difference;
2014              //change reduced cost on this
2015              dj_[iPivot] = -difference;
2016              which[newNumber++]=iRow;
2017            }
2018          }
2019        }
2020      }
2021    }
2022  } else {
2023    // values pass so look at all
2024    for (iIndex=0;iIndex<number;iIndex++) {
2025     
2026      int iRow = which[iIndex];
2027      double alpha = work[iIndex];
2028      work[iIndex]=0.0;
2029      int iPivot=pivotVariable_[iRow];
2030      double change = theta*alpha;
2031      double value = solution_[iPivot] - change;
2032      solution_[iPivot]=value;
2033      clearActive(iRow);
2034      // But make sure one going out is feasible
2035      if (change>0.0) {
2036        // going down
2037        if (value<=lower_[iPivot]+primalTolerance_) {
2038          if (iPivot==sequenceOut_&&value>lower_[iPivot]-relaxedTolerance)
2039            value=lower_[iPivot];
2040          double difference = nonLinearCost_->setOne(iPivot,value);
2041          if (difference) {
2042            if (iRow==pivotRow_)
2043              pivotPosition=newNumber;
2044            work[newNumber] = difference;
2045            //change reduced cost on this
2046            dj_[iPivot] = -difference;
2047            which[newNumber++]=iRow;
2048          }
2049        }
2050      } else {
2051        // going up
2052        if (value>=upper_[iPivot]-primalTolerance_) {
2053          if (iPivot==sequenceOut_&&value<upper_[iPivot]+relaxedTolerance)
2054            value=upper_[iPivot];
2055          double difference = nonLinearCost_->setOne(iPivot,value);
2056          if (difference) {
2057            if (iRow==pivotRow_)
2058              pivotPosition=newNumber;
2059            work[newNumber] = difference;
2060            //change reduced cost on this
2061            dj_[iPivot] = -difference;
2062            which[newNumber++]=iRow;
2063          }
2064        }
2065      }
2066    }
2067  }
2068  objectiveChange += nonLinearCost_->changeInCost();
2069  rowArray->setPacked();
2070#if 0
2071  rowArray->setNumElements(newNumber);
2072  rowArray->expand();
2073  if (pivotRow_>=0) {
2074    dualIn_ += (oldCost-cost_[sequenceOut_]);
2075    // update change vector to include pivot
2076    rowArray->add(pivotRow_,-dualIn_);
2077    // and convert to packed
2078    rowArray->scanAndPack();
2079  } else {
2080    // and convert to packed
2081    rowArray->scanAndPack();
2082  }
2083#else
2084  if (pivotRow_>=0) {
2085    double dualIn = dualIn_+(oldCost-cost_[sequenceOut_]);
2086    // update change vector to include pivot
2087    if (pivotPosition>=0) {
2088      work[pivotPosition] -= dualIn;
2089    } else {
2090      work[newNumber]=-dualIn;
2091      which[newNumber++]=pivotRow_;
2092    }
2093  }
2094  rowArray->setNumElements(newNumber);
2095#endif
2096  return 0;
2097}
2098// Perturbs problem
2099void 
2100ClpSimplexPrimal::perturb(int type)
2101{
2102  if (perturbation_>100)
2103    return; //perturbed already
2104  if (perturbation_==100)
2105    perturbation_=50; // treat as normal
2106  int savePerturbation = perturbation_;
2107  int i;
2108  if (!numberIterations_)
2109    cleanStatus(); // make sure status okay
2110  // Make sure feasible bounds
2111  if (nonLinearCost_)
2112    nonLinearCost_->feasibleBounds();
2113  // look at element range
2114  double smallestNegative;
2115  double largestNegative;
2116  double smallestPositive;
2117  double largestPositive;
2118  matrix_->rangeOfElements(smallestNegative, largestNegative,
2119                           smallestPositive, largestPositive);
2120  smallestPositive = CoinMin(fabs(smallestNegative),smallestPositive);
2121  largestPositive = CoinMax(fabs(largestNegative),largestPositive);
2122  double elementRatio = largestPositive/smallestPositive;
2123  if (!numberIterations_&&perturbation_==50) {
2124    // See if we need to perturb
2125    double * sort = new double[numberRows_];
2126    for (i=0;i<numberRows_;i++) {
2127      double lo = fabs(lower_[i]);
2128      double up = fabs(upper_[i]);
2129      double value=0.0;
2130      if (lo&&lo<1.0e20) {
2131        if (up&&up<1.0e20)
2132          value = 0.5*(lo+up);
2133        else
2134          value=lo;
2135      } else {
2136        if (up&&up<1.0e20)
2137          value = up;
2138      }
2139      sort[i]=value;
2140    }
2141    std::sort(sort,sort+numberRows_);
2142    int number=1;
2143    double last = sort[0];
2144    for (i=1;i<numberRows_;i++) {
2145      if (last!=sort[i])
2146        number++;
2147      last=sort[i];
2148    }
2149    delete [] sort;
2150    //printf("ratio number diff rhs %g, element ratio %g\n",((double)number)/((double) numberRows_),
2151    //                                                                elementRatio);
2152    if (number*4>numberRows_||elementRatio>1.0e12) {
2153      perturbation_=100;
2154      return; // good enough
2155    }
2156  }
2157  // primal perturbation
2158  double perturbation=1.0e-20;
2159  int numberNonZero=0;
2160  // maximum fraction of rhs/bounds to perturb
2161  double maximumFraction = 1.0e-5;
2162  if (perturbation_>=50) {
2163    perturbation = 1.0e-4;
2164    for (i=0;i<numberColumns_+numberRows_;i++) {
2165      if (upper_[i]>lower_[i]+primalTolerance_) {
2166        double lowerValue, upperValue;
2167        if (lower_[i]>-1.0e20)
2168          lowerValue = fabs(lower_[i]);
2169        else
2170          lowerValue=0.0;
2171        if (upper_[i]<1.0e20)
2172          upperValue = fabs(upper_[i]);
2173        else
2174          upperValue=0.0;
2175        double value = CoinMax(fabs(lowerValue),fabs(upperValue));
2176        value = CoinMin(value,upper_[i]-lower_[i]);
2177#if 1
2178        if (value) {
2179          perturbation += value;
2180          numberNonZero++;
2181        }
2182#else
2183        perturbation = CoinMax(perturbation,value);
2184#endif
2185      }
2186    }
2187    if (numberNonZero) 
2188      perturbation /= (double) numberNonZero;
2189    else
2190      perturbation = 1.0e-1;
2191  } else if (perturbation_<100) {
2192    perturbation = pow(10.0,perturbation_);
2193    // user is in charge
2194    maximumFraction = 1.0;
2195  }
2196  double largestZero=0.0;
2197  double largest=0.0;
2198  double largestPerCent=0.0;
2199  bool printOut=(handler_->logLevel()==63);
2200  printOut=false; //off
2201  // Check if all slack
2202  int number=0;
2203  int iSequence;
2204  for (iSequence=0;iSequence<numberRows_;iSequence++) {
2205    if (getRowStatus(iSequence)==basic) 
2206      number++;
2207  }
2208  if (rhsScale_>100.0) {
2209    // tone down perturbation
2210    maximumFraction *= 0.1;
2211  }
2212  if (number!=numberRows_)
2213    type=1;
2214  // modify bounds
2215  // Change so at least 1.0e-5 and no more than 0.1
2216  // For now just no more than 0.1
2217  // printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
2218  if (type==1) {
2219    double tolerance = 100.0*primalTolerance_;
2220    //double multiplier = perturbation*maximumFraction;
2221    for (iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
2222      if (getStatus(iSequence)==basic) {
2223        double lowerValue = lower_[iSequence];
2224        double upperValue = upper_[iSequence];
2225        if (upperValue>lowerValue+tolerance) {
2226          double solutionValue = solution_[iSequence];
2227          double difference = upperValue-lowerValue;
2228          difference = CoinMin(difference,perturbation);
2229          difference = CoinMin(difference,fabs(solutionValue)+1.0);
2230          double value = maximumFraction*(difference+1.0);
2231          value = CoinMin(value,0.1);
2232          value *= CoinDrand48();
2233          if (solutionValue-lowerValue<=primalTolerance_) {
2234            lower_[iSequence] -= value;
2235          } else if (upperValue-solutionValue<=primalTolerance_) {
2236            upper_[iSequence] += value;
2237          } else {
2238#if 0
2239            if (iSequence>=numberColumns_) {
2240              // may not be at bound - but still perturb (unless free)
2241              if (upperValue>1.0e30&&lowerValue<-1.0e30)
2242                value=0.0;
2243              else
2244                value = - value; // as -1.0 in matrix
2245            } else {
2246              value = 0.0;
2247            }
2248#else
2249            value=0.0;
2250#endif
2251          }
2252          if (value) {
2253            if (printOut)
2254              printf("col %d lower from %g to %g, upper from %g to %g\n",
2255                     iSequence,lower_[iSequence],lowerValue,upper_[iSequence],upperValue);
2256            if (solutionValue) {
2257              largest = CoinMax(largest,value);
2258              if (value>(fabs(solutionValue)+1.0)*largestPerCent)
2259                largestPerCent=value/(fabs(solutionValue)+1.0);
2260            } else {
2261              largestZero = CoinMax(largestZero,value);
2262            } 
2263          }
2264        }
2265      }
2266    }
2267  } else {
2268    double tolerance = 100.0*primalTolerance_;
2269    for (i=0;i<numberColumns_;i++) {
2270      double lowerValue=lower_[i], upperValue=upper_[i];
2271      if (upperValue>lowerValue+primalTolerance_) {
2272        double value = perturbation*maximumFraction;
2273        value = CoinMin(value,0.1);
2274        value *= CoinDrand48();
2275        if (savePerturbation!=50) {
2276          if (fabs(value)<=primalTolerance_)
2277            value=0.0;
2278          if (lowerValue>-1.0e20&&lowerValue)
2279            lowerValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2280          if (upperValue<1.0e20&&upperValue)
2281            upperValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(upperValue))); 
2282        } else if (value) {
2283          double valueL =value *(CoinMax(1.0e-2,1.0e-5*fabs(lowerValue)));
2284          // get in range
2285          if (valueL<=tolerance) {
2286            valueL *= 10.0;
2287            while (valueL<=tolerance) 
2288              valueL *= 10.0;
2289          } else if (valueL>1.0) {
2290            valueL *= 0.1;
2291            while (valueL>1.0) 
2292              valueL *= 0.1;
2293          }
2294          if (lowerValue>-1.0e20&&lowerValue)
2295            lowerValue -= valueL;
2296          double valueU =value *(CoinMax(1.0e-2,1.0e-5*fabs(upperValue)));
2297          // get in range
2298          if (valueU<=tolerance) {
2299            valueU *= 10.0;
2300            while (valueU<=tolerance) 
2301              valueU *= 10.0;
2302          } else if (valueU>1.0) {
2303            valueU *= 0.1;
2304            while (valueU>1.0) 
2305              valueU *= 0.1;
2306          }
2307          if (upperValue<1.0e20&&upperValue)
2308            upperValue += valueU;
2309        }
2310        if (lowerValue!=lower_[i]) {
2311          double difference = fabs(lowerValue-lower_[i]);
2312          largest = CoinMax(largest,difference);
2313          if (difference>fabs(lower_[i])*largestPerCent)
2314            largestPerCent=fabs(difference/lower_[i]);
2315        } 
2316        if (upperValue!=upper_[i]) {
2317          double difference = fabs(upperValue-upper_[i]);
2318          largest = CoinMax(largest,difference);
2319          if (difference>fabs(upper_[i])*largestPerCent)
2320            largestPerCent=fabs(difference/upper_[i]);
2321        } 
2322        if (printOut)
2323          printf("col %d lower from %g to %g, upper from %g to %g\n",
2324                 i,lower_[i],lowerValue,upper_[i],upperValue);
2325      }
2326      lower_[i]=lowerValue;
2327      upper_[i]=upperValue;
2328    }
2329    for (;i<numberColumns_+numberRows_;i++) {
2330      double lowerValue=lower_[i], upperValue=upper_[i];
2331      double value = perturbation*maximumFraction;
2332      value = CoinMin(value,0.1);
2333      value *= CoinDrand48();
2334      if (upperValue>lowerValue+tolerance) {
2335        if (savePerturbation!=50) {
2336          if (fabs(value)<=primalTolerance_)
2337            value=0.0;
2338          if (lowerValue>-1.0e20&&lowerValue)
2339            lowerValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2340          if (upperValue<1.0e20&&upperValue)
2341            upperValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(upperValue))); 
2342        } else if (value) {
2343          double valueL =value *(CoinMax(1.0e-2,1.0e-5*fabs(lowerValue)));
2344          // get in range
2345          if (valueL<=tolerance) {
2346            valueL *= 10.0;
2347            while (valueL<=tolerance) 
2348              valueL *= 10.0;
2349          } else if (valueL>1.0) {
2350            valueL *= 0.1;
2351            while (valueL>1.0) 
2352              valueL *= 0.1;
2353          }
2354          if (lowerValue>-1.0e20&&lowerValue)
2355            lowerValue -= valueL;
2356          double valueU =value *(CoinMax(1.0e-2,1.0e-5*fabs(upperValue)));
2357          // get in range
2358          if (valueU<=tolerance) {
2359            valueU *= 10.0;
2360            while (valueU<=tolerance) 
2361              valueU *= 10.0;
2362          } else if (valueU>1.0) {
2363            valueU *= 0.1;
2364            while (valueU>1.0) 
2365              valueU *= 0.1;
2366          }
2367          if (upperValue<1.0e20&&upperValue)
2368            upperValue += valueU;
2369        }
2370      } else if (upperValue>0.0) {
2371        upperValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2372        lowerValue -= value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2373      } else if (upperValue<0.0) {
2374        upperValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2375        lowerValue += value * (CoinMax(1.0e-2,1.0e-5*fabs(lowerValue))); 
2376      } else {
2377      }
2378      if (lowerValue!=lower_[i]) {
2379        double difference = fabs(lowerValue-lower_[i]);
2380        largest = CoinMax(largest,difference);
2381        if (difference>fabs(lower_[i])*largestPerCent)
2382          largestPerCent=fabs(difference/lower_[i]);
2383      } 
2384      if (upperValue!=upper_[i]) {
2385        double difference = fabs(upperValue-upper_[i]);
2386        largest = CoinMax(largest,difference);
2387        if (difference>fabs(upper_[i])*largestPerCent)
2388          largestPerCent=fabs(difference/upper_[i]);
2389      } 
2390      if (printOut)
2391        printf("row %d lower from %g to %g, upper from %g to %g\n",
2392               i-numberColumns_,lower_[i],lowerValue,upper_[i],upperValue);
2393      lower_[i]=lowerValue;
2394      upper_[i]=upperValue;
2395    }
2396  }
2397  // Clean up
2398  for (i=0;i<numberColumns_+numberRows_;i++) {
2399    switch(getStatus(i)) {
2400     
2401    case basic:
2402      break;
2403    case atUpperBound:
2404      solution_[i]=upper_[i];
2405      break;
2406    case isFixed:
2407    case atLowerBound:
2408      solution_[i]=lower_[i];
2409      break;
2410    case isFree:
2411      break;
2412    case superBasic:
2413      break;
2414    }
2415  }
2416  handler_->message(CLP_SIMPLEX_PERTURB,messages_)
2417    <<100.0*maximumFraction<<perturbation<<largest<<100.0*largestPerCent<<largestZero
2418    <<CoinMessageEol;
2419  // redo nonlinear costs
2420  // say perturbed
2421  perturbation_=101;
2422}
2423// un perturb
2424bool
2425ClpSimplexPrimal::unPerturb()
2426{
2427  if (perturbation_!=101)
2428    return false;
2429  // put back original bounds and costs
2430  createRim(1+4);
2431  sanityCheck();
2432  // unflag
2433  unflag();
2434  // get a valid nonlinear cost function
2435  delete nonLinearCost_;
2436  nonLinearCost_= new ClpNonLinearCost(this);
2437  perturbation_ = 102; // stop any further perturbation
2438  // move non basic variables to new bounds
2439  nonLinearCost_->checkInfeasibilities(0.0);
2440#if 1
2441  // Try using dual
2442  return true;
2443#else
2444  gutsOfSolution(NULL,NULL,ifValuesPass!=0);
2445  return false;
2446#endif
2447 
2448}
2449// Unflag all variables and return number unflagged
2450int 
2451ClpSimplexPrimal::unflag()
2452{
2453  int i;
2454  int number = numberRows_+numberColumns_;
2455  int numberFlagged=0;
2456  // we can't really trust infeasibilities if there is dual error
2457  // allow tolerance bigger than standard to check on duals
2458  double relaxedToleranceD=dualTolerance_ + CoinMin(1.0e-2,10.0*largestDualError_);
2459  for (i=0;i<number;i++) {
2460    if (flagged(i)) {
2461      clearFlagged(i);
2462      // only say if reasonable dj
2463      if (fabs(dj_[i])>relaxedToleranceD)
2464        numberFlagged++;
2465    }
2466  }
2467  numberFlagged += matrix_->generalExpanded(this,8,i);
2468  if (handler_->logLevel()>2&&numberFlagged&&objective_->type()>1)
2469    printf("%d unflagged\n",numberFlagged);
2470  return numberFlagged;
2471}
2472// Do not change infeasibility cost and always say optimal
2473void 
2474ClpSimplexPrimal::alwaysOptimal(bool onOff)
2475{
2476  if (onOff)
2477    specialOptions_ |= 1;
2478  else
2479    specialOptions_ &= ~1;
2480}
2481bool 
2482ClpSimplexPrimal::alwaysOptimal() const
2483{
2484  return (specialOptions_&1)!=0;
2485}
2486// Flatten outgoing variables i.e. - always to exact bound
2487void 
2488ClpSimplexPrimal::exactOutgoing(bool onOff)
2489{
2490  if (onOff)
2491    specialOptions_ |= 4;
2492  else
2493    specialOptions_ &= ~4;
2494}
2495bool 
2496ClpSimplexPrimal::exactOutgoing() const
2497{
2498  return (specialOptions_&4)!=0;
2499}
2500/*
2501  Reasons to come out (normal mode/user mode):
2502  -1 normal
2503  -2 factorize now - good iteration/ NA
2504  -3 slight inaccuracy - refactorize - iteration done/ same but factor done
2505  -4 inaccuracy - refactorize - no iteration/ NA
2506  -5 something flagged - go round again/ pivot not possible
2507  +2 looks unbounded
2508  +3 max iterations (iteration done)
2509*/
2510int
2511ClpSimplexPrimal::pivotResult(int ifValuesPass)
2512{
2513
2514  bool roundAgain=true;
2515  int returnCode=-1;
2516
2517  // loop round if user setting and doing refactorization
2518  while (roundAgain) {
2519    roundAgain=false;
2520    returnCode=-1;
2521    pivotRow_=-1;
2522    sequenceOut_=-1;
2523    rowArray_[1]->clear();
2524#if 0
2525    {
2526      int seq[]={612,643};
2527      int k;
2528      for (k=0;k<sizeof(seq)/sizeof(int);k++) {
2529        int iSeq=seq[k];
2530        if (getColumnStatus(iSeq)!=basic) {
2531          double djval;
2532          double * work;
2533          int number;
2534          int * which;
2535         
2536          int iIndex;
2537          unpack(rowArray_[1],iSeq);
2538          factorization_->updateColumn(rowArray_[2],rowArray_[1]);
2539          djval = cost_[iSeq];
2540          work=rowArray_[1]->denseVector();
2541          number=rowArray_[1]->getNumElements();
2542          which=rowArray_[1]->getIndices();
2543         
2544          for (iIndex=0;iIndex<number;iIndex++) {
2545           
2546            int iRow = which[iIndex];
2547            double alpha = work[iRow];
2548            int iPivot=pivotVariable_[iRow];
2549            djval -= alpha*cost_[iPivot];
2550          }
2551          double comp = 1.0e-8 + 1.0e-7*(CoinMax(fabs(dj_[iSeq]),fabs(djval)));
2552          if (fabs(djval-dj_[iSeq])>comp)
2553            printf("Bad dj %g for %d - true is %g\n",
2554                   dj_[iSeq],iSeq,djval);
2555          assert (fabs(djval)<1.0e-3||djval*dj_[iSeq]>0.0);
2556          rowArray_[1]->clear();
2557        }
2558      }
2559    }
2560#endif
2561       
2562    // we found a pivot column
2563    // update the incoming column
2564    unpackPacked(rowArray_[1]);
2565    // save reduced cost
2566    double saveDj = dualIn_;
2567    factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
2568    // Get extra rows
2569    matrix_->extendUpdated(this,rowArray_[1],0);
2570    // do ratio test and re-compute dj
2571    primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
2572              ifValuesPass);
2573    if (ifValuesPass) {
2574      saveDj=dualIn_;
2575      //assert (fabs(alpha_)>=1.0e-5||(objective_->type()<2||!objective_->activated())||pivotRow_==-2);
2576      if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
2577        if(fabs(dualIn_)<1.0e2*dualTolerance_&&objective_->type()<2) {
2578          // try other way
2579          directionIn_=-directionIn_;
2580          primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
2581                    0);
2582        }
2583        if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
2584          if (solveType_==1) {
2585            // reject it
2586            char x = isColumn(sequenceIn_) ? 'C' :'R';
2587            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2588              <<x<<sequenceWithin(sequenceIn_)
2589              <<CoinMessageEol;
2590            setFlagged(sequenceIn_);
2591            progress_->clearBadTimes();
2592            lastBadIteration_ = numberIterations_; // say be more cautious
2593            clearAll();
2594            pivotRow_=-1;
2595          }
2596          returnCode=-5;
2597          break;
2598        }
2599      }
2600    }
2601    // need to clear toIndex_ in gub
2602    // ? when can I clear stuff
2603    // Clean up any gub stuff
2604    matrix_->extendUpdated(this,rowArray_[1],1);
2605    double checkValue=1.0e-2;
2606    if (largestDualError_>1.0e-5)
2607      checkValue=1.0e-1;
2608    if (!ifValuesPass&&solveType_==1&&(saveDj*dualIn_<1.0e-20||
2609        fabs(saveDj-dualIn_)>checkValue*(1.0+fabs(saveDj))||
2610                        fabs(dualIn_)<dualTolerance_)) {
2611      char x = isColumn(sequenceIn_) ? 'C' :'R';
2612      handler_->message(CLP_PRIMAL_DJ,messages_)
2613        <<x<<sequenceWithin(sequenceIn_)
2614        <<saveDj<<dualIn_
2615        <<CoinMessageEol;
2616      if(lastGoodIteration_ != numberIterations_) {
2617        clearAll();
2618        pivotRow_=-1; // say no weights update
2619        returnCode=-4;
2620        if(lastGoodIteration_+1 == numberIterations_) {
2621          // not looking wonderful - try cleaning bounds
2622          // put non-basics to bounds in case tolerance moved
2623          nonLinearCost_->checkInfeasibilities(0.0);
2624        }
2625        sequenceOut_=-1;
2626        break;
2627      } else {
2628        // take on more relaxed criterion
2629        if (saveDj*dualIn_<1.0e-20||
2630            fabs(saveDj-dualIn_)>2.0e-1*(1.0+fabs(dualIn_))||
2631            fabs(dualIn_)<dualTolerance_) {
2632          // need to reject something
2633          char x = isColumn(sequenceIn_) ? 'C' :'R';
2634          handler_->message(CLP_SIMPLEX_FLAG,messages_)
2635            <<x<<sequenceWithin(sequenceIn_)
2636            <<CoinMessageEol;
2637          setFlagged(sequenceIn_);
2638          progress_->clearBadTimes();
2639          lastBadIteration_ = numberIterations_; // say be more cautious
2640          clearAll();
2641          pivotRow_=-1;
2642          returnCode=-5;
2643          sequenceOut_=-1;
2644          break;
2645        }
2646      }
2647    }
2648    if (pivotRow_>=0) {
2649      if (solveType_==2) {
2650        // **** Coding for user interface
2651        // do ray
2652        primalRay(rowArray_[1]);
2653        // update duals
2654        // as packed need to find pivot row
2655        //assert (rowArray_[1]->packedMode());
2656        //int i;
2657       
2658        //alpha_ = rowArray_[1]->denseVector()[pivotRow_];
2659        CoinAssert (fabs(alpha_)>1.0e-8);
2660        double multiplier = dualIn_/alpha_;
2661        rowArray_[0]->insert(pivotRow_,multiplier);
2662        factorization_->updateColumnTranspose(rowArray_[2],rowArray_[0]);
2663        // put row of tableau in rowArray[0] and columnArray[0]
2664        matrix_->transposeTimes(this,-1.0,
2665                                rowArray_[0],columnArray_[1],columnArray_[0]);
2666        // update column djs
2667        int i;
2668        int * index = columnArray_[0]->getIndices();
2669        int number = columnArray_[0]->getNumElements();
2670        double * element = columnArray_[0]->denseVector();
2671        for (i=0;i<number;i++) {
2672          int ii = index[i];
2673          dj_[ii] += element[ii];
2674          reducedCost_[ii] = dj_[ii];
2675          element[ii]=0.0;
2676        }
2677        columnArray_[0]->setNumElements(0);
2678        // and row djs
2679        index = rowArray_[0]->getIndices();
2680        number = rowArray_[0]->getNumElements();
2681        element = rowArray_[0]->denseVector();
2682        for (i=0;i<number;i++) {
2683          int ii = index[i];
2684          dj_[ii+numberColumns_] += element[ii];
2685          dual_[ii] = dj_[ii+numberColumns_];
2686          element[ii]=0.0;
2687        }
2688        rowArray_[0]->setNumElements(0);
2689        // check incoming
2690        CoinAssert (fabs(dj_[sequenceIn_])<1.0e-1);
2691      }
2692      // if stable replace in basis
2693      // If gub or odd then alpha and pivotRow may change
2694      int updateType=0;
2695      int updateStatus = matrix_->generalExpanded(this,3,updateType);
2696      if (updateType>=0)
2697        updateStatus = factorization_->replaceColumn(this,
2698                                                     rowArray_[2],
2699                                                     rowArray_[1],
2700                                                     pivotRow_,
2701                                                     alpha_);
2702
2703      // if no pivots, bad update but reasonable alpha - take and invert
2704      if (updateStatus==2&&
2705          lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
2706        updateStatus=4;
2707      if (updateStatus==1||updateStatus==4) {
2708        // slight error
2709        if (factorization_->pivots()>5||updateStatus==4) {
2710          returnCode=-3;
2711        }
2712      } else if (updateStatus==2) {
2713        // major error
2714        // better to have small tolerance even if slower
2715        factorization_->zeroTolerance(1.0e-15);
2716        int maxFactor = factorization_->maximumPivots();
2717        if (maxFactor>10) {
2718          if (forceFactorization_<0)
2719            forceFactorization_= maxFactor;
2720          forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
2721        } 
2722        // later we may need to unwind more e.g. fake bounds
2723        if(lastGoodIteration_ != numberIterations_) {
2724          clearAll();
2725          pivotRow_=-1;
2726          if (solveType_==1) {
2727            returnCode=-4;
2728            break;
2729          } else {
2730            // user in charge - re-factorize
2731            int lastCleaned;
2732            ClpSimplexProgress dummyProgress;
2733            if (saveStatus_)
2734              statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true,ifValuesPass);
2735            else
2736              statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true,ifValuesPass);
2737            roundAgain=true;
2738            continue;
2739          }
2740        } else {
2741          // need to reject something
2742          if (solveType_==1) {
2743            char x = isColumn(sequenceIn_) ? 'C' :'R';
2744            handler_->message(CLP_SIMPLEX_FLAG,messages_)
2745              <<x<<sequenceWithin(sequenceIn_)
2746              <<CoinMessageEol;
2747            setFlagged(sequenceIn_);
2748            progress_->clearBadTimes();
2749          }
2750          lastBadIteration_ = numberIterations_; // say be more cautious
2751          clearAll();
2752          pivotRow_=-1;
2753          sequenceOut_=-1;
2754          returnCode = -5;
2755          break;
2756
2757        }
2758      } else if (updateStatus==3) {
2759        // out of memory
2760        // increase space if not many iterations
2761        if (factorization_->pivots()<
2762            0.5*factorization_->maximumPivots()&&
2763            factorization_->pivots()<200)
2764          factorization_->areaFactor(
2765                                     factorization_->areaFactor() * 1.1);
2766        returnCode =-2; // factorize now
2767      } else if (updateStatus==5) {
2768        problemStatus_=-2; // factorize now
2769      }
2770      // here do part of steepest - ready for next iteration
2771      if (!ifValuesPass)
2772        primalColumnPivot_->updateWeights(rowArray_[1]);
2773    } else {
2774      if (pivotRow_==-1) {
2775        // no outgoing row is valid
2776        if (valueOut_!=COIN_DBL_MAX) {
2777          double objectiveChange=0.0;
2778          theta_=valueOut_-valueIn_;
2779          updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,ifValuesPass);
2780          solution_[sequenceIn_] += theta_;
2781        }
2782        rowArray_[0]->clear();
2783        if (!factorization_->pivots()&&acceptablePivot_<=1.0e-8) {
2784          returnCode = 2; //say looks unbounded
2785          // do ray
2786          primalRay(rowArray_[1]);
2787        } else if (solveType_==2) {
2788          // refactorize
2789          int lastCleaned;
2790          ClpSimplexProgress dummyProgress;
2791          if (saveStatus_)
2792            statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true,ifValuesPass);
2793          else
2794            statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true,ifValuesPass);
2795          roundAgain=true;
2796          continue;
2797        } else {
2798          acceptablePivot_=1.0e-8;
2799          returnCode = 4; //say looks unbounded but has iterated
2800        }
2801        break;
2802      } else {
2803        // flipping from bound to bound
2804      }
2805    }
2806
2807    double oldCost = 0.0;
2808    if (sequenceOut_>=0)
2809      oldCost=cost_[sequenceOut_];
2810    // update primal solution
2811   
2812    double objectiveChange=0.0;
2813    // after this rowArray_[1] is not empty - used to update djs
2814    // If pivot row >= numberRows then may be gub
2815    int savePivot = pivotRow_;
2816    if (pivotRow_>=numberRows_)
2817      pivotRow_=-1;
2818    updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange,ifValuesPass);
2819    pivotRow_=savePivot;
2820   
2821    double oldValue = valueIn_;
2822    if (directionIn_==-1) {
2823      // as if from upper bound
2824      if (sequenceIn_!=sequenceOut_) {
2825        // variable becoming basic
2826        valueIn_ -= fabs(theta_);
2827      } else {
2828        valueIn_=lowerIn_;
2829      }
2830    } else {
2831      // as if from lower bound
2832      if (sequenceIn_!=sequenceOut_) {
2833        // variable becoming basic
2834        valueIn_ += fabs(theta_);
2835      } else {
2836        valueIn_=upperIn_;
2837      }
2838    }
2839    objectiveChange += dualIn_*(valueIn_-oldValue);
2840    // outgoing
2841    if (sequenceIn_!=sequenceOut_) {
2842      if (directionOut_>0) {
2843        valueOut_ = lowerOut_;
2844      } else {
2845        valueOut_ = upperOut_;
2846      }
2847      if(valueOut_<lower_[sequenceOut_]-primalTolerance_)
2848        valueOut_=lower_[sequenceOut_]-0.9*primalTolerance_;
2849      else if (valueOut_>upper_[sequenceOut_]+primalTolerance_)
2850        valueOut_=upper_[sequenceOut_]+0.9*primalTolerance_;
2851      // may not be exactly at bound and bounds may have changed
2852      // Make sure outgoing looks feasible
2853      directionOut_=nonLinearCost_->setOneOutgoing(sequenceOut_,valueOut_);
2854      // May have got inaccurate
2855      //if (oldCost!=cost_[sequenceOut_])
2856      //printf("costchange on %d from %g to %g\n",sequenceOut_,
2857      //       oldCost,cost_[sequenceOut_]);
2858      if (solveType_!=2)
2859        dj_[sequenceOut_]=cost_[sequenceOut_]-oldCost; // normally updated next iteration
2860      solution_[sequenceOut_]=valueOut_;
2861    }
2862    // change cost and bounds on incoming if primal
2863    nonLinearCost_->setOne(sequenceIn_,valueIn_); 
2864    int whatNext=housekeeping(objectiveChange);
2865    //nonLinearCost_->validate();
2866#if CLP_DEBUG >1
2867    {
2868      double sum;
2869      int ninf= matrix_->checkFeasible(this,sum);
2870      if (ninf)
2871        printf("infeas %d\n",ninf);
2872    }
2873#endif
2874    if (whatNext==1) {
2875        returnCode =-2; // refactorize
2876    } else if (whatNext==2) {
2877      // maximum iterations or equivalent
2878      returnCode=3;
2879    } else if(numberIterations_ == lastGoodIteration_
2880              + 2 * factorization_->maximumPivots()) {
2881      // done a lot of flips - be safe
2882      returnCode =-2; // refactorize
2883    }
2884    // Check event
2885    {
2886      int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2887      if (status>=0) {
2888        problemStatus_=5;
2889        secondaryStatus_=ClpEventHandler::endOfIteration;
2890        returnCode=4;
2891      }
2892    }
2893  }
2894  if (solveType_==2&&(returnCode == -2||returnCode==-3)) {
2895    // refactorize here
2896    int lastCleaned;
2897    ClpSimplexProgress dummyProgress;
2898    if (saveStatus_)
2899      statusOfProblemInPrimal(lastCleaned,1,&dummyProgress,true,ifValuesPass);
2900    else
2901      statusOfProblemInPrimal(lastCleaned,0,&dummyProgress,true,ifValuesPass);
2902    if (problemStatus_==5) {
2903      printf("Singular basis\n");
2904      problemStatus_=-1;
2905      returnCode=5;
2906    }
2907  }
2908#ifdef CLP_DEBUG
2909  {
2910    int i;
2911    // not [1] as may have information
2912    for (i=0;i<4;i++) {
2913      if (i!=1)
2914        rowArray_[i]->checkClear();
2915    }   
2916    for (i=0;i<2;i++) {
2917      columnArray_[i]->checkClear();
2918    }   
2919  }     
2920#endif
2921  return returnCode;
2922}
2923// Create primal ray
2924void 
2925ClpSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
2926{
2927  delete [] ray_;
2928  ray_ = new double [numberColumns_];
2929  CoinZeroN(ray_,numberColumns_);
2930  int number=rowArray->getNumElements();
2931  int * index = rowArray->getIndices();
2932  double * array = rowArray->denseVector();
2933  double way=-directionIn_;
2934  int i;
2935  double zeroTolerance=1.0e-12;
2936  if (sequenceIn_<numberColumns_)
2937    ray_[sequenceIn_]=directionIn_;
2938  if (!rowArray->packedMode()) {
2939    for (i=0;i<number;i++) {
2940      int iRow=index[i];
2941      int iPivot=pivotVariable_[iRow];
2942      double arrayValue = array[iRow];
2943      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2944        ray_[iPivot] = way* arrayValue;
2945    }
2946  } else {
2947    for (i=0;i<number;i++) {
2948      int iRow=index[i];
2949      int iPivot=pivotVariable_[iRow];
2950      double arrayValue = array[i];
2951      if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
2952        ray_[iPivot] = way* arrayValue;
2953    }
2954  }
2955}
2956/* Get next superbasic -1 if none,
2957   Normal type is 1
2958   If type is 3 then initializes sorted list
2959   if 2 uses list.
2960*/
2961int 
2962ClpSimplexPrimal::nextSuperBasic(int superBasicType,CoinIndexedVector * columnArray)
2963{
2964  if (firstFree_>=0&&superBasicType) {
2965    int returnValue=-1;
2966    bool finished=false;
2967    while (!finished) {
2968      returnValue=firstFree_;
2969      int iColumn=firstFree_+1;
2970      if (superBasicType>1) {
2971        if (superBasicType>2) {
2972          // Initialize list
2973          // Wild guess that lower bound more natural than upper
2974          int number=0;
2975          double * work=columnArray->denseVector();
2976          int * which=columnArray->getIndices();
2977          for (iColumn=0;iColumn<numberRows_+numberColumns_;iColumn++) {
2978            if (!flagged(iColumn)) {
2979              if (getStatus(iColumn)==superBasic) {
2980                if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
2981                  solution_[iColumn]=lower_[iColumn];
2982                  setStatus(iColumn,atLowerBound);
2983                } else if (fabs(solution_[iColumn]-upper_[iColumn])
2984                           <=primalTolerance_) {
2985                  solution_[iColumn]=upper_[iColumn];
2986                  setStatus(iColumn,atUpperBound);
2987                } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
2988                  setStatus(iColumn,isFree);
2989                  break;
2990                } else if (!flagged(iColumn)) {
2991                  // put ones near bounds at end after sorting
2992                  work[number]= - CoinMin(0.1*(solution_[iColumn]-lower_[iColumn]),
2993                                      upper_[iColumn]-solution_[iColumn]);
2994                  which[number++] = iColumn;
2995                }
2996              }
2997            }
2998          }
2999          CoinSort_2(work,work+number,which);
3000          columnArray->setNumElements(number);
3001          CoinZeroN(work,number);
3002        }
3003        int * which=columnArray->getIndices();
3004        int number = columnArray->getNumElements();
3005        if (!number) {
3006          // finished
3007          iColumn = numberRows_+numberColumns_;
3008          returnValue=-1;
3009        } else {
3010          number--;
3011          returnValue=which[number];
3012          iColumn=returnValue;
3013          columnArray->setNumElements(number);
3014        }     
3015      } else {
3016        for (;iColumn<numberRows_+numberColumns_;iColumn++) {
3017          if (!flagged(iColumn)) {
3018            if (getStatus(iColumn)==superBasic) {
3019              if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
3020                solution_[iColumn]=lower_[iColumn];
3021                setStatus(iColumn,atLowerBound);
3022              } else if (fabs(solution_[iColumn]-upper_[iColumn])
3023                         <=primalTolerance_) {
3024                solution_[iColumn]=upper_[iColumn];
3025                setStatus(iColumn,atUpperBound);
3026              } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
3027                setStatus(iColumn,isFree);
3028                break;
3029              } else {
3030                break;
3031              }
3032            }
3033          }
3034        }
3035      }
3036      firstFree_ = iColumn;
3037      finished=true;
3038      if (firstFree_==numberRows_+numberColumns_)
3039        firstFree_=-1;
3040      if (returnValue>=0&&getStatus(returnValue)!=superBasic&&getStatus(returnValue)!=isFree)
3041        finished=false; // somehow picked up odd one
3042    }
3043    return returnValue;
3044  } else {
3045    return -1;
3046  }
3047}
3048void
3049ClpSimplexPrimal::clearAll()
3050{
3051  // Clean up any gub stuff
3052  matrix_->extendUpdated(this,rowArray_[1],1);
3053  int number=rowArray_[1]->getNumElements();
3054  int * which=rowArray_[1]->getIndices();
3055 
3056  int iIndex;
3057  for (iIndex=0;iIndex<number;iIndex++) {
3058   
3059    int iRow = which[iIndex];
3060    clearActive(iRow);
3061  }
3062  rowArray_[1]->clear();
3063  // make sure any gub sets are clean
3064  matrix_->generalExpanded(this,11,sequenceIn_);
3065 
3066}
Note: See TracBrowser for help on using the repository browser.