source: branches/devel-1/ClpSimplexPrimal.cpp @ 7

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

For infeasible problems

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