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

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

Synchronizing

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