# source:branches/devel-1/ClpSimplexPrimal.cpp@3

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

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Author Date Id Revision`
File size: 41.6 KB
Line
3
4
5/* Notes on implementation of primal simplex algorithm.
6
7   When primal feasible(A):
8
9   If dual feasible, we are optimal.  Otherwise choose an infeasible
10   basic variable to enter basis from a bound (B).  We now need to find an
11   outgoing variable which will leave problem primal feasible so we get
12   the column of the tableau corresponding to the incoming variable
13   (with the correct sign depending if variable will go up or down).
14
15   We now perform a ratio test to determine which outgoing variable will
16   preserve primal feasibility (C).  If no variable found then problem
17   is unbounded (in primal sense).  If there is a variable, we then
18   perform pivot and repeat.  Trivial?
19
20   -------------------------------------------
21
22   A) How do we get primal feasible?  All variables have fake costs
23   outside their feasible region so it is trivial to declare problem
24   feasible.  OSL did not have a phase 1/phase 2 approach but
25   instead effectively put an extra cost on infeasible basic variables
26   I am taking the same approach here, although it is generalized
27   to allow for non-linear costs and dual information.
28
29   In OSL, this weight was changed heuristically, here at present
30   it is only increased if problem looks finished.  If problem is
31   feasible I check for unboundedness.  If not unbounded we
32   could play with going into dual.  As long as weights increase
33   any algorithm would be finite.
34
35   B) Which incoming variable to choose is a virtual base class.
36   For difficult problems steepest edge is preferred while for
37   very easy (large) problems we will need partial scan.
38
39   C) Sounds easy, but this is hardest part of algorithm.
40      1) Instead of stopping at first choice, we may be able
41      to allow that variable to go through bound and if objective
42      still improving choose again.  These mini iterations can
43      increase speed by orders of magnitude but we may need to
44      go to more of a bucket choice of variable rather than looking
45      at them one by one (for speed).
46      2) Accuracy.  Basic infeasibilities may be less than
47      tolerance.  Pivoting on these makes objective go backwards.
48      OSL modified cost so a zero move was made, Gill et al
49      modified so a strictly positive move was made.
50      The two problems are that re-factorizations can
51      change rinfeasibilities above and below tolerances and that when
52      finished we need to reset costs and try again.
53      3) Degeneracy.  Gill et al helps but may not be enough.  We
54      may need more.  Also it can improve speed a lot if we perturb
55      the costs significantly.
56
57  References:
58     Forrest and Goldfarb, Steepest-edge simplex algorithms for
59       linear programming - Mathematical Programming 1992
60     Forrest and Tomlin, Implementing the simplex method for
61       the Optimization Subroutine Library - IBM Systems Journal 1992
62     Gill, Murray, Saunders, Wright A Practical Anti-Cycling
63       Procedure for Linear and Nonlinear Programming SOL report 1988
64
65
66  TODO:
67
68  a) Better recovery procedures.  At present I never check on forward
70     re-factorization frequency, but this is only on singular
71     factorizations.
72  b) Fast methods for large easy problems (and also the option for
73     the code to automatically choose which method).
74  c) We need to be able to stop in various ways for OSI - this
75     is fairly easy.
76
77 */
78
79#if defined(_MSC_VER)
80// Turn off compiler warning about long names
81#  pragma warning(disable:4786)
82#endif
83
84#include <math.h>
85
86#include "CoinHelperFunctions.hpp"
87#include "ClpSimplexPrimal.hpp"
88#include "ClpFactorization.hpp"
89#include "ClpNonLinearCost.hpp"
90#include "OsiPackedMatrix.hpp"
91#include "OsiIndexedVector.hpp"
92#include "OsiWarmStartBasis.hpp"
93#include "ClpPrimalColumnPivot.hpp"
94#include "ClpMessage.hpp"
95#include <cfloat>
96#include <cassert>
97#include <string>
98#include <stdio.h>
99#include <iostream>
100// This returns a non const array filled with input from scalar
101// or actual array
102template <class T> inline T*
103copyOfArray( const T * array, const int size, T value)
104{
105  T * arrayNew = new T[size];
106  if (array)
107    CoinDisjointCopyN(array,size,arrayNew);
108  else
109    CoinFillN ( arrayNew, size,value);
110  return arrayNew;
111}
112
113// This returns a non const array filled with actual array (or NULL)
114template <class T> inline T*
115copyOfArray( const T * array, const int size)
116{
117  if (array) {
118    T * arrayNew = new T[size];
119    CoinDisjointCopyN(array,size,arrayNew);
120    return arrayNew;
121  } else {
122    return NULL;
123  }
124}
125// primal
126int ClpSimplexPrimal::primal (int ifValuesPass )
127{
128
129  /*
130      Method
131
132     It tries to be a single phase approach with a weight of 1.0 being
133     given to getting optimal and a weight of infeasibilityCost_ being
134     given to getting primal feasible.  In this version I have tried to
135     be clever in a stupid way.  The idea of fake bounds in dual
136     seems to work so the primal analogue would be that of getting
137     bounds on reduced costs (by a presolve approach) and using
138     these for being above or below feasible region.  I decided to waste
139     memory and keep these explicitly.  This allows for non-linear
140     costs!
141
142     The code is designed to take advantage of sparsity so arrays are
143     seldom zeroed out from scratch or gone over in their entirety.
144     The only exception is a full scan to find incoming variable for
145     Dantzig row choice.  For steepest edge we keep an updated list
146     of dual infeasibilities (actually squares).
147     On easy problems we don't need full scan - just
148     pick first reasonable.
149
150     One problem is how to tackle degeneracy and accuracy.  At present
151     I am using the modification of costs which I put in OSL and which was
152     extended by Gill et al.  I am still not sure of the exact details.
153
154     The flow of primal is three while loops as follows:
155
156     while (not finished) {
157
158       while (not clean solution) {
159
160          Factorize and/or clean up solution by changing bounds so
161          primal feasible.  If looks finished check fake primal bounds.
162          Repeat until status is iterating (-1) or finished (0,1,2)
163
164       }
165
166       while (status==-1) {
167
168         Iterate until no pivot in or out or time to re-factorize.
169
170         Flow is:
171
172         choose pivot column (incoming variable).  if none then
173         we are primal feasible so looks as if done but we need to
174         break and check bounds etc.
175
176         Get pivot column in tableau
177
178         Choose outgoing row.  If we don't find one then we look
179         primal unbounded so break and check bounds etc.  (Also the
180         pivot tolerance is larger after any iterations so that may be
181         reason)
182
183         If we do find outgoing row, we may have to adjust costs to
184         keep going forwards (anti-degeneracy).  Check pivot will be stable
185         and if unstable throw away iteration and break to re-factorize.
186         If minor error re-factorize after iteration.
187
188         Update everything (this may involve changing bounds on
189         variables to stay primal feasible.
190
191       }
192
193     }
194
195     At present we never check we are going forwards.  I overdid that in
196     OSL so will try and make a last resort.
197
198     Needs partial scan pivot in option.
199
200     May need other anti-degeneracy measures, especially if we try and use
201     loose tolerances as a way to solve in fewer iterations.
202
203     I like idea of dynamic scaling.  This gives opportunity to decouple
204     different implications of scaling for accuracy, iteration count and
205     feasibility tolerance.
206
207  */
208
209  // sanity check
210  assert (numberRows_==matrix_->getNumRows());
211  assert (numberColumns_==matrix_->getNumCols());
212  // for moment all arrays must exist
213  assert(columnLower_);
214  assert(columnUpper_);
215  assert(rowLower_);
216  assert(rowUpper_);
217
218#ifdef CLP_DEBUG
219  int debugIteration=-1;
220#endif
221
222  algorithm_ = +1;
223  primalTolerance_=dblParam_[OsiPrimalTolerance];
224  dualTolerance_=dblParam_[OsiDualTolerance];
225
226  // put in standard form (and make row copy)
227  // create modifiable copies of model rim and do optional scaling
228  createRim(7+8+16,true);
229
230  // save infeasibility cost
231  double saveInfeasibilityCost = infeasibilityCost_;
232
233  int iRow,iColumn;
234  // Do initial factorization
235  // and set certain stuff
236  // We can either set increasing rows so ...IsBasic gives pivot row
237  // or we can just increment iBasic one by one
238  // for now let ...iBasic give pivot row
239  factorization_->increasingRows(2);
240  // row activities have negative sign
241  factorization_->slackValue(-1.0);
242  factorization_->zeroTolerance(1.0e-13);
243
244
245  // If user asked for perturbation - do it
246  int savePerturbation = perturbation_;
247
248  if (perturbation_<100)
249    perturb();
250
251  double objectiveChange;
252  // for primal we will change bounds using infeasibilityCost_
253  if (nonLinearCost_==NULL) {
254    // get a valid nonlinear cost function
255    delete nonLinearCost_;
256    nonLinearCost_= new ClpNonLinearCost(this);
257  }
258
259  // loop round to clean up solution if values pass
260  int numberThrownOut = -1;
261  int firstSuperBasic=numberRows_+numberColumns_;
262  while(numberThrownOut) {
263    if (internalFactorize(0+10*ifValuesPass))
264      return 1; // some error
265
266    // for this we need clean basis so it is after factorize
267    numberThrownOut=gutsOfSolution(rowActivityWork_,columnActivityWork_,
268                                   ifValuesPass);
269
270    // find first superbasic - columns, then rows
271    if (ifValuesPass) {
272      nextSuperBasic(firstSuperBasic);
273      if (firstSuperBasic==numberRows_+numberColumns_)
274        ifValuesPass=0; // signal no values pass
275    }
276  }
277
278  problemStatus_ = -1;
279  numberIterations_=0;
280
281  int lastCleaned=0; // last time objective or bounds cleaned up
282
283  // number of times we have declared optimality
284  numberTimesOptimal_=0;
285
286  // Say no pivot has occurred (for steepest edge and updates)
287  pivotRow_=-2;
288
289  // This says whether to restore things etc
290  int factorType=0;
291  /*
292    Status of problem:
293    0 - optimal
294    1 - infeasible
295    2 - unbounded
296    -1 - iterating
297    -2 - factorization wanted
298    -3 - redo checking without factorization
299    -4 - looks infeasible
300    -5 - looks unbounded
301  */
302  while (problemStatus_<0) {
303    // clear
304    for (iRow=0;iRow<4;iRow++) {
305      rowArray_[iRow]->clear();
306    }
307
308    for (iColumn=0;iColumn<2;iColumn++) {
309      columnArray_[iColumn]->clear();
310    }
311
312    // give matrix (and model costs and bounds a chance to be
313    // refreshed (normally null)
314    matrix_->refresh(this);
315    // If getting nowhere - why not give it a kick
316#if 0
317    // primal perturbation not coded yet
318    if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
319      perturb();
320#endif
321    // may factorize, checks if problem finished
322    statusOfProblemInPrimal(lastCleaned,factorType);
323
324    // Say good factorization
325    factorType=1;
326
327    // Say no pivot has occurred (for steepest edge and updates)
328    pivotRow_=-2;
329
330    // Save iteration number
331    int saveNumber = numberIterations_;
332
333    // status stays at -1 while iterating, >=0 finished, -2 to invert
334    // status -3 to go to top without an invert
335    while (problemStatus_==-1) {
336#ifdef CLP_DEBUG
337      {
338        int i;
339        // not [1] as has information
340        for (i=0;i<4;i++) {
341          if (i!=1)
342            rowArray_[i]->checkClear();
343        }
344        for (i=0;i<2;i++) {
345          columnArray_[i]->checkClear();
346        }
347      }
348#endif
349#if CLP_DEBUG>2
350      // very expensive
351      if (numberIterations_>0&&numberIterations_<-2534) {
352        handler_->setLogLevel(63);
353        double saveValue = objectiveValue_;
354        double * saveRow1 = new double[numberRows_];
355        double * saveRow2 = new double[numberRows_];
356        memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
357        memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
358        double * saveColumn1 = new double[numberColumns_];
359        double * saveColumn2 = new double[numberColumns_];
360        memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
361        memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
362        createRim(7);
363        gutsOfSolution(rowActivityWork_,columnActivityWork_);
364        printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
365               numberIterations_,
366               saveValue,objectiveValue_,sumPrimalInfeasibilities_);
367        memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
368        memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
369        memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
370        memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
371        delete [] saveRow1;
372        delete [] saveRow2;
373        delete [] saveColumn1;
374        delete [] saveColumn2;
375        objectiveValue_=saveValue;
376      }
377#endif
378#ifdef CLP_DEBUG
379      if(numberIterations_==debugIteration) {
380        printf("dodgy iteration coming up\n");
381      }
382#endif
383      if (!ifValuesPass) {
384        // choose column to come in
385        // can use pivotRow_ to update weights
386        // pass in list of cost changes so can do row updates (rowArray_[1])
387        // NOTE rowArray_[0] is used by computeDuals which is a
388        // slow way of getting duals but might be used
389        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
390                     columnArray_[0],columnArray_[1]);
391      } else {
392        // in values pass
393        if (ifValuesPass>0) {
394          nextSuperBasic(firstSuperBasic);
395          if (firstSuperBasic==numberRows_+numberColumns_)
396            ifValuesPass=-1; // signal end of values pass after this
397        } else {
398          // end of values pass - initialize weights etc
399          primalColumnPivot_->saveWeights(this,5);
400          ifValuesPass=0;
401          if(saveNumber != numberIterations_) {
402            problemStatus_=-2; // factorize now
403            pivotRow_=-1; // say no weights update
404            break;
405          }
406
407          // and get variable
408          primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
409                     columnArray_[0],columnArray_[1]);
410        }
411      }
412      pivotRow_=-1;
413      sequenceOut_=-1;
414      rowArray_[1]->clear();
415      if (sequenceIn_>=0) {
416        // we found a pivot column
417#ifdef CLP_DEBUG
418        if ((handler_->logLevel()&32)) {
419          char x = isColumn(sequenceIn_) ? 'C' :'R';
420          std::cout<<"pivot column "<<
421                  x<<sequenceWithin(sequenceIn_)<<std::endl;
422        }
423#endif
424        // update the incoming column
425        unpack(rowArray_[1]);
426        // save reduced cost
427        double saveDj = dualIn_;
428        factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
429        // do ratio test and re-compute dj
430        primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
431                  ifValuesPass);
432        if (ifValuesPass) {
433          saveDj=dualIn_;
434          if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
435            if(fabs(dualIn_)<1.0e2*dualTolerance_) {
436              // try other way
437              directionIn_=-directionIn_;
438              primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
439                        0);
440            }
441            if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
442              // reject it
443              char x = isColumn(sequenceIn_) ? 'C' :'R';
444              handler_->message(CLP_SIMPLEX_FLAG,messages_)
445                <<x<<sequenceWithin(sequenceIn_)
446                <<OsiMessageEol;
447              setFlagged(sequenceIn_);
448              rowArray_[1]->clear();
449              pivotRow_=-1;
450              continue;
451            }
452          }
453        }
454
455#ifdef CLP_DEBUG
456        if ((handler_->logLevel()&32))
457          printf("btran dj %g, ftran dj %g\n",saveDj,dualIn_);
458#endif
459        if (saveDj*dualIn_<1.0e-20||
460            fabs(saveDj-dualIn_)>1.0e-7*(1.0+fabs(dualIn_))) {
461          handler_->message(CLP_PRIMAL_DJ,messages_)
462            <<saveDj<<dualIn_
463            <<OsiMessageEol;
464          if(saveNumber != numberIterations_) {
465            problemStatus_=-2; // factorize now
466            rowArray_[1]->clear();
467            pivotRow_=-1; // say no weights update
468            break;
469          } else {
470            // take on more relaxed criterion
471            if (saveDj*dualIn_<1.0e-20||
472                fabs(saveDj-dualIn_)>1.0e-4*(1.0+fabs(dualIn_))) {
473              // need to reject something
474              char x = isColumn(sequenceIn_) ? 'C' :'R';
475              handler_->message(CLP_SIMPLEX_FLAG,messages_)
476                <<x<<sequenceWithin(sequenceIn_)
477                <<OsiMessageEol;
478              setFlagged(sequenceIn_);
479              rowArray_[1]->clear();
480              pivotRow_=-1;
481              continue;
482            }
483          }
484        }
485        if (pivotRow_>=0) {
486          // if stable replace in basis
488                                                           pivotRow_,
489                                                           alpha_);
491            // slight error
492            if (factorization_->pivots()>5)
493              problemStatus_=-2; // factorize now
494          } else if (updateStatus==2) {
495            // major error
496            // later we may need to unwind more e.g. fake bounds
497            if(saveNumber != numberIterations_) {
498              problemStatus_=-2; // factorize now
499              break;
500            } else {
501              // need to reject something
502              char x = isColumn(sequenceIn_) ? 'C' :'R';
503              handler_->message(CLP_SIMPLEX_FLAG,messages_)
504                <<x<<sequenceWithin(sequenceIn_)
505                <<OsiMessageEol;
506              setFlagged(sequenceIn_);
507              rowArray_[1]->clear();
508              pivotRow_=-1;
509              continue;
510            }
511          } else if (updateStatus==3) {
512            // out of memory
513            // increase space if not many iterations
514            if (factorization_->pivots()<
515                0.5*factorization_->maximumPivots()&&
516                factorization_->pivots()<200)
517              factorization_->areaFactor(
518                                         factorization_->areaFactor() * 1.1);
519            problemStatus_=-2; // factorize now
520          }
521          // here do part of steepest - ready for next iteration
522          primalColumnPivot_->updateWeights(rowArray_[1]);
523        } else {
524          if (pivotRow_==-1) {
525            // no outgoing row is valid
526#ifdef CLP_DEBUG
527            if (handler_->logLevel()&32)
528              printf("** no row pivot\n");
529#endif
530            if (!factorization_->pivots()) {
531              problemStatus_=-5; //say looks unbounded
532              // do ray
533              delete [] ray_;
534              ray_ = new double [numberColumns_];
535              CoinFillN(ray_,numberColumns_,0.0);
536              int number=rowArray_[1]->getNumElements();
537              int * index = rowArray_[1]->getIndices();
538              double * array = rowArray_[1]->denseVector();
539              double way=-directionIn_;
540              int i;
541              double zeroTolerance=1.0e-12;
542              if (sequenceIn_<numberColumns_)
543                ray_[sequenceIn_]=directionIn_;
544              for (i=0;i<number;i++) {
545                int iRow=index[i];
546                int iPivot=pivotVariable_[iRow];
547                double arrayValue = array[iRow];
548                if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
549                  ray_[iPivot] = way* array[iRow];
550              }
551            }
552            rowArray_[0]->clear();
553            break;
554          } else {
555            // flipping from bound to bound
556          }
557        }
558
559        // update primal solution
560
561        objectiveChange=0.0;
562        // Cost on pivot row may change - may need to change dualIn
563        double oldCost=0.0;
564        if (pivotRow_>=0)
565          oldCost = cost(pivotVariable_[pivotRow_]);
566        // rowArray_[1] is not empty - used to update djs
567        updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
568        if (pivotRow_>=0)
569          dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
570
571        int whatNext=housekeeping(objectiveChange);
572
573        if (whatNext==1) {
574          problemStatus_ =-2; // refactorize
575        } else if (whatNext==2) {
576          // maximum iterations or equivalent
577          problemStatus_= 3;
578          break;
579        }
580      } else {
581        // no pivot column
582#ifdef CLP_DEBUG
583        if (handler_->logLevel()&32)
584          printf("** no column pivot\n");
585#endif
586        break;
587      }
588    }
589#if CLP_DEBUG>2
590    if (numberIterations_>620&&numberIterations_<-2534) {
591      handler_->setLogLevel(63);
592      double saveValue = objectiveValue_;
593      double * saveRow1 = new double[numberRows_];
594      double * saveRow2 = new double[numberRows_];
595      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
596      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
597      double * saveColumn1 = new double[numberColumns_];
598      double * saveColumn2 = new double[numberColumns_];
599      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
600      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
601      createRim(7);
602      gutsOfSolution(rowActivityWork_,columnActivityWork_);
603      printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
604             numberIterations_,
605             saveValue,objectiveValue_,sumPrimalInfeasibilities_);
606      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
607      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
608      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
609      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
610      delete [] saveRow1;
611      delete [] saveRow2;
612      delete [] saveColumn1;
613      delete [] saveColumn2;
614      objectiveValue_=saveValue;
615    }
616#endif
617  }
618
619  // at present we are leaving factorization around
620  // maybe we should empty it
621  deleteRim();
622  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
623    <<objectiveValue()
624    <<OsiMessageEol;
625  // Restore any saved stuff
626  perturbation_ = savePerturbation;
627  infeasibilityCost_ = saveInfeasibilityCost;
628  return problemStatus_;
629}
630/* Checks if finished.  Updates status */
631void
632ClpSimplexPrimal::statusOfProblemInPrimal(int & lastCleaned,int type)
633{
634  if (type==2) {
635    // trouble - restore solution
636    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
637    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
638           numberRows_*sizeof(double));
639    memcpy(columnActivityWork_,savedSolution_ ,
640           numberColumns_*sizeof(double));
641    forceFactorization_=1; // a bit drastic but ..
642    pivotRow_=-1; // say no weights update
644  }
645  int tentativeStatus = problemStatus_;
646
647  if (problemStatus_>-3) {
648    // factorize
649    // later on we will need to recover from singularities
650    // also we could skip if first time
651    // do weights
652    // This may save pivotRow_ for use
653    primalColumnPivot_->saveWeights(this,1);
654    // is factorization okay?
655    if (internalFactorize(1)) {
656      // no - restore previous basis
657      assert (type==1);
658      memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
659      memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
660             numberRows_*sizeof(double));
661      memcpy(columnActivityWork_,savedSolution_ ,
662             numberColumns_*sizeof(double));
663      forceFactorization_=1; // a bit drastic but ..
664      type = 2;
665      assert (internalFactorize(1)==0);
667    }
668    problemStatus_=-3;
669  }
670  // at this stage status is -3 or -5 if looks unbounded
671  // get primal and dual solutions
672  // put back original bounds and then check
673  createRim(7);
674  gutsOfSolution(rowActivityWork_, columnActivityWork_);
675  handler_->message(CLP_SIMPLEX_STATUS,messages_)
676    <<numberIterations_<<objectiveValue();
677  handler_->printing(sumPrimalInfeasibilities_>0.0)
678    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
679  handler_->printing(sumDualInfeasibilities_>0.0)
680    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
681  handler_->printing(numberDualInfeasibilitiesWithoutFree_
682                     <numberDualInfeasibilities_)
683                       <<numberDualInfeasibilities_-
684    numberDualInfeasibilitiesWithoutFree_;
685  handler_->message()<<OsiMessageEol;
686  assert (primalFeasible());
687  // we may wish to say it is optimal even if infeasible
688  bool alwaysOptimal = (specialOptions_&1)!=0;
689  if (dualFeasible()) {
690    if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
691      //may need infeasiblity cost changed
692      // we can see if we can construct a ray
693      // make up a new objective
694      double saveWeight = infeasibilityCost_;
695      // save nonlinear cost as we are going to switch off costs
696      ClpNonLinearCost * nonLinear = nonLinearCost_;
697      infeasibilityCost_=1.0e100;
698      // put back original bounds
699      createRim(7);
700      nonLinearCost_->checkInfeasibilities(true);
701      nonLinearCost_=NULL;
702      // scale
703      int i;
704      for (i=0;i<numberRows_+numberColumns_;i++)
705        cost_[i] *= 1.0e-100;
706      gutsOfSolution(rowActivityWork_, columnActivityWork_);
707      nonLinearCost_=nonLinear;
708      infeasibilityCost_=saveWeight;
709      if (infeasibilityCost_>=1.0e20||
710          numberDualInfeasibilities_==0) {
711        // we are infeasible - use as ray
712        delete [] ray_;
713        ray_ = new double [numberRows_];
714        memcpy(ray_,dual_,numberRows_*sizeof(double));
715        // and get feasible duals
716        infeasibilityCost_=0.0;
717        createRim(7);
718        nonLinearCost_->checkInfeasibilities(true);
719        gutsOfSolution(rowActivityWork_, columnActivityWork_);
720        // so will exit
721        infeasibilityCost_=1.0e30;
722      }
723
724      if (infeasibilityCost_<1.0e20) {
725        infeasibilityCost_ *= 5.0;
727        handler_->message(CLP_PRIMAL_WEIGHT,messages_)
728          <<infeasibilityCost_
729          <<OsiMessageEol;
730        // put back original bounds and then check
731        createRim(7);
732        nonLinearCost_->checkInfeasibilities(true);
733        gutsOfSolution(rowActivityWork_, columnActivityWork_);
734        problemStatus_=-1; //continue
735      } else {
736        // say infeasible
737        problemStatus_ = 1;
738      }
739    } else {
740      // may be optimal
741      if ( lastCleaned!=numberIterations_) {
742        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
743          <<primalTolerance_
744          <<OsiMessageEol;
745        if (numberTimesOptimal_<4) {
746          numberTimesOptimal_++;
748          if (numberTimesOptimal_==1) {
749            // better to have small tolerance even if slower
750            factorization_->zeroTolerance(1.0e-15);
751          }
752          lastCleaned=numberIterations_;
753          handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
754            <<OsiMessageEol;
755          primalTolerance_=dblParam_[OsiPrimalTolerance];
756
757          // put back original bounds and then check
758          createRim(7);
759          nonLinearCost_->checkInfeasibilities(true);
760          gutsOfSolution(rowActivityWork_, columnActivityWork_);
761          problemStatus_ = -1;
762        } else {
763          problemStatus_=0; // optimal
764          if (lastCleaned<numberIterations_) {
765            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
766              <<OsiMessageEol;
767          }
768        }
769      } else {
770        problemStatus_=0; // optimal
771      }
772    }
773  } else {
774    // see if looks unbounded
775    if (problemStatus_==-5) {
776      if (nonLinearCost_->numberInfeasibilities()) {
777        //we need infeasiblity cost changed
778        if (infeasibilityCost_<1.0e20) {
779          infeasibilityCost_ *= 5.0;
781          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
782            <<infeasibilityCost_
783            <<OsiMessageEol;
784          // put back original bounds and then check
785          createRim(7);
786          gutsOfSolution(rowActivityWork_, columnActivityWork_);
787          problemStatus_=-1; //continue
788        } else {
789          // say unbounded
790          problemStatus_ = 2;
791        }
792      } else {
793        // say unbounded
794        problemStatus_ = 2;
795      }
796    } else {
797      // carry on
798      problemStatus_ = -1;
799    }
800  }
801  if (type==0||type==1) {
802    if (!type) {
803      // create save arrays
804      delete [] saveStatus_;
805      delete [] savedSolution_;
806      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
807      savedSolution_ = new double [numberRows_+numberColumns_];
808    }
809    // save arrays
810    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
811    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
812           numberRows_*sizeof(double));
813    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
814  }
815  // restore weights (if saved) - also recompute infeasibility list
816  if (tentativeStatus>-3)
817    primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
818  else
819    primalColumnPivot_->saveWeights(this,3);
821    problemStatus_=4; // unknown
822  }
823}
824/*
825   Row array has pivot column
826   This chooses pivot row.
827   For speed, we may need to go to a bucket approach when many
828   variables go through bounds
829   On exit rhsArray will have changes in costs of basic variables
830*/
831void
832ClpSimplexPrimal::primalRow(OsiIndexedVector * rowArray,
833                            OsiIndexedVector * rhsArray,
834                            OsiIndexedVector * spareArray,
835                            OsiIndexedVector * spareArray2,
836                            int valuesPass)
837{
838  if (valuesPass) {
839    dualIn_ = cost_[sequenceIn_];
840
841    double * work=rowArray->denseVector();
842    int number=rowArray->getNumElements();
843    int * which=rowArray->getIndices();
844
845    int iIndex;
846
847    for (iIndex=0;iIndex<number;iIndex++) {
848
849      int iRow = which[iIndex];
850      double alpha = work[iRow];
851      int iPivot=pivotVariable_[iRow];
852      dualIn_ -= alpha*cost(iPivot);
853    }
854    // determine direction here
855    if (dualIn_<-dualTolerance_) {
856      directionIn_=1;
857    } else if (dualIn_>dualTolerance_) {
858      directionIn_=-1;
859    } else {
860      // towards nearest bound
861      if (valueIn_-lowerIn_<upperIn_-valueIn_) {
862        directionIn_=-1;
863        dualIn_=dualTolerance_;
864      } else {
865        directionIn_=1;
866        dualIn_=-dualTolerance_;
867      }
868    }
869  }
870
871  // sequence stays as row number until end
872  pivotRow_=-1;
873  int numberSwapped=0;
874  int numberRemaining=0;
875
876  int numberThru =0; // number gone thru a barrier
877  int lastThru =0; // number gone thru a barrier on last time
878
879  double totalThru=0.0; // for when variables flip
880  double acceptablePivot=1.0e-7;
881  if (factorization_->pivots())
882    acceptablePivot=1.0e-5; // if we have iterated be more strict
883  double bestEverPivot=acceptablePivot;
884  int lastPivotRow = -1;
885  double lastPivot=0.0;
886  double lastTheta=1.0e50;
887  int lastNumberSwapped=0;
888
889  // use spareArrays to put ones looked at in
890  // First one is list of candidates
891  // We could compress if we really know we won't need any more
892  // Second array has current set of pivot candidates
893  // with a backup list saved in double * part of indexed vector
894
895  // for zeroing out arrays after
896  int maximumSwapped=0;
897  // pivot elements
898  double * spare;
899  // indices
900  int * index, * indexSwapped;
901  int * saveSwapped;
902  spareArray->clear();
903  spareArray2->clear();
904  spare = spareArray->denseVector();
905  index = spareArray->getIndices();
906  saveSwapped = (int *) spareArray2->denseVector();
907  indexSwapped = spareArray2->getIndices();
908
909  // we also need somewhere for effective rhs
910  double * rhs=rhsArray->denseVector();
911
912  /*
913    First we get a list of possible pivots.  We can also see if the
914    problem looks unbounded.
915
916    At first we increase theta and see what happens.  We start
917    theta at a reasonable guess.  If in right area then we do bit by bit.
918    We save possible pivot candidates
919
920   */
921
922  // do first pass to get possibles
923  // We can also see if unbounded
924  // We also re-compute reduced cost
925
926  dualIn_ = cost_[sequenceIn_];
927
928  double * work=rowArray->denseVector();
929  int number=rowArray->getNumElements();
930  int * which=rowArray->getIndices();
931
932  // we need to swap sign if coming in from ub
933  double way = directionIn_;
934  double maximumMovement;
935  if (way>0.0)
936    maximumMovement = min(1.0e30,upperIn_-valueIn_);
937  else
938    maximumMovement = min(1.0e30,valueIn_-lowerIn_);
939
940  double tentativeTheta = maximumMovement;
941  double upperTheta = maximumMovement;
942
943  int iIndex;
944
945  for (iIndex=0;iIndex<number;iIndex++) {
946
947    int iRow = which[iIndex];
948    double alpha = work[iRow];
949    int iPivot=pivotVariable_[iRow];
950    dualIn_ -= alpha*cost(iPivot);
951    alpha *= way;
952    double oldValue = solution(iPivot);
953    // get where in bound sequence
954    if (alpha>0.0) {
955      // basic variable going towards lower bound
956      double bound = lower(iPivot);
957      oldValue -= bound;
958    } else if (alpha<0.0) {
959      // basic variable going towards upper bound
960      double bound = upper(iPivot);
961      oldValue = bound-oldValue;
962    }
963    double value = oldValue-tentativeTheta*fabs(alpha);
964    assert (oldValue>=-primalTolerance_*1.0001);
965    if (value<-primalTolerance_) {
967      spare[numberRemaining]=alpha;
968      rhs[iRow]=oldValue;
969      index[numberRemaining++]=iRow;
970      double value=oldValue-upperTheta*fabs(alpha);
971      if (value<-primalTolerance_)
972        upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
973    }
974  }
975
976  // we need to keep where rhs non-zeros are
977  int numberInRhs=numberRemaining;
978  memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
979  rhsArray->setNumElements(numberInRhs);
980
981  theta_=maximumMovement;
982
983  double dualCheck = fabs(dualIn_);
984  // but make a bit more pessimistic
985  dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
986
987  bool goBackOne = false;
988
989  if (numberRemaining) {
990
991    // looks like pivoting
992    // now try until reasonable theta
993    tentativeTheta = max(10.0*upperTheta,1.0e-7);
994    tentativeTheta = min(tentativeTheta,maximumMovement);
995
996    // loops increasing tentative theta until can't go through
997
998    while (tentativeTheta <= maximumMovement) {
999      double thruThis = 0.0;
1000
1001      double bestPivot=acceptablePivot;
1002      pivotRow_ = -1;
1003
1004      numberSwapped = 0;
1005
1006      upperTheta = maximumMovement;
1007
1008      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1009
1010        int iRow = index[iIndex];
1011        double alpha = spare[iIndex];
1012        double oldValue = rhs[iRow];
1013        double value = oldValue-tentativeTheta*fabs(alpha);
1014
1015        if (value<-primalTolerance_) {
1016          // how much would it cost to go thru
1017          thruThis += alpha*
1018            nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1019          // goes on swapped list (also means candidates if too many)
1020          indexSwapped[numberSwapped++]=iRow;
1021          if (fabs(alpha)>bestPivot) {
1022            bestPivot=fabs(alpha);
1023            pivotRow_ = iRow;
1024            theta_ = oldValue/bestPivot;
1025          }
1026        } else {
1027          value = oldValue-upperTheta*fabs(alpha);
1028          if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
1029            upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1030        }
1031      }
1032
1033      maximumSwapped = max(maximumSwapped,numberSwapped);
1034
1035      if (totalThru+thruThis>=dualCheck) {
1036        // We should be pivoting in this batch
1037        // so compress down to this lot
1038
1039        int saveNumber = numberRemaining;
1040        numberRemaining=0;
1041        for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1042          int iRow = indexSwapped[iIndex];
1043          spare[numberRemaining]=way*work[iRow];
1044          index[numberRemaining++]=iRow;
1045        }
1046        memset(spare+numberRemaining,0,
1047               (saveNumber-numberRemaining)*sizeof(double));
1048        int iTry;
1049#define MAXTRY 100
1050        // first get ratio with tolerance
1051        for (iTry=0;iTry<MAXTRY;iTry++) {
1052
1053          upperTheta=maximumMovement;
1054          numberSwapped = 0;
1055
1056          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1057
1058            int iRow = index[iIndex];
1059            double alpha = fabs(spare[iIndex]);
1060            double oldValue = rhs[iRow];
1061            double value = oldValue-upperTheta*alpha;
1062
1063            if (value<-primalTolerance_ && alpha>=acceptablePivot)
1064              upperTheta = (oldValue+primalTolerance_)/alpha;
1065
1066          }
1067
1068          // now look at best in this lot
1069          bestPivot=acceptablePivot;
1070          pivotRow_=-1;
1071          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1072
1073            int iRow = index[iIndex];
1074            double alpha = spare[iIndex];
1075            double oldValue = rhs[iRow];
1076            double value = oldValue-upperTheta*fabs(alpha);
1077
1078            if (value<=0) {
1079              // how much would it cost to go thru
1080              totalThru += alpha*
1081                nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1082              // goes on swapped list (also means candidates if too many)
1083              indexSwapped[numberSwapped++]=iRow;
1084              if (fabs(alpha)>bestPivot) {
1085                bestPivot=fabs(alpha);
1086                theta_ = oldValue/bestPivot;
1087                pivotRow_=iRow;
1088              }
1089            } else {
1090              value = oldValue-upperTheta*fabs(alpha);
1091              if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
1092                upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1093            }
1094          }
1095
1096          maximumSwapped = max(maximumSwapped,numberSwapped);
1097          if (bestPivot<0.1*bestEverPivot&&
1098              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1099            // back to previous one
1100            goBackOne = true;
1101            break;
1102          } else if (pivotRow_==-1&&upperTheta>largeValue_) {
1103            if (lastPivot>acceptablePivot) {
1104              // back to previous one
1105              goBackOne = true;
1106            } else {
1107              // can only get here if all pivots too small
1108            }
1109            break;
1110          } else if (totalThru>=dualCheck) {
1111            break; // no point trying another loop
1112          } else {
1113            // skip this lot
1114            nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1115            lastPivotRow=pivotRow_;
1116            lastTheta = theta_;
1117            lastThru = numberThru;
1118            numberThru += numberSwapped;
1119            lastNumberSwapped = numberSwapped;
1120            memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1121            if (bestPivot>bestEverPivot)
1122              bestEverPivot=bestPivot;
1123          }
1124        }
1125        break;
1126      } else {
1127        // skip this lot
1128        nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1129        lastPivotRow=pivotRow_;
1130        lastTheta = theta_;
1131        lastThru = numberThru;
1132        numberThru += numberSwapped;
1133        lastNumberSwapped = numberSwapped;
1134        memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1135        if (bestPivot>bestEverPivot)
1136          bestEverPivot=bestPivot;
1137        totalThru += thruThis;
1138        tentativeTheta = 2.0*upperTheta;
1139      }
1140    }
1141    // can get here without pivotRow_ set but with lastPivotRow
1142    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
1143      // back to previous one
1144      pivotRow_=lastPivotRow;
1145      theta_ = lastTheta;
1146            // undo this lot
1147      nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs);
1148      memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int));
1149      numberSwapped = lastNumberSwapped;
1150    }
1151  }
1152
1153  if (pivotRow_>=0) {
1154
1155#define MINIMUMTHETA 1.0e-12
1156    // will we need to increase tolerance
1157#ifdef CLP_DEBUG
1158    bool found=false;
1159#endif
1160    double largestInfeasibility = primalTolerance_;
1161    if (theta_<MINIMUMTHETA) {
1162      theta_=MINIMUMTHETA;
1163      for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1164        int iRow = indexSwapped[iIndex];
1165#ifdef CLP_DEBUG
1166        if (iRow==pivotRow_)
1167          found=true;
1168#endif
1169        largestInfeasibility = max (largestInfeasibility,
1170                                    -(rhs[iRow]-fabs(work[iRow])*theta_));
1171      }
1172#ifdef CLP_DEBUG
1173      assert(found);
1174      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32))
1175        printf("Primal tolerance increased from %g to %g\n",
1176               primalTolerance_,largestInfeasibility);
1177#endif
1178      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
1179    }
1180    alpha_ = work[pivotRow_];
1181    // translate to sequence
1182    sequenceOut_ = pivotVariable_[pivotRow_];
1183    valueOut_ = solution(sequenceOut_);
1184    if (way<0.0)
1185      theta_ = - theta_;
1186    double newValue = valueOut_ - theta_*alpha_;
1187    if (alpha_*way<0.0) {
1188      directionOut_=-1;      // to upper bound
1189      if (fabs(theta_)>0.1)
1190        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1191      else
1192        upperOut_ = newValue;
1193    } else {
1194      directionOut_=1;      // to lower bound
1195      if (fabs(theta_)>0.1)
1196        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1197      else
1198        lowerOut_ = newValue;
1199    }
1200    dualOut_ = reducedCost(sequenceOut_);
1201  } else if (maximumMovement<1.0e20) {
1202    // flip
1203    pivotRow_ = -2; // so we can tell its a flip
1204    sequenceOut_ = sequenceIn_;
1205    valueOut_ = valueIn_;
1206    dualOut_ = dualIn_;
1207    lowerOut_ = lowerIn_;
1208    upperOut_ = upperIn_;
1209    alpha_ = 0.0;
1210    if (way<0.0) {
1211      directionOut_=1;      // to lower bound
1212      theta_ = lowerOut_ - valueOut_;
1213    } else {
1214      directionOut_=-1;      // to upper bound
1215      theta_ = upperOut_ - valueOut_;
1216    }
1217  }
1218
1219  // clear arrays
1220
1221  memset(spare,0,numberRemaining*sizeof(double));
1222  memset(saveSwapped,0,maximumSwapped*sizeof(int));
1223
1224  // put back original bounds etc
1225  nonLinearCost_->goBackAll(rhsArray);
1226
1227  rhsArray->clear();
1228
1229}
1230/*
1231   Chooses primal pivot column
1232   updateArray has cost updates (also use pivotRow_ from last iteration)
1233   Would be faster with separate region to scan
1234   and will have this (with square of infeasibility) when steepest
1235   For easy problems we can just choose one of the first columns we look at
1236*/
1237void
1239                               OsiIndexedVector * spareRow1,
1240                               OsiIndexedVector * spareRow2,
1241                               OsiIndexedVector * spareColumn1,
1242                               OsiIndexedVector * spareColumn2)
1243{
1245                                               spareRow2,spareColumn1,
1246                                               spareColumn2);
1247  if (sequenceIn_>=0) {
1248    valueIn_=solution_[sequenceIn_];
1249    lowerIn_=lower_[sequenceIn_];
1250    upperIn_=upper_[sequenceIn_];
1251    dualIn_=dj_[sequenceIn_];
1252    if (dualIn_>0.0)
1253      directionIn_ = -1;
1254    else
1255      directionIn_ = 1;
1256  } else {
1257    sequenceIn_ = -1;
1258  }
1259}
1260/* The primals are updated by the given array.
1261   Returns number of infeasibilities.
1262   After rowArray will have list of cost changes
1263*/
1264int
1265ClpSimplexPrimal::updatePrimalsInPrimal(OsiIndexedVector * rowArray,
1266                  double theta,
1267                  double & objectiveChange)
1268{
1269  double * work=rowArray->denseVector();
1270  int number=rowArray->getNumElements();
1271  int * which=rowArray->getIndices();
1272
1273  int newNumber = 0;
1274
1275  nonLinearCost_->setChangeInCost(0.0);
1276  int iIndex;
1277
1278  for (iIndex=0;iIndex<number;iIndex++) {
1279
1280    int iRow = which[iIndex];
1281    double alpha = work[iRow];
1282    int iPivot=pivotVariable_[iRow];
1283    double & value = solutionAddress(iPivot);
1284    double change = theta*alpha;
1285    value -= change;
1286
1287    if (change>0.0) {
1288      // going down
1289      if (value<=lower(iPivot)+primalTolerance_) {
1290        double difference = nonLinearCost_->setOne(iPivot,value);
1291        work[iRow] = difference;
1292        if (difference) {
1293          //change reduced cost on this
1295          which[newNumber++]=iRow;
1296        }
1297      } else {
1298        work[iRow]=0.0;
1299      }
1300    } else {
1301      // going up
1302      if (value>=upper(iPivot)-primalTolerance_) {
1303        double difference = nonLinearCost_->setOne(iPivot,value);
1304        work[iRow] = difference;
1305        if (difference) {
1306          //change reduced cost on this
1308          which[newNumber++]=iRow;
1309        }
1310      } else {
1311        work[iRow]=0.0;
1312      }
1313    }
1314  }
1315  objectiveChange += nonLinearCost_->changeInCost();
1316  rowArray->setNumElements(newNumber);
1317  return 0;
1318}
1319void
1320ClpSimplexPrimal::nextSuperBasic(int & firstSuperBasic)
1321{
1322  int iColumn;
1323  if (firstSuperBasic==numberRows_+numberColumns_) {
1324    // initialization
1325    iColumn=0;
1326  } else {
1327    // normal
1328    sequenceIn_=firstSuperBasic;
1329    valueIn_=solution_[sequenceIn_];
1330    lowerIn_=lower_[sequenceIn_];
1331    upperIn_=upper_[sequenceIn_];
1332    dualIn_=dj_[sequenceIn_];
1333    iColumn=firstSuperBasic+1;
1334  }
1335  for (;iColumn<numberRows_+numberColumns_;iColumn++) {
1336    if (getStatus(iColumn)==ClpSimplex::superBasic) {
1337      // is it really super basic
1338      if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
1339        solution_[iColumn]=lower_[iColumn];
1340        setStatus(iColumn,ClpSimplex::atLowerBound);
1341      } else if (fabs(solution_[iColumn]-upper_[iColumn])
1342                 <=primalTolerance_) {
1343        solution_[iColumn]=upper_[iColumn];
1344        setStatus(iColumn,ClpSimplex::atUpperBound);
1345      } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
1346        setStatus(iColumn,ClpSimplex::isFree);
1347      } else {
1348        break;
1349      }
1350    }
1351  }
1352  firstSuperBasic = iColumn;
1353}
1354// Perturbs problem
1355void
1356ClpSimplexPrimal::perturb()
1357{
1358  if (perturbation_>100)