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

Last change on this file since 19 was 19, checked in by ladanyi, 17 years ago

reordering Clp

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