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

Last change on this file since 15 was 15, checked in by forrest, 17 years ago

Hope this works from wincvs

Fix error in values pass

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