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

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

Breaking out whileIterating

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