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

Last change on this file since 10 was 8, checked in by forrest, 18 years ago

Improving reliability using IBM Burlington problems (Primal)

• Property svn:eol-style set to `native`
• Property svn:keywords set to `Author Date Id Revision`
File size: 42.4 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    // for this we need clean basis so it is after factorize
266    numberThrownOut=gutsOfSolution(rowActivityWork_,columnActivityWork_,
267                                   ifValuesPass);
268
269    // find first superbasic - columns, then rows
270    if (ifValuesPass) {
271      nextSuperBasic(firstSuperBasic);
272      if (firstSuperBasic==numberRows_+numberColumns_)
273        ifValuesPass=0; // signal no values pass
274    }
275  }
276
277  problemStatus_ = -1;
278  numberIterations_=0;
279
280  int lastCleaned=0; // last time objective or bounds cleaned up
281
282  // number of times we have declared optimality
283  numberTimesOptimal_=0;
284
285  // Say no pivot has occurred (for steepest edge and updates)
286  pivotRow_=-2;
287
288  // This says whether to restore things etc
289  int factorType=0;
290  // Save iteration number
291  int saveNumber = -1;
292  /*
293    Status of problem:
294    0 - optimal
295    1 - infeasible
296    2 - unbounded
297    -1 - iterating
298    -2 - factorization wanted
299    -3 - redo checking without factorization
300    -4 - looks infeasible
301    -5 - looks unbounded
302  */
303  while (problemStatus_<0) {
304    // clear
305    for (iRow=0;iRow<4;iRow++) {
306      rowArray_[iRow]->clear();
307    }
308
309    for (iColumn=0;iColumn<2;iColumn++) {
310      columnArray_[iColumn]->clear();
311    }
312
313    // give matrix (and model costs and bounds a chance to be
314    // refreshed (normally null)
315    matrix_->refresh(this);
316    // If getting nowhere - why not give it a kick
317#if 0
318    // primal perturbation not coded yet
319    if (perturbation_<101&&numberIterations_>2*(numberRows_+numberColumns_))
320      perturb();
321#endif
322    // If we have done no iterations - special
323    if (saveNumber==numberIterations_)
324      factorType=3;
325    // may factorize, checks if problem finished
326    statusOfProblemInPrimal(lastCleaned,factorType);
327
328    // Say good factorization
329    factorType=1;
330
331    // Say no pivot has occurred (for steepest edge and updates)
332    pivotRow_=-2;
333
334    // Save iteration number
335    saveNumber = numberIterations_;
336
337    // status stays at -1 while iterating, >=0 finished, -2 to invert
338    // status -3 to go to top without an invert
339    while (problemStatus_==-1) {
340#ifdef CLP_DEBUG
341      {
342        int i;
343        // not [1] as has information
344        for (i=0;i<4;i++) {
345          if (i!=1)
346            rowArray_[i]->checkClear();
347        }
348        for (i=0;i<2;i++) {
349          columnArray_[i]->checkClear();
350        }
351      }
352#endif
353#if CLP_DEBUG>2
354      // very expensive
355      if (numberIterations_>0&&numberIterations_<-2534) {
356        handler_->setLogLevel(63);
357        double saveValue = objectiveValue_;
358        double * saveRow1 = new double[numberRows_];
359        double * saveRow2 = new double[numberRows_];
360        memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
361        memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
362        double * saveColumn1 = new double[numberColumns_];
363        double * saveColumn2 = new double[numberColumns_];
364        memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
365        memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
366        createRim(7);
367        gutsOfSolution(rowActivityWork_,columnActivityWork_);
368        printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
369               numberIterations_,
370               saveValue,objectiveValue_,sumPrimalInfeasibilities_);
371        memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
372        memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
373        memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
374        memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
375        delete [] saveRow1;
376        delete [] saveRow2;
377        delete [] saveColumn1;
378        delete [] saveColumn2;
379        objectiveValue_=saveValue;
380      }
381#endif
382#ifdef CLP_DEBUG
383      if(numberIterations_==debugIteration) {
384        printf("dodgy iteration coming up\n");
385      }
386#endif
387      if (!ifValuesPass) {
388        // choose column to come in
389        // can use pivotRow_ to update weights
390        // pass in list of cost changes so can do row updates (rowArray_[1])
391        // NOTE rowArray_[0] is used by computeDuals which is a
392        // slow way of getting duals but might be used
393        primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
394                     columnArray_[0],columnArray_[1]);
395      } else {
396        // in values pass
397        if (ifValuesPass>0) {
398          nextSuperBasic(firstSuperBasic);
399          if (firstSuperBasic==numberRows_+numberColumns_)
400            ifValuesPass=-1; // signal end of values pass after this
401        } else {
402          // end of values pass - initialize weights etc
403          primalColumnPivot_->saveWeights(this,5);
404          ifValuesPass=0;
405          if(saveNumber != numberIterations_) {
406            problemStatus_=-2; // factorize now
407            pivotRow_=-1; // say no weights update
408            break;
409          }
410
411          // and get variable
412          primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
413                     columnArray_[0],columnArray_[1]);
414        }
415      }
416      pivotRow_=-1;
417      sequenceOut_=-1;
418      rowArray_[1]->clear();
419      if (sequenceIn_>=0) {
420        // we found a pivot column
421#ifdef CLP_DEBUG
422        if ((handler_->logLevel()&32)) {
423          char x = isColumn(sequenceIn_) ? 'C' :'R';
424          std::cout<<"pivot column "<<
425                  x<<sequenceWithin(sequenceIn_)<<std::endl;
426        }
427#endif
428        // update the incoming column
429        unpack(rowArray_[1]);
430        // save reduced cost
431        double saveDj = dualIn_;
432        factorization_->updateColumn(rowArray_[2],rowArray_[1],true);
433        // do ratio test and re-compute dj
434        primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
435                  ifValuesPass);
436        if (ifValuesPass) {
437          saveDj=dualIn_;
438          if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
439            if(fabs(dualIn_)<1.0e2*dualTolerance_) {
440              // try other way
441              directionIn_=-directionIn_;
442              primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
443                        0);
444            }
445            if (pivotRow_==-1||(pivotRow_>=0&&fabs(alpha_)<1.0e-5)) {
446              // reject it
447              char x = isColumn(sequenceIn_) ? 'C' :'R';
448              handler_->message(CLP_SIMPLEX_FLAG,messages_)
449                <<x<<sequenceWithin(sequenceIn_)
450                <<OsiMessageEol;
451              setFlagged(sequenceIn_);
452              lastBadIteration_ = numberIterations_; // say be more cautious
453              rowArray_[1]->clear();
454              pivotRow_=-1;
455              continue;
456            }
457          }
458        }
459
460#ifdef CLP_DEBUG
461        if ((handler_->logLevel()&32))
462          printf("btran dj %g, ftran dj %g\n",saveDj,dualIn_);
463#endif
464        if (saveDj*dualIn_<1.0e-20||
465            fabs(saveDj-dualIn_)>1.0e-5*(1.0+fabs(saveDj))) {
466          handler_->message(CLP_PRIMAL_DJ,messages_)
467            <<saveDj<<dualIn_
468            <<OsiMessageEol;
469          if(saveNumber != numberIterations_) {
470            problemStatus_=-2; // factorize now
471            rowArray_[1]->clear();
472            pivotRow_=-1; // say no weights update
473            break;
474          } else {
475            // take on more relaxed criterion
476            if (saveDj*dualIn_<1.0e-20||
477                fabs(saveDj-dualIn_)>1.0e-4*(1.0+fabs(dualIn_))) {
478              // need to reject something
479              char x = isColumn(sequenceIn_) ? 'C' :'R';
480              handler_->message(CLP_SIMPLEX_FLAG,messages_)
481                <<x<<sequenceWithin(sequenceIn_)
482                <<OsiMessageEol;
483              setFlagged(sequenceIn_);
484              lastBadIteration_ = numberIterations_; // say be more cautious
485              rowArray_[1]->clear();
486              pivotRow_=-1;
487              continue;
488            }
489          }
490        }
491        if (pivotRow_>=0) {
492          // if stable replace in basis
494                                                           pivotRow_,
495                                                           alpha_);
497            // slight error
498            if (factorization_->pivots()>5)
499              problemStatus_=-2; // factorize now
500          } else if (updateStatus==2) {
501            // major error
502            // later we may need to unwind more e.g. fake bounds
503            if(saveNumber != numberIterations_) {
504              problemStatus_=-2; // factorize now
505              break;
506            } else {
507              // need to reject something
508              char x = isColumn(sequenceIn_) ? 'C' :'R';
509              handler_->message(CLP_SIMPLEX_FLAG,messages_)
510                <<x<<sequenceWithin(sequenceIn_)
511                <<OsiMessageEol;
512              setFlagged(sequenceIn_);
513              lastBadIteration_ = numberIterations_; // say be more cautious
514              rowArray_[1]->clear();
515              pivotRow_=-1;
516              continue;
517            }
518          } else if (updateStatus==3) {
519            // out of memory
520            // increase space if not many iterations
521            if (factorization_->pivots()<
522                0.5*factorization_->maximumPivots()&&
523                factorization_->pivots()<200)
524              factorization_->areaFactor(
525                                         factorization_->areaFactor() * 1.1);
526            problemStatus_=-2; // factorize now
527          }
528          // here do part of steepest - ready for next iteration
529          primalColumnPivot_->updateWeights(rowArray_[1]);
530        } else {
531          if (pivotRow_==-1) {
532            // no outgoing row is valid
533#ifdef CLP_DEBUG
534            if (handler_->logLevel()&32)
535              printf("** no row pivot\n");
536#endif
537            if (!factorization_->pivots()) {
538              problemStatus_=-5; //say looks unbounded
539              // do ray
540              delete [] ray_;
541              ray_ = new double [numberColumns_];
542              CoinFillN(ray_,numberColumns_,0.0);
543              int number=rowArray_[1]->getNumElements();
544              int * index = rowArray_[1]->getIndices();
545              double * array = rowArray_[1]->denseVector();
546              double way=-directionIn_;
547              int i;
548              double zeroTolerance=1.0e-12;
549              if (sequenceIn_<numberColumns_)
550                ray_[sequenceIn_]=directionIn_;
551              for (i=0;i<number;i++) {
552                int iRow=index[i];
553                int iPivot=pivotVariable_[iRow];
554                double arrayValue = array[iRow];
555                if (iPivot<numberColumns_&&fabs(arrayValue)>=zeroTolerance)
556                  ray_[iPivot] = way* array[iRow];
557              }
558            }
559            rowArray_[0]->clear();
560            break;
561          } else {
562            // flipping from bound to bound
563          }
564        }
565
566        // update primal solution
567
568        objectiveChange=0.0;
569        // Cost on pivot row may change - may need to change dualIn
570        double oldCost=0.0;
571        if (pivotRow_>=0)
572          oldCost = cost(pivotVariable_[pivotRow_]);
573        // rowArray_[1] is not empty - used to update djs
574        updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
575        if (pivotRow_>=0)
576          dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
577
578        int whatNext=housekeeping(objectiveChange);
579
580        if (whatNext==1) {
581          problemStatus_ =-2; // refactorize
582        } else if (whatNext==2) {
583          // maximum iterations or equivalent
584          problemStatus_= 3;
585          break;
586        }
587      } else {
588        // no pivot column
589#ifdef CLP_DEBUG
590        if (handler_->logLevel()&32)
591          printf("** no column pivot\n");
592#endif
593        if (nonLinearCost_->numberInfeasibilities())
594          problemStatus_=-4; // might be infeasible
595        break;
596      }
597    }
598#if CLP_DEBUG>2
599    if (numberIterations_>620&&numberIterations_<-2534) {
600      handler_->setLogLevel(63);
601      double saveValue = objectiveValue_;
602      double * saveRow1 = new double[numberRows_];
603      double * saveRow2 = new double[numberRows_];
604      memcpy(saveRow1,rowReducedCost_,numberRows_*sizeof(double));
605      memcpy(saveRow2,rowActivityWork_,numberRows_*sizeof(double));
606      double * saveColumn1 = new double[numberColumns_];
607      double * saveColumn2 = new double[numberColumns_];
608      memcpy(saveColumn1,reducedCostWork_,numberColumns_*sizeof(double));
609      memcpy(saveColumn2,columnActivityWork_,numberColumns_*sizeof(double));
610      createRim(7);
611      gutsOfSolution(rowActivityWork_,columnActivityWork_);
612      printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
613             numberIterations_,
614             saveValue,objectiveValue_,sumPrimalInfeasibilities_);
615      memcpy(rowReducedCost_,saveRow1,numberRows_*sizeof(double));
616      memcpy(rowActivityWork_,saveRow2,numberRows_*sizeof(double));
617      memcpy(reducedCostWork_,saveColumn1,numberColumns_*sizeof(double));
618      memcpy(columnActivityWork_,saveColumn2,numberColumns_*sizeof(double));
619      delete [] saveRow1;
620      delete [] saveRow2;
621      delete [] saveColumn1;
622      delete [] saveColumn2;
623      objectiveValue_=saveValue;
624    }
625#endif
626  }
627
628  // if infeasible get real values
629  if (problemStatus_) {
630    infeasibilityCost_=0.0;
631    createRim(7);
632    nonLinearCost_->checkInfeasibilities(true);
633    sumPrimalInfeasibilities_=nonLinearCost_->sumInfeasibilities();
634    numberPrimalInfeasibilities_= nonLinearCost_->numberInfeasibilities();
635  }
636  // at present we are leaving factorization around
637  // maybe we should empty it
638  deleteRim();
639  handler_->message(CLP_SIMPLEX_FINISHED+problemStatus_,messages_)
640    <<objectiveValue()
641    <<OsiMessageEol;
642  // Restore any saved stuff
643  perturbation_ = savePerturbation;
644  infeasibilityCost_ = saveInfeasibilityCost;
645  return problemStatus_;
646}
647/* Checks if finished.  Updates status */
648void
649ClpSimplexPrimal::statusOfProblemInPrimal(int & lastCleaned,int type)
650{
651  if (type==2) {
652    // trouble - restore solution
653    memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
654    memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
655           numberRows_*sizeof(double));
656    memcpy(columnActivityWork_,savedSolution_ ,
657           numberColumns_*sizeof(double));
658    forceFactorization_=1; // a bit drastic but ..
659    pivotRow_=-1; // say no weights update
661  }
662  int tentativeStatus = problemStatus_;
663
664  if (problemStatus_>-3||problemStatus_==-4) {
665    // factorize
666    // later on we will need to recover from singularities
667    // also we could skip if first time
668    // do weights
669    // This may save pivotRow_ for use
670    primalColumnPivot_->saveWeights(this,1);
671    // is factorization okay?
672    if (internalFactorize(1)) {
673      // no - restore previous basis
674      assert (type==1);
675      memcpy(status_ ,saveStatus_,(numberColumns_+numberRows_)*sizeof(char));
676      memcpy(rowActivityWork_,savedSolution_+numberColumns_ ,
677             numberRows_*sizeof(double));
678      memcpy(columnActivityWork_,savedSolution_ ,
679             numberColumns_*sizeof(double));
680      forceFactorization_=1; // a bit drastic but ..
681      type = 2;
682      assert (internalFactorize(1)==0);
684    }
685    if (problemStatus_!=-4)
686      problemStatus_=-3;
687  }
688  // at this stage status is -3 or -5 if looks unbounded
689  // get primal and dual solutions
690  // put back original bounds and then check
691  createRim(7);
692  gutsOfSolution(rowActivityWork_, columnActivityWork_);
693  handler_->message(CLP_SIMPLEX_STATUS,messages_)
694    <<numberIterations_<<objectiveValue();
695  handler_->printing(sumPrimalInfeasibilities_>0.0)
696    <<sumPrimalInfeasibilities_<<numberPrimalInfeasibilities_;
697  handler_->printing(sumDualInfeasibilities_>0.0)
698    <<sumDualInfeasibilities_<<numberDualInfeasibilities_;
699  handler_->printing(numberDualInfeasibilitiesWithoutFree_
700                     <numberDualInfeasibilities_)
701                       <<numberDualInfeasibilities_-
702    numberDualInfeasibilitiesWithoutFree_;
703  handler_->message()<<OsiMessageEol;
704  assert (primalFeasible());
705  // we may wish to say it is optimal even if infeasible
706  bool alwaysOptimal = (specialOptions_&1)!=0;
707  if (dualFeasible()||problemStatus_==-4||type==3) {
708    if (nonLinearCost_->numberInfeasibilities()&&!alwaysOptimal) {
709      //may need infeasiblity cost changed
710      // we can see if we can construct a ray
711      // make up a new objective
712      double saveWeight = infeasibilityCost_;
713      // save nonlinear cost as we are going to switch off costs
714      ClpNonLinearCost * nonLinear = nonLinearCost_;
715      infeasibilityCost_=1.0e100;
716      // put back original bounds
717      createRim(7);
718      nonLinearCost_->checkInfeasibilities(true);
719      nonLinearCost_=NULL;
720      // scale
721      int i;
722      for (i=0;i<numberRows_+numberColumns_;i++)
723        cost_[i] *= 1.0e-100;
724      gutsOfSolution(rowActivityWork_, columnActivityWork_);
725      nonLinearCost_=nonLinear;
726      infeasibilityCost_=saveWeight;
727      if (infeasibilityCost_>=1.0e20||
728          numberDualInfeasibilities_==0) {
729        // we are infeasible - use as ray
730        delete [] ray_;
731        ray_ = new double [numberRows_];
732        memcpy(ray_,dual_,numberRows_*sizeof(double));
733        // and get feasible duals
734        infeasibilityCost_=0.0;
735        createRim(7);
736        nonLinearCost_->checkInfeasibilities(true);
737        gutsOfSolution(rowActivityWork_, columnActivityWork_);
738        // so will exit
739        infeasibilityCost_=1.0e30;
740      }
741
742      if (infeasibilityCost_<1.0e20) {
743        infeasibilityCost_ *= 5.0;
745        handler_->message(CLP_PRIMAL_WEIGHT,messages_)
746          <<infeasibilityCost_
747          <<OsiMessageEol;
748        // put back original bounds and then check
749        createRim(7);
750        nonLinearCost_->checkInfeasibilities(true);
751        gutsOfSolution(rowActivityWork_, columnActivityWork_);
752        problemStatus_=-1; //continue
753      } else {
754        // say infeasible
755        problemStatus_ = 1;
756      }
757    } else {
758      // may be optimal
759      if ( lastCleaned!=numberIterations_) {
760        handler_->message(CLP_PRIMAL_OPTIMAL,messages_)
761          <<primalTolerance_
762          <<OsiMessageEol;
763        if (numberTimesOptimal_<4) {
764          numberTimesOptimal_++;
766          if (numberTimesOptimal_==1) {
767            // better to have small tolerance even if slower
768            factorization_->zeroTolerance(1.0e-15);
769          }
770          lastCleaned=numberIterations_;
771          handler_->message(CLP_PRIMAL_ORIGINAL,messages_)
772            <<OsiMessageEol;
773          primalTolerance_=dblParam_[OsiPrimalTolerance];
774
775          // put back original bounds and then check
776          createRim(7);
777          nonLinearCost_->checkInfeasibilities(true);
778          gutsOfSolution(rowActivityWork_, columnActivityWork_);
779          problemStatus_ = -1;
780        } else {
781          problemStatus_=0; // optimal
782          if (lastCleaned<numberIterations_) {
783            handler_->message(CLP_SIMPLEX_GIVINGUP,messages_)
784              <<OsiMessageEol;
785          }
786        }
787      } else {
788        problemStatus_=0; // optimal
789      }
790    }
791  } else {
792    // see if looks unbounded
793    if (problemStatus_==-5) {
794      if (nonLinearCost_->numberInfeasibilities()) {
795        //we need infeasiblity cost changed
796        if (infeasibilityCost_<1.0e20) {
797          infeasibilityCost_ *= 5.0;
799          handler_->message(CLP_PRIMAL_WEIGHT,messages_)
800            <<infeasibilityCost_
801            <<OsiMessageEol;
802          // put back original bounds and then check
803          createRim(7);
804          gutsOfSolution(rowActivityWork_, columnActivityWork_);
805          problemStatus_=-1; //continue
806        } else {
807          // say unbounded
808          problemStatus_ = 2;
809        }
810      } else {
811        // say unbounded
812        problemStatus_ = 2;
813      }
814    } else {
815      // carry on
816      problemStatus_ = -1;
817    }
818  }
819  if (type==0||type==1) {
820    if (!type) {
821      // create save arrays
822      delete [] saveStatus_;
823      delete [] savedSolution_;
824      saveStatus_ = new unsigned char [numberRows_+numberColumns_];
825      savedSolution_ = new double [numberRows_+numberColumns_];
826    }
827    // save arrays
828    memcpy(saveStatus_,status_,(numberColumns_+numberRows_)*sizeof(char));
829    memcpy(savedSolution_+numberColumns_ ,rowActivityWork_,
830           numberRows_*sizeof(double));
831    memcpy(savedSolution_ ,columnActivityWork_,numberColumns_*sizeof(double));
832  }
833  // restore weights (if saved) - also recompute infeasibility list
834  if (tentativeStatus>-3)
835    primalColumnPivot_->saveWeights(this,(type <2) ? 2 : 4);
836  else
837    primalColumnPivot_->saveWeights(this,3);
839    problemStatus_=4; // unknown
840  }
841}
842/*
843   Row array has pivot column
844   This chooses pivot row.
845   For speed, we may need to go to a bucket approach when many
846   variables go through bounds
847   On exit rhsArray will have changes in costs of basic variables
848*/
849void
850ClpSimplexPrimal::primalRow(OsiIndexedVector * rowArray,
851                            OsiIndexedVector * rhsArray,
852                            OsiIndexedVector * spareArray,
853                            OsiIndexedVector * spareArray2,
854                            int valuesPass)
855{
856  if (valuesPass) {
857    dualIn_ = cost_[sequenceIn_];
858
859    double * work=rowArray->denseVector();
860    int number=rowArray->getNumElements();
861    int * which=rowArray->getIndices();
862
863    int iIndex;
864
865    for (iIndex=0;iIndex<number;iIndex++) {
866
867      int iRow = which[iIndex];
868      double alpha = work[iRow];
869      int iPivot=pivotVariable_[iRow];
870      dualIn_ -= alpha*cost(iPivot);
871    }
872    // determine direction here
873    if (dualIn_<-dualTolerance_) {
874      directionIn_=1;
875    } else if (dualIn_>dualTolerance_) {
876      directionIn_=-1;
877    } else {
878      // towards nearest bound
879      if (valueIn_-lowerIn_<upperIn_-valueIn_) {
880        directionIn_=-1;
881        dualIn_=dualTolerance_;
882      } else {
883        directionIn_=1;
884        dualIn_=-dualTolerance_;
885      }
886    }
887  }
888
889  // sequence stays as row number until end
890  pivotRow_=-1;
891  int numberSwapped=0;
892  int numberRemaining=0;
893
894  int numberThru =0; // number gone thru a barrier
895  int lastThru =0; // number gone thru a barrier on last time
896
897  double totalThru=0.0; // for when variables flip
898  double acceptablePivot=1.0e-7;
899  if (factorization_->pivots())
900    acceptablePivot=1.0e-5; // if we have iterated be more strict
901  double bestEverPivot=acceptablePivot;
902  int lastPivotRow = -1;
903  double lastPivot=0.0;
904  double lastTheta=1.0e50;
905  int lastNumberSwapped=0;
906
907  // use spareArrays to put ones looked at in
908  // First one is list of candidates
909  // We could compress if we really know we won't need any more
910  // Second array has current set of pivot candidates
911  // with a backup list saved in double * part of indexed vector
912
913  // for zeroing out arrays after
914  int maximumSwapped=0;
915  // pivot elements
916  double * spare;
917  // indices
918  int * index, * indexSwapped;
919  int * saveSwapped;
920  spareArray->clear();
921  spareArray2->clear();
922  spare = spareArray->denseVector();
923  index = spareArray->getIndices();
924  saveSwapped = (int *) spareArray2->denseVector();
925  indexSwapped = spareArray2->getIndices();
926
927  // we also need somewhere for effective rhs
928  double * rhs=rhsArray->denseVector();
929
930  /*
931    First we get a list of possible pivots.  We can also see if the
932    problem looks unbounded.
933
934    At first we increase theta and see what happens.  We start
935    theta at a reasonable guess.  If in right area then we do bit by bit.
936    We save possible pivot candidates
937
938   */
939
940  // do first pass to get possibles
941  // We can also see if unbounded
942  // We also re-compute reduced cost
943
944  dualIn_ = cost_[sequenceIn_];
945
946  double * work=rowArray->denseVector();
947  int number=rowArray->getNumElements();
948  int * which=rowArray->getIndices();
949
950  // we need to swap sign if coming in from ub
951  double way = directionIn_;
952  double maximumMovement;
953  if (way>0.0)
954    maximumMovement = min(1.0e30,upperIn_-valueIn_);
955  else
956    maximumMovement = min(1.0e30,valueIn_-lowerIn_);
957
958  double tentativeTheta = maximumMovement;
959  double upperTheta = maximumMovement;
960
961  int iIndex;
962
963  for (iIndex=0;iIndex<number;iIndex++) {
964
965    int iRow = which[iIndex];
966    double alpha = work[iRow];
967    int iPivot=pivotVariable_[iRow];
968    dualIn_ -= alpha*cost(iPivot);
969    alpha *= way;
970    double oldValue = solution(iPivot);
971    // get where in bound sequence
972    if (alpha>0.0) {
973      // basic variable going towards lower bound
974      double bound = lower(iPivot);
975      oldValue -= bound;
976    } else if (alpha<0.0) {
977      // basic variable going towards upper bound
978      double bound = upper(iPivot);
979      oldValue = bound-oldValue;
980    }
981    double value = oldValue-tentativeTheta*fabs(alpha);
982    assert (oldValue>=-primalTolerance_*1.0001);
983    if (value<-primalTolerance_) {
985      spare[numberRemaining]=alpha;
986      rhs[iRow]=oldValue;
987      index[numberRemaining++]=iRow;
988      double value=oldValue-upperTheta*fabs(alpha);
989      if (value<-primalTolerance_)
990        upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
991    }
992  }
993
994  // we need to keep where rhs non-zeros are
995  int numberInRhs=numberRemaining;
996  memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
997  rhsArray->setNumElements(numberInRhs);
998
999  theta_=maximumMovement;
1000
1001  double dualCheck = fabs(dualIn_);
1002  // but make a bit more pessimistic
1003  dualCheck=max(dualCheck-100.0*dualTolerance_,0.99*dualCheck);
1004
1005  bool goBackOne = false;
1006
1007  if (numberRemaining) {
1008
1009    // looks like pivoting
1010    // now try until reasonable theta
1011    tentativeTheta = max(10.0*upperTheta,1.0e-7);
1012    tentativeTheta = min(tentativeTheta,maximumMovement);
1013
1014    // loops increasing tentative theta until can't go through
1015
1016    while (tentativeTheta <= maximumMovement) {
1017      double thruThis = 0.0;
1018
1019      double bestPivot=acceptablePivot;
1020      pivotRow_ = -1;
1021
1022      numberSwapped = 0;
1023
1024      upperTheta = maximumMovement;
1025
1026      for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1027
1028        int iRow = index[iIndex];
1029        double alpha = spare[iIndex];
1030        double oldValue = rhs[iRow];
1031        double value = oldValue-tentativeTheta*fabs(alpha);
1032
1033        if (value<-primalTolerance_) {
1034          // how much would it cost to go thru
1035          thruThis += alpha*
1036            nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1037          // goes on swapped list (also means candidates if too many)
1038          indexSwapped[numberSwapped++]=iRow;
1039          if (fabs(alpha)>bestPivot) {
1040            bestPivot=fabs(alpha);
1041            pivotRow_ = iRow;
1042            theta_ = oldValue/bestPivot;
1043          }
1044        } else {
1045          value = oldValue-upperTheta*fabs(alpha);
1046          if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
1047            upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1048        }
1049      }
1050
1051      maximumSwapped = max(maximumSwapped,numberSwapped);
1052
1053      if (totalThru+thruThis>=dualCheck) {
1054        // We should be pivoting in this batch
1055        // so compress down to this lot
1056
1057        int saveNumber = numberRemaining;
1058        numberRemaining=0;
1059        for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1060          int iRow = indexSwapped[iIndex];
1061          spare[numberRemaining]=way*work[iRow];
1062          index[numberRemaining++]=iRow;
1063        }
1064        memset(spare+numberRemaining,0,
1065               (saveNumber-numberRemaining)*sizeof(double));
1066        int iTry;
1067#define MAXTRY 100
1068        // first get ratio with tolerance
1069        for (iTry=0;iTry<MAXTRY;iTry++) {
1070
1071          upperTheta=maximumMovement;
1072          numberSwapped = 0;
1073
1074          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1075
1076            int iRow = index[iIndex];
1077            double alpha = fabs(spare[iIndex]);
1078            double oldValue = rhs[iRow];
1079            double value = oldValue-upperTheta*alpha;
1080
1081            if (value<-primalTolerance_ && alpha>=acceptablePivot)
1082              upperTheta = (oldValue+primalTolerance_)/alpha;
1083
1084          }
1085
1086          // now look at best in this lot
1087          bestPivot=acceptablePivot;
1088          pivotRow_=-1;
1089          for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1090
1091            int iRow = index[iIndex];
1092            double alpha = spare[iIndex];
1093            double oldValue = rhs[iRow];
1094            double value = oldValue-upperTheta*fabs(alpha);
1095
1096            if (value<=0) {
1097              // how much would it cost to go thru
1098              totalThru += alpha*
1099                nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1100              // goes on swapped list (also means candidates if too many)
1101              indexSwapped[numberSwapped++]=iRow;
1102              if (fabs(alpha)>bestPivot) {
1103                bestPivot=fabs(alpha);
1104                theta_ = oldValue/bestPivot;
1105                pivotRow_=iRow;
1106              }
1107            } else {
1108              value = oldValue-upperTheta*fabs(alpha);
1109              if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
1110                upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1111            }
1112          }
1113
1114          maximumSwapped = max(maximumSwapped,numberSwapped);
1115          if (bestPivot<0.1*bestEverPivot&&
1116              bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1117            // back to previous one
1118            goBackOne = true;
1119            break;
1120          } else if (pivotRow_==-1&&upperTheta>largeValue_) {
1121            if (lastPivot>acceptablePivot) {
1122              // back to previous one
1123              goBackOne = true;
1124            } else {
1125              // can only get here if all pivots too small
1126            }
1127            break;
1128          } else if (totalThru>=dualCheck) {
1129            break; // no point trying another loop
1130          } else {
1131            // skip this lot
1132            nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1133            lastPivotRow=pivotRow_;
1134            lastTheta = theta_;
1135            lastThru = numberThru;
1136            numberThru += numberSwapped;
1137            lastNumberSwapped = numberSwapped;
1138            memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1139            if (bestPivot>bestEverPivot)
1140              bestEverPivot=bestPivot;
1141          }
1142        }
1143        break;
1144      } else {
1145        // skip this lot
1146        nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1147        lastPivotRow=pivotRow_;
1148        lastTheta = theta_;
1149        lastThru = numberThru;
1150        numberThru += numberSwapped;
1151        lastNumberSwapped = numberSwapped;
1152        memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1153        if (bestPivot>bestEverPivot)
1154          bestEverPivot=bestPivot;
1155        totalThru += thruThis;
1156        tentativeTheta = 2.0*upperTheta;
1157      }
1158    }
1159    // can get here without pivotRow_ set but with lastPivotRow
1160    if (goBackOne||(pivotRow_<0&&lastPivotRow>=0)) {
1161      // back to previous one
1162      pivotRow_=lastPivotRow;
1163      theta_ = lastTheta;
1164            // undo this lot
1165      nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs);
1166      memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int));
1167      numberSwapped = lastNumberSwapped;
1168    }
1169  }
1170
1171  if (pivotRow_>=0) {
1172
1173#define MINIMUMTHETA 1.0e-12
1174    // will we need to increase tolerance
1175#ifdef CLP_DEBUG
1176    bool found=false;
1177#endif
1178    double largestInfeasibility = primalTolerance_;
1179    if (theta_<MINIMUMTHETA) {
1180      theta_=MINIMUMTHETA;
1181      for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1182        int iRow = indexSwapped[iIndex];
1183#ifdef CLP_DEBUG
1184        if (iRow==pivotRow_)
1185          found=true;
1186#endif
1187        largestInfeasibility = max (largestInfeasibility,
1188                                    -(rhs[iRow]-fabs(work[iRow])*theta_));
1189      }
1190#ifdef CLP_DEBUG
1191      assert(found);
1192      if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32))
1193        printf("Primal tolerance increased from %g to %g\n",
1194               primalTolerance_,largestInfeasibility);
1195#endif
1196      primalTolerance_ = max(primalTolerance_,largestInfeasibility);
1197    }
1198    alpha_ = work[pivotRow_];
1199    // translate to sequence
1200    sequenceOut_ = pivotVariable_[pivotRow_];
1201    valueOut_ = solution(sequenceOut_);
1202    if (way<0.0)
1203      theta_ = - theta_;
1204    double newValue = valueOut_ - theta_*alpha_;
1205    if (alpha_*way<0.0) {
1206      directionOut_=-1;      // to upper bound
1207      if (fabs(theta_)>0.1)
1208        upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1209      else
1210        upperOut_ = newValue;
1211    } else {
1212      directionOut_=1;      // to lower bound
1213      if (fabs(theta_)>0.1)
1214        lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1215      else
1216        lowerOut_ = newValue;
1217    }
1218    dualOut_ = reducedCost(sequenceOut_);
1219  } else if (maximumMovement<1.0e20) {
1220    // flip
1221    pivotRow_ = -2; // so we can tell its a flip
1222    sequenceOut_ = sequenceIn_;
1223    valueOut_ = valueIn_;
1224    dualOut_ = dualIn_;
1225    lowerOut_ = lowerIn_;
1226    upperOut_ = upperIn_;
1227    alpha_ = 0.0;
1228    if (way<0.0) {
1229      directionOut_=1;      // to lower bound
1230      theta_ = lowerOut_ - valueOut_;
1231    } else {
1232      directionOut_=-1;      // to upper bound
1233      theta_ = upperOut_ - valueOut_;
1234    }
1235  }
1236
1237  // clear arrays
1238
1239  memset(spare,0,numberRemaining*sizeof(double));
1240  memset(saveSwapped,0,maximumSwapped*sizeof(int));
1241
1242  // put back original bounds etc
1243  nonLinearCost_->goBackAll(rhsArray);
1244
1245  rhsArray->clear();
1246
1247}
1248/*
1249   Chooses primal pivot column
1250   updateArray has cost updates (also use pivotRow_ from last iteration)
1251   Would be faster with separate region to scan
1252   and will have this (with square of infeasibility) when steepest
1253   For easy problems we can just choose one of the first columns we look at
1254*/
1255void
1257                               OsiIndexedVector * spareRow1,
1258                               OsiIndexedVector * spareRow2,
1259                               OsiIndexedVector * spareColumn1,
1260                               OsiIndexedVector * spareColumn2)
1261{
1263                                               spareRow2,spareColumn1,
1264                                               spareColumn2);
1265  if (sequenceIn_>=0) {
1266    valueIn_=solution_[sequenceIn_];
1267    lowerIn_=lower_[sequenceIn_];
1268    upperIn_=upper_[sequenceIn_];
1269    dualIn_=dj_[sequenceIn_];
1270    if (dualIn_>0.0)
1271      directionIn_ = -1;
1272    else
1273      directionIn_ = 1;
1274  } else {
1275    sequenceIn_ = -1;
1276  }
1277}
1278/* The primals are updated by the given array.
1279   Returns number of infeasibilities.
1280   After rowArray will have list of cost changes
1281*/
1282int
1283ClpSimplexPrimal::updatePrimalsInPrimal(OsiIndexedVector * rowArray,
1284                  double theta,
1285                  double & objectiveChange)
1286{
1287  double * work=rowArray->denseVector();
1288  int number=rowArray->getNumElements();
1289  int * which=rowArray->getIndices();
1290
1291  int newNumber = 0;
1292
1293  nonLinearCost_->setChangeInCost(0.0);
1294  int iIndex;
1295
1296  for (iIndex=0;iIndex<number;iIndex++) {
1297
1298    int iRow = which[iIndex];
1299    double alpha = work[iRow];
1300    int iPivot=pivotVariable_[iRow];
1301    double & value = solutionAddress(iPivot);
1302    double change = theta*alpha;
1303    value -= change;
1304
1305    if (change>0.0) {
1306      // going down
1307      if (value<=lower(iPivot)+primalTolerance_) {
1308        double difference = nonLinearCost_->setOne(iPivot,value);
1309        work[iRow] = difference;
1310        if (difference) {
1311          //change reduced cost on this
1313          which[newNumber++]=iRow;
1314        }
1315      } else {
1316        work[iRow]=0.0;
1317      }
1318    } else {
1319      // going up
1320      if (value>=upper(iPivot)-primalTolerance_) {
1321        double difference = nonLinearCost_->setOne(iPivot,value);
1322        work[iRow] = difference;
1323        if (difference) {
1324          //change reduced cost on this
1326          which[newNumber++]=iRow;
1327        }
1328      } else {
1329        work[iRow]=0.0;
1330      }
1331    }
1332  }
1333  objectiveChange += nonLinearCost_->changeInCost();
1334  rowArray->setNumElements(newNumber);
1335  return 0;
1336}
1337void
1338ClpSimplexPrimal::nextSuperBasic(int & firstSuperBasic)
1339{
1340  int iColumn;
1341  if (firstSuperBasic==numberRows_+numberColumns_) {
1342    // initialization
1343    iColumn=0;
1344  } else {
1345    // normal
1346    sequenceIn_=firstSuperBasic;
1347    valueIn_=solution_[sequenceIn_];
1348    lowerIn_=lower_[sequenceIn_];
1349    upperIn_=upper_[sequenceIn_];
1350    dualIn_=dj_[sequenceIn_];
1351    iColumn=firstSuperBasic+1;
1352  }
1353  for (;iColumn<numberRows_+numberColumns_;iColumn++) {
1354    if (getStatus(iColumn)==ClpSimplex::superBasic) {
1355      // is it really super basic
1356      if (fabs(solution_[iColumn]-lower_[iColumn])<=primalTolerance_) {
1357        solution_[iColumn]=lower_[iColumn];
1358        setStatus(iColumn,ClpSimplex::atLowerBound);
1359      } else if (fabs(solution_[iColumn]-upper_[iColumn])
1360                 <=primalTolerance_) {
1361        solution_[iColumn]=upper_[iColumn];
1362        setStatus(iColumn,ClpSimplex::atUpperBound);
1363      } else if (lower_[iColumn]<-1.0e20&&upper_[iColumn]>1.0e20) {
1364        setStatus(iColumn,ClpSimplex::isFree);
1365      } else {
1366        break;
1367      }
1368    }
1369  }
1370  firstSuperBasic = iColumn;
1371}
1372// Perturbs problem
1373void
1374ClpSimplexPrimal::perturb()
1375{
1376  if (perturbation_>100)