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

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

Silly me

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