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

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

Changes to make more reliable if problem size changes

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